首页 > 代码库 > 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模板
声明:以上内容来自用户投稿及互联网公开渠道收集整理发布,本网站不拥有所有权,未作人工编辑处理,也不承担相关法律责任,若内容有误或涉及侵权可进行投诉: 投诉/举报 工作人员会在5个工作日内联系你,一经查实,本站将立刻删除涉嫌侵权内容。