「学习笔记」二次剩余

更好的阅读体验

问题引入

对于方程 (x^2 equiv n pmod P)，求解 (x) 的值。

本文只讨论 (P) 为奇素数的情况。

二次剩余 / 非二次剩余

对于方程 (x^2 equiv n pmod P)，如果有解则称 (n) 是模 (P) 的二次剩余，反之则称 (n) 是模 (P) 的非二次剩余。

定理 1

模 (P) 意义下的二次剩余的数量为 (frac{P-1}{2})（不包括 (0)），非二次剩余的数量为 (frac{P-1}{2}) 个。

证明

由于 ((x - P)^2 = x^2 - 2xP + P^2 = x^2 - P(2x - P) equiv x^2 pmod P)，那么当 (x in [1,frac{P-1}{2}]) 我们可以取到所有解。

因此我们只需要证明当 (x in [1,frac{P-1}{2}]) 时 (x^2 mod P) 两两不同。

运用反证法，设存在两个数 (x,y (x eq y)) 满足 (x,y in [1,frac{P-1}{2}]) 且 (x^2 equiv y^2 pmod P)，即 ((x + y)(x - y) equiv 0 pmod P)。

因为 (x eq y)、(x,y in [1,frac{P-1}{2}])，则 (P mid x - y) 且 (x - y eq 0)，因此 (Pmid (x + y))。

由于 (x,y in [1,frac{P-1}{2}]) ，所以当且仅当 (x = y = frac{P-1}{2}) 时，(Pmid (x + y))，与 (x eq y) 不符，故假设不成立。

所以模 (P) 意义下的二次剩余的数量为 (frac{P-1}{2})（不包括 (0)），非二次剩余的数量为 (frac{P-1}{2}) 个。

阶

由欧拉定理可知若 (a⊥P)，则 (a^{varphi(P)} equiv 1 pmod P)。

那么我们称满足 (a^n equiv 1 pmod P) 的最小正整数 (n) 为 (a) 模 (P) 的阶，记作 (delta_P(a))。

这也被称为阶的最小性。

定理 2

若 (a^n equiv 1 pmod P)，则 (delta_P(a)mid n)。

证明

运用反证法，设 (n=delta_P(a)k+b, 0 < b < delta_P(a))，则

[a^nequiv a^{delta_P(a)k+b} = a^{delta_P(a)k} cdot a^b equiv a^b equiv 1 pmod P ]

这与阶的最小性矛盾，故假设不成立，即 (b = 0)，故 (delta_P(a)mid n)。

原根

若 (a⊥P)，且 (delta_P(a) = varphi(P))，则称 (a) 为模 (P) 的原根。

定理 3

设 (a) 为模素数 (P) 的原根，则 (a, a^2, a^3,..., a^{varphi(P)}) 各数对 (P) 取模可以得到一个 ([1, P-1]) 的排列。

证明

我们先证明 (a, a^2, a^3,..., a^{varphi(P)}) 各数对 (P) 取模互不相等。

运用反证法，假设存在两个数 (i eq j (1 le i,j le varphi(P)))，且 (a^i equiv a^j pmod P)，则 (a^{mid i-jmid} equiv 1 pmod P)。

因为 (a) 为模 (P) 的原根，所以 (delta_P(a) = varphi(P))，则 (1 le i,j le delta_P(a))。

所以 (mid i-jmid < delta_P(a))，但这与阶的最小性矛盾，故假设不成立，则 (a, a^2, a^3,..., a^{varphi(P)}) 各数对 (P) 取模互不相等。

又因为 (P) 为素数，则 (delta_P(a) = varphi(P) = P - 1)，且 (a, a^2, a^3,..., a^{varphi(P)}) 这 (P-1) 个数对 (P) 取模的结果都在 ([1,P-1]) 中，因而 (a, a^2, a^3,..., a^{varphi(P)}) 这 (P-1) 个数对 (P) 取模后就是一个 ([1, P-1]) 的排列。

勒让德符号

[left(frac{n}{p} ight)=egin{cases} 1,\,&p mid n ext{ 且 }n ext{ 是 }p ext{ 的二次剩余}\ -1,\,&p mid n ext{ 且 }n ext{ 不是 }p ext{ 的二次剩余}\ 0,\,&pmid n end{cases} ]

欧拉判别准则

[left(frac{n}{P} ight)equiv n^{frac{P-1}{2}}pmod P ]

若 (n) 是二次剩余，当且仅当 $n^{frac{P-1}{2}}equiv 1pmod P $。

若 (n) 是非二次剩余，当且仅当 $n^{frac{P-1}{2}}equiv -1pmod P $。

证明

由费马小定理得 (n^{P-1} equiv 1 pmod P)，则 ((n^{frac{P-1}{2}}+1)(n^{frac{P-1}{2}}-1) equiv 0 pmod P)，所以 (n^{frac{P-1}{2}}) 在模 (P) 意义下只能为 (1) 或 (-1)。

先证明 $n^{frac{P-1}{2}}equiv 1pmod P $ 的情况，由于只有两种情况故一种情况对应着 (n) 是二次剩余则另一种情况必定对应着 (n) 是非二次剩余。

设 (f) 为模 (P) 下的原根，由定理 3 得 (n) 一定可以表示为 (f^k (1 le k le P - 1)) 的形式，于是可得

[f^{kcdotfrac{P-1}{2}} equiv 1 ]

由定理 2 可得 (delta_P(a)mid (kcdotfrac{P-1}{2}))，由于 (P) 为素数，则 ((P - 1)mid kcdotfrac{P-1}{2})，由此可得 (frac{k}{2}) 为整数。

令 (i = frac{k}{2})，则有 ((a^i)^2 = a^{2i} equiv n pmod P)，则 (n) 是模 (P) 的二次剩余。

Cipolla 算法

对于任意的数 (a)，令 (i^2 = a^2 - n)，若 (i^2) 为非二次剩余，则方程 (x^2 equiv n pmod P) 的解为 (x_0 = (a + i)^{frac{P + 1}{2}}, x_1 = P - (a + i)^{frac{P + 1}{2}})。

算法实现时随机一个 (a) 判断 (i^2) 是否为非二次剩余，若是则可以解出 (x)。

证明

引理 1：((a+b)^Pequiv a^P + b^P pmod P)
证明：

[egin{aligned}(a+b)^P &equiv sum_{i=0}^{P} ext{C}_p^i a^{P-i}b^i \ &equiv sum_{i=0}^{P}frac{P!}{(P-i)!i!}a^{P-i}b^i \ &equiv a^P+b^Ppmod Pend{aligned} ]

发现式子在 (i=0) 和 (i=P) 时没有含 (P)，因而没有被模去，留下 (a^P + b^P)。

引理 2：(i^P equiv -i pmod P)
证明：

[egin{aligned} i^P &equiv i^{P-1} cdot i \ &equiv (i^2)^{frac{P-1}{2}}cdot i \ &equiv (a^2-n)^{frac{P-1}{2}}cdot i \ &equiv -i pmod P end{aligned} ]

由上面两条引理可以开始推导：

[egin{aligned} x_0 &equiv (a+i)^{frac{P+1}{2}} \ &equiv ((a+i)^{P+1})^{frac{1}{2}} \ &equiv ((a+i)^Pcdot (a+i))^{frac{1}{2}} \ &equiv ((a^P+i^P)cdot(a+i))^{frac{1}{2}} \ &equiv ((a-i)cdot(a+i))^{frac{1}{2}} \ &equiv (a^2-i^2)^{frac{1}{2}} \ &equiv (a^2-(a^2-n))^{frac{1}{2}} \ &equiv n^{frac{1}{2}}pmod P end{aligned} ]

由于 ((x - P)^2 = x^2 - 2xP + P^2 = x^2 - P(2x - P) equiv x^2 pmod P)，因此 (x_1 = P - x_0 = P - (a + i)^{frac{P + 1}{2}})。

因而该算法的正确性得以证明。

但是，既然 (a^2 - n) 为非二次剩余，又怎么能开根呢？

这里采用扩张数域的方式，由实数扩大到复数，即 (i) 是一个复数（与实际的复数并不完全相同，与实际中的复数比起来虚数单位由 (sqrt{-1}) 变为了 (sqrt{a^2 - n})），这样就可以解决无法开根的问题，但于此同来的是另一个问题：(x) 的式子中含有 (i) ，那么得到的 (x) 是实数吗？

让我们令 (x_0 = (a + i)^{frac{P + 1}{2}} = A + Bi)，则 ((A + Bi)^2 equiv n pmod P)，假设 (x) 是一个复数，那么 (B eq 0)，可得

[(A + Bi)^2 = A^2 + 2ABi + B^2i^2 equiv n pmod P ]

[A^2 + B^2(a^2 - n) - n equiv -2ABi pmod P ]

由于式子左边不存在虚部，那么式子右边也不存在虚部，即 (AB equiv 0 pmod P)。

由于 (B eq 0)，因此 (A equiv 0 pmod P)，可得

[B^2(a^2 - n) = B^2i^2 equiv n pmod P ]

[i^2 equiv n(B^{-1})^2 pmod P ]

由于 ((B^{-1})^2) 为二次剩余，乘上另一个二次剩余 (n) 之后得到的 (n(B^{-1})^2) 也为二次剩余，这与 (i^2) 是一个非二次剩余相矛盾，故 (B eq 0) 不成立，因此 (B = 0)，即 (x) 为实数。

至此 Cipolla 算法的正确性得到完全证明。

算法复杂度分析

由定理 1 得模 (P) 意义下的非二次剩余的数量为 (frac{P-1}{2}) 个，可以近似看为 (frac{P}{2})，因此在 ([0, P - 1]) 中随机选取一个 (a) 得到非二次剩余的次数期望约为 (2) 次。

因此随机化部分复杂度为 (Theta(1))，而运用欧拉判别准则判断 (a^2 - n) 是否为非二次剩余的复杂度为快速幂的时间复杂度即 (Theta(log P))。

所以 Cipolla 算法单次的时间复杂度为 (Theta(log P))。

例题

洛谷 P5491 【模板】二次剩余，二次剩余模板题，记得特判 (n=0) 的情况和用欧拉判别准则判断是否有解。

代码如下：

/*
    I will never forget it.
*/

// 392699

#include <bits/stdc++.h>

#define int long long

using namespace std;

void fr(int &a, bool f = 0, char ch = getchar()) {
    for (a = 0; ch < '0' || ch > '9'; ch = getchar()) ch == '-' ? f = 1 : 0;
    for (; ch >= '0' && ch <= '9'; ch = getchar()) a = a * 10 + ch - '0';
    a = f ? -a : a;
}
int fr() {
    int a;
    return fr(a), a;
}

int n, P, i2;

struct Complex {
    int real, imag;
    Complex(int real = 0, int imag = 0) : real(real), imag(imag) {}
    bool operator == (const Complex &rhs) const { return real == rhs.real && imag == rhs.imag; }
    Complex operator * (const Complex &rhs) const {
        return Complex((real * rhs.real % P + i2 * imag % P * rhs.imag % P) % P, (imag * rhs.real % P + real * rhs.imag % P) % P);
    }
};

Complex qpow(Complex a, int b) {
    Complex ret(1, 0);
    for (; b; b >>= 1, a = a * a) b & 1 ? (ret = ret * a, 0) : 0;
    return ret;
}

int qpow(int a, int b) {
    int ret = 1;
    for (; b; b >>= 1, a = a * a % P) b & 1 ? ret = ret * a % P : 0;
    return ret;
}

bool check(int x) { return qpow(x, (P - 1) >> 1) == 1; }

struct OI {
    int RP, score;
} CSPS2021, NOIP2021;

signed main() {
    CSPS2021.RP++, CSPS2021.score++, NOIP2021.RP++, NOIP2021.score++, 392699;
    srand(999392699);
    for (int T = fr(), a, x0, x1; T; T--) {
        fr(n), fr(P);
        if (n == 0) {
            puts("0");
            continue;
        } else if (check(n) == 0) {
            puts("Hola!");
            continue;
        }
        a = rand() % P;
        while (a == 0 || check((a * a % P - n + P) % P)) a = rand() % P;
        i2 = (a * a % P - n + P) % P, x0 = qpow(Complex(a, 1), (P + 1) >> 1).real, x1 = P - x0;
        if (x0 > x1) swap(x0, x1);
        if (x0 == x1) printf("%lld
", x0);
        else printf("%lld %lld
", x0, x1);
    }
    return 0;
}