本文主要记录一些不是那么熟悉的高级数论算法的推导与应用。
1、exBSGS算法
解决模数、底数不互质的离散对数问题。
（1）为何$BSGS$算法不再适用：$A$不一定存在逆元，而且无法保证解的循环性。
（2）无解的结论：
设方程为$A^x=B \pmod{P}$
当 $(A,P) \nmid B$且$B\ne 1$ 时，原方程无自然数解。
还有就是$A=0,B≠0$这种。
（3）算法流程：
先判无解。
然后若$B=1$，显然$x=0$，特判之。
否则令$G=(A,P)$
若$G>1$，
两边同除以$G$,得$A^{x-1}=\frac{B}{G}*(\frac{A}{G})^{-1} \pmod{\frac{P}{G}}$
迭代至$(A,P)==1$时即可套用$BSGS$算法了。
1 int BSGS(int a,int b,int P) 2 { 3 int s,o; 4 S.clear(); 5 int m=(int)ceil(sqrt(P)); 6 o=a; 7 s=b; 8 for(int i=1; i<=m; i++) 9 { 10 s=(ll)s*o%P; 11 S[s]=i; 12 } 13 o=fpow(a,m,P); 14 s=1; 15 for(int i=1; i<=m; i++) 16 { 17 s=(ll)s*o%P; 18 if((it=S.find(s))!=S.end()) 19 return i*m-it->second; 20 } 21 return -1; 22 } 23 24 int exbsgs(int a,int b,int c) 25 { 26 a%=c; b%=c; 27 if(b==1||c==1)return 0; 28 if(!a)return b?-1:1; 29 int d=gcd(a,c); 30 if(d==1)return BSGS(a,b,c); 31 if(b%d!=0)return -1; 32 int o=exbsgs(a,(ll)b/d*Inv(a/d,c/d)%(c/d),c/d); 33 return o==-1?-1:o+1; 34 }
2、快速阶乘算法
即快速计算：
对于质数$p$，把$N!$去掉所有$p$因子的部分对$p^e$取模，$p^e\le 10^5$
（1）为了便于统计出现了多少个p的次幂，先将阶乘中所有p的倍数提出来。可以简单算出：共有
$\displaystyle \lfloor \frac{n}{p} \rfloor $ 个。
这中间每一项都除去p，
可以得到 $\displaystyle \lfloor \frac{n}{p} \rfloor !$ 。该部分可以选择递归求解。
（2）那么接下来只剩下非p的倍数的几项了。通过简单观察可以知道，剩余几项具有循环节，循环节长度小于$p^e$ 。
原因是剩余项关于p具有循环节，而
$x+p^e \equiv x \pmod{p^e}$，
所以可以一起计算。
结果会剩下几项凑不齐一个循环节，但是这几项长度已经小于$p^e$了，可以选择暴力乘法求解。
有必要举个例子模拟一下：
$N=22,p=3,e=2,$
$22!=1*2*3*......*22$
$=(1*2*4*5*7*8*10*11*13*14*16*17*19*20*22)*3^7*7!$
$7!$递归处理；
而
$(1*2*4*5*7*8*10*11*13*14*16*17*19*20*22)$
$\equiv (1*2*4*5*7*8)^2*(19*20*22) \pmod {3^2}$
两个括号里的暴力算；
时间复杂度为似乎是$O(p^elog(log_p(n-k)))$
以上部分参考洛谷Great_Influence dalao的题解。
参考代码：
1 ll fac(ll n,ll p,ll P) 2 { 3 if(!n) 4 return 1; 5 ll s=1; 6 for(ll i=1;i<=P;i++) 7 if(i%p)s=s*i%P; 8 s=fpow(s,n/P,P); 9 for(ll i=1;i<=n%P;i++) 10 if(i%p)s=s*i%P; 11 s=s*fac(n/p,p,P)%P; 12 return s; 13 }
然后如果P是固定数，那么P以内的阶乘就可以预处理，该算法就飞快。
3、扩展卢卡斯算法
即快速计算$\dbinom{N}{M} \pmod{P}$，$N,M\le 10^9,P不一定是质数,P\le 10^9$,但$P$中任何一个质因子的总积不超过$10^5$
考虑中国剩余定理：
令$Ans=\dbinom{N}{M}$
对质因子分开计算有：
$Ans \equiv a_1 \pmod{p_1^{q_1}}$
$Ans \equiv a_2 \pmod{p_2^{q_2}}$
$......$
对于一个式子：
不含$p_i$的部分使用上述快速阶乘方法求解。
含有$p_i$的重数使用$\displaystyle \sum_{i=1} \lfloor \frac{n}{p^i} \rfloor$计算。
然后两部分分开算，乘起来就是对应的$a_i$。
最后CRT合并起来就可以了。
4、二次剩余
即解方程$x^2=a \pmod{P}$，P是质数。
（1）勒让德符号的定义
$(\frac{a}{p})=1$：a在模p意义下有二次剩余。
$(\frac{a}{p})=-1$：a在模p意义下无二次剩余。
计算公式：$(\frac{a}{p}) = a^{\frac{p-1}{2}} \pmod p$
（2）求解二次剩余
随机一个$a$使得$b=\sqrt{a^2-n}$（此时 $a^2-n$ 无二次剩余，就直接写成根号形式），然后把$b$作为二次域的单位元。
那么$x=(a+b)^{\frac{p+1}{2}}$（证明需要不少性质以及二项式定理啥的，具体参考这里 Acdreamer）
写一个扩域乘法类就好了。
注意这里$p$得是奇质数，但$p=2$也是很容易特判的。
代码：
1 namespace Cipolla 2 { 3 struct II { 4 LL x,y ; 5 II(LL _x=0,LL _y=0) 6 { x=_x;y=_y; } 7 }; 8 9 LL w; 10 11 II Mut (const II& a,const II& b,const LL p) 12 { 13 return II((a.x*b.x%p+a.y*b.y%p*w%p)%p , (a.x*b.y%p+a.y*b.x%p)%p); 14 } 15 16 II Ipow(II a,LL b,LL c) 17 { 18 II r(1,0); 19 while(b) 20 { 21 if(b&1)r=Mut(r,a,c); 22 a=Mut(a,a,c); 23 b>>=1; 24 } 25 return r; 26 } 27 28 LL Le(LL x,LL p) 29 { 30 return fpow(x,(p-1)>>1,p); 31 } 32 33 // x^2 == n (mod p) 34 PP work(LL n,LL p) 35 { 36 LL a; 37 while(1) 38 { 39 a=rand() % p; 40 w=((a*a-n)%p+p)%p; 41 if(Le(w,p) == p-1) 42 break; 43 } 44 II b(a,1); 45 LL r = Ipow(b,(p+1)>>1,p).x % p; 46 if(r > p-r) r = p-r; 47 return mp(r,p-r); 48 } 49 }
（3）顺便提一下：
在模$1000000009$意义下，$5$存在一个二次剩余是$383008016$（也就是说$\sqrt{5} \equiv 383008016 \pmod{1000000009}$）
这样式子里的$\sqrt{5}$就能化成最简单的代数运算，大多在斐波那契数列相关的题目里有应用。
5、原根
定义性质啥的看这里 原根
这里说下用途：
（1）乘法变成指数的加法（N次剩余 --> 线性同余方程）
（2）NTT多项式乘法（大家都会）
（3）单位根反演（......）
6、快速乘法取模
某些情况下机器不支持int128，或者毒瘤出题人卡int128常数，就需要写这种 $O(1)$ 快速乘。
1 inline ll mul(ll a,ll b,ll p){ 2 ll res=a*b-((ll)((long double)a*b/p+0.5))*p; 3 return res<0?res+p:res; 4 }
注意那个 $+0.5$ 不写/写成 $0.1^9$ 之类的都是有可能错误的，随便随一组数据都能卡掉。
至少以上的写法通过了洛谷的 $\text{Pollard_Rho}$ 模板题。
7、合并同余方程组
$x \equiv now \pmod m$
$x \equiv a \pmod b$
由第一式得
$x=now+k*m$
代入二式
$now+km \equiv a \pmod b$
$mk \equiv a-now \pmod b$
拿个扩欧就能解出来最小的$k$
然后$now=now+k\times m$
有一个坑点就是你的$k$不要对$m$取模，
归纳易证你这个$now$一定是最小的正整数解，
只要题目保证了啥啥，$now$就肯定不会溢出。
但是好多题目中$m$都会溢出，变成负数，再取模就会出bug。
其实如果数据保证答案合法，每次都会得到 $c=0$，不影响结果。
1 void excrt(ll a,ll b) 2 { 3 ll c=((a-now)%b+b)%b; 4 ll x,y; 5 ll d=exgcd(M,b,x,y); 6 if(c%d!=0) err(); 7 c/=d; 8 b/=d; 9 x=(x%b+b)%b; 10 x=(u128)x*c%b; 11 now+=x*M; 12 M*=b; 13 }
某些极限情况下不一定能通过，还需要一定特技，如random_shuffle。
8、MillerRabin和Pollard-rho
人类的智慧。
都是使用了玄学方法（其实是不懂233）来增加准确性。
（1）素数测试 MillerRabin
二次探测原理：
当p是质数时，则一定有
若$x^2 \equiv 1 \pmod p$，则$x=1 \pmod p$或$x=-1 \pmod p$。
这样我们就考虑反过来想办法验证$n$是质数。
令$n=2^k \times m + 1$（m是奇数），
拿个小质数来进行二次探测，进行$k$次迭代即可。
取10-12个质数就足够稳了。
代码21行是二次探测，25行是费马小定理。
1 int ispr(ll n) 2 { 3 if(n<2) 4 return 0; 5 for(int i=0; i<12; i++) 6 if(n%Pr[i]==0) 7 return n==Pr[i]; 8 ll m=n-1,x,y; 9 int k=0; 10 while(~m&1) 11 { 12 m>>=1; 13 k++; 14 } 15 for(int i=0; i<12; i++) 16 { 17 x=fpow(Pr[i],m,n); 18 for(int j=0; j<k; j++) 19 { 20 y=Mul(x,x,n); 21 if(y==1&&x!=1&&x!=n-1) 22 return 0; 23 x=y; 24 } 25 if(x!=1) 26 return 0; 27 } 28 return 1; 29 }
（2）质因数分解 Pollard-rho
首先说一下，
@whzzt 增加的几组强力数据，目前网上的很多Pollard-rho都无法通过，然后洛谷时限就2secs，然后那个模板题都成黑题了233
先说一下怎么为大整数 $N(N\le 10^{18})$ 分解质因子：
（1）弄出来一个 $N$ 的非平凡因子 $x$。
（2）递归分解 $x$ 和 $N/x$。
1 x=find(n); 2 while(x==n)x=find(n); 3 sol(x); 4 sol(n/x);
Pollard-rho就是高效的寻找非平凡因子的算法。
Pollard-rho基于随机化实现，然后算法的效率保证是“生日悖论”，这里的含义如下：
如果随机一个数 $x$，那么 $x$ 是 $N$ 约数的概率极其微小。
但是我们如果随机出 $x,y$，那么 $|x-y|$ 是 $N$ 约数的概率将显著提高。
（1）不断随机数字$x$，假如$gcd(x,n)>1$直接返回这个$gcd$
额样例都跑不出来，没有任何前途。
（2）按照上面讲的这么写，并不断调整伪随机数：
1 x=y=0; 2 for(int k=2; ; k<<=1) 3 { 4 for(int i=1; i<=k; i++) 5 { 6 x=(x*x+sd)%n; 7 v=gcd(abs(x-y),n); 8 if(v!=1 && v!=n) 9 return v; 10 if(x==y) 11 return n; 12 } 13 y=x; 14 }
效率玄学地大幅提升，但是依然TLE飞起。。。
（3）然而你发现这个gcd常数巨大，我们把改成用一个模 $N$ 意义下的变量记一下，这样每隔一会查一次就行！
static int Size(1<<7); x=y=0; for(int k=2; ;k<<=1) { v=1; y=x; for(int i=1; i<=k; i++) { x=((u128)x*x+sd)%n; v=((u128)v*_abs(x-y))%n; if(i%Size==0) { z=gcd(v,n); if(z>1) return z; } } z=gcd(v,n); if(z>1) return z; }
就能过了，常数极小。
Pollard-rho时间复杂度是$O(n^{0.25})$
实测每次找质因子内层乘法能稳定在$10^5$之内。
感觉特别暴力。
9、任意模数N次剩余
题目链接点这里 51nod
鸽着
扫一扫
1392
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
