简介

二次剩余即模意义开根，求解形如 x2≡n(modp)x^2\equiv n\pmod p 的方程，若方程存在非 00 解，称 nn 是模 pp 的二次剩余，欧拉准则可以判断 nn 是否为模 pp 的二次剩余，本文只讨论 pp 为奇素数的情况。

欧拉准则

当 pp 为奇素数时，根据费马小定理：

n2p−12≡1(modp)(np−12−1)(np−12+1)≡0(modp)\begin{aligned} n^{2{\frac{p-1}{2}}}&\equiv 1\pmod p\\ (n^{\frac{p-1}{2}}-1)(n^{\frac{p-1}{2}}+1)&\equiv 0\pmod p\\ \end{aligned}

为方便表达引入勒让德符号：

(np)≡np−12(modp)(\frac{n}{p})\equiv n^{\frac{p-1}{2}}\pmod p

下面证明 nn 为模 pp 的二次剩余是 (np)≡1(modp)(\frac{n}{p})\equiv 1\pmod p 的充要条件。

  • 充分性：若存在 x2≡n(modp)x^2\equiv n\pmod p 那么 np−12≡xp−1≡1(modp)n^{\frac{p-1}{2}}\equiv x^{p-1}\equiv 1\pmod p。
  • 必要性：设 gg 为模 pp 的原根，由于 pp 为奇素数所以一定存在，令 gk≡n(modp)g^k\equiv n\pmod p，那么 gkp−12≡1(modp)g^{k\frac{p-1}{2}}\equiv 1\pmod p，那么就有 kp−12≡0(modφ(p))k\frac{p-1}{2}\equiv 0\pmod {\varphi(p)}，也就是说 (p−1)∣k2(p−1)(p-1)\mid \frac{k}{2}(p-1)，所以 kk 是偶数，那么 x=gk2(modp)x=g^{\frac{k}{2}}\pmod p 就是原方程的一个解。

关于解的数量，当 x1x_1 是原方程的一个解时，若存在 x2≢x1(modp)x_2\not\equiv x_1\pmod p 并且 x12≡x22x_1^2\equiv x_2^2 可以推出：

(x1−x2)(x1+x2)≡0(modp)(x_1-x_2)(x_1+x_2)\equiv 0\pmod p

所以只有可能是 x1≡−x2(modp)x_1\equiv -x_2\pmod p，方程有两个解，当 n=0n=0 时 x1≡x2≡0(modp)x_1\equiv x_2\equiv0\pmod p。

Cipolla

首先找到一个 aa 使得 (a2−np)≡−1(modp)(\frac{a^2-n}{p})\equiv -1\pmod p，设 gg 为模 pp 的原根，且 a2−n≡gk(modp)a^2-n\equiv g^k\pmod p 那么 (a2−n)p−12≡gkp−12(modp)(a^2-n)^{\frac{p-1}{2}}\equiv g^{k\frac{p-1}{2}}\pmod p，kk 的奇偶性决定 (a2−n)p−12(a^2-n)^{\frac{p-1}{2}} 的值，因此平均两次就可以找到满足条件的 aa。

类比虚数的定义，令 i2≡a2−n(modp)i^2\equiv a^2-n\pmod p，那么有一个很好的性质：

ip≡i⋅ip−1≡i(a2−n)p−12≡−i(modp)i^p\equiv i\cdot i^{p-1}\equiv i(a^2-n)^{\frac{p-1}{2}}\equiv -i\pmod p

下面证明 (a+i)p+1≡n(modp)(a+i)^{p+1}\equiv n\pmod p，其中 (a+i)p≡(ap+ip)(modp)(a+i)^p\equiv (a^p+i^p)\pmod p 原理是二项式定理：

(a+i)p+1≡(a+i)p(a+i)≡(ap+ip)(a+i)≡(a−i)(a+i)≡a2−i2≡n(modp)\begin{aligned} (a+i)^{p+1}&\equiv(a+i)^p(a+i)\\ &\equiv(a^p+i^p)(a+i)\\ &\equiv(a-i)(a+i)\\ &\equiv a^2-i^2\\ &\equiv n\pmod p \end{aligned}

所以 x=(a+i)p+12x=(a+i)^{\frac{p+1}{2}} 就是原方程的一个解，下面用反证法证明 xx 的虚部为 00：

设存在 c,dc,d 使得 (c+di)2≡n(modp)(c+di)^2\equiv n\pmod p 且 d≠0d\ne 0，那么有：

c2+2cdi+d2(a2−n)≡n(modp)c^2+2cdi+d^2(a^2-n)\equiv n\pmod p\\

既然 d≠0d\ne 0 一定是 c=0c=0 否则无法满足等式，此时：

a2−n≡nd−2(modp)a^2-n\equiv nd^{-2}\pmod p

由于 nd−2nd^{-2} 一定是二次剩余，所以 a2−na^2-n 也是二次剩余，与假设矛盾，所以 xx 的虚部一定为 00。

因此算法的时间复杂度主要在求幂上，用快速幂做到 O(log⁡p)O(\log p)，题目是 P5491。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
#include <bits/stdc++.h>
using namespace std;

namespace IO {}
using IO::read;
using IO::write;

#define int long long
int T, n, p, i2;

struct Complex {
    int real, imag;
    Complex operator*(const Complex& c) const {
        // 这里三个大数相乘注意取模
        return {(real * c.real + imag * c.imag % p * i2) % p, (real * c.imag + imag * c.real) % p};
    }
};

mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());

int qmi(int a, int k) {
    int res = 1;
    while (k) {
        if (k & 1) res = res * a % p;
        a = a * a % p;
        k >>= 1;
    }
    return res;
}

Complex qmi(Complex a, int k) {
    Complex res = {1, 0};
    while (k) {
        if (k & 1) res = res * a;
        a = a * a;
        k >>= 1;
    }
    return res;
}

int cipolla() {
    if (!n) return 0;
    if (qmi(n, (p-1) / 2) != 1) return -1;

    int a;
    do {
        a = uniform_int_distribution<>(1, p-1)(rng);
        i2 = (a * a - n) % p;
        if (i2 < 0) i2 += p;
    } while (qmi(i2, (p-1) / 2) == 1);

    return qmi(Complex{a, 1}, (p+1) / 2).real;
}

signed main() {
    read(T);
    while (T--) {
        read(n, p);
        int x = cipolla(), y = p-x;
        if (x == 0) write(0);
        else if (x == -1) write("Hola!");
        else x > y && (swap(x, y), 0), write(x, ' '), write(y);
    }
    return 0;
}

扩域

Cipolla 算法用到了扩域的思想，那么这个思想其实可以应用到诸多递推数列上，如果存在通项公式，就能通过扩域的方法求解。

众所周知，斐波那契的通项公式是：

fn=15((1+52)n−(1−52)n)f_n=\frac{1}{\sqrt 5}\left(\Big(\frac{1+\sqrt 5}{2}\Big)^n-\Big(\frac{1-\sqrt 5}{2}\Big)^n\right)

如果模数 p=109+7p=10^9+7，那么 55 是非二次剩余，但是可以这样定义：

i2≡5(modp)i^2\equiv 5\pmod p

此时可以改写通项公式了：

fn≡i5((1+i2)n−(1−i2)n)(modp)f_n\equiv \frac{i}{5}\left(\Big(\frac{1+i}{2}\Big)^n-\Big(\frac{1-i}{2}\Big)^n\right)\pmod p

实现了 P1962：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
#include <bits/stdc++.h>
using namespace std;

namespace IO {}
using IO::read;
using IO::write;

#define int long long
const int P = 1000000007, I5 = 400000003, I2 = 500000004;
int n;

struct Complex {
    int real, imag;

    Complex operator*(const Complex& c) const {
        return {(real * c.real + imag * c.imag * 5) % P, (real * c.imag + imag * c.real) % P};
    }

    Complex operator-(const Complex& c) const {
        int a = (real - c.real) % P, b = (imag - c.imag) % P;
        if (a < 0) a += P;
        if (b < 0) b += P;
        return {a, b};
    }
};

Complex qmi(Complex a, int k) {
    Complex res = {1, 0};
    while (k) {
        if (k & 1) res = res * a;
        a = a * a;
        k >>= 1;
    }
    return res;
}

signed main() {
    read(n);
    Complex a = {0, I5}, b = qmi({I2, I2}, n) - qmi({I2, P-I2}, n);
    return write((a * b).real), 0;
}