模板汇总3.2_数学2
1.(EX)BSGS
①什么是BSGS
用来在$O(sqrt(n))$的时间内求关于$x$的方程 $a^x \equiv b (mod$ $p)$的解
②BSGS的原理
这里分块思想又出现了！
以$m=sqrt(p)$(注意向上取整)为标准分块，打出一张$a^k(k∈[0,m)$ $)$的表，然后开始捣鼓原式
把左边变成$a^{i*m-j}$这样的形式，然后移项之后就变成了$a^{i*m}\equiv a^j+b(mod$ $p)$这样子，于是就可以从枚举$i$然后查表验证了，这样不计算查表的复杂度的话就是$O(sqrt(p))$的复杂度了
然后发现这个只能解决$a,p$是互质的情况，不互质怎么办？
设$g=gcd(a,p)$，然后将原式
$a^x \equiv b(mod$ $p)$
变形为(注意此时若$b$不被$g$整除则无解)
$a^{x-1}*\frac{a}{g} \equiv \frac{b}{g}(mod$ $\frac{p}{g})$
不断按照这种方式规约，直到$gcd(a,p')=1$，这样规约复杂度为$sqrt(p)$，不影响整体的复杂度
③BSGS的复杂度
时间复杂度$O(sqrt(n))$(如果不算表的复杂度的话，因为可以哈希表/unordered_map(C++11) )
④BSGS的具体实现
(不知道为什么SPOJ过不了(WA)，然后懒得写哈希了就用map+O2滚过去了)
1 // luogu-judger-enable-o2 2 #include<map> 3 #include<cmath> 4 #include<cstdio> 5 #include<cstring> 6 #include<algorithm> 7 using namespace std; 8 const int inf=2147483647; 9 int T,x,y,t1,t2,t3; 10 map<long long,int> mmp; 11 int gcd(int a,int b) 12 { 13 return b?gcd(b,a%b):a; 14 } 15 int EXGCD(int a,int b,int &x,int &y) 16 { 17 if(!b) {x=1,y=0; return a;} 18 int tmp=EXGCD(b,a%b,x,y),tep=x; 19 x=y,y=tep-a/b*x; return tmp; 20 } 21 int Inv(int a,int b) 22 { 23 EXGCD(a,b,x,y); 24 return (x%b+b)%b; 25 } 26 long long qpow(long long x,long long k,long long m) 27 { 28 if(!k) return 1%m; 29 if(k==1) return x%m; 30 long long tmp=qpow(x,k/2,m); 31 return k%2?tmp*tmp%m*x%m:tmp*tmp%m; 32 } 33 int BSGS(int bas,long long res,int mod) 34 { 35 long long tmp,tep,sqr=ceil(sqrt(mod)); 36 mmp.clear(),mmp[tmp=res%mod]=0; 37 for(int i=1;i<=sqr;i++) 38 tmp=tmp*bas%mod,mmp[tmp]=i; 39 tmp=1,tep=qpow(bas,sqr,mod); 40 for(int i=1;i<=sqr;i++) 41 { 42 tmp=tmp*tep%mod; 43 if(mmp.count(tmp)) 44 return sqr*i-mmp[tmp]; 45 } 46 return -inf; 47 } 48 int EXBSGS(int bas,long long res,int mod) 49 { 50 int gnm,knm=1,ret=0; 51 while((gnm=gcd(bas,mod))!=1) 52 { 53 if(res%gnm) return -inf; 54 ret++,res/=gnm,mod/=gnm; 55 knm=1ll*bas/gnm*knm%mod; 56 } 57 res*=Inv(knm,mod); 58 if(1%mod==res) return ret; 59 if(bas%mod==0) return res?-inf:ret+1; 60 return BSGS(bas,res,mod)+ret; 61 } 62 int main() 63 { 64 while(true) 65 { 66 scanf("%d%d%d",&t1,&t2,&t3); 67 if(t1||t2||t3) 68 { 69 int ans=EXBSGS(t1,t3,t2); 70 if(ans<0) puts("No Solution"); 71 else printf("%d\n",ans); 72 } 73 else break; 74 } 75 return 0; 76 }
2.Miller-Rabin素性测试
①什么是素性测试
就是判素数，大家都会的$O(sqrt(n))$的那个
②Miller-Rabin素性测试的原理
首先我们需要了解一个东西
费马小定理：若正整数$a$小于$p$且$p$是素数，则有$a^(p-1)$ $mod$ $p=1$，不证。
那么它的逆定理就应该长这样：若所有小于$p$的正整数$a$满足$a^(p-1)$ $mod$ $p=1$，则$p$是素数— —看起来这个逆定理用来判素数很有用，我们只需要枚举所有小于$p$的素数判就可以了，而一般判不了几个素数就能发现$p$不是素数了。或者为了提高效率可以牺牲一下复杂度— —只用几个素数来判断，然后把那些判不出来的打个表......这个东西大概叫做“Fermat素性测试”
遗憾的是它的逆定理是假的，不过好消息是不满足它的逆定理的数一定不是素数，只是满足的也不一定是素数，人们把这些数叫做“伪素数”(于是我们知道打表也是行不通的，在$10^{10}$范围内就有$600$个这样的数— —所有小于它的素数都无法探测出它是一个合数，其中最小的是$561$)，于是有了Miller-Rabin素性测试。
我们又需要了解一个东西：
二次探测定理：对于一个奇素数$p$，在小于它的所有正整数$b$中，只有$1$和$p-1$满足$b^2 \equiv 1(mod$ $p)$，不证
当一个数通过以某数为底的“Fermat素性测试时”，再使用二次探测定理探测一次(下面说怎么做)，当然这样还是可能有一些“强伪素数”通过了探测，不过能通过以$2$和$3$为底的第一个强伪素数已经达到了$1373653$，.所以多搞几个素数或者找个素数表从里面多rand几次基本不会出错，事实上：
Miller-Rabin素性探测失败的概率为$\frac{1}{4^k}$，其中$k$是选取的底数的个数
然而我们还没说如何用二次探测定理探测，具体来说是这样的：
1.将$p−1$提取出表示为$x*2^k$的形式
2.选取一个底数$b$并求出$b'=b^{x}(mod$ $p)$
3.将$b'$再平方$k$次，如果结果中未出现$p-1$，则这个数未通过二次探测
③Miller-Rabin素性测试的复杂度
时间复杂度$O(klog$ $p)$
④Miller-Rabin素性测试的具体实现
中间压得有点丑=。=
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 using namespace std; 5 const int pri[8]={0,2,3,7,23,59,101,233}; 6 long long n,T,rd; 7 long long bmul(long long a,long long b,long long m) 8 { 9 long long tmp=((long double)a/m*b+1e-7); 10 long long ret=a*b-tmp*m; 11 return ret<0?ret+m:ret; 12 } 13 long long qpow(long long x,long long k,long long p) 14 { 15 if(k==1) return x%p; 16 long long tmp=qpow(x,k/2,p); 17 return k%2?bmul(bmul(tmp,tmp,p),x,p):bmul(tmp,tmp,p); 18 } 19 bool Miller_Rabin(long long x) 20 { 21 for(int i=1;i<=7;i++) 22 if(x==pri[i]) return true; 23 if(x%2==0||x==1) return false; 24 long long nm=x-1,pw=0; 25 while(nm%2==0) nm/=2,pw++; 26 for(int i=1;i<=7;i++) 27 { 28 long long fm=qpow(pri[i],nm,x); 29 if(fm==1||fm==x-1) continue; 30 bool mr=false; 31 for(int j=1;j<=pw;j++) 32 { 33 fm=bmul(fm,fm,x); 34 if(fm==x-1) {mr=true; break;} 35 } 36 if(!mr) return false; 37 } 38 return true; 39 } 40 int main() 41 { 42 scanf("%lld%lld",&n,&T); 43 while(T--) 44 scanf("%lld",&rd),Miller_Rabin(rd)?puts("Yes"):puts("No"); 45 return 0; 46 }
3.Pollard-Rho因数分解
①Pollard-Rho是干什么的
在期望复杂度$O(x^{\frac{1}{4}})$的复杂度下找到$x$的一个因数
②Pollard-Rho的具体实现
首先Miller-Rabin判素，不谈
然后考虑我们最朴素的分解方法：我们从$1$到$sqrt(n)$枚举然后试除
然后发现这样很慢，考虑优化这个过程，我们随机一个数$x$，然后求$gcd(x,n)$，如果大于$1$就递归分解。
不是咋办？再rand一个显然很睿智，我们用一种类似于在模$x$剩余系下“行走”的方法，用“距离”表示这个用来尝试的数：
找两个数$stp1$和$stp2$，每次$stp1$平方后加上一个$[1,x-1]$范围内的“步长”，然后用$abs(stp1-stp2)$乘上当前的“距离”，这样“走”2的整次幂步后看看距离是否和$x$有公因数，没有就用$stp1$代替$stp2$，然后“步数”加倍再继续走。首先这样一定是能找到$x$的因数的，然后根据生日悖论这样的期望复杂度是$O(x^{\frac{1}{4}})$，显然我不会证明=。=
③Pollard-Rho的复杂度
上面说了。。。
④Pollard-Rho的具体实现
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 using namespace std; 5 const int pri[8]={0,2,3,7,23,59,101,233}; 6 long long T,rd,maxf; 7 long long bmul(long long a,long long b,long long m) 8 { 9 long long tmp=((long double)a/m*b+1e-7); 10 long long ret=a*b-tmp*m; 11 return ret<0?ret+m:ret; 12 } 13 long long qpow(long long x,long long k,long long p) 14 { 15 if(k==1) return x%p; 16 long long tmp=qpow(x,k/2,p); 17 return k%2?bmul(bmul(tmp,tmp,p),x,p):bmul(tmp,tmp,p); 18 } 19 long long gcd(long long a,long long b) 20 { 21 return b?gcd(b,a%b):a; 22 } 23 bool Miller_Rabin(long long x) 24 { 25 for(int i=1;i<=7;i++) 26 if(x==pri[i]) return true; 27 if(x%2==0||x==1) return false; 28 long long nm=x-1,pw=0; 29 while(nm%2==0) nm/=2,pw++; 30 for(int i=1;i<=7;i++) 31 { 32 long long fm=qpow(pri[i],nm,x); 33 if(fm==1||fm==x-1) continue; 34 bool mr=false; 35 for(int j=1;j<=pw;j++) 36 { 37 fm=bmul(fm,fm,x); 38 if(fm==x-1) {mr=true; break;} 39 } 40 if(!mr) return false; 41 } 42 return true; 43 } 44 long long Pollard_Rho(long long x) 45 { 46 long long stp1=0,stp2=0,ret=1; 47 long long dist=1,cnt=2,walk=rand()%(x-1)+1; 48 while(true) 49 { 50 for(int i=1;i<=cnt;i++) 51 { 52 stp1=bmul(stp1,stp1,x)+walk; 53 dist=bmul(dist,abs(stp1-stp2),x); 54 } 55 if(ret>1||(ret=gcd(dist,x))>1) break; 56 cnt*=2,stp2=stp1,dist=1; 57 } 58 if(ret==x) 59 { 60 ret=1; 61 while(ret==1) 62 { 63 stp1=bmul(stp1,stp1,x)+walk; 64 ret=gcd(abs(stp1-stp2),x); 65 } 66 } 67 return ret; 68 } 69 void Decompose(long long x) 70 { 71 if(x==1||x<=maxf) return; 72 if(Miller_Rabin(x)) 73 { 74 maxf=max(x,maxf); 75 return ; 76 } 77 long long p=x; 78 while(p==x) p=Pollard_Rho(x); 79 Decompose(p),Decompose(x/p); 80 } 81 int main() 82 { 83 scanf("%lld",&T); 84 while(T--) 85 { 86 scanf("%lld",&rd); 87 maxf=0,Decompose(rd); 88 if(maxf==rd) puts("Prime"); 89 else printf("%lld\n",maxf); 90 } 91 return 0; 92 }
4.杜教筛
①这是啥
一种亚线性筛，可以$O(n^{\frac{2}{3}})$筛出一些特定的数论函数的前缀和
②实现原理
构造一个函数和要求的数论函数卷起来，对卷积求前缀和
$\sum\limits_{i=1}^n(f×g)(i)$
$=\sum\limits_{i=1}^n\sum_{d|i}f(d)g(\frac{i}{d})$
$=\sum\limits_{i=1}^ng(i)\sum\limits_{j=1}^{\left\lfloor\frac{n}{i}\right\rfloor}f(j)$
下面为了方便设$f$的前缀和为$S$
$=\sum\limits_{i=1}^ng(i)S(\left\lfloor\frac{n}{i}\right\rfloor)$
那怎么搞出来我们要求的前缀和呢？直接作差
$ans=\sum\limits_{i=1}^ng(i)S(\left\lfloor\frac{n}{i}\right\rfloor)-\sum\limits_{i=2}^ng(i)S(\left\lfloor\frac{n}{i}\right\rfloor)$
前面那一项换回刚才的$\sum\limits_{i=1}^n(f×g)(i)$，那么我们只需要快速算$g$和$f×g$的前缀和，然后递归下去就可以了。根据(我不会的)积分知识，当我们预处理出前$O(n^{\frac{2}{3}})$的$f$时复杂度最好，为$O(n^{\frac{2}{3}})$
顺便说一句常见的两个函数φ和μ的配法：都是卷常函数，卷出来分别是单位函数和狄利克雷单位函数，很好求
③时间复杂度
$O(n^{\frac{2}{3}})$
④具体实现
1 #include<map> 2 #include<cstdio> 3 #include<cstring> 4 #include<algorithm> 5 using namespace std; 6 const int N=5e6+50,Edge=3e6,mod=1e9+7; 7 int T,n,cnt,pri[N],npr[N],mul[N]; 8 long long phi[N]; map<int,long long> mp1,mp2; 9 void Prework() 10 { 11 mul[1]=phi[1]=1,npr[1]=true; 12 for(int i=2;i<=Edge;i++) 13 { 14 if(!npr[i]) pri[++cnt]=i,mul[i]=-1,phi[i]=i-1; 15 for(int j=1;j<=cnt&&1ll*i*pri[j]<=Edge;j++) 16 { 17 npr[i*pri[j]]=true; 18 phi[i*pri[j]]=phi[i]*pri[j]; 19 if(i%pri[j]) 20 mul[i*pri[j]]=-mul[i], 21 phi[i*pri[j]]-=phi[i]; 22 else 23 break; 24 } 25 phi[i]+=phi[i-1],mul[i]+=mul[i-1]; 26 } 27 } 28 long long Dusieve_phi(int x) 29 { 30 if(x<=Edge) return phi[x]; 31 if(mp1.count(x)) return mp1[x]; 32 long long ret=1ll*(x+1)*x/2; 33 for(int i=2,j;i<=x;i=j+1) 34 j=x/(x/i),ret-=1ll*Dusieve_phi(x/i)*(j-i+1); 35 return mp1[x]=ret; 36 } 37 int Dusieve_mul(int x) 38 { 39 if(x<=Edge) return mul[x]; 40 if(mp2.count(x)) return mp2[x]; 41 int ret=1; 42 for(int i=2,j;i<=x;i=j+1) 43 j=x/(x/i),ret-=Dusieve_mul(x/i)*(j-i+1); 44 return mp2[x]=ret; 45 } 46 int main() 47 { 48 Prework(),scanf("%d",&T); 49 while(T--) 50 { 51 scanf("%d",&n); 52 long long ans1=Dusieve_phi(n); 53 int ans2=Dusieve_mul(n); 54 printf("%lld %d\n",ans1,ans2); 55 } 56 return 0; 57 }
浙公网安备 33010602011771号