首页 > 代码库 > xudyh的gcd模板

xudyh的gcd模板

hdu 5019

  1 #include <cstdlib>  2 #include <cctype>  3 #include <cstring>  4 #include <cstdio>  5 #include <cmath>  6 #include <algorithm>  7 #include <vector>  8 #include <string>  9 #include <iostream> 10 #include <map> 11 #include <set> 12 #include <queue> 13 #include <stack> 14 #include <bitset> 15 #include <list> 16 #include <cassert> 17 #include <complex> 18 using namespace std; 19 #define rep(i,a,n) for (int i=a;i<n;i++) 20 #define per(i,a,n) for (int i=n-1;i>=a;i--) 21 #define all(x) (x).begin(),(x).end() 22 //#define fi first 23 #define se second 24 #define SZ(x) ((int)(x).size()) 25 #define TWO(x) (1<<(x)) 26 #define TWOL(x) (1ll<<(x)) 27 #define clr(a) memset(a,0,sizeof(a)) 28 #define POSIN(x,y) (0<=(x)&&(x)<n&&0<=(y)&&(y)<m) 29 typedef vector<int> VI; 30 typedef vector<string> VS; 31 typedef vector<double> VD; 32 typedef long long ll; 33 typedef long double LD; 34 typedef pair<int,int> PII; 35 typedef pair<ll,ll> PLL; 36 typedef vector<ll> VL; 37 typedef vector<PII> VPII; 38 typedef complex<double> CD; 39 const int inf=0x20202020; 40 const ll mod=1000000007; 41 const double eps=1e-9; 42  43 ll powmod(ll a,ll b)             //return (a*b)%mod 44 {ll res=1;a%=mod;for(;b;b>>=1){if(b&1)res=res*a%mod;a=a*a%mod;}return res;} 45 ll powmod(ll a,ll b,ll mod)     //return (a*b)%mod 46 {ll res=1;a%=mod;for(;b;b>>=1){if(b&1)res=res*a%mod;a=a*a%mod;}return res;} 47 ll gcd(ll a,ll b)                 //return gcd(a,b) 48 { return b?gcd(b,a%b):a;} 49 // head 50  51 namespace Factor { 52     const int N=1010000; 53     ll C,fac[10010],n,mut,a[1001000]; 54     int T,cnt,i,l,prime[N],p[N],psize,_cnt; 55     ll _e[100],_pr[100]; 56     vector<ll> d; 57  58     inline ll mul(ll a,ll b,ll p) {     //return (a*b)%p 59         if (p<=1000000000) return a*b%p; 60         else if (p<=1000000000000ll) return (((a*(b>>20)%p)<<20)+(a*(b&((1<<20)-1))))%p; 61         else { 62             ll d=(ll)floor(a*(long double)b/p+0.5); 63             ll ret=(a*b-d*p)%p; 64             if (ret<0) ret+=p; 65             return ret; 66         } 67     } 68  69     void prime_table(){     //prime[1..tot]: prime[i]=ith prime 70         int i,j,tot,t1; 71         for (i=1;i<=psize;i++) p[i]=i; 72         for (i=2,tot=0;i<=psize;i++){ 73             if (p[i]==i) prime[++tot]=i; 74             for (j=1;j<=tot && (t1=prime[j]*i)<=psize;j++){ 75                 p[t1]=prime[j]; 76                 if (i%prime[j]==0) break; 77             } 78         } 79     } 80  81     void init(int ps) {                 //initial 82         psize=ps; 83         prime_table(); 84     } 85  86     ll powl(ll a,ll n,ll p) {           //return (a^n)%p 87         ll ans=1; 88         for (;n;n>>=1) { 89             if (n&1) ans=mul(ans,a,p); 90             a=mul(a,a,p); 91         } 92         return ans; 93     } 94  95     bool witness(ll a,ll n) { 96         int t=0; 97         ll u=n-1; 98         for (;~u&1;u>>=1) t++; 99         ll x=powl(a,u,n),_x=0;100         for (;t;t--) {101             _x=mul(x,x,n);102             if (_x==1 && x!=1 && x!=n-1) return 1;103             x=_x;104         }105         return _x!=1;106     }107 108     bool miller(ll n) {109         if (n<2) return 0;110         if (n<=psize) return p[n]==n;111         if (~n&1) return 0;112         for (int j=0;j<=7;j++) if (witness(rand()%(n-1)+1,n)) return 0;113         return 1;114     }115 116     ll gcd(ll a,ll b) {117         ll ret=1;118         while (a!=0) {119             if ((~a&1) && (~b&1)) ret<<=1,a>>=1,b>>=1;120             else if (~a&1) a>>=1; else if (~b&1) b>>=1;121             else {122                 if (a<b) swap(a,b);123                 a-=b;124             }125         }126         return ret*b;127     }128 129     ll rho(ll n) {130         for (;;) {131             ll X=rand()%n,Y,Z,T=1,*lY=a,*lX=lY;132             int tmp=20;133             C=rand()%10+3;134             X=mul(X,X,n)+C;*(lY++)=X;lX++;135             Y=mul(X,X,n)+C;*(lY++)=Y;136             for(;X!=Y;) {137                 ll t=X-Y+n;138                 Z=mul(T,t,n);139                 if(Z==0) return gcd(T,n);140                 tmp--;141                 if (tmp==0) {142                     tmp=20;143                     Z=gcd(Z,n);144                     if (Z!=1 && Z!=n) return Z;145                 }146                 T=Z;147                 Y=*(lY++)=mul(Y,Y,n)+C;148                 Y=*(lY++)=mul(Y,Y,n)+C;149                 X=*(lX++);150             }151         }152     }153 154     void _factor(ll n) {155         for (int i=0;i<cnt;i++) {156             if (n%fac[i]==0) n/=fac[i],fac[cnt++]=fac[i];}157         if (n<=psize) {158             for (;n!=1;n/=p[n]) fac[cnt++]=p[n];159             return;160         }161         if (miller(n)) fac[cnt++]=n;162         else {163             ll x=rho(n);164             _factor(x);_factor(n/x);165         }166     }167 168     void dfs(ll x,int dep) {169         if (dep==_cnt) d.push_back(x);170         else {171             dfs(x,dep+1);172             for (int i=1;i<=_e[dep];i++) dfs(x*=_pr[dep],dep+1);173         }174     }175 176     void norm() {177         sort(fac,fac+cnt);178         _cnt=0;179         rep(i,0,cnt) if (i==0||fac[i]!=fac[i-1]) _pr[_cnt]=fac[i],_e[_cnt++]=1;180             else _e[_cnt-1]++;181     }182 183     vector<ll> getd() {184         d.clear();185         dfs(1,0);186         return d;187     }188 189     vector<ll> factor(ll n) {       //return all factors of n        cnt:the number of factors190         cnt=0;191         _factor(n);192         norm();193         return getd();194     }195 196     vector<PLL> factorG(ll n) {197         cnt=0;198         _factor(n);199         norm();200         vector<PLL> d;201         rep(i,0,_cnt) d.push_back(make_pair(_pr[i],_e[i]));202         return d;203     }204 205     bool is_primitive(ll a,ll p) {206         //assert(miller(p));207         vector<PLL> D=factorG(p-1);208         rep(i,0,SZ(D)) if (powmod(a,(p-1)/D[i].first,p)==1) return 0;209         return 1;210     }211 }212 213 ll x,y,k;214 int _;215 int main() {216     Factor::init(200000);217     for (scanf("%d",&_);_;_--) {218         scanf("%I64d%I64d%I64d",&x,&y,&k);219         vector<ll> c=Factor::factor(gcd(x,y));      //c:all factors of gcd(x,y)220         sort(all(c));                               // =all common factors of x and y221         reverse(all(c));222         if (SZ(c)<k) puts("-1"); else printf("%I64d\n",c[k-1]);223     }224 }

 

xudyh的gcd模板