Algorithms with Python

はじめに

今回は位数と原始根について簡単に説明し、位数と原始根を求めるプログラムを作ってみましょう。

●位数と原始根

互いに素となる整数 a と n に対して、合同式 \(a^k \equiv 1 \pmod n\) を満たす最小の正整数 k を、n を法とする a の位数 (order of a modulo n) といい、\(ord_n(a)\) や \(O_n(a)\) と記述します。オイラーの定理 \(a^{\varphi(n)} \equiv 1 \pmod n\) より、位数 k は \(\varphi(n)\) 以下になります。そして、\(\varphi(n)\) と位数が等しいとき、a を n の原始根 (primitive root of n) と呼びます。

簡単な例を示しましょう。2 ** k % 11 と 3 ** k % 11 (k = 1, 2, ..., 10) を計算します。

>>> for k in range(1, 11): print(2 ** k % 11, end=" ")
...
2 4 8 5 10 9 7 3 6 1
>>> for k in range(1, 11): print(3 ** k % 11, end=" ")
...
3 9 5 4 1 3 9 5 4 1

mod 11 の世界では、2 の位数は 10 で、3 の位数は 5 になりました。n が素数のとき \(\varphi(n) = n - 1\) になるので、\(\varphi(11)\) は 10 になります。したがって、2 は 11 の原始根ですが、3 は 11 の原始根ではありません。それから、2 ** k % 11 (k = 1, ..., 10) の計算結果を見ると、1 から 10 までの値がひとつずつ現れていますね。これは原始根で一般的に成立する性質 (定理) です。

• 定理１
• a が n の原始根であるとき、\(a, a^2, \cdots, a^{n-1}\) を n で割った余りはすべて異なり、1 から n - 1 までの値がひとつずつ現れる

• 証明 (背理法)
• \(a^i \equiv a^j \pmod n\) となる自然数 \(i, j \ (n - 1 \geq i \gt j)\) が存在すると仮定する
• 式を変形する
\(\begin{array}{l} a^i \equiv a^j \\ a^i - a^j \equiv 0 \\ a^j \times (a^{i-j} - 1) \equiv 0 \end{array}\)
• a と n は互いに素なので、a^j で除算することができる
\(\begin{array}{l} a^{i-j} - 1 \equiv 0 \pmod n \\ a^{i-j} \equiv 1 \pmod n \end{array}\)
• n - 1 よりも小さな整数 i - j が位数となるため、原始根の定義に反する
• \(a^i \equiv a^j \pmod n\) となる自然数 i, j は存在しない (n で割った余りはすべて異なる)

n が素数のとき、n の原始根 a が必ず存在することが知られています。これを「原始根の存在定理」とか「原始根定理」といいます。n の原始根が存在するとき、その個数は \(\varphi(n - 1)\) になります。

それから、位数には次の性質 (定理) があります。

• 定理２
• mod n における a の位数を s とする
• \(a^k \equiv 1 \pmod n\) を満たすならば、k は s の倍数になる

たとえば、mod 11 における 3 の位数は 5 ですが、3^k % 11 (k = 1, 2, ..., 10) を計算すると、k = 5 * 2 = 10 のときも 1 になります。

• 証明
• \(k = p \times s + q\) とおく
\(\begin{array}{l} a^k \equiv 1 \\ a^{p \times s + q} \equiv 1 \\ (a^s)^p \times a^q \equiv 1 \end{array}\)
• s は位数なので、\(a^s \equiv 1 \pmod n\) より \(a^q \equiv 1 \pmod n\) となる
• q は余りなので \(0 \leq q \lt s\) だが、\(q \ne 0\) とすると位数の定義に反する
• よって \(q = 0\) となり、k は s の倍数 (\(k = p \times s\)) であることが示された

●ウィルソンの定理

位数と原始根を使うと、ウィルソンの定理 (Wilson's theorem) を簡単に証明することができます。

• ウィルソンの定理
• \(n\) が素数ならば、\((n - 1)! \equiv -1 \pmod n\) が成立する
• 逆に、整数 \(n \gt 1\) に対し、\((n - 1)! \equiv -1 \pmod n\) ならば、\(n\) は素数である

実際に試してみましょう。

>>> for x in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]:
...     print(x, fact(x - 1) % x)
...
2 1
3 2
5 4
7 6
11 10
13 12
17 16
19 18
23 22
29 28
31 30
37 36
41 40
43 42
47 46

確かに (n - 1)! % n は n - 1 (-1 % n = n - 1) になります。n = 2 の場合、定理が成り立つのは自明なので、奇素数の場合について証明します。

• 証明
• 素数 n の原始根を a とすると以下の式が成り立つ
\(\begin{array}{l} a \times a^2 \times \cdots \times a^{n-1} \equiv 1 \times 2 \times \cdots \times (n - 1) \pmod n \\ a^m \equiv (n - 1)! \pmod n, \quad \left(m = 1 + 2 + \cdots + (n - 1) = \dfrac{n(n - 1)}{2}\right) \end{array}\)
• \(a^{2m}\) にフェルマーの小定理 \(a^{n-1} \equiv 1 \pmod n\) を適用する
\(\begin{array}{l} a^{2m} = a^{n(n-1)} = (a^{n-1})^n \\ a^{2m} \equiv 1^n \pmod n \\ (a^m + 1) \times (a^m - 1) \equiv 0 \pmod n \end{array}\)
• \(a^m \equiv 1\) と仮定すると、m は n - 1 の倍数になる (m は n - 1 で割り切れる)
\( \dfrac{m}{p - 1} = \dfrac{n(n - 1)}{2(n - 1)} = \dfrac{n}{2} \)
• n は奇素数だから n / 2 は整数にならない、よって \(a^m \equiv -1 \pmod n\)
• \(a^m \equiv (n - 1)!\) より \((n - 1)! \equiv -1 \pmod n\) が成り立つ
• 逆に、\(n\) を合成数とする
• \(n\) を割り切る素数 \(q\) (\(2 \leq q \lt n)\) が存在するので \((n - 1)! \equiv 0 \pmod q\) になる
• \((n - 1)! \equiv -1 \pmod n\) であれば \((n - 1)! \equiv -1 \pmod q\) になるので矛盾する (\(n\) は素数である)

●プログラムの作成

位数を求めるプログラムは、多少遅くてもよければ totient() と divisors() を使って簡単に作ることができます。次のリストを見てください。

リスト : 位数を求める

def n_order(a, n):
    if math.gcd(a, n) == 1:
        for d in divisors(totient(n)):
            if pow(a, d, n) == 1: return d
    return None

関数名 n_order は SymPy より拝借しました。引数 a と n が互いに素でなければ None を返します。あとは、totient(n) で \(\varphi(n)\) を求め、divisors() でその約数を求めます。divsors() は昇順にソートした結果を返すので、先頭から pow(a, d, n) を計算して 1 になる約数 d が位数となります。

簡単な実行例を示します。

>>> n_order(2, 11)
10
>>> n_order(3, 11)
5
>>> n_order(2, 13)
12
>>> n_order(3, 13)
3
>>> n_order(4, 13)
6
>>> n_order(5, 13)
4
>>> n_order(6, 13)
12
>>> n_order(7, 13)
12

a が n の原始根か判定する関数 is_primitive_root() も簡単です。次のリストを見てください。

リスト :  原始根の判定

def is_primitive_root(a, n):
    if math.gcd(a, n) == 1:
        phi = totient(n)
        if phi == n - 1:
            for d in divisors(phi):
                if pow(a, d, n) == 1: return phi == d
    return False

関数名は SymPy より拝借しました。totient() で \(\varphi(n)\) の値を求め、変数 phi にセットします。phi が n - 1 と等しくなければ原始根ではありません。return で False を返します。n - 1 と等しい場合、phi の約数 d を取り出して、a^d % n を計算します。その値が 1 であれば d は位数になります。d が phi と等しければ a は原始根なので True を返し、そうでなければ False を返します。

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

>>> is_primitive_root(2, 11)
True
>>> is_primitive_root(3, 11)
False
>>> is_primitive_root(2, 13)
True
>>> is_primitive_root(3, 13)
False
>>> is_primitive_root(4, 13)
False
>>> is_primitive_root(5, 13)
False
>>> is_primitive_root(6, 13)
True
>>> is_primitive_root(7, 13)
True

●プログラムの高速化

位数を求める n_order() や原始根を判定する is_primitive_root() は、次の定理を使うともっと速くすることができます。

• 定理３
• a が n の原始根で n - 1 の素因数分解が \({p_1}^a \times {p_2}^b \times \cdots \times {p_m}^k\) と表されるとき、以下の式が成り立つ
\( a^{\frac{n-1}{p_i}} \not\equiv 1 \pmod n, \quad (i = 1, 2, \cdots, m) \)

• 証明
• n - 1 は \(p^i\) で割り切れるので、\(1 \leq \frac{n - 1}{p_i} \lt n - 1\) が成立する
• a は n の原始根なので、\(a^k \equiv 1 \pmod n\) を満たす最小の \(k\) は \(n - 1\)
• \(\frac{n - 1}{p_i}\) は \(n - 1\) よりも小さいので、\(a^{\frac{n - 1}{p_i}} \equiv 1 \pmod n\) にはならない

• 定理３の対偶
• a が n の原始根ではなく、n - 1 の素因数分解が \({p_1}^a \times {p_2}^b \times \cdots \times {p_m}^k\) と表されるとき、以下の式が成り立つ \(p_i\) が存在する
\( a^{\frac{n-1}{p_i}} \equiv 1 \pmod n \)

• 証明
• a が n の原始根でなければ、位数 s は n - 1 未満で、n - 1 の約数になる (n - 1 は s の倍数)
• 位数 s は n - 1 と同様に \(p_1, \ p_2, \ \cdots, \ p_m\) で素因数分解できる
• ただし、指数は a, b, ..., k よりも小さいものが少なくとも一つある
\( s = {p_1}^{a'} \times {p_2}^{b'} \times \cdots \times {p_m}^{k'}, \quad (a' \leq a, \ b' \leq b, \cdots, \ k' \leq k) \)
• それを \({p_i}^{e'}, \ (e' \lt e)\) とすると、次式が成り立つ
\(\begin{eqnarray} \dfrac{n-1}{p_i} &=& \left(\dfrac{1}{p_i}\right) \times {p_1}^a \times {p_2}^b \times \cdots \times {p_i}^e \times \cdots \times {p_m}^k \\ &=& \left(\dfrac{1}{p_i}\right) \times {p_1}^{a-a'} \times {p_2}^{b-b'} \times \cdots \times {p_i}^{e-e'} \times \cdots \times {p_m}^{k-k'} \times s \\ &=& {p_i}^{e-e'-1} \times s \end{eqnarray}\)

\(e - e' - 1 \geq 0\) より \(\dfrac{n - 1}{p_i} = (数) * s\)

\(\begin{eqnarray} a^{\frac{n-1}{p_i}} &=& a^{(数) \times s} \\ &=& (a^s)^{(数)} \end{eqnarray}\)
• \(a^s \equiv 1 \pmod n\) より \(a^{\frac{n-1}{p_i}} \equiv 1 \pmod n\) が成り立つ

is_primitive_root() は次のようになります。

リスト : 原始根の判定 (2)

def is_primitive_root1(a, n):
    if math.gcd(a, n) != 1: return False
    for p, _ in factorization(n - 1):
        if pow(a, (n - 1) // p, n) == 1: return False
    return True

factorization() で n - 1 を素因数分解し、for ループで因数 p を順番に取り出します。pow(a, (n - 1) // p, n) が 1 になる場合、a は原始根ではありません。return で False を返します。すべての因数で pow() が 1 にならなければ a は原始根です。最後に True を返します。

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

>>> for x in range(2, 21): print(is_primitive_root(x, 2**31 - 1), is_primitive_root1(x, 2**31 - 1))
...
False False
False False
False False
False False
False False
True True
False False
False False
False False
True True
False False
False False
True True
False False
False False
False False
False False
False False
False False

>>> s = time.time(); print([is_primitive_root(x, 2**31 - 1) for x in range(2, 30)]); time.time() - s
[False, False, False, False, False, True, False, False, False, True, False, False, True, False, False, 
False, False, False, False, False, True, False, False, False, False, False, True, False]
0.15649771690368652

>>> s = time.time(); print([is_primitive_root1(x, 2**31 - 1) for x in range(2, 30)]); time.time() - s
[False, False, False, False, False, True, False, False, False, True, False, False, True, False, False, 
False, False, False, False, False, True, False, False, False, False, False, True, False]
0.0010783672332763672

2 ** 31 - 1 (= 2147483647) はメルセンヌ素数です。結果を見ればおわかりのように、is_primitive_root1() のほうが高速になりました。

次は位数を求める n_order() の高速版を作ります。

リスト : 位数を求める (2)

def n_order1(a, n):
    if math.gcd(a, n) != 1: return None
    d = n - 1
    xs = factorization(d)
    while True:
        for p, _ in xs:
            if d % p == 0:
                d1 = d // p
                if pow(a, d1, n) == 1:
                    d = d1
                    break
        else:
            return d

変数 d を n - 1 に初期化して、factorization() で素因数分解します。d // p を d1 にセットし、pow(a, d1, n) を計算します。結果が 1 になる場合、d よりも小さい d1 が位数になる可能性があります。d1 を d にセットし、d 未満の約数をチェックします。d を割り切る因数 p がないか、pow(a, d1, n) が 1 になる因数がない場合、d が位数になります。for ループの else 節で d を返します。

それでは実行していましょう。

>>> for x in range(2, 21): print(n_order(x, 2**31 - 1), n_order1(x, 2**31 - 1))
...
31 31
715827882 715827882
31 31
195225786 195225786
715827882 715827882
2147483646 2147483646
31 31
357913941 357913941
195225786 195225786
2147483646 2147483646
715827882 715827882
1073741823 1073741823
2147483646 2147483646
1073741823 1073741823
31 31
1073741823 1073741823
357913941 357913941
34636833 34636833
195225786 195225786

>>> s = time.time(); print([n_order(x, 2**31 - 1) for x in range(2, 30)]); time.time() - s
[31, 715827882, 31, 195225786, 715827882, 2147483646, 31, 357913941, 195225786, 2147483646, 
715827882, 1073741823, 2147483646, 1073741823, 31, 1073741823, 357913941, 34636833, 195225786, 
153391689, 2147483646, 1073741823, 715827882, 97612893, 1073741823, 238609294, 2147483646, 357913941]
0.1576385498046875

>>> s = time.time(); print([n_order1(x, 2**31 - 1) for x in range(2, 30)]); time.time() - s
[31, 715827882, 31, 195225786, 715827882, 2147483646, 31, 357913941, 195225786, 2147483646, 
715827882, 1073741823, 2147483646, 1073741823, 31, 1073741823, 357913941, 34636833, 195225786, 
153391689, 2147483646, 1073741823, 715827882, 97612893, 1073741823, 238609294, 2147483646, 357913941]
0.0017619132995605469

結果を見ればおわかりのように、n_order1() のほうが高速になりました。

●原始根を求める

定理３を使うと、mod n の原始根を求める関数 primitive_root() を簡単に作ることができます。

リスト :  原始根を求める (n は素数であること)

def primitive_root(n):
    xs = factorization(n - 1)
    for a in range(2, n):
        for p, _ in xs:
            if pow(a, (n - 1) // p, n) == 1: break
        else:
            yield a

primitive_root() は mod n の原始根を返すジェネレータを生成します。n - 1 を素因数分解して、結果を変数 xs にセットします。あとは for ループで a を +1 しながら、a が原始根かチェックするだけです。2 番目の for ループで、pow(a, (n - 1) // p, n) が 1 ならば、a は原始根ではありません。break で for ループを脱出します。for ループを途中で break しなければ、for ループの else 節が実行されるので、yield a で a を出力します。

簡単な実行例を示します。

>>> for x in [3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]:
...     print(x, list(primitive_root(x)))
...
3 [2]
5 [2, 3]
7 [3, 5]
11 [2, 6, 7, 8]
13 [2, 6, 7, 11]
17 [3, 5, 6, 7, 10, 11, 12, 14]
19 [2, 3, 10, 13, 14, 15]
23 [5, 7, 10, 11, 14, 15, 17, 19, 20, 21]
29 [2, 3, 8, 10, 11, 14, 15, 18, 19, 21, 26, 27]
31 [3, 11, 12, 13, 17, 21, 22, 24]
37 [2, 5, 13, 15, 17, 18, 19, 20, 22, 24, 32, 35]
41 [6, 7, 11, 12, 13, 15, 17, 19, 22, 24, 26, 28, 29, 30, 34, 35]
43 [3, 5, 12, 18, 19, 20, 26, 28, 29, 30, 33, 34]
47 [5, 10, 11, 13, 15, 19, 20, 22, 23, 26, 29, 30, 31, 33, 35, 38, 39, 40, 41, 43, 44, 45]

>>> for x in [3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]:
...     print(x, totient(x - 1))
...
3 1
5 2
7 2
11 4
13 4
17 8
19 6
23 10
29 12
31 8
37 12
41 16
43 12
47 22

●平方剰余

合同式 \(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 が平方非剰余になることが分かります。

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

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

ルジャンドル記号は (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)\) となる

ルジャンドル記号を Python でプログラムすると次のようになります。

リスト : ルジャンドル記号

def legendre_symbol(a, n):
    if not miller_rabin_test(n) or n == 2:
        raise ValueError("n should be an odd prime")
    if a % n == 0:
        return 0
    elif pow(a, (n - 1) // 2, n) == 1:
        return 1
    return -1

miller_rabin_test() は 素数 で作成した素数を判定 (ミラーラビン法) する関数です。引数 n が奇素数でなければエラーを送出します。引数 a が n で割り切れる場合は 0 を返します。次に、pow(a, (n - 1) // 2, n) を計算し、それが 1 ならば a は平方剰余なので 1 を返します。そうでなければ平方非剰余なので -1 を返します。

簡単な実行例を示します。

>>> for a in range(1, 20): print(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(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(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(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(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)\) 型の奇素数ならば、次の公式で求めることができます。

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

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

>>> 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\) の場合
1. \(x \equiv a^{\frac{p+3}{8}}\)
2. \(if \ x^2 \not\equiv a\) ならば
1. \(x \equiv 2^{\frac{p-1}{4}} \times x\)
3. \(return \ x\)
• \(p \equiv 1 \pmod 8\) の場合
1. d として平方非剰余を選ぶ
2. \(p - 1 = 2^s \times t\) として s, t を定める、ただし t は奇数
3. \(A \equiv a^t\)
4. \(D \equiv d^t\)
5. \(m = 0\)
6. \(For \ i \ in \ [0, s)\)
1. \(if \ (AD^m)^{2^{s-i-1}} \equiv -1\)
1. \(m = m + 2^i\)
7. \(x \equiv a^{\frac{t+1}{2}} \times D^{\frac{m}{2}}\)
8. \( return \ x\)

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

これをそのままプログラムすると次のようになります。

リスト :  平方剰余を求める

def modsqrt(a, n):
    if legendre_symbol(a, n) == -1: return None
    if n % 4 == 3:
        return pow(a, (n + 1) // 4, n)
    elif n % 8 == 5:
        x = pow(a, (n + 3) // 8, n)
        if pow(x, 2, n) != a:
            x = x * pow(2, (n - 1) // 4, n) % n
        return x
    else:
        # Tonelli–Shanks のアルゴリズム
        # 平方非剰余を選ぶ
        for i in range(2, n):
            if legendre_symbol(i, n) == -1:
                d = i
                break
        s, t = factor2(n - 1)  # n - 1 = 2**s * t
        a1 = pow(a, t, n)
        d1 = pow(d, t, n)
        m = 0
        for i in range(s):
            if pow(a1 * (d1 ** m), 2 ** (s - i - 1), n) == n - 1:
                m += 2 ** i
        return (pow(a, (t + 1) // 2, n) * pow(d1, m // 2, n)) % n

関数 modsqrt(a, n) は合同式 x² ≡ a (mod n) の解 x を求めます。最初に legendre_symbol() を呼び出して、a が平方剰余であることを確認します。解が無い場合は None を返します。あとはアルゴリズムをそのままプログラムしただけなので、難しいところはないと思います。

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

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

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

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

>>> 97 % 8
1
>>> for a in range(1, 97): print((a, modsqrt(a, 97)), end=" ")
...
(1, 1) (2, 83) (3, 87) (4, 95) (5, None) (6, 54) (7, None) (8, 69) (9, 94) (10, None) (11, 37) 
(12, 77) (13, None) (14, None) (15, None) (16, 93) (17, None) (18, 55) (19, None) (20, None) 
(21, None) (22, 33) (23, None) (24, 11) (25, 92) (26, None) (27, 67) (28, None) (29, None) 
(30, None) (31, 82) (32, 41) (33, 79) (34, None) (35, 61) (36, 91) (37, None) (38, None) (39, None) 
(40, None) (41, None) (42, None) (43, 72) (44, 74) (45, None) (46, None) (47, 85) (48, 57) (49, 7) 
(50, 27) (51, None) (52, None) (53, 76) (54, 65) (55, None) (56, None) (57, None) (58, None) 
(59, None) (60, None) (61, 35) (62, 81) (63, None) (64, 89) (65, 68) (66, 39) (67, None) (68, None) 
(69, None) (70, 19) (71, None) (72, 13) (73, 49) (74, None) (75, 47) (76, None) (77, None) (78, None) 
(79, 51) (80, None) (81, 88) (82, None) (83, None) (84, None) (85, 52) (86, 38) (87, None) (88, 66) 
(89, 34) (90, None) (91, 73) (92, None) (93, 44) (94, 26) (95, 17) (96, 22)

正常に動作しているようです。興味のある方はいろいろ試してみてください。

●ヤコビ記号

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

上式の左辺がヤコビ記号、右辺がルジャンドル記号です。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) の値になります。実際に計算すると次のようになります。

このように、ヤコブ記号 (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 となります。

これをプログラムすると次のようになります。

リスト : ヤコビ記号

def jacobi_sym(a, n):
    if a == 0: return 0
    elif a == 1: return 1
    elif a % 2 == 0:
        k = -1 if n % 8 in [3, 5] else 1
        return k * jacobi_sym(a // 2, n)
    else:
        k = -1 if a % 4 == 3 and n % 4 == 3 else 1
        return k * jacobi_sym(n % a, a)

def jacobi_symbol0(a, n):
    if n <= 0 or n % 2 == 0:
        raise ValueError("n should be an odd positive integer")
    return jacobi_sym(a % n, n)

関数 jacobi_symbol0(a, n) はヤコビ記号 (a / n) の値を求めます。引数 n が正の奇数でなければエラーを送出します。実際の処理は関数 jacobi_sym() で行います。a > n の場合もあるので、jacobi_sym() の第 1 引数には a % n を渡します。

jacobi_sym(a, n) は簡単です。最初に (0 / n) と (1 / n) をチェックします。a == 0 であれば 0 を、a == 1 であれば 1 を返します。次に a が偶数の場合、第 2 補充法則を適用します。符号を変数 k に求め、jacobi_sym(a // 2, n) の返り値に k を乗算して返します。

最後は相互法則を適用します。符号を変数 k に求め、jacobi_sym(n % a, a) の返り値に k を乗算して返します。jacobi_sym() を再帰呼び出しするとき、a と n を逆に渡していることに注意してください。

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

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

n が奇素数の場合、jacobi_symbol0() の値はルジャンドル記号と同じ値になっていますね。n が奇素数ではない場合、gcd(a, n) != 1 のとき jacobi_symbol0() の値は 0 になります。プログラムは正常に動作しているようです。興味のある方はいろいろ試してみてください。

このように、再帰定義を使うと簡単にプログラムできますが、実用的には再帰定義を繰り返しに変換したほうがよいでしょう。プログラムは次のようになります。参考 URL 11 のプログラムを参考にさせていただきました。作者様に感謝いたします。

リスト : ヤコビ記号 (繰り返し版)

def jacobi_symbol(a, n):
    if n <= 0 or n % 2 == 0:
        raise ValueError("n should be an odd positive integer")
    a %= n
    k = 1
    while a > 0:
        # 第 2 補充法則
        while a % 2 == 0:
            a //= 2
            if n % 8 in [3, 5]:
                k = -k
        # 相互法則
        a, n = n, a
        if a % 4 == 3 and n % 4 == 3:
            k = -k
        a %= n
    return k if n == 1 else 0

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

●参考 URL

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

●プログラムリスト

#
# modulo.py : 位数と原始根
#
#             Copyright (C) 2023 Makoto Hiroi
#
import math

#
# ミラーラビンテスト
#

# n を (2 ** c) * d に分解する
def factor2(n):
    c = 0
    while n % 2 == 0:
        n //= 2
        c += 1
    return c, n

def miller_rabin_sub(n, s, d, ks):
    for a in ks:
        if a >= n: break
        if pow(a, d, n) != 1:
            r = 0
            while r < s:
                if pow(a, (2**r) * d, n) == n - 1: break
                r += 1
            else:
                return False
    return True

def miller_rabin_test(n, k = 20):
    if n < 2: return False
    elif n == 2: return True
    elif n % 2 == 0: return False
    else:
        s, d = factor2(n - 1)
        if n < 4759123141:
            ks = [2, 7, 61]
        elif n <= 0x10000000000000000:
            ks = [2, 325, 9375, 28178, 450775, 9780504, 1795265022]
        else:
            ks = [random.randrange(1, n) for _ in range(k)]
        return miller_rabin_sub(n, s, d, ks)

# 素因数分解
def factor_sub(n, m):
    c = 0
    while n % m == 0:
        c += 1
        n //= m
    return c, n

# wheel factorization
def factorization(n):
    buff = []
    wheel = [1, 2, 2, 4, 2, 4, 2, 4, 6, 2, 6]
    s = 3
    i, x, m = 0, 2, n
    while m >= x * x:
        c, m = factor_sub(m, x)
        if c > 0: buff.append((x, c))
        x += wheel[i]
        i += 1
        if i >= len(wheel): i = s
    if m > 1: buff.append((m, 1))
    return buff

# 約数
def divisor_sub(p, q):
    return [p ** i for i in range(q + 1)]

def divisors(n):
    xs = factorization(n)
    ys = divisor_sub(xs[0][0], xs[0][1])
    for p, q in xs[1:]:
        ys = [x * y for x in divisor_sub(p, q) for y in ys]
    return sorted(ys)

# オイラーのトーシェント関数
def totient(n):
    d = n
    for p, _ in factorization(n):
        d = d * (p - 1) // p
    return d

# 位数を求める
def n_order(a, n):
    if math.gcd(a, n) == 1:
        for d in divisors(totient(n)):
            if pow(a, d, n) == 1: return d
    return None

# 高速化
def n_order1(a, n):
    if math.gcd(a, n) != 1: return None
    d = n - 1
    xs = factorization(d)
    while True:
        for p, _ in xs:
            if d % p == 0:
                d1 = d // p
                if pow(a, d1, n) == 1:
                    d = d1
                    break
        else:
            return d

# 原始根の判定
def is_primitive_root(a, n):
    if math.gcd(a, n) == 1:
        phi = totient(n)
        if phi == n - 1:
            for d in divisors(phi):
                if pow(a, d, n) == 1: return phi == d
    return False

# 高速化
def is_primitive_root1(a, n):
    if math.gcd(a, n) != 1: return False
    for p, _ in factorization(n - 1):
        if pow(a, (n - 1) // p, n) == 1: return False
    return True

# 原始根を求める (n は素数であること)
def primitive_root(n):
    xs = factorization(n - 1)
    for a in range(2, n):
        for p, _ in xs:
            if pow(a, (n - 1) // p, n) == 1: break
        else:
            yield a

# ルジャンドル記号
def legendre_symbol(a, n):
    if not miller_rabin_test(n) or n == 2:
        raise ValueError("n should be an odd prime")
    if a % n == 0:
        return 0
    elif pow(a, (n - 1) // 2, n) == 1:
        return 1
    return -1

# ヤコビ記号
def jacobi_sym(a, n):
    if a == 0: return 0
    elif a == 1: return 1
    elif a % 2 == 0:
        k = -1 if n % 8 in [3, 5] else 1
        return k * jacobi_sym(a // 2, n)
    else:
        k = -1 if a % 4 == 3 and n % 4 == 3 else 1
        return k * jacobi_sym(n % a, a)

def jacobi_symbol0(a, n):
    if n <= 0 or n % 2 == 0:
        raise ValueError("n should be an odd positive integer")
    return jacobi_sym(a % n, n)

# 繰り返し版
def jacobi_symbol(a, n):
    if n <= 0 or n % 2 == 0:
        raise ValueError("n should be an odd positive integer")
    a %= n
    k = 1
    while a > 0:
        # 第２法則
        while a % 2 == 0:
            a //= 2
            if n % 8 in [3, 5]:
                k = -k
        # 相互法則
        a, n = n, a
        if a % 4 == 3 and n % 4 == 3:
            k = -k
        a %= n
    return k if n == 1 else 0

# 合同式を解く
def modsqrt(a, n):
    if legendre_symbol(a, n) == -1: return None
    if n % 4 == 3:
        return pow(a, (n + 1) // 4, n)
    elif n % 8 == 5:
        x = pow(a, (n + 3) // 8, n)
        if pow(x, 2, n) != a:
            x = x * pow(2, (n - 1) // 4, n) % n
        return x
    else:
        # Tonelli–Shanks のアルゴリズム
        # 平方非剰余を選ぶ
        for i in range(2, n):
            if legendre_symbol(i, n) == -1:
                d = i
                break
        s, t = factor2(n - 1)  # n - 1 = 2**s * t
        a1 = pow(a, t, n)
        d1 = pow(d, t, n)
        m = 0
        for i in range(0, s):
            if pow(a1 * (d1 ** m), 2 ** (s - i - 1), n) == n - 1:
                m += 2 ** i
        return (pow(a, (t + 1) // 2, n) * pow(d1, m // 2, n)) % n