對於同餘式
\[x^2 \equiv n \pmod p\]
若對於給定的\(n, P\)，存在\(x\)滿足上面的式子，則乘\(n\)在模\(p\)意義下是二次剩餘，否則為非二次剩餘
我們需要計算的是在給定範圍內所有滿足條件的\(x\)，同時為了方便，我們只討論\(p\)是奇質數的情況
前置定理
- \(x^2 \equiv (x+p)^2 \pmod p\)
證明：\(x^2 \equiv x^2 + 2xp + p^2 \pmod p\)顯然成立
- 對於\(x^2 \equiv n \pmod p\)，除\(n=0\)外，總共有\(\frac{p-1}{2}\)個\(n\)使得該方程有解
我局的參考資料裡對於這條性質的證明漏洞很大，所以下面的是自己yy的
根據第一個前置定理的式子，我們只需討論\(x \in [1, p - 1]\)即可(當\(x=0\)時對應了\(n=0\)的特殊情況)
一個顯然的性質是
\[x^2 \equiv (p - x)^2 \pmod p\]
那麼當\(x \in [1, \frac{p - 1}{2}]\)我們可以取到所有解。
接下來我們只需要證明當\(x\in[1, \frac{p-1}{2}]\)時\(x^2 \bmod p\)均兩兩不同
可以用反證法，若存在不同的\(u, v\)滿足\(u^2 \equiv v^2 \pmod p\)
那麼有\((u + v)(u - v) \equiv 0 \pmod p\)
顯然\(-p < u + v < p\)且\(-p < u - v < p\)且\(u + v \not = 0, u - v \not = 0\)，故該假設不成立，故原命題成立。
Q.E.D
- 勒讓德符號(Legender symbol)
記
\[ (\frac{a}{p}) = \begin{cases} 1 , &\text{a在模$p$意義下是二次剩餘}\\ -1, &\text{a在模$p$意義下是非二次剩餘}\\ 0, &\text{a mod p = 0} \end{cases} \]
這個東西的分佈大概是這個樣子
計算公式
我局的這個公式就是構造出來的
\[ (\frac{a}{p}) = a^{\frac{p - 1}{2}} \pmod p \]
證明：
費馬小定理：對於任意互質的\(x, p\)，有\(x^{p - 1} = 1 \pmod p\)
一條同餘式的性質：若\(a^k \equiv b^k \pmod p\)，那麼\(a^{kx} \equiv b^{xk} \pmod p\)
然後直接把這玩意兒帶到\(x^2 \equiv a \pmod p\)裡就行了
這裡簡單的寫一下:
首先要明確我們的目的，我們現在要驗證這個公式的正確性，也就是說我們要證明當\(a^{\frac{p-1}{2}}=1 \pmod p\)時滿足條件的\(x\)存在，當\(a^{\frac{p-1}{2}}= -1 \pmod p\)時\(x\)不存在，當\(a^{\frac{p-1}{2}}= 0 \pmod p\)時\(a\mod p = 0\)
- 當\(a^{\frac{p-1}{2}}=1 \pmod p\)時
我們假設有\(x^2 \equiv a \pmod p\)
\[x^{2\frac{p-1}{2}} \equiv a^{\frac{p-1}{2}} \pmod p\]
\[x^{p-1} \equiv 1 \pmod p\]
根據費馬小定理\(x\)顯然存在，因此\(a\)是模\(p\)意義下的二次剩餘
- 當\(a^{\frac{p-1}{2}}= -1 \pmod p\)時
假設有\(x^2 \equiv a \pmod p\)
同理可知
\[x^{p-1} \equiv -1 \pmod p\]
顯然\(x\)不存在，因此\(a\)不是模\(p\)意義下的二次剩餘
- 當\(a^{\frac{p-1}{2}}= 0 \pmod p\)時
顯然有\(a \bmod p = 0\)
Cipolla演算法
演算法流程
這個演算法其實用兩句話就能說完，但是背後的理論卻非常高深(對於我這種菜雞而言)。
首先使用隨機的方法找到一個\((\frac{a^2 - n}{p}) = -1\)，記\(\omega = \sqrt{a^2-n}\)
那麼\(x \equiv (a + w)^{\frac{p+1}{2}} \pmod p\)
做完了。。。期望複雜度\(O(\log^2 n)\)
但是實際上實現起來並沒有這麼簡單，因為要自定義類似於虛數的乘法/冪運算
演算法理論
首先要有一點抽代基礎(群/環/域什麼的要知道定義)
我們來逐步分析這個演算法(按照我的敘述風格應該是從發明者的角度出發一步一步推出這玩意兒來，但是十分抱歉我實在是搞不明白他當時的腦回路qwq)
對於第一步，根據前面的定理，如果在\([1, p]\)內隨機，每次有\(\frac{1}{p}*\frac{p-1}{2}\)的概率找到一個解，那麼期望步數大約為兩次，因此複雜度是可以保證的。
但是找到這個東西有什麼用呢？。如果我們把之前的數域記做\(\mathbf F_p\)，\(\omega\)在這個數域下是不能開根的，但是我們可以構造一個新的數域\(\mathbf F_p\)，使得\(\omega\)在\(\mathbf F_{p2}\)下能夠開根。類比於\(-1\)在複數域下能夠表示為\(\sqrt{-1}\)一樣。
這樣的話\(\mathbf F_{p2}\)內的數都可以寫作\(a + k\omega\)的形式。可以證明這玩意兒確實是個合法的域，證明過程，同時也可以證明在\(\mathbf F_{p2}\)下得到的解在\(\mathbf F_{p1}\)下也成立，同時最後的答案中\(\omega\)的係數一定為\(0\)
現在來簡單說明一下為什麼\(x \equiv (a+\omega)^{\frac{p+1}{2}}\)
先來了解兩個性質
- \(\omega^p \equiv -\omega \pmod p\)
證明：
\[
\begin{aligned}
\omega^p &= (a^2-n)^{\frac{p}{2}}\\
&= (a^2 - n)^{\frac{p - 1}{2}} (a^2 - n)^{\frac{1}{2}}\\
&= -\omega
\end{aligned}
\]
- \((a + b)^p \equiv a^p + b^p \pmod p\)
證明就直接考慮二項式定理中的組合數展開，發現除了第一項和最後一項之外都無法把\(n!\)消掉。
那麼要證明\(x \equiv (a+\omega)^{\frac{p+1}{2}}\)，實際上我們只需要證明\((a+\omega)^{p+1}\equiv n \pmod p\)就行了
\[ \begin{aligned} &(a + \omega)^{p + 1}\\ =&(a + \omega)^p(a + \omega)\\ =&(a - \omega)(a + \omega)(\text{根據費馬小定理$a^p \equiv p \pmod p$})\\ =&(a^2 - \omega^2)\\ =&(a^2 - (a^2 - n))\\ =&n \end{aligned} \]
演算法的大概思想就講完了，下面煮個栗子~。
對於\(x^2 \equiv n \pmod p\)
假設此時\(p=13, n = 10\)。
首先要找到一個\(a\)滿足\((\frac{a^ - 10}{13}) = -1\)，然後臉黑的attack在經過1e9 +7次嘗試後終於找到了一個\(a =2\)它滿足條件，因為\((\frac{7}{13}) = -1\)此時\(\omega = \sqrt{a^2 - n} = \sqrt{-6}\)
按照老祖宗講給我們的
\(x \equiv (2 + \sqrt{-6})^{7} \pmod {13}\)
\[ \begin{aligned} &\left(2+{\sqrt {-6}}\right)^{2}=4+4{\sqrt {-6}}-6=-2+4{\sqrt {-6}}\\ &\left(2+{\sqrt {-6}}\right)^{4}=\left(-2+4{\sqrt {-6}}\right)^{2}=-1-3{\sqrt {-6}}\\ &\left(2+{\sqrt {-6}}\right)^{6}=\left(-2+4{\sqrt {-6}}\right)\left(-1-3{\sqrt {-6}}\right)=9+2{\sqrt {-6}}\\ &\left(2+{\sqrt {-6}}\right)^{7}=\left(9+2{\sqrt {-6}}\right)\left(2+{\sqrt {-6}}\right)=6. \end{aligned} \]
然後不難發現\(36 \equiv 10 \pmod {13}\)
同時因為平方的性質，\(-x\)也是一個合法解，因此\(-6 + 13 = 7\)也是合法的
最後有一個小問題就是為什麼最後\(\omega\)的係數一定是\(0\)，參考資料中給出的解釋我實在是不能理解，如果有看得懂的大佬歡迎給本菜雞講一下qwq
程式碼模板
#include<bits/stdc++.h>
using namespace std;
const int mod = 13;
namespace TwoRemain {
template <typename A, typename B> inline int add(A x, B y) {if(x + y < 0) return x + y + mod; return x + y >= mod ? x + y - mod : x + y;}
template <typename A, typename B> inline void add2(A &x, B y) {if(x + y < 0) x = x + y + mod; else x = (x + y >= mod ? x + y - mod : x + y);}
template <typename A, typename B> inline int mul(A x, B y) {return 1ll * x * y % mod;}
template <typename A, typename B> inline void mul2(A &x, B y) {x = (1ll * x * y % mod + mod) % mod;}
int fmul(int a, int p, int Mod = mod) {
int base = 0;
while(p) {
if(p & 1) base = (base + a) % Mod;
a = (a + a) % Mod; p >>= 1;
}
return base;
}
int fp(int a, int p, int Mod = mod) {
int base = 1;
while(p) {
if(p & 1) base = fmul(base, a, Mod);
p >>= 1; a = fmul(a, a, Mod);
}
return base;
}
int f(int x) {
return fp(x, (mod - 1) >> 1);
}
struct MyComplex {
int a, b;
int cn;
MyComplex operator * (const MyComplex &rhs) {
return {
add(fmul(a, rhs.a), fmul(cn, fmul(b, rhs.b, mod))),
add(fmul(a, rhs.b), fmul(b, rhs.a)),
cn
};
}
};
MyComplex fp(MyComplex a, int p) {
MyComplex base = {1, 0, a.cn};
while(p) {
if(p & 1) base = base * a;
a = a * a; p >>= 1;
}
return base;
}
int TwoSqrt(int n) {
if(f(n) == mod - 1) return -1;
if(f(n) == 0) return 0;
int a = -1, val = -1;
while(val == -1) {
a = rand() << 15 | rand();
val = add(mul(a, a), -n);
if(f(val) != mod - 1) val = -1;
}
return fp({a, 1, val}, (mod + 1) / 2).a;
}
}
using namespace TwoRemain;
signed main() {
cout << TwoSqrt(10);
return 0;
}