好像有不少更新：）
本文主要記錄一些不是那么熟悉的高級數論算法的推導與應用。
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\)算法了。
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^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}\)
兩個括號里的暴力計算；
考慮對 \(p^e\) 以內做預處理，查詢復雜度可以做到 \(O(\log_p n)\) 左右。
代碼在跟下面一起貼
組合數取模（擴展盧卡斯）
即快速計算\(\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合並起來就可以了。
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(x\in [0,p-2])\) 與整數 \(y\in [1,p-1]\) 一一對應。
我們一般用暴力枚舉的方法求原根，原理類似試除法，判斷是否存在 \(d|p-1\) 使得 \(g^d\equiv 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^k\equiv 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}\)
單位根反演可以應用於 \(\text{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,b\lt 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\times 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 \times 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(N\le 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\)之內。
其實感覺特別暴力。