$\therefore x\equiv (a+i)^{\frac{p+1}{2}} \equiv n^{\frac{1}{2}}\pmod p$
-### 参考实现
-
-```c++
-#include <bits/stdc++.h>
-using namespace std;
-
-typedef long long ll;
-int t;
-ll n, p;
-ll w;
-
-struct num { //建立一个复数域
-
- ll x, y;
-};
-
-num mul(num a, num b, ll p) { //复数乘法
- num ans = {0, 0};
- ans.x = ((a.x * b.x % p + a.y * b.y % p * w % p) % p + p) % p;
- ans.y = ((a.x * b.y % p + a.y * b.x % p) % p + p) % p;
- return ans;
-}
-
-ll binpow_real(ll a, ll b, ll p) { //实部快速幂
- ll ans = 1;
- while (b) {
- if (b & 1) ans = ans * a % p;
- a = a * a % p;
- b >>= 1;
- }
- return ans % p;
-}
-
-ll binpow_imag(num a, ll b, ll p) { //虚部快速幂
- num ans = {1, 0};
- while (b) {
- if (b & 1) ans = mul(ans, a, p);
- a = mul(a, a, p);
- b >>= 1;
- }
- return ans.x % p;
-}
-
-ll cipolla(ll n, ll p) {
- n %= p;
- if (p == 2) return n;
- if (binpow_real(n, (p - 1) / 2, p) == p - 1) return -1;
- ll a;
- while (1) { //生成随机数再检验找到满足非二次剩余的a
- a = rand() % p;
- w = ((a * a % p - n) % p + p) % p;
- if (binpow_real(w, (p - 1) / 2, p) == p - 1) break;
- }
- num x = {a, 1};
- return binpow_imag(x, (p + 1) / 2, p);
-}
-```
+??? note "参考实现"
+ ```c++
+ #include <bits/stdc++.h>
+ using namespace std;
+
+ typedef long long ll;
+ int t;
+ ll n, p;
+ ll w;
+
+ struct num { //建立一个复数域
+ ll x, y;
+ };
+
+ num mul(num a, num b, ll p) { //复数乘法
+ num ans = {0, 0};
+ ans.x = ((a.x * b.x % p + a.y * b.y % p * w % p) % p + p) % p;
+ ans.y = ((a.x * b.y % p + a.y * b.x % p) % p + p) % p;
+ return ans;
+ }
+
+ ll binpow_real(ll a, ll b, ll p) { //实部快速幂
+ ll ans = 1;
+ while (b) {
+ if (b & 1) ans = ans * a % p;
+ a = a * a % p;
+ b >>= 1;
+ }
+ return ans % p;
+ }
+
+ ll binpow_imag(num a, ll b, ll p) { //虚部快速幂
+ num ans = {1, 0};
+ while (b) {
+ if (b & 1) ans = mul(ans, a, p);
+ a = mul(a, a, p);
+ b >>= 1;
+ }
+ return ans.x % p;
+ }
+
+ ll cipolla(ll n, ll p) {
+ n %= p;
+ if (p == 2) return n;
+ if (binpow_real(n, (p - 1) / 2, p) == p - 1) return -1;
+ ll a;
+ while (1) { //生成随机数再检验找到满足非二次剩余的a
+ a = rand() % p;
+ w = ((a * a % p - n) % p + p) % p;
+ if (binpow_real(w, (p - 1) / 2, p) == p - 1) break;
+ }
+ num x = {a, 1};
+ return binpow_imag(x, (p + 1) / 2, p);
+ }
+ ```
## Tonelli-Shanks 算法
-前提条件： $p$ 为质数， $a$ 是 $p$ 的二次剩余
-
-## 证明
-
-- 令 $p-1 = 2^s*t$ , 则 $t$ 为奇数
-- 若 $s = 1$ 则：
- 1. $x^2 = a\ (mod\ p)$ , $x = \sqrt{a}\ (mod\ p)$
- 2. a 是 p 的二次剩余，所以 $a^{\frac{p-1}{2}} = 1 (mod\ p)$ ，所以 $\sqrt{a^{\frac{p-1}{2}}} = 1 (mod\ p)$
- 3. 由 1 和 2，所以 $x = \sqrt{a^{\frac{p-1}{2}}*a}\ (mod\ p)$
- 4. 注意到 $s = 1$ , $x = \sqrt{a^{t}*a}\ (mod\ p)$
- 5. $x = a^{\frac{t+1}{2}}\ (mod\ p)$
-- 因此 $s = 1$ 时，容易求得 $x$
-- 对于 $s > 1$ ，设 $x_{s-1} = a^{\frac{t+1}{2}} (\ mod\ p)$
- 1. $a^{\frac{p-1}{2}} = 1 (mod\ p)$
- 2. $a^{2^{(s-1)}*t} = 1 (mod\ p)$
- 3. $\left(a^{-1} \times\left(a^{\frac{t+1}{2}}\right)^{2}\right)^{2^{s-1}} = 1 (mod\ p)$
- 4. $\left(a^{-1} \times x_{s-1}^{2}\right)^{2^{s-1}} = 1 (mod\ p)$
-- 所以 $a^{-1} \times x_{s-1}^{2}$ 是 1 的 $2^{s-1}$ 次根（模 $p$ 下）
-- 设 $e_k$ 为 $2^k$ 次单位根（模 $p$ 下）
- - $e_{s-1} = a^{-1}*x^{2}_{s-1}$
-- 假设已知 $e_{s-k}, x_{s-k}$
- - $e_{s-k}^{2^{s-k}} = 1 (mod\ p)$
- - $sqrt(e_{s-k}^{2^{s-k}}) = \pm 1 (mod\ p)$
- - $e_{s-k}^{2^{s-k-1}} = \pm 1 (mod\ p)$
-- 计算 $e_{s-k}^{2^{s-k-1}} (mod\ p)$
- - $e_{s-k}^{2^{s-k-1}} = 1 (mod\ p)$
- - $e_{s-k-1} = e_{s-k} (mod\ p)$
- - $x_{s-k-1} = x_{s-k} (mod\ p)$
- - $e_{s-k}^{2^{s-k-1}} = -1 (mod\ p)$
- - $\left(a^{-1} \times x_{s-k}^{2}\right)^{2^{s-k-1}} = -1 (mod\ p)$
- - 我们此时需要找到一个 $q$ ，使得 $\left(a^{-1} \times (qx_{s-k})^{2}\right)^{2^{s-k-1}} = 1 (mod\ p)$
- - 寻找 $q$
- - $\left(a^{-1} \times (qx_{s-k})^{2}\right)^{2^{s-k-1}} = 1 (mod\ p)$
- - $q^{2^{s-k}} = -1 (mod\ p)$
- - 寻找一个非二次剩余 $b$
- - $b^{\frac{p-1}{2}} = -1 (mod\ p)$
- - $b^{t2^{s-1}} = -1 (mod\ p)$
- - $b^{t*2^{s-k}*2^{k-1}} = -1 (mod\ p)$
- - $(b^{t2^{k-1}})^{2^{s-k}} = -1 (mod\ p)$
- - $q = b^{t2^{k-1}}(mod \ p)$
- - $x_{s-k-1} = x_{s-k} * q (mod\ p)$
-- 迭代至最初即可获得答案
-- 复杂度约为 $O(log_2p)$
-
-## 参考实现
-
-```python3
-import random
-# example a, p
-a = 8479994658316772151941616510097127087554541274812435112009425778595495359700244470400642403747058566807127814165396640215844192327900454116257979487432016769329970767046735091249898678088061634796559556704959846424131820416048436501387617211770124292793308079214153179977624440438616958575058361193975686620046439877308339989295604537867493683872778843921771307305602776398786978353866231661453376056771972069776398999013769588936194859344941268223184197231368887060609212875507518936172060702209557124430477137421847130682601666968691651447236917018634902407704797328509461854842432015009878011354022108661461024768
-p = 30531851861994333252675935111487950694414332763909083514133769861350960895076504687261369815735742549428789138300843082086550059082835141454526618160634109969195486322015775943030060449557090064811940139431735209185996454739163555910726493597222646855506445602953689527405362207926990442391705014604777038685880527537489845359101552442292804398472642356609304810680731556542002301547846635101455995732584071355903010856718680732337369128498655255277003643669031694516851390505923416710601212618443109844041514942401969629158975457079026906304328749039997262960301209158175920051890620947063936347307238412281568760161
-
-
-b = random.randint(1, p)
-while (pow(b, (p-1)//2, p) == 1):
+大致思路如下：
+
+先令 $p-1 = 2^s \times t$ , 则 $t$ 为奇数。
+
+针对 $s$ 的值分两种情况 $s = 1$ 和 $s > 1$ 进行讨论。
+
+如果 $s = 1$：
+
+1. 因为 $a$ 是 $p$ 的二次剩余，所以 $a^{\frac{p-1}{2}} \equiv 1 \pmod p$ ，即 $\sqrt{a^{\frac{p-1}{2}}} \equiv 1 \pmod p$；
+2. 由 1 可知，$x \equiv \sqrt{a^{\frac{p-1}{2}} \times a} \pmod p$；
+3. 注意到 $s = 1$，可以得出 $x \equiv \sqrt{a^{t} \times a} \pmod p$，进一步得到 $x \equiv a^{\frac{t+1}{2}} \pmod p$。
+
+因此 $s = 1$ 时，$x$ 的值容易求出。
+
+对于 $s > 1$ 的情况，设 $x_{s-1} \equiv a^{\frac{t+1}{2}} \pmod p$。
+
+1. 由欧拉判别准则可知 $a^{\frac{p-1}{2}} \equiv 1 \pmod p$，进一步有 $a^{2^{(s-1)} \times t} \equiv 1 \pmod p$；
+2. 稍作变形得到 $\left(a^{-1} \times \left(a^{\frac{t+1}{2}}\right)^{2}\right)^{2^{s-1}} \equiv 1 \pmod p)$，即 $\left(a^{-1} \times x_{s-1}^{2}\right)^{2^{s-1}} \equiv 1 \pmod p)$。
+
+所以 $a^{-1} \times x_{s-1}^{2}$ 是模 $p$ 意义下的 1 的 $2^{s-1}$ 次根。
+
+接下来设 $e_k$ 为模 $p$ 意义下的 $2^k$ 次单位根。容易发现 $e_{s-1} = a^{-1} \times x^{2}_{s-1}$。
+
+假设我们已经知道 $e_{s-k}, x_{s-k}$：因为 $e_{s-k}^{2^{s-k}} \equiv 1 \pmod p$，所以有 $\sqrt(e_{s-k}^{2^{s-k}}) \equiv \pm 1 \pmod p$，即 $e_{s-k}^{2^{s-k-1}} \equiv \pm 1 \pmod p$；
+
+现在的任务变成了计算 $e_{s-k}^{2^{s-k-1}} \pmod p$。这时候又可以分为两种情况：$e_{s-k}^{2^{s-k-1}} \equiv 1 \pmod p$ 和 $e_{s-k}^{2^{s-k-1}} \equiv -1 \pmod p$。
+
+对于第一种情况，容易得到：$e_{s-k-1} \equiv e_{s-k} \pmod p$，$x_{s-k-1} \equiv x_{s-k} \pmod p$。
+
+对于第二种情况，则有：$\left(a^{-1} \times x_{s-k}^{2}\right)^{2^{s-k-1}} \equiv -1 \pmod p$。我们此时需要找到一个 $q$ ，使得 $\left(a^{-1} \times (qx_{s-k})^{2}\right)^{2^{s-k-1}} \equiv 1 \pmod p$。
+
+寻找 $q$ 的思路如下：因为 $\left(a^{-1} \times (qx_{s-k})^{2}\right)^{2^{s-k-1}} \equiv 1 \pmod p$，所以 $q^{2^{s-k}} \equiv -1 \pmod p$。
+
+接下来寻找一个非二次剩余 $b$：
+
+- 因为 $b$ 是非二次剩余，所以 $b^{\frac{p-1}{2}} \equiv -1 \pmod p$，即 $b^{t2^{s-1}} \equiv -1 \pmod\ p$；
+- 稍微变个形得到 $b^{t \times 2^{s-k} \times 2^{k-1}} \equiv -1 \pmod p$，即 $(b^{t \times 2^{k-1}})^{2^{s-k}} \equiv -1 \pmod p$；
+- $q \equiv b^{t \times 2^{k-1}} \pmod p$。
+
+最后得到 $x_{s-k-1} \equiv x_{s-k} \times q \pmod p$。
+
+不断迭代即可得到答案。该做法的时间复杂度约为 $O(\log p)$。
+
+??? note "参考实现"
+ ```python3
+ import random
+ # example a, p
+ a = 8479994658316772151941616510097127087554541274812435112009425778595495359700244470400642403747058566807127814165396640215844192327900454116257979487432016769329970767046735091249898678088061634796559556704959846424131820416048436501387617211770124292793308079214153179977624440438616958575058361193975686620046439877308339989295604537867493683872778843921771307305602776398786978353866231661453376056771972069776398999013769588936194859344941268223184197231368887060609212875507518936172060702209557124430477137421847130682601666968691651447236917018634902407704797328509461854842432015009878011354022108661461024768
+ p = 30531851861994333252675935111487950694414332763909083514133769861350960895076504687261369815735742549428789138300843082086550059082835141454526618160634109969195486322015775943030060449557090064811940139431735209185996454739163555910726493597222646855506445602953689527405362207926990442391705014604777038685880527537489845359101552442292804398472642356609304810680731556542002301547846635101455995732584071355903010856718680732337369128498655255277003643669031694516851390505923416710601212618443109844041514942401969629158975457079026906304328749039997262960301209158175920051890620947063936347307238412281568760161
+
b = random.randint(1, p)
+ while (pow(b, (p-1)//2, p) == 1):
+ b = random.randint(1, p)
-t = p - 1
-s = 0
-while (t % 2 == 0):
- s += 1
- t = t // 2
-x = pow(a, (t + 1) // 2, p)
-e = pow(a, t, p)
-k = 1
-while (k < s):
- if (pow(e, 1 << (s - k - 1), p) != 1):
- x = x * pow(b, (1 << (k - 1)) * t, p) % p
- e = pow(a, p-2, p) * x % p * x % p
- k += 1
-print(x)
-```
+ t = p - 1
+ s = 0
+ while (t % 2 == 0):
+ s += 1
+ t = t // 2
+ x = pow(a, (t + 1) // 2, p)
+ e = pow(a, t, p)
+ k = 1
+ while (k < s):
+ if (pow(e, 1 << (s - k - 1), p) != 1):
+ x = x * pow(b, (1 << (k - 1)) * t, p) % p
+ e = pow(a, p-2, p) * x % p * x % p
+ k += 1
+ print(x)
+ ```
## 习题
OSDN Git Service
