咱们来安利一个黑科技。(实际上是Claris安利来的ios
好比我如今有一坨询问,每次询问两个不超过n的数的gcd。算法
n大概1kw,询问大概300w(怎么输入就不是个人事了,大不了交互库ide
http://mimuw.edu.pl/~kociumaka/files/stacs2013_slides.pdf测试
http://drops.dagstuhl.de/opus/volltexte/2013/3938/pdf/26.pdf spa
咱们定义一个数k的一种因式分解k=k1*k2*k3为“迷之分解”当且仅当k一、k二、k3为质数或小于等于$\sqrt{k}$ 。3d
咱们发现线筛的时候对于一个数x,设x最小的质因子为p,x/p=g,那么x的“迷之分解”能够经过g的“迷之分解”中三个数最小的一个乘上p获得。code
证实彷佛能够用数学概括法证(然而我证不出来啊blog
而后对于每两个小于等于$\sqrt{n}$ 的数咱们能够打一张gcd表出来。ci
最后若是咱们要询问gcd(x,y),咱们找到x的“迷之分解”,而后若是分解的一部分小于等于$\sqrt{n}$ 那就查表,不然那就是一个质数,分类讨论一下就好了。get
伪代码:
UPD:实际测试了一下随机数据跑得并无沙茶gcd快。多是我实现的姿式不够优越(雾
你们能够测试一下跑gcd(5702887,9227465)这个算法比沙茶gcd不知道快到哪里去了
//跑得比谁都快的gcd? #include <iostream> #include <stdio.h> #include <stdlib.h> #include <algorithm> #include <string.h> #include <vector> #include <math.h> #include <time.h> #include <limits> #include <set> #include <map> using namespace std; const int N=10000000; const int sn=sqrt(N); bool np[N+5]; int ps[N+5],pn=0; int cs[N+5][3]; void xs() { np[1]=cs[1][0]=cs[1][1]=cs[1][2]=1; for(int i=2;i<=N;i++) { if(!np[i]) {cs[i][0]=cs[i][1]=1; cs[i][2]=i; ps[++pn]=i;} for(int j=1;j<=pn&&i*ps[j]<=N;j++) { np[i*ps[j]]=1; int cm=cs[i][0]*ps[j]; if(cm<cs[i][1]) { cs[i*ps[j]][0]=cm; cs[i*ps[j]][1]=cs[i][1]; cs[i*ps[j]][2]=cs[i][2]; } else if(cm<cs[i][2]) { cs[i*ps[j]][0]=cs[i][1]; cs[i*ps[j]][1]=cm; cs[i*ps[j]][2]=cs[i][2]; } else { cs[i*ps[j]][0]=cs[i][1]; cs[i*ps[j]][1]=cs[i][2]; cs[i*ps[j]][2]=cm; } if(i%ps[j]);else break; } } } int gcdd[sn+3][sn+3]; void smgcd() { for(int i=0;i<=sn;i++) gcdd[i][0]=gcdd[0][i]=i; for(int i=1;i<=sn;i++) { for(int j=1;j<=i;j++) gcdd[i][j]=gcdd[j][i]=gcdd[i-j][j]; } } void pre_gcd() {xs(); smgcd();} int gcd(int a,int b) { if(a>N||b>N) { puts("Fuck You\n"); return -1; } int *x=cs[a],g=1; for(int i=0;i<3;i++) { int d; if(x[i]<=sn) d=gcdd[x[i]][b%x[i]]; else if(b%x[i]) d=1; else d=x[i]; g*=d; b/=d; } return g; } int euclid_gcd(int x,int y) { while(y) { int t=x%y; x=y; y=t; } return x; } int tmd=-1; void gc() { if(tmd==-1) tmd=clock(); else { printf("Passed: %dms\n",clock()-tmd); tmd=-1; } } int main() { int seed=time(0); //1kw个随机数测试 int ans; printf("Euclid gcd...\n"); srand(seed); gc(); ans=0; for(int i=1;i<=10000000;i++) { int a=(rand()*32768+rand())%N+1,b=(rand()*32768+rand())%N+1; ans^=euclid_gcd(a,b); } printf("Ans = %d\n",ans); gc(); printf("New gcd...\n"); srand(seed); gc(); pre_gcd(); ans=0; for(int i=1;i<=10000000;i++) { int a=(rand()*32768+rand())%N+1,b=(rand()*32768+rand())%N+1; ans^=gcd(a,b); } printf("Ans = %d\n",ans); gc(); }