同餘系基本知識
炒冷飯
求逆元
exgcd(a,m,x,y); ans=(x%m+m)%m;，萬能。qpow(a,m-2,m)，這個僅限於m是一個質數。inv[i]=(m-m/i)*inv[m%i]%m，線性求逆元。
中國剩餘定理(CRT)
現在已知一個同餘方程組:
\[\begin{cases} x\equiv a_1\mod m_1\\ x\equiv a_2\mod m_2\\ ...\\ x\equiv a_n\mod m_n \end{cases} \]
滿足\(m_i\)兩兩互質。
解法：
首先計算出下面的各個值：
\[M=\prod\limits_{i=1}^m m_i\\ P_i={M\over m_i}\\ Q_i=P_i^{-1}\mod m_i \]
那麼最終的答案就是：
\[ans\equiv\sum\limits_{i=1}^n a_iP_iQ_i \mod M \]
我不會證明
擴充套件中國剩餘定理(EXCRT)
現在解決的是\(m_i\)可能不兩兩互質的情況。
符號約定：
\(\mathrm{inv}(a,m)\)表示\(a\)在模\(m\)下的逆元。
\((a,b)=\gcd(a,b)\)
考慮將兩條方程進行合併。
現在有兩條方程：
\[\begin{cases} x\equiv a_1\mod m_1\\ x\equiv a_2\mod m_2 \end{cases} \]
轉化為不定方程的形式：
\[\begin{cases} x=a_1+m_1 y_1\\ x=a_2+m_2 y_2 \end{cases} \Rightarrow \begin{cases} y_1=\frac{x-a_1}{m_1}\\ y_1=\frac{x-a_2}{m_2} \end{cases} \]
那麼可以得到：
\[\begin{aligned} a_1+m_1 y_1&=a_2+m_2 y_2\\ m_1 y_1-m_2 y_2&=a_2-a_1\\ \end{aligned}\\ \]
可以發現，這個方程有解的條件是\((m_1, m_2)|(a_2-a_1)\)，方程兩邊同時除以\((m_1, m_2)\),然後化為同餘式：
\[\frac{m_1}{(m_1,m_2)}y_1\equiv\frac{a_2-a_1}{(m_1, m_2)}\mod{\frac{m_2}{(m_1,m_2)}} \]
可以發現此時\((\frac{m_1}{(m_1, m_2)}, \frac{m_2}{(m_1,m_2)})=1\)
\[\Rightarrow y_1\equiv \frac{a_2-a_1}{(m_1, m_2)}\times \mathrm{inv}(\frac{m_1}{(m_1, m_2)},\frac{m_2}{(m_1,m_2)}) \mod \frac{m_2}{(m_1,m_2)} \]
接著將其再次轉化為一個不定方程，然後帶入\(x\)：
\[\Rightarrow y_1= \frac{a_2-a_1}{(m_1, m_2)} \times \mathrm{inv}(\frac{m_1} {(m_1,m_2)},\frac{m_2}{(m_1,m_2)}) \mod \frac{m_2}{(m_1,m_2)} + k \times \frac{m_2}{(m_1,m_2)}\\ \Rightarrow \frac{x-a_1}{m_1}= \frac{a_2-a_1}{(m_1, m_2)} \times \mathrm{inv}(\frac{m_1}{(m_1, m_2)},\frac{m_2}{(m_1,m_2)}) \mod \frac{m_2}{(m_1,m_2)} + k \times \frac{m_2}{(m_1,m_2)}\\ \Rightarrow x= \frac{a_2-a_1} {(m_1, m_2)} \times \mathrm{inv}(\frac{m_1}{(m_1, m_2)},\frac{m_2}{(m_1,m_2)}) \mod \frac{m_2}{(m_1,m_2)} \times m_1 + a_1 + k \times \frac{m_1 m_2}{(m_1,m_2)} \]
然後轉化為同餘式：
\[\Rightarrow x\equiv \frac{a_2-a_1} {(m_1, m_2)} \times \mathrm{inv}(\frac{m_1}{(m_1, m_2)},\frac{m_2}{(m_1,m_2)}) \mod \frac{m_2}{(m_1,m_2)} \times m_1 + a_1 ~~~(\mathrm{mod} \frac{m_1 m_2}{(m_1,m_2)}) \]
這樣就合併成功了。
同餘式的轉化
在實際問題中，可能給出的同餘式是類似下面這個樣子的：
\[ax\equiv b\mod M \]
我們需要將其轉化為(EX)CRT的形式的形式。
將其轉化為不定方程的形式：
\[ax+My=b \]
可以發現，這個方程有解的條件是\((a,M)|b\)。
首先使用exgcd求出方程\(ax+My=(a,M)\)的特解，得到的一組解為\(x',y'\)。
那麼原方程的一組解為
\[Sx=\frac{b}{(a,M)}x'\\ Sy=\frac{b}{(a,M)}y' \]
那麼原方程的通解就是
\[x=Sx+k\times \frac{b}{(a,M)}\\ y=Sy-k\times \frac{b}{(a,M)} \]
我們取x的通解表示式，將其轉化為同餘式：
\[x\equiv Sx~~~~(\mathrm{mod}~~\frac{b}{(a,M)}) \]
新內容
基本概念
階:當\((a,p)=1\)時，滿足\(a^x\equiv 1~~~~(\mathrm{mod}~~p)\)的最小\(x\)，計作\(\delta_p(a)\)。
原根:如果\(\delta_p(a)=\varphi(p)\)，那麼\(a\)就是模\(p\)下的原根。一個模數由原根，就必須先滿足\(p=2,4,p^n,2p^n\)。
如果\(p\)是一個質數，那麼可以用\(a^x\%p\)表示\([1,p-1]\)中的所有數，其中\(1\leq x\leq p-1\)
二次剩餘
形式化的，就是找到滿足下面的同餘式的\(x\):
\[x^2\equiv a~~~~(\mathrm{mod}~~p) \]
通俗點來講，就是將\(a\)在模\(p\)以一下進行開根。
如果\(a\)能夠開根（上面的同餘方程有解），那麼我們就可以認為“\(a\)是模\(p\)意義下的二次剩餘”，否則認為“\(a\)不是模\(p\)意義下的二次剩餘”。
上面同餘方程要麼有兩個解（可能相同），要麼沒有解。
如果\(p=2\)，那麼顯然的，x的取值只有0/1，而且很容易找到，所以我們只討論\(p\)屬於奇素數的情況。
引理1:對於一個模數\(p\)，是二次剩餘的數共有\(\frac{p-1}{2}\)（不包括\(0\)）,非二次剩餘的數也有\(\frac{p-1}{2}\)個，各佔一半。
雖然我們很難在較短時間內準確知道是二次剩餘的數的分佈，但是沒有關係。
引理2:對於一個數\(a\)和一個模數\(p\),可以得到：
\[v\equiv a^{\frac{p-1}{2}}~~~~(\mathrm{mod}~~p)\\ \begin{cases} v=1&\Leftrightarrow a是模p的二次剩餘\\ v=-1&\Leftrightarrow a不是模p的二次剩餘 \end{cases} \]
接下來考慮演算法（忽略各種特判）。
首先可以想到，第一步一定是判斷一個數\(a\)是否是二次剩餘。
如果是一個二次剩餘，那麼繼續執行演算法，找到\(x\)，否則直接返回一個表示不是二次剩餘的結果。
接下來就是一系列神仙操作。
隨機一個數\(a\)，範圍是\([1,p-1]\)，滿足\(a^2-n\)不是模\(p\)的二次剩餘。
接下來使用類似虛數這種拓展數域的方式，令\(\omega=\sqrt{a^2-n}\)，這個東西在模\(p\)的數域下是不存在的。用二元組\((a,b)\)表示\(a+b\omega\)，那麼可以得到：
\[\begin{aligned} (a+\omega)^{p-1}&=(a+\omega)^p(a+\omega)\\ &=(\sum\limits_{i=0}^p C(p,i)a^i\omega^{p-i})(a+\omega)\\ &=(\sum\limits_{i=1}^p C(\lfloor{p/p}\rfloor,\lfloor{i/p}\rfloor) C(p\%p,i\%p)a^i\omega^{p-i})(a+\omega)\\ \end{aligned} \]
很明顯，\(C(\lfloor{p/p}\rfloor,\lfloor{i/p}\rfloor)\)的值永遠為\(1\),而\(C(p\%p,i\%p)\)只有在\(i=0或p\)的時候才會是\(1\)，其他時刻都為\(0\)，那麼式子可以繼續化簡：
\[\begin{aligned} &=(a^p+\omega^p)(a+\omega)\\ &=(a^{p-1}a+\omega^{p-1}\omega)(a+\omega)\\ &=(a^{p-1}a+(a^2+n)^{\frac{p-1}{2}}\omega)(a+\omega) \end{aligned} \]
因為\(a^{p-1}\equiv 1\), \((a^2+n)^{\frac{p-1}{2}}\equiv -1\)，所以
\[\begin{aligned} &=(a-\omega)(a+\omega)\\ &=a^2-\omega^2\\ &=n \end{aligned} \]
所以\((a+\omega)^{\frac{p+1}{2}}\)就是其中一個\(x\)，用\(p-x\)就能得到另一個\(x\)。
struct Complex {
Complex(LL x, LL y) : x(x), y(y) {}
LL x, y;
};
inline Complex mul(Complex a, Complex b) {
return Complex(a.x*b.x%M + a.y*b.y%M*w%M, a.y*b.x*M+a.x*b.y%M);
}
LL qpow(LL a, LL b, LL M) {
LL ret = 1;
while (b) {
if (b & 1) ret = ret * a % M;
a = a * a % M, b >>= 1;
}
return ret;
}
Complex qpow(Complex a, LL b) {
Complex ret = Complex(1, 0);
while (b) {
if (b & 1) ret = mul(ret, a);
a = mul(a, a), b >>= 1;
}
return ret;
}
pair<LL, LL> solve(LL n, LL M) {
n %= M;
if (!n) return make_pair(1, 1);
if (M == 2) return make_pair(1, 1);
if (qpow(n, (M - 1) >> 1, M) == M - 1) return make_pair(-1,-1);
LL ans, ans1, ans2;
a = rand() % M;
while (!a || qpow((a * a - n + M) % M, (M - 1) >> 1, M)) a = rand() % M;
w = (a * a - n + M) % M;
ans1 = qpow(Complex(a, 1), (M + 1) >> 1).x;
ans2 = M - ans1;
if (ans1 > ans2) swap(ans1, ans2);
return make_pair(ans1, ans2);
}
大步小步(BSGS)
就是在模\(p\)的情況下求出\(log_ab\)，非常暴力。
現在考慮的是\(a,p\)互質的特殊情況，
首先使用Hash表或者std::map儲存\(b\times a^{i}(0\leq i\leq \lceil\sqrt{b}\rceil)\)，以及對應的\(i\)，然後列舉\(a^{i\lceil{\sqrt{p}}\rceil}(0\leq i\lceil{\sqrt{p}}\rceil\leq p)\)，如果得到一個值存在於表中，那麼直接將指數相減就好了。
namespace BSGS {
map<LL, LL> val_map;
LL log(LL a, LL b, LL M) { // a^x = b (mod M)
val_map.clear();
LL sqrt_m = ceil(sqrt(M));
b %= M;
for (int i = 1; i <= sqrt_m; i++) {
b = b * a % M;
val_map[b] = i;
}
LL vec = qpow(a, sqrt_m, M); b = 1;
for (int i = 1; i <= sqrt_m; i++) {
b = b * vec % M;
if (val_map.count(b)) return (i * sqrt_m - val_map[b] + M) % M;
}
return -1;
}
}
拓展大步小步(EX_BSGS)
現在考慮\(a,p\)不互質的情況。
\[設g=(a,p)\\ a^x\equiv b~~~~(\mathrm{mod}~~p)\\ (\frac{a}{g})a^{x-1}\equiv \frac{b}{g}~~~~(\mathrm{mod}~~\frac{p}{g}) \]
如果\(b\)不能被\(g\)整除，那麼式子無解，因為\(\frac{a}{g},\frac{p}{g}\)互質，所以一定有\(\frac{a}{g}\)在模\(\frac{p}{g}\)下的逆元，對式子進行變化：
\[a^{x-1}\equiv \frac{b}{g}(\frac{a}{g})^{-1}~~~~(\mathrm{mod}~~\frac{p}{g}) \]
這樣子就得到另一個形式相同的式子，但是\(a\)和\(\frac{p}{g}\)可能依然不互質，需要繼續變化。如果最終互質了，就用BSGS進行求解。
namespace EX_BSGS {
LL log(LL a, LL b, LL M) {
if (b == 1 || M == 1) return 0;
LL g = gcd(a, M), cnt = 0, k = 1;
while (g > 1) {
if (b % g) return -1;
cnt++, b /= g, M /= g, k = k * (a / g) % M;
if (k == b) return cnt;
g = gcd(a, M);
}
LL ans2 = BSGS::log(a, b * inv(k, M) % M, M);
if (ans2 == -1) return -1;
else return ans2 + cnt;
}
}