From 0a0181c85ddba8231ee6415a477e66a434dcefa6 Mon Sep 17 00:00:00 2001 From: Shuhao Zhang Date: Wed, 30 Dec 2020 22:27:29 +0800 Subject: [PATCH] fix(quad-residue): fix format --- docs/math/quad-residue.md | 249 +++++++++++++++++++++++----------------------- 1 file changed, 122 insertions(+), 127 deletions(-) diff --git a/docs/math/quad-residue.md b/docs/math/quad-residue.md index 34acec1b..1fccc827 100644 --- a/docs/math/quad-residue.md +++ b/docs/math/quad-residue.md @@ -128,138 +128,133 @@ $$ $\therefore x\equiv (a+i)^{\frac{p+1}{2}} \equiv n^{\frac{1}{2}}\pmod p$ -### 参考实现 - -```c++ -#include -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 + 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) + ``` ## 习题 -- 2.11.0