mydcwfy's Blog

二次剩余与 Cipolla 算法

Posted on
Symbols count in article: 2.3k Reading time ≈ 5 mins.

虽说比较困难，但是代码实现十分简单。

1. 定义

我们求解 $x ^ 2 \equiv n \pmod p$，其中 $n, p$ 给定，$p$ 是质数。就是二次剩余，可以记为 $x = \sqrt n$。

一般使用 Cipolla 算法，时间复杂度是玄学，期望是 $O(\log n)$，一般也比较稳定。

2. 推导过程

首先，证明一下在 $\bmod p$ 意义下的有 $\sqrt n$ 的 $n$ 的个数有 $\dfrac{p - 1}{2}$ 个。

首先，我们容易发现，$x$ 与 $p - x$ 在 $\bmod p$ 意义下是相等的。于是我们只需要判断 $x\in[1, \dfrac{p - 1}{2}]$ 的 $x ^ 2$ 是不相同的。这样的话，对于 $\dfrac{p - 1}{2}$ 个 $n$ 都有两个不同的二次剩余对应它。

使用反证法，我们假设 $x_1,x_2\in[1, \dfrac{p - 1}{2}]$ 的平方是一样的，那么可以得到 $(x_2 +x_1)(x_2 - x_1)\equiv 0\pmod p$。但是 $x_2 + x_1$ 和 $x_2 - x_1$ 都是不可能相等的，所以不可能存在这样的 $x_1$ 和 $x_2$。

第二个，我们有一个定理：如果 $n^{\frac{p - 1}{2}}\equiv 1\pmod p$，那么 $n$ 在 $\bmod p$ 下是有二次剩余的。

这一个比较难证，这里就略过了。但是一个可以明确的地方是 $\forall n, n^{\frac{p - 1}2}\bmod p\in \{1, p - 1\}$。因为 $n^{p - 1}\bmod p = 1$。

第三个，我们找到一个 $a$，使得 $a^2 - n$ 是一个非二次剩余。找到 $a$ 的话，我们使用随机化算法，在 $[0, p - 1]$ 中随机，然后判断 $(a ^ 2 - n)^{\frac{p - 1}2}$ 是否等于 $-1\bmod p$。由于有 $\dfrac{p - 1}{2}$ 个数不是二次剩余，所以我们期望 2 次就可以找到一个 $a ^ 2 - n$。这一步的期望是 $O(\log n)$。

设 $w = a ^ 2 - n$。

最后一步，我们定义一个二维向量 $(a, b)$，表示 $a + b\sqrt w$。

我们可以得到第三个结论：将二维向量当作一个数来乘，$\sqrt n = (a, 1) ^{\frac{p - 1}{2}}$。

其中，$(a, 1)$ 表示 $a + \sqrt w$ 这个数，我们将这个数 $\dfrac{p - 1}{2}$ 次方，实部就是 $\sqrt n$。

证明也太难了，我们就不讲了。

3. 代码实现

首先，我们将向量直接看做复数封装好。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
struct Complex{
    ll x, y;
};
ll w;

Complex mul(Complex t1, Complex t2)
{
    return {(t1.x * t2.x % p + t1.y * t2.y % p * w % p) % p, (t1.x * t2.y % p + t1.y * t2.x % p) % p};
}

Complex qpow_com(Complex a, ll k)
{
    Complex res = {1, 0};
    while (k)
    {
        if (k & 1) res = mul(res, a);
        a = mul(a, a);
        k >>= 1;
    }
    return res;
}

其中 qpow_com 表示复数快速幂，其实是比较简单的。

然后就可以直接实现了。

1
2
3
4
5
6
7
8
9
10
11
ll calc(ll n)
{
    ll a;
    while (1)
    {
        a = rand() % p;
        w = (a * a % p - n + p) % p;
        if (qpow(w, (p - 1) / 2) == p - 1) break;
    }
    return qpow_com({a, 1}, (p + 1) / 2).x;
}

完全版。

详细代码 
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
class Cipolla{
private:
struct Complex{
    LL x, y;
};
LL w, Mod;

Complex mul(Complex a, Complex b)
{
    return {(a.x * b.x % Mod + a.y * b.y % Mod * w % Mod) % Mod, (a.x * b.y + a.y * b.x) % Mod};
}

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

LL qpow(LL a, LL k)
{
    LL res = 1;
    while (k)
    {
        if (k & 1) res = res * a % Mod;
        a = a * a % Mod;
        k >>= 1;
    }
    return res;
}

public:
Cipolla(int _Mod = 998244353) : Mod(_Mod) {}

LL Mod_Sqrt(LL n)
{
    LL a;
    while (1)
    {
        a = rand() % Mod;
        w = (a * a - n + Mod) % Mod;
        if (qpow(w, (Mod - 1) >> 1) == Mod - 1) break;
    }
    return qpow_com({a, 1}, (Mod + 1) >> 1).x;
}
};