好像有不少更新：）
本文主要记录一些不是那么熟悉的高级数论算法的推导与应用。
exBSGS算法
解决模数、底数不互质的离散对数问题。
（1）为何(BSGS)算法不再适用：(A)不一定存在逆元，而且无法保证解的循环性。
（2）无解的结论：
设方程为(A^x=B pmod{P})
当 ((A,P)
mid B)且(B
e 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)算法了。
int BSGS(int a,int b,int P)
{
int s,o;
S.clear();
int m=(int)ceil(sqrt(P));
o=a;
s=b;
for(int i=1; i<=m; i++)
{
s=(ll)s*o%P;
S[s]=i;
}
o=fpow(a,m,P);
s=1;
for(int i=1; i<=m; i++)
{
s=(ll)s*o%P;
if((it=S.find(s))!=S.end())
return i*m-it->second;
}
return -1;
}
int exbsgs(int a,int b,int c)
{
a%=c; b%=c;
if(b==1||c==1)return 0;
if(!a)return b?-1:1;
int d=gcd(a,c);
if(d==1)return BSGS(a,b,c);
if(b%d!=0)return -1;
int o=exbsgs(a,(ll)b/d*Inv(a/d,c/d)%(c/d),c/d);
return o==-1?-1:o+1;
}
快速计算阶乘
即快速计算：
对于质数(p)，把(N!)去掉所有(p)因子的部分对(p^e)取模，(p^ele 10^5)
（1）为了便于统计出现了多少个p的次幂，先将阶乘中所有p的倍数提出来。可以简单算出：共有
$displaystyle lfloor frac{n}{p}
floor $ 个。
这中间每一项都除去p，
可以得到 (displaystyle lfloor frac{n}{p}
floor !) 。该部分可以选择递归求解。
（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})
两个括号里的暴力计算；
考虑对 (p^e) 以内做预处理，查询复杂度可以做到 (O(log_p n)) 左右。
代码在跟下面一起贴
组合数取模（扩展卢卡斯）
即快速计算(dbinom{N}{M} pmod{P})，(N,Mle 10^9,P不一定是质数,Ple 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}
floor)计算。
然后两部分分开算，乘起来就是对应的(a_i)。
最后CRT合并起来就可以了。
struct bnm {
int p, Mod, fac[1 << 20];
void init(int x, int y)
{
p = x;
Mod = y;
fac[0] = 1;
for(int s = 1; s != Mod; ++s)
if(s % p)
fac[s] = 1LL * fac[s - 1] * s % Mod;
else
fac[s] = fac[s - 1];
}
int Fac(ll x)
{
if(!x)
return 1;
return 1LL * Fac(x / p) * fexp(fac[Mod - 1], x / Mod, Mod) % Mod * fac[x % Mod] % Mod;
}
int operator()(ll n, ll m) {
if(n < m || m < 0) return 0;
int res = Fac(n) * Inv(Fac(m), Mod) % Mod * Inv(Fac(n - m), Mod) % Mod;
int z = 0;
for(ll d = n; d; z += d /= p);
for(ll d = m; d; z -= d /= p);
for(ll d = n - m; d; z -= d /= p);
res = 1LL * res * fexp(p, z, Mod) % Mod;
return res;
}
}Solve;
普通二次剩余
即解方程(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)也是很容易特判的。
二次剩余有一正一负两解，下面代码随便取了一个。
namespace Cipolla
{
inline int Le(int x){return fpow(x,(P-1)/2);}
int w;
struct Cp
{
int x,y;
Cp(int a=0,int b=0){x=a;y=b;}
};
Cp operator+(Cp a,Cp b){return Cp((a.x+b.x)%P,(a.y+b.y)%P);}
Cp operator-(Cp a,Cp b){return Cp((a.x-b.x+P)%P,(a.y-b.y+P)%P);}
Cp operator*(Cp a,Cp b){return Cp(((ll)a.x*b.x+(ll)w*a.y%P*b.y)%P,((ll)a.x*b.y+(ll)a.y*b.x)%P);}
Cp operator^(Cp a,int b){
Cp o(1,0);
for(;b;b>>=1) {
if(b&1)o=o*a;
a=a*a;
}
return o;
}
int calc(int x)
{
x%=P;
int a;
while(1) {
a=rand();
w=((ll)a*a-x+P)%P;
if(Le(w)==P-1) break;
}
Cp s=Cp(a,1)^((P+1)/2);
return min(s.x,P-s.x);
}
}
原根和阶
奇素数与奇素数的方幂有原根，原根 (g) 是阶恰好等于 (p-1) 的数，也就是 (g^x(xin [0,p-2])) 与整数 (yin [1,p-1]) 一一对应。
我们一般用暴力枚举的方法求原根，原理类似试除法，判断是否存在 (d|p-1) 使得 (g^dequiv 1 pmod p) 即可。
bool judge(int x)
{
for(int j = 0; j < tot; j ++) //只需对质因子做
if(fexp(x, (P - 1) / c[j], P) == 1)
return false;
return true;
}
对一个任意模数求阶也不是很难，只要对 (varphi (m)) 不断试除保证 (x^kequiv 1 pmod m)即可。
原根可以转换为单位根，即若 (n|(varphi(p)-1)) 则 (omega_n=g^{frac{varphi(p)-1}{n}})。
常用于单位根反演 ([n|i]=dfrac{1}{n}sum_{k=0}^{n-1} omega_n^{ki})
单位根反演可以应用于 ( ext{NTT}) 以及结合“二项式定理”优化复杂度。
快速乘法取模
某些情况下机器不支持int128，或者毒瘤出题人卡int128常数，就需要写这种 (O(1)) 快速乘。
inline ll mul(ll a,ll b,ll p){
ll res=a*b-((ll)((long double)a*b/p+0.5))*p;
return res<0?res+p:res;
}
一定要保证 (a,blt p) 的时候计算，因为它是由这个保证精度的。
合并同余方程组
(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 imes m)
有一个坑点就是你的(k)不要对(m)取模，
归纳易证你这个(now)一定是最小的正整数解，只要题目保证了啥啥，(now)就肯定不会溢出。
但是好多题目中(m)都会溢出，变成负数，再取模就会出bug。
当然还是不能排除求出假答案的情况，最后可以 check 一遍。
void excrt(ll a,ll b)
{
ll c=((a-now)%b+b)%b;
ll x,y;
ll d=exgcd(M,b,x,y);
if(c%d!=0) err();
c/=d;
b/=d;
x=(x%b+b)%b;
x=(__int128)x*c%b;
now+=x*M;
M*=b;
}
Miller-Rabin素数测试
二次探测原理：
当p是质数时，则一定有
若(x^2 equiv 1 pmod p)，则(x=1 pmod p)或(x=-1 pmod p)。
这样我们就考虑反过来想办法验证(n)是质数。
令(n=2^k imes m + 1)（m是奇数），
拿个小质数来进行二次探测，进行(k)次迭代即可。
取10-12个质数就足够稳了。
代码21行是二次探测，25行是费马小定理。
int ispr(ll n)
{
if(n<2)
return 0;
for(int i=0; i<12; i++)
if(n%Pr[i]==0)
return n==Pr[i];
ll m=n-1,x,y;
int k=0;
while(~m&1)
{
m>>=1;
k++;
}
for(int i=0; i<12; i++)
{
x=fpow(Pr[i],m,n);
for(int j=0; j<k; j++)
{
y=Mul(x,x,n);
if(y==1&&x!=1&&x!=n-1)
return 0;
x=y;
}
if(x!=1)
return 0;
}
return 1;
}
质因数分解 Pollard-rho
先说一下为大整数 (N(Nle 10^{18})) 分解质因子的流程：
（1）弄出来一个 (N) 的非平凡因子 (x)。
（2）递归分解 (x) 和 (N/x)。
而Pollard-rho就是高效地寻找非平凡因子的算法。
Pollard-rho基于随机化实现，然后算法的效率保证是“生日悖论”，这里的含义如下：
如果随机一个数 (x)，那么 (x) 是 (N) 约数的概率极其微小。
但是我们如果随机出 (x,y)，那么 (|x-y|) 是 (N) 约数的概率将显著提高。
（1）不断随机数字(x)，假如(gcd(x,n)>1)直接返回这个(gcd)
额样例都跑不出来，没有任何前途。
（2）按照上面讲的这么写，并按照倍增的形式不断调整：
x=y=0;
for(int k=2; ; k<<=1)
{
for(int i=1; i<=k; i++)
{
x=(x*x+sd)%n;
v=gcd(abs(x-y),n);
if(v!=1 && v!=n)
return v;
if(x==y)
return n;
}
y=x;
}
效率玄学地大幅提升，但是依然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)之内。
其实感觉特别暴力。