M.Hiroi's Home Page / Memorandum 数学編

初等整数論 [4]

●平方剰余

合同式 \(x^2 \equiv a \pmod n\) を満たす整数 \(x\) が存在するとき、\(a\) は \(n\) を法とする「平方剰余」といいます。整数 \(x\) が存在しない場合、\(a\) は \(n\) を法とする「平方非剰余」といいます。簡単な例を示しましょう。次の表を見てください。

                      mod n
 x : x*x :  2  3  4  5  6  7  8  9 10 11
---+-----+-------------------------------
 0 :   0 :  0  0  0  0  0  0  0  0  0  0
 1 :   1 :  1  1  1  1  1  1  1  1  1  1
 2 :   4 :     1  0  4  4  4  4  4  4  4
 3 :   9 :        1  4  3  2  1  0  9  9
 4 :  16 :           1  4  2  0  7  6  5
 5 :  25 :              1  4  1  7  5  3
 6 :  36 :                 1  4  0  6  3
 7 :  49 :                    1  4  9  5
 8 :  64 :                       1  4  9
 9 :  81 :                          1  4
10 : 100 :                             1

       :    平方剰余      :  平方非剰余
-------+------------------+---------------
mod  2 : 0, 1             : None
mod  3 : 0, 1             : 2
mod  4 : 0, 1             : 2, 3
mod  5 : 0, 1, 4          : 2, 3
mod  6 : 0, 1, 3, 4       : 2, 5
mod  7 : 0, 1, 2, 4       : 3, 5, 6
mod  8 : 0, 1, 4          : 2, 3, 5, 6, 7
mod  9 : 0, 1, 4, 7       : 2, 3, 5, 6, 8
mod 10 : 0, 1, 4, 5, 6, 9 : 2, 3, 7, 8
mod 11 : 0, 1, 3, 4, 5, 9 : 2, 6, 7, 8, 10

上の表では x * x % n を計算しています。0 から x 未満の数で、式 x * x % n で現れる数が平方剰余、現れない数が平方非剰余となります。たとえば、mod 2 の世界では、数は 0 と 1 しかありませんが、0 * 0 % 2 = 0, 1 * 1 % 2 = 1 となるので、0, 1 は平方剰余、平方非剰余の数はありません。

同様に計算していくと、たとえば mod 7 の世界では、0, 1, 2, 4 が平方剰余、3, 5, 6 が平方非剰余となります。つまり、平方数を 7 で割り算した場合、余りが 3, 5, 6 になることはない、ということです。平方数を 3 で割り算すると余りが 0, 1 だけで 2 になることはない、平方数を 4 で割り算すると余りは 0, 1 だけで 2, 3 になることはない、などはよく知られている性質です。

\(m = 3k\) (k は整数) とすると、\(m^2 = 9k^2\) となるので、3 で割り切れる (余り 0)
\(m = 3k + 1\) とすると、\(9k^2 + 6k + 1\) となるので、3 で割ると 1 余る
\(m = 3k + 2\) とすると、\(9k^2 + 12k + 4\) となるので、3 で割ると 1 余る
余りは 0, 1 だけで 2 になることはない

同様に、\(4k, \ 4k + 1, \ 4k + 2, \ 4k + 3\) の二乗を計算すると、0, 1 が平方剰余で 2, 3 が平方非剰余になることが分かります。

\(\begin{array}{l} (4k)^2 = 16k^2 \equiv 0 \pmod 4 \\ (4k + 1)^2 = 16k^2 + 8k + 1 \equiv 1 \pmod 4 \\ (4k + 2)^2 = 16k^2 + 16k + 4 \equiv 0 \pmod 4 \\ (4k + 3)^2 = 16k^2 + 24k + 9 \equiv 1 \pmod 4 \end{array}\)

●ルジャンドル記号とオイラーの基準

数学の整数論では、\(a\) が \(n\) の平方剰余か否かを表す「ルジャンドル記号」が定義されています。

a が n の平方剰余であるとき

\(\left(\dfrac{a}{n}\right) = 1\)

a が n の平方非剰余であるとき

\(\left(\dfrac{a}{n}\right) = -1\)

ルジャンドル記号は (a | n) や (a / n) と表記されることがあります。また、\(a \equiv 0 \pmod n\) ならば (a / n) = 0 とする定義もあります。

平方剰余を判定する定理に「オイラーの基準」があります。

オイラーの基準
n を奇素数、a と n は互いに素な 0 ではない整数とすると、以下の式が成り立つ

\( \left(\dfrac{a}{n}\right) \equiv a^{\frac{n-1}{2}} \pmod n \)

簡単な例を示しましょう。

>>> for x in [3,5,7,11,13]:
...     for a in range(1, x):
...         print(a ** ((x - 1) // 2) % x, end=" ")
...     print()
...
1 2
1 4 4 1
1 1 6 1 6 6
1 10 1 1 1 10 10 10 1 10
1 12 1 1 12 12 12 12 1 1 12 1

このように、オイラーの基準で平方剰余を判定することができます。

証明
ルジャンドル記号 (a / n) = 1 の場合
\(x^2 \equiv a \pmod n\) を満たす整数 x が存在する

\( a^{\frac{n-1}{2}} = (x^2)^{\frac{n-1}{2}} = x^{n-1} \)

フェルマーの小定理により \(x^{n-1} \equiv 1 \pmod n\), よって (a / n) = 1 となる

ルジャンドル記号 (a / n) = -1 の場合
次の定理が正しいと認める

[定理] a と n が互いに素のとき、a, 2a, ..., (n - 1)a を n で割った余りはすべて異なる

\(x = (1, 2, \cdots, n - 1)\) に要素 \(x_i\) を選んで乗算する
n で割った余りはすべて異なるので、その中で \(x_i \times x_j \equiv a \pmod n\) となる \(x_j\) が必ずひとつ存在する
そのペアの総数は \(\frac{n - 1}{2}\) 個になるので、次式が成り立つ

\( (n - 1)! \equiv a^{\frac{n-1}{2}} \pmod n \)

ウィルソンの定理 \((n - 1)! \equiv -1 \pmod n\) より、\(a^{\frac{n-1}{2}} \equiv -1 = \left(\frac{a}{n}\right)\) となる

SymPy にはルジャンドル記号を求める関数 legendre_symbol(a, n) が用意されています。簡単な実行例を示します。

>>> for a in range(1, 20): print(sy.legendre_symbol(a, 3), end=" ")
...
1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1
>>> for a in range(1, 20): print(sy.legendre_symbol(a, 5), end=" ")
...
1 -1 -1 1 0 1 -1 -1 1 0 1 -1 -1 1 0 1 -1 -1 1
>>> for a in range(1, 20): print(sy.legendre_symbol(a, 7), end=" ")
...
1 1 -1 1 -1 -1 0 1 1 -1 1 -1 -1 0 1 1 -1 1 -1
>>> for a in range(1, 20): print(sy.legendre_symbol(a, 11), end=" ")
...
1 -1 1 1 1 -1 -1 -1 1 -1 0 1 -1 1 1 1 -1 -1 -1
>>> for a in range(1, 20): print(sy.legendre_symbol(a, 13), end=" ")
...
1 -1 1 1 -1 -1 -1 -1 1 1 -1 1 0 1 -1 1 1 -1 -1

●平方剰余の相互法則

平方剰余 \(x^2 \equiv a \pmod p\) の解の有無を判定するとき、「平方剰余の相互法則」を使うと簡単に判定できることがあります。

[平方剰余の相互法則]
相異なる奇素数 p, q に対して、次式が成立する

\( \left(\dfrac{q}{p}\right) \left(\dfrac{p}{q}\right) = (-1)^{\frac{(p-1)(q-1)}{4}} \)

\(x^2 \equiv q \pmod p\) の解の有無がわかれば、\(x^2 \equiv p \pmod q\) の解の有無もわかる
その逆もしかり

\(\begin{eqnarray} (7 / 11) &=& (-1)^{6 \times 10/4} \times (11 / 7) \\ &=& (-1)^{15} \times (4 / 7) \\ &=& - (4 / 7) \end{eqnarray}\)

>>> - legendre_symbol(4, 7)
-1
>>> legendre_symbol(7, 11)
-1

\(\begin{eqnarray} (7 / 13) &=& (-1)^{6 \times 12/4} \times (13 / 7) \\ &=& (-1)^{18} \times (6 / 7) \\ &=& (6 / 7) \end{eqnarray}\)

>>> legendre_symbol(6, 7)
-1
>>> legendre_symbol(7, 13)
-1

第一補充法則

\( \left(\dfrac{-1}{p}\right) = (-1)^{\frac{p-1}{2}} \)

\(\frac{p - 1}{2}\) が偶数であれば、-1 が法 p の平方剰余になる
\(\frac{p - 1}{2} = 2k\) => \(p = 4k + 1\)
つまり、p が \(4k + 1\) 型の素数ならば、平方剰余は \(p - 1 \ (\equiv -1 \pmod p)\)

\(\begin{array}{l} x^2 \equiv 4 \pmod 5 \\ x^2 \equiv 12 \pmod {13} \\ x^2 \equiv 16 \pmod {17} \\ x^2 \equiv 28 \pmod {29} \end{array}\)

>>> legendre_symbol(4, 5)
1
>>> legendre_symbol(12, 13)
1
>>> legendre_symbol(16, 17)
1
>>> legendre_symbol(28, 29)
1

第二補充法則

\( \left(\dfrac{2}{p}\right) = (-1)^{\frac{p^2 - 1}{8}} \)

>>> legendre_symbol(2, 3)
-1
>>> (-1)**((3**2 - 1) // 8)
-1
>>> legendre_symbol(2, 5)
-1
>>> (-1)**((5**2 - 1) // 8)
-1
>>> legendre_symbol(2, 7)
1
>>> (-1)**((7**2 - 1) // 8)
1
>>> legendre_symbol(2, 11)
-1
>>> (-1)**((11**2 - 1) // 8)
-1

整数 \(a, b\) が \(p\) と互いに素の場合、次式が成り立つ

\( \left(\dfrac{ab}{p}\right) = \left(\dfrac{a}{p}\right) \times \left(\dfrac{b}{p}\right) \)

(4 / 7) = (2 / 7) * (2 / 7) = 1 * 1 = 1

>>> legendre_symbol(4, 7)
1

(6 / 11) = (2 / 11) * (3 / 11) = (-1) * (3 / 11)

>>> - legendre_symbol(3, 11)
-1
>>> legendre_symbol(6, 11)
-1

(15 / 11) = (3 / 11) * (5 / 11) = (5 / 11)

>>> legendre_symbol(5, 11)
1
>>> legendre_symbol(15, 11)
1

●合同式 x² ≡ a (mod n) の解き方

合同式 \(x^2 \equiv a \pmod n\) の解 \(x\) は、\(n\) が \(4k + 3 \ (k = 0, 1, 2, \cdots)\) 型の奇素数ならば、次の公式で求めることができます。

\( x \equiv \pm a^{\frac{n+1}{4}} \pmod n \)

上の公式はオイラーの基準から簡単に導出することができます。

\(\begin{array}{l} x^2 \equiv a \pmod n\\ a = a^{\frac{n+1}{2}} \times a^{\frac{n-1}{2}} \end{array}\)

a は平方剰余なので (a / n) = 1, よって \(a^{\frac{n-1}{2}} \equiv 1 \pmod n\)

\(\begin{array}{l} x^2 \equiv a^{\frac{n+1}{2}} \pmod n \\ x \equiv \pm a^{\frac{n+1}{4}} \pmod n \end{array}\)

ただし \(\dfrac{n + 1}{4} = k \)は整数であること

\(n = 4k - 1 \ (k = 1, 2, \cdots) \ or \ n = 4k + 3 \ (k = 0, 1, 2, \cdots)\)

それでは実際に試してみましょう。

>>> pow(1, (3 + 1) // 4, 3)
1
>>> 1 * 1 % 3
1

>>> pow(1, (7 + 1) // 4, 7)
1
>>> 1 * 1 % 7
1
>>> pow(2, (7 + 1) // 4, 7)
4
>>> 4 * 4 % 7
2
>>> pow(4, (7 + 1) // 4, 7)
2
>>> 2 * 2 % 7
4

>>> pow(1, (11 + 1) // 4, 11)
1
>>> 1 * 1 % 11
1
>>> pow(3, (11 + 1) // 4, 11)
5
>>> 5 * 5 % 11
3
>>> pow(4, (11 + 1) // 4, 11)
9
>>> 9 * 9 % 11
4
>>> pow(5, (11 + 1) // 4, 11)
4
>>> 4 * 4 % 11
5
>>> pow(9, (11 + 1) // 4, 11)
3
>>> 3 * 3 % 11
9

このように \(4k + 3\) 型の素数であれば、公式を使って簡単に解を求めることができます。

\(4k + 1\) 型の素数の場合、mod 8 で \(8l + 1\) と \(8l + 5\) に分けると、前者は Tonelli–Shanks のアルゴリズムで、後者は公式で求めることができます。アルゴリズムの概要を参考 URL 7 より引用します。

\(p \equiv 5 \pmod 8\) の場合

\(x \equiv a^{\frac{p+3}{8}}\)
\(if \ x^2 \not\equiv a\) ならば

\(x \equiv 2^{\frac{p-1}{4}} \times x\)

\(return \ x\)

\(p \equiv 1 \pmod 8\) の場合

d として平方非剰余を選ぶ
\(p - 1 = 2^s \times t\) として s, t を定める、ただし t は奇数
\(A \equiv a^t\)
\(D \equiv d^t\)
\(m = 0\)
\(For \ i \ in \ [0, s)\)

\(if \ (AD^m)^{2^{s-i-1}} \equiv -1\)

\(m = m + 2^i\)

\(x \equiv a^{\frac{t+1}{2}} \times D^{\frac{m}{2}}\)
\( return \ x\)

これらの公式やアルゴリズムについて、M.Hiroi はほとんど理解しておりません (ごめんなさい)。詳しい説明は参考 URL 7, 8, 13 をお読みくださいませ。有益な情報を公開されている作者様に感謝いたします。

SymPy には合同式 \(x^2 \equiv a \pmod n\) を解く関数 sqrt_mod(a, n) が用意されています。それでは実際に試してみましょう。

>>> for a in range(1, 11): print((a, sy.sqrt_mod(a, 11)))
...
(1, 1)
(2, None)
(3, 5)
(4, 2)
(5, 4)
(6, None)
(7, None)
(8, None)
(9, 3)
(10, None)

>>> for a in range(1, 13): print((a, sy.sqrt_mod(a, 13)))
...
(1, 1)
(2, None)
(3, 4)
(4, 2)
(5, None)
(6, None)
(7, None)
(8, None)
(9, 3)
(10, 6)
(11, None)
(12, 5)

>>> for a in range(1, 17): print((a, sy.sqrt_mod(a, 17)))
...
(1, 1)
(2, 6)
(3, None)
(4, 2)
(5, None)
(6, None)
(7, None)
(8, 5)
(9, 3)
(10, None)
(11, None)
(12, None)
(13, 8)
(14, None)
(15, 7)
(16, 4)

●ヤコビ記号

ルジャンドル記号 (a / n) の n は奇素数でなければいけませんが、それを奇数に拡張したものが「ヤコビ記号 (jacobi symbol)」です。ヤコビ記号はルジャンドル記号と同じ記号 (a / n) で表します。ヤコビ記号の定義を示します。

\(n = {p_1}^{e1} \times {p_2}^{e2} \times \cdots \times {p_i}^{ei}\) とすると

\( \left(\dfrac{a}{n}\right) = \left(\dfrac{a}{p_1}\right)^{e1} \times \left(\dfrac{a}{p_2}\right)^{e2} \times \cdots \times \left(\dfrac{a}{p_i}\right)^{ei} \)

上式の左辺がヤコビ記号、右辺がルジャンドル記号です。n が奇素数の場合、ヤコビ記号はルジャンドル記号と同じになります。ヤコビ記号 (a / n) が -1 の場合、ルジャンドル記号と同様に a は平方非剰余ですが、(a / n) が 1 になったからといって、a が平方剰余になるとは限りません。たとえば、n = 15 の場合を考えてみましょう。

>>> for a in range(1, 15): print(a, a*a % 15)
...
1 1
2 4
3 9
4 1
5 10
6 6
7 4
8 4
9 6
10 10
11 1
12 9
13 4
14 1

n = 15 の場合、平方剰余は 1, 4, 6, 9, 10 で、それ以外は平方非剰余になります。たとえば、(2 / 15) を考えてみましょう。ヤコビ記号の定義によると、ルジャンドル記号 (2 / 3) * (2 / 5) の値がヤコブ記号 (2 / 15) の値になります。実際に計算すると次のようになります。

\( \left(\dfrac{2}{15}\right) = \left(\dfrac{2}{3}\right) \times \left(\dfrac{2}{5}\right) = (-1) \times (-1) = 1 \)

このように、ヤコブ記号 (2 / 15) の値は 1 になりますが、2 は平方剰余ではありません。ヤコブ記号 (a / n) が 1 になるからといって、a が平方剰余であるとは確定できないことに注意してください。

●ヤコビ記号の性質

ヤコビ記号はルジャンドル記号と同様の性質や相互法則があります。

\(\gcd(a, n) \ne 0\) のとき (a / n) = 0, とくに (0 / n) = 0
\(\gcd(a, n) = 1\) のとき (a / n) = 1 or -1, とくに (1 / n) = 1
\(a \equiv b \pmod n\) のとき (a / n) = (b / n)

ex) \(16 \equiv 5 \pmod {11}, \ \left(\dfrac{16}{11}\right) = \left(\dfrac{5}{11}\right)\)

ようするに (a / n) = ((a % n) / n) が成り立つ

a, b を整数とすると、次式が成り立つ

\( \left(\dfrac{ab}{p}\right) = \left(\dfrac{a}{p}\right) \times\left(\dfrac{b}{p}\right) \)

[相互法則]
p, q が奇数で互いに素のとき、次式が成立する

\( \left(\dfrac{q}{p}\right) \left(\dfrac{p}{q}\right) = (-1)^{\frac{(p-1)(q-1)}{4}} \)

\((p / q) = (q / p), \ if \ p \equiv 1 \pmod 4 \ or \ q \equiv 1 \pmod 4\)
\((p / q) = - (q / p), \ if \ p \equiv 3 \pmod 4 \ and \ q \equiv 3 \pmod 4\)
p, q は奇数だから、p mod 4, q mod 4 の取りうる値は 1 と 3 だけ
1 の場合:

(4*i + 1 - 1)*(4*j + 1 - 1)/4 = 4*i*j = 偶数, (-1)^偶数 = 1
(4*i + 1 - 1)*(4*j + 3 - 1)/4 = 2*i*(2*j + 1) = 偶数 * 奇数 = 偶数, (-1)^偶数 = 1

3 の場合: (4*i + 3 - 1)*(4*j + 3 - 1)/4 = (2*i + 1)*(2*j + 1) = 奇数, (-1)^奇数 = -1

第一補充法則

\( \left(\dfrac{-1}{p}\right) = (-1)^{\frac{p-1}{2}} \)

\((-1 / p) = 1 \ if \ p \equiv 1 \pmod 4\)
\((-1 / p) = -1 \ if \ p \equiv 3 \pmod 4\)
(4*k + 1 - 1) / 2 = 2*k = 偶数, (-1)^偶数 = 1
(4*k + 3 - 1) / 2 = 2*k + 1= 奇数, (-1)^奇数 = -1

第二補充法則

\( \left(\dfrac{2}{p}\right) = (-1)^{\frac{p^2 - 1}{8}} \)

\((2 / p) = 1 \ if \ p \equiv 1 \pmod 8 \ or \ p \equiv 7 \pmod 8\)
\((2 / p) = -1 \ if \ p \equiv 3 \pmod 8 \ or \ p \equiv 5 \pmod 8\)
(8*k + 1 + 1)*(8*k + 1 - 1) / 8 = 2*k*(4*k + 1) = 偶数, (-1)^偶数 = 1
(8*k + 3 + 1)*(8*k + 3 - 1) / 8 = (2*k + 1)*(4*k + 1) = 奇数, (-1)^奇数 = -1
(8*k + 5 + 1)*(8*k + 5 - 1) / 8 = (4*k + 3)*(2*k + 1) = 奇数, (-1)^奇数 = -1
(8*k + 7 + 1)*(8*k + 7 - 1) / 8 = 2*k*(4*k + 3) = 偶数, (-1)^偶数 = 1

●ヤコビ記号の計算

これらの性質を使うと、素因数分解しなくてもヤコビ記号を計算することができます。次の例を見てください。

(1001 / 12345) = 1 * (12345 / 1001)  [相互法則] 
               = (333 / 1001)        [(a / n) = ((a % n) / n]
               = 1 * (1001 / 333)    [相互法則]
               = (2 / 333)           [(a / n) = ((a % n) / n]
               = (2 / 333)*(1 / 333) [(a*b / n) = (a / n)*(b / n)
               = (-1) * (1 / 333)    [第 2 補充法則]
               = -1

ヤコビ記号 (1001 / 12345) を求めます。便宜上、(a / n) の a を分子、n を分母と呼ぶことにします。分母 n は奇数でなので、分子 a が奇数であれば、相互法則を適用することができます。1001 % 4 == 1 or 12345 % 4 == 1 より 1 * (123345 / 1001) になるので、分子 % 分母を計算して、ヤコビ記号は (333 / 1001) になります。このように、分子が奇数であれば相互法則と (a / n) = ((a % n) / n) を適用することで、分子と分母の値を小さくしていくことができます。

同様に、(333 / 1001) は (2 / 333) になります。分子が偶数の場合、((a * b) / n) = (a / n) * (b / n) を適用して、式を (2 / n) *((a / 2) / n) に変形し、そこに第 2 補充法則を適用します。333 % 8 == 5 より、(2 / 333) * (1 / 333) = (-1) * (1 / 333) になります。(1 / 333) の値は 1 なので、ヤコビ記号 (1001 / 12345) の値は -1 となります。

SymPy にはヤコビ記号を求める関数 jacobi_symbol(a, n) が用意されています。それでは実際に試してみましょう。

>>> for a in range(1, 3): print(sy.jacobi_symbol(a, 3), end=" ")
...
1 -1
>>> for a in range(1, 5): print(sy.jacobi_symbol(a, 5), end=" ")
...
1 -1 -1 1
>>> for a in range(1, 7): print(sy.jacobi_symbol(a, 7), end=" ")
...
1 1 -1 1 -1 -1
>>> for a in range(1, 9): print(sy.jacobi_symbol(a, 9), end=" ")
...
1 1 0 1 1 0 1 1
>>> for a in range(1, 11): print(sy.jacobi_symbol(a, 11), end=" ")
...
1 -1 1 1 1 -1 -1 -1 1 -1
>>> for a in range(1, 13): print(sy.jacobi_symbol(a, 13), end=" ")
...
1 -1 1 1 -1 -1 -1 -1 1 1 -1 1
>>> for a in range(1, 15): print(sy.jacobi_symbol(a, 15), end=" ")
...
1 1 0 1 0 0 -1 1 0 0 -1 0 -1 -1

n が奇素数の場合、jacobi_symbol() の値はルジャンドル記号と同じ値になっていますね。n が奇素数ではない場合、gcd(a, n) != 1 のとき jacobi_symbol() の値は 0 になります。

●参考 URL

原始根の定義と具体例（高校生向け）, (マスオさん)
位数の性質と原始根の応用, (マスオさん)
平方剰余と基本的な問題, (マスオさん)
ルジャンドル記号とオイラーの規準, (マスオさん)
平方剰余の相互法則の意味と応用, (マスオさん)
【internal_math編②】AtCoder Library 解読〜Pythonでの実装まで〜, (Akari さん)
mod pでの平方根, (Peria さん)
平方剰余の根, (kenken さん)
ウィルソンの定理 - Wikipedia
ルジャンドル記号 - Wikipedia
ヤコビ記号 - Wikipedia
平方剰余 - Wikipedia
Tonelli–Shanks algorithm - Wikipedia, (en)

Memorandum

数学編