一种优美而神奇的数论算法。

二次剩余

        考虑下面方程的求解：

        这里的$p$为奇质数。（可以是其它数，这里只考虑奇质数。）
        当解存在时，称$n$为二次剩余，x为它的解。
        倘若存在两个解$x_1$与$x_2$有$x1\not \equiv x2\pmod p$，但有$x_1^2\equiv x_2^2\pmod p$，则：

        于是$x_1\equiv -x_2\pmod p$，说明当有多解时，它们一定互为相反数，且当$x_1\equiv 0\pmod p$时，必有$x_1\not \equiv -x_1=x_2\pmod p$。
        这里证明了解最多有两个。

有解性的欧拉准则

        考虑$n$是否为模$p$意义下的二次剩余。
        当$n\equiv 0\pmod p$时，显然有解，这里不再考虑。
        由费马小定理得：

        即：

        $1$显然是有解的，这样可得$n^{\frac {p-1}{2}}$为$1$或$-1$。
        下证$n$是二次剩余与$n^{\frac {p-1}{2}}\equiv 1\pmod p$等价。
        若$n$是二次剩余，则有：

        反过来，将$n$表示成原根的形式$g^k$，$g$为$p$的原根。则：

        由于$g$为原根，它的指数为阶（即欧拉函数$\phi(p)=p-1$）的倍数，有：

        于是必有$2|k$，那么$g^{\frac {k}{2}}$就是$x^2\equiv n\pmod p$的根。证毕。
        同时可以知道$n$不是二次剩余与$n^{\frac {p-1}{2}}\equiv -1\pmod p$等价。
        定义勒让德符号：

        得到欧拉准则：

        

Cipolla算法

        现在来求$x^2\equiv n\pmod p$的解。
        先找一个数$a$，使得$a^a-n$为非二次剩余。由于每一对相反数都对应一个二次剩余，并且它们互不相同，这样就会有$\frac{p-1}{2}$个二次剩余，找$a$的期望次数为2。
        接下来定义$i^2\equiv a^a-n\pmod p$，这里的$i$定义了一个不存在的数，相当于虚数。它具有下面的性质：

  1. $i^p\equiv -i\pmod p$
  2. $(a+b)^p\equiv a^p+b^p\pmod p$
            展开二项式：        对于$0<i<p$，分子阶乘上的$p$消不掉（因为$p$是质数），因此只有$a^p+b^p$剩下来。

        由上面两条性质可得：

        这样$(a+i)^{\frac {p+1}{2}}$就是解，另一个是它的相反数。可以证明这个数的“虚部”为0，答案就是其实部的整数。
        模板。

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
#include <bits/stdc++.h>
using namespace std;
int n, p, a2mn;
struct Complex
{
    int a, b;
    Complex(int aa, int bb)
        : a(aa), b(bb) {}
    void operator*=(Complex x)
    {
        int ta = (1ll * a * x.a % p + 1ll * b * x.b % p * a2mn % p) % p;
        int tb = (1ll * a * x.b % p + 1ll * b * x.a % p) % p;
        a = ta, b = tb;
    }
};

inline int qPow(int x, int y)
{
    int ans = 1, sta = x;
    while (y) {
        if (y & 1) ans = 1ll * ans * sta % p;
        sta = 1ll * sta * sta % p, y >>= 1;
    }
    return (ans % p + p) % p;
}
int main()
{
    srand(time(nullptr));
    int t;
    scanf("%d", &t);
    while (t--) {
        scanf("%d%d", &n, &p);
        if (n == 0) {
            printf("0\n");
            continue;
        }
        if (qPow(n, (p - 1) >> 1) == p - 1) {  //无解性判定
            printf("Hola!\n");
        }
        else {
            int a, x1, x2;
            do
                a = rand() % (p - 1) + 1;
            while (qPow((1ll * a * a - n) % p, (p - 1) >> 1) == 1);
            a2mn = (1ll * a * a - n) % p;
            Complex s(a, 1), ans(1, 0);
            a = (p + 1) >> 1;
            while (a) {
                if (a & 1) ans *= s;
                s *= s, a >>= 1;
            }
            x1 = (ans.a + p) % p, x2 = p - x1;
            printf("%d %d\n", min(x1, x2), max(x1, x2));
        }
    }
    return 0;
}