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

二次剩余

        考虑下面方程的求解:

        这里的$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;
}