显然，我们只需考虑 $0 \leq a < n$ 的情形。这个问题应该很早就被人考虑好了，不过无所谓吧（反正只是个博文而已）。以下内容都是独立完成的。并可看作 二次剩余和 Gauss 互反律 这篇博文的延续。最后再给出方程的一个解（如果存在的话）。

记 $A_n = \{0 < x < n \mid \gcd(x,n) = 1\}$，则 $A_n$ 关于 $\mod n$ 的乘法构成一个群（根据 Euler 定理 $|A_n| = \psi(n)$）

所以若 $\gcd(a,n)= 1$，则任意 $x \in A_n$，存在唯一 $y \in A_n$ 使得 $xy \equiv a \mod n$

若 $x^2 \equiv a \mod n$ 无解，则 $\prod_{x \in A_n} x = a^{\frac{\psi(n)}{2}} \mod n$

若 $x^2 \equiv a \mod n$ 有解，设 $0 < a_1 < a_2 < \cdots < a_r < n$ 是全部解。则 $\prod_{x \in A_n} x = a^{\frac{\psi(n)-r}{2}} a_1 \cdots a_r \mod n$，注意到 $r$ 只与 $n$ 有关，与 $a$ 无关。

用 $[\frac{a}{n}] = 1$ 表示有解，用 $[\frac{a}{n}] = 0$ 表示无解。

为了考虑这个问题，我们分情况讨论

$x^2 \equiv a \mod 2^m, \gcd(a,p^m) = 1$ 何时有解

当 $m=1,2$ 时，有解当且仅当 $a=1$

当 $m \geq 3$ 时，$2^m | (a_2 -a_1)(a_2+a_1)$, 若 $(a_2-a_1)$ 和 $a_2+a_1$ 都是 4 的倍数，则 $a_1,a_2$ 也是 2 的倍数矛盾于 $\gcd(a,2^m) = 1$。所以 $2^{m-1} \mid a_2+a_1)$ 或者 $2^{m-1} \mid (a_2-a_1)$。所以 $0 < a_1 < 2^{m-2}$，$a_2 = 2^{m-1} - a_1, a_3 = 2^{m-1}+a_1, a_4 = 2^m-a_1$。

但是，若 $x$ 是奇数，则 $x^2 \equiv 1\mod 8$ ，令一方面 $\gcd(a,2^m) = 1$ 当且仅当 $a \equiv \mod 2$，而且每一个解对应 4 个 $x$。

即 $1^2,2^2,\cdots,(2^m-1)^2 \mod 2^m$，只有 $2^{m-3}$ 个数，但是小于 $2^m$ 且模 8 为 1 的数也只有，$2^{m-3}$ 个。因此，

有解 当且仅当 $a \equiv 1 \mod 8$

$x^2 \equiv a \mod p^m,\;p>2, \gcd(a,p^m) = 1$ 何时有解

因为 $p^m \mid (a_2-a_1)(a_2+a_1)$，若$(a_2-a_1)$ 和 $a_2+a_1$ 都是 $p$ 的倍数，则 $a_1, a_2$ 也是 $p$ 的倍数矛盾于 $\gcd(a,p^m) = 1$，所以 $a_2 = p^m - a_1$。从而 $r=2, \; a_1a_2 \equiv -a \mod p^m$，所以

  • 若 $x^2 \equiv a \mod p^m$ 无解，则 $\prod_{x \in A_n} x = a^{\frac{p^{m-1}(p-1)}{2}} \mod p^m$

  • 若 $x^2 \equiv a \mod p^m$ 有解，则 $\prod_{x \in A_n} x = -a^{\frac{p^{m-1}(p-1)}{2}} \mod p^m$

取 $a=1$，则 $\prod_{x \in A_n} x = -1 \mod p^m$，即

$x^2 \equiv a p^i \mod p^m, \; i<m, \; \gcd(a,p^m) = 1$ 何时有解

若 $i=2k-1$ 为奇数，则 $p^k \mid x$，从而 $p^{i+1} \mid ap^i$，矛盾于 $\gcd(a,p^m) = 1$

所以 $i$ 为偶数，显然此时上式有解当且仅当 $x^2 \equiv a \mod p^{m-i},\gcd(a,p^{m-i}) = 1$ 有解

综上

$x^2 \equiv a \mod m_1m_2, \gcd(m_1,m_2)=1$ 何时有解

考虑

则，存在 $t_1,t_2$ 使得，$t_1m_1 + t_2m_2 = 1$，所以 $t_1^2 m_1^2 + t_2^2 m_2^2 + 2 t_1 m_1 t_2 m_2=1$, 所以 $x \equiv (t_2m_2 a_1)^2$ 且 $x \equiv (t_1 m_1 a_2)^2$，所以

令一方面若 $x \equiv a^2 \mod m_1 m_2$，则 $x \equiv (a \mod m_i)^2 \mod m_i$

即 $x^2 = a \mod m_1 m_2$ 有解当且仅当 $[\frac{a}{m_1}] = [\frac{a}{m_2}] = 1$

最终方案：$n = p_1^{m_1} \cdots p_r^{m_r}, \; p_1< \cdots <p_r$, $[\frac{a}{m}] = \prod [\frac{a}{p_i^{m_i}}]$

判断是否有解的 sagemath 代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# modSqrtSymbol.ipynb
def modSqrtSymbol(an,n):
    for p,m in factor(n):
        a = an%p^m
        if a == 0:  continue
        flag = True
        while a%p == 0:
            a//=p
            m-=1
            flag = not flag
        if not flag:  return False
        if p == 2:
            if a%8 != 1:  return  False
        else:
            if power_mod(a,p^(m-1)*(p-1)//2,p^m) != 1:
                return False
    return True

因为 factor 复杂度是 $O(\sqrt{n})$， 所以整体复杂度 $O(\sqrt{n})$

求出 $x^2 \equiv a \mod n$ 的一个解

本来是想求全部解的，但是考虑到如果 $\gcd(a,n) \neq 1$ 则解的个数实在太多，就不想求全部解了，当然如果真的有这个需要，求全部解也不难。

另外，无论是是否有解还是求全部解都可以直接暴力做，只是复杂度为 $O(n)$ ，效率太低。

求 $x^2 \equiv a \mod 2^m$ 的全部解，不妨设 $a \equiv 1 \mod 8$

若 $m=1$，则 $x=1$，若 $m=2$，则 $x = 1, 3$，若 $m=3$，则 $x=1,3,5,7$，若 $m \geq 4$，不妨设 $x_{m-1}$ 是满足方程 $x^2 \equiv a \mod 2^{m-1}$的最小正整数解，即 $x_{m-1}^2 = a + 2^{m-1}k$

  • 若 $k$ 是偶数，则必然是 $x_{m-1}^2 \equiv a \mod 2^m$ 的解，最小的解就是 $\min(x_{m-1},2^{m-2}-x_{m-1})$，但由于 $x_{m-1} \leq 2^{m-3}$，所以 $x_m=x_{m-1}$

  • 若 $k$ 为奇数，则 $(2^{m-2} - x_{m-1})^2 = 2^{2m-4} - 2^{m-1} + x_{m-1}^2 \equiv a \mod 2^m$，且 $0<2^{m-2}-x_{m-1}<2^{m-2}$，所以 $x_m =2^{m-2} - x_{m-1}$

所以经过 $m$ 步就可找到所有解就是 $x_m, 2^{m-1} - x_m,2^{m-1} + x_m, 2^m - x_m$

求 $x^2 \equiv a \mod p^m, \; p>2$ 的全部解，其中 $\gcd(a,p^m) = 1$

令 $q = \frac{\psi(p^m)}{2} = p^{m-1}\frac{p-1}{2}$，若 $q$ 是奇数，则 $(a^{\frac{q+1}{2}})^2 = a^q a \equiv a\mod p^m$，从而解就是 $a^{\frac{q+1}{2}},p^m-a^{\frac{q+1}{2}}$

否则，记 $q = 2^{i}q’$（想模仿上面 $q$ 是奇数的情况），我们找一个非二次剩余 $b$，即 $b^q \equiv -1 \mod p^m$。我们记 $x = a^{\frac{q’+1}{2}},\;y = b^{q’}, \; t= a^{q’}$。则

$x^2 = ta, t^{2^i} = 1, y^{2^i}=-1$。我们想保持 $x^2 = ta$，让 $t$ 变成 1。就打到了我们的目的。

若 $t^{2^{s}} = -1, t^{2^{s+1}} = 1$, 则 $s<i$，那么 $(y^{2^{i-s}}t)^{2^{s}} = y^{2^i}t^{2^{s}} = 1$，并且 $(y^{2^{i-s-1}}x)^2 = (y^{2^{i-s}}t)a$，我们更新 $t = y^{2^{i-s}}t$，更新 $x=y^{2^{i-s-1}}x$。重复上述过程直到 $t=1$。

我们也可以类似 $p=2$ 的情况，经过 $m$ 步得到全部解，只是 $m=1$ 的时候是困难的，并且和上述分析一致，所以就没用 $p=2$ 时的方法。若 $x_{m-1}$ 满足 $x_{m-1}^2 = a + p^{m-1}k$

  • 若 $k$ 是偶数，$(x_{m-1} - \frac{k}{2}p^{m-1})^2 \equiv a \mod p^m$

  • 若 $k$ 是奇数，$(x_{m-1} + \frac{p-k}{2}p^{m-1})^2 \equiv a \mod p^m$

对一般的 $n = p_1^{m_1} \cdots p_r ^{m_r}$，我们对每个 $(a, p_i ^m)$ 给出一个解，再用 中国剩余定理 拼出最终解。

$x^2 \equiv a \mod n$ 的一个解的 Sagemath 代码

需要用到上面的 modSqrtSymbol 函数

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
66
# modSqrt.ipynb

# find one x s.t. x^2 = a mod 2^m, gcd(a,2) = 1 and 0 < a < 2^m and must have an answer
def modSqrt2Core(a,m):
    x = 1
    for i in range(4,m+1):
        if (x*x-a)%(1<<i) != 0:
            x = (1<<i-2)-x
    return x

# find one x s.t. x^2 = a mod p^m, gcd(a,p)=1, p>2 and  0 < a < p^m and must have an answer
def modSqrtPCore(a,p,m):
    n = p^m
    q = p^(m-1)*(p-1)>>1
    if q%2 == 1:  return power_mod(a,(q+1)>>1,n)
    b=randint(0,n-1)
    while(power_mod(b,q,n)==1):  b=randint(0,n-1)
    ni = q&(-q)
    q //= ni
    x = power_mod(a,(q+1)>>1,n)
    y = power_mod(b,q,n)
    t = power_mod(a,q,n)
    # x^2 = ta mod n, y^ni = -1 mod n
    while t != 1:
        ns = 1;tt=t*t;
        while tt!=1:
            tt=(tt*tt)%n
            ns<<=1
        # t^ns = -1 mod n
        t = power_mod(y,ni//ns,n)*t%n
        x = power_mod(y,ni//(ns*2),n)*x%n
    return x

# find one x s.t. x^2 = a mod p^m and must have an answer
def modPowerSqrt(a,p,m):
    n = p^m
    a %= n
    if a == 0:  return 0
    time = 0
    while a%p==0:
        a//=p
        m-=1
        time+=1
    if a == 0:  return 0
    return 2**(time//2)*(modSqrt2Core(a,m) if p==2 else modSqrtPCore(a,p,m))

def modSqrt(a,n):
    if not modSqrtSymbol(a,n):  return []
    fact = list(factor(n))
    x = [modPowerSqrt(a,p,m) for p,m in fact]
    m = [p^m for p,m in fact]
    return crt(x,m)


# 测试几个样例
a,m = 17,2**8
x =  modSqrt(a,m)
print(x,x*x%m,a,m)

a,m = 15,7**8
x =  modSqrt(a,m)
print(x,x*x%m,a,m)

a,m = 16,3*5*5*7*7*7
x =  modSqrt(a,m)
print(x,x*x%m,a,m)