$x^2 \equiv a \mod n$ 何时有解

最后更新于：2022年6月25日 下午

显然，我们只需考虑 \(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\)，即 \[ [\frac{a}{p^m}] = (1 \equiv a^{\frac{p^{m-1}(p-1)}{2}} \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\) 有解

综上 \[ [ \frac{a}{p^m} ] = \left\{ \begin{array}{cc} 1 & a=0 \\ 1 & p=2,\; a \equiv 1 \mod 8 \\ -1 & p=2,\; a \not \equiv 1 \mod 8 \\ (1 \equiv a^{\frac{p^{m-1}(p-1)}{2}} \mod p^m)? & p>2,\; \gcd(a,p^m) = 1\\ -1 & \gcd(a,p^m) = p^i, i \equiv 1 \mod 2 \\ [\frac{ap^{-i}}{p^{m-i}}] & \gcd(a,p^m) = p^i, i \equiv 0 \mod 2 \\ \end{array} \right. \]

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

考虑 \[ \begin{aligned} x \equiv a_1^2 \mod m_1 \\ x \equiv a_2^2 \mod m_2 \end{aligned} \]

则，存在 \(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 (t_2 m_2 a_1)^2 + (t_1 m_1 a_2)^2 \equiv (t_2 m_2 a_1 + t_1 m_1 a_2)^2 \mod m_1 m_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
67
68
69
70
# 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)