0. はじめに

2020年9月7日にAtCoder公式のアルゴリズム集 AtCoder Library (ACL)が公開されました。
私はACLに収録されているアルゴリズムのほとんどが初見だったのでいい機会だと思い、アルゴリズムの勉強からPythonでの実装までを行いました。

この記事では internal_math をみていきます。

internal_mathはACLの内部で使われる数論的アルゴリズムの詰め合わせで内容は以下の通りです。

名称	概要
safe_mod	整数 $x$ の正整数 $m$ による剰余($x \% m$)。ただし $0 \leq x \% m < m $ を満たす。
barrett	高速な剰余演算。
pow_mod	$x^n \pmod{m}$ の計算。
is_prime	高速な素数判定。
inv_gcd	整数 $a$ と正整数 $b$ の最大公約数 $g$ および $xa \equiv g \pmod{b}$ となる $x$ の計算。ただし $0 \leq x < \frac{b}{g} $ を満たす。
primitive_root	素数 $m$ の原始根。
floor_sum_unsigned	$0$ 以上の整数 $n, a, b$ と自然数 $m$ について $\sum_{i=0}^{n-1}{\left\lfloor\frac{ai+b}{m}\right\rfloor}$ の計算。

本記事ではこれらの内、

is_prime
primitive_root

を扱います。なお、constexpr(定数式)自体については触れません。

本記事で扱わない

safe_mod
pow_mod
inv_gcd
barrett

については以下の記事で扱っています。よろしければそちらもご覧ください。

【internal_math編①】AtCoder Library 解読〜Pythonでの実装まで〜

なお、floor_sum_unsigned は math 内の floor_sum の補助的な役割なので説明は

【math編】AtCoder Library 解読〜Pythonでの実装まで〜

で行なっています。

対象としている読者

ACLの内部で使われている数学系アルゴリズムを知りたい方。
ACLのコードを見てみたけど何をしているのかわからない方。
C++はわからないのでPythonで読み進めたい方。

参考にしたもの

is_primeに関連するwikipediaのページです。

強い擬素数についての論文です。

GERHARD JAESCHKE (1993) "ON STRONG PSEUDOPRIMES TO SEVERAL BASES"

ミラー-ラビン素数判定法について書かれている@srtk86さんの記事です。わかりやすいです。

素数判定いろいろ - フェルマーテスト・ミラーラビン素数判定法

原始根について説明されています。

Wikibooksへの寄稿者ら, "初等整数論/原始根と指数," Wikibooks, フリー教科書

1. is_prime

正整数 $n$ が素数であるか判定することを考えます。

1.1. 決定的素数判定法

素数の定義は、「正の約数が1と自分自身のみである、1より大きい自然数」ですから、２から $n - 1$ までの自然数で $n$ を割って、割り切れなかったら $n$ は素数だといえます。いわゆる試し割りという方法です。Pythonで実装すると以下のようになります。

def deterministicPT(n):
    if n <= 1: return False
    if n == 2: return True
    for i in range(2, n):
        if n % i == 0: return False
    return True

これを使って素数判定を行うと以下のようになります。

for n in range(1, 10):
    if deterministicPT(n):
        print(f'{n}は素数である')
    else:
        print(f'{n}は素数でない')

# 1は素数でない
# 2は素数である
# 3は素数である
# 4は素数でない
# 5は素数である
# 6は素数でない
# 7は素数である
# 8は素数でない
# 9は素数でない

この方法では定義に忠実に確かめたので、$n$ が素数であるか否かを確実に判定することができます。このような方法を決定的素数判定法といいます。
決定的素数判定法は言い換えれば次のような性質を持つテストを行うことです。

素数ならば必ず合格になる
素数でないならば必ず不合格になる

1.2. 確率的素数判定法

「素数である」「素数でない」のどちらかを判定する決定的素数判定法に対して、「素数かもしれない」「素数でない」のどちらかを判定するアルゴリズムを確率的素数判定法と言います。
確率的素数判定法に用いられるテストは次の性質を持っています。

素数ならば必ず合格する
素数でないならば合格不合格のどちらにもなりうる
自然数 $a$ が組み込まれている

そして、このようなテストに合格する自然数を「底 $a$ の確率的素数」と言います。

具体例を見てみましょう。判定したい自然数 $n$ に対して次のようなテストを考えます。

$n=1$ ならば不合格
$2 \leq n \leq a$ ならば合格
$n > a$ ならば $a$ で割り切れなければ合格、割り切れたら不合格

このテストは上で挙げた３つの性質を満たしているので確率的素数判定法に用いることができます。Pythonで実装すると次のようになるでしょう。

class ProbabilisticPT:
    def __init__(self, a=2):
        self._a = a
    
    def change_base(self, a):
        self._a = a
    
    def test(self, n):
        if n == 1: return False
        if n <= self._a: return True
        if n % self._a != 0:
            return True
        else:
            return False

$a = 2$ の場合に素数判定をしてみましょう。

a = 2
ppt = ProbabilisticPT(a)
for n in range(1, 10):
    if ppt.test(n):
        print(f'{n}は底{a}の確率的素数')
    else:
        print(f'{n}は素数でない')

# 1は素数でない
# 2は底２の確率的素数
# 3は底２の確率的素数
# 4は素数でない
# 5は底２の確率的素数
# 6は素数でない
# 7は底２の確率的素数
# 8は素数でない
# 9は底２の確率的素数

確率的素数判定において「素数でない」という判定は信用できます。なぜならテストの持つ性質の一つ、「素数ならば必ず合格になる」の対偶は「不合格ならば必ず素数でない」だからです。よって4, 6, 8は素数でないことが確定しました。しかし、「確率的素数」と判定された場合、そこからはあまり情報が得られませんので注意が必要です。

確率的素数判定のメリットはやはりその計算速度でしょう。今回の例で言えばたった１回の割り算で「確率的素数」か「素数でない」かが判定できます。しかしあくまでも得られる結果は「素数の可能性あり」であり、本当に素数かはわかりません。実際、合成数である９も確率的素数と判定されています。

1.3. 精度を向上させるために

確率的素数判定法では判定の精度を向上するために複数の底でテストを行います。もしある一つの底で「素数でない」と判定されたならその自然数は素数でないことが確定するので判定精度の向上が期待できます。
実際に底が2, 3, 5の場合に３０未満の自然数についてテストを行ってみます。

ppt = ProbabilisticPT()
p = {}
for a in [2, 3, 5]:
    p[a] = set()
    ppt.change_base(a)  # 底をaに設定
    for n in range(1, 30):
        if ppt.test(n):
            p[a].add(n)
for k, v in p.items(): print(f'底{k}の確率的素数 : {v}')

# 各底の確率的素数の集合の共通部分を求める
print('全ての底での確率的素数 :', p[2] & p[3] & p[5])



# 底2の確率的素数 : {2, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29}
# 底3の確率的素数 : {2, 3, 4, 5, 7, 8, 10, 11, 13, 14, 16, 17, 19, 20, 22, 23, 25, 26, 28, 29}
# 底5の確率的素数 : {2, 3, 4, 5, 6, 7, 8, 9, 11, 12, 13, 14, 16, 17, 18, 19, 21, 22, 23, 24, 26, 27, 28, 29}
# 全ての底での確率的素数 : {2, 3, 5, 7, 11, 13, 17, 19, 23, 29}

底が2の場合に確率的素数と判定された９は底が3のテストで合成数であることが確定しました。
このように、複数の底を組み合わせることで全ての底での確率的素数が素数である確率を十分な精度まで向上させるのが確率的素数判定法の考え方です。今回用いたテストでは３つの底（2, 3, 5）を組み合わせることで30未満の自然数を100%の精度で素数判定することができました。

1.4. フェルマーテスト

フェルマーの小定理をテストに用いるものをフェルマーテストと言います。
フェルマーの小定理とは次のものです。

$p$ を素数とし、$a$ を $p$ の倍数でない整数（$a$ と $p$ は互いに素）とするときに、

a^{p-1} \equiv 1 \pmod{p}

が成り立つ。

つまり、フェルマーテストとは

$n$ と互いに素な自然数 $a$ に対し $a^{n-1} \equiv 1 \pmod{n}$ を満たせば合格
満たさなければ不合格

というものです。
フェルマーテストは先ほど用いたテストに比べて非常に強力で、$25\times 10^9$ までの奇数の合成数11,408,012,595個の内、底が2のフェルマーテストに合格できるのは21,853個しかないそうです。(「日本語版Wikipedia-確率的素数」より)

しかし、フェルマーテストにも欠点があります。カーマイケル数と呼ばれる合成数は任意の底のフェルマーテストに合格します。よって、いくらたくさんの底を組み合わせてもフェルマーテストでは誤判定の可能性が無視できない大きさで残ります。

1.5. modの世界における1の平方根

フェルマーテストを改良する準備としてまず、modの世界における1の平方根について考えます。
modの世界における1の平方根とは２以上の自然数 $n$ に対して次の合同式を満たす $x$ のことです。

x^2 \equiv 1 \pmod{n}

明らかに $x = 1, n-1$ はこの式を満たすので、これらを自明な平方根と言い、これら以外を非自明な平方根と言います。

非自明な平方根について考えてみましょう。つまり $x$ は $１$ でも $n-1$ でもない場合を考えます。例えば $n=15$ のとき、$x = 4, 11$ は非自明な平方根となります。

\begin{aligned}
(4)^2 &= 16 \equiv 1 \pmod{15}\\
(11)^2 &= 121 \equiv 1 \pmod{15}
\end{aligned}

では $n$ が素数の場合はどうでしょうか。実はこのとき非自明な平方根は存在しません。これを背理法で示しましょう。
いま素数 $p$ を法とする非自明な平方根が存在し、これを $x$ とします。つまり、$x$ は $１$ でも $p-1$ でもなく、

x^2 \equiv 1 \pmod{p}

を満たします。このとき、

\begin{aligned}
x^2 &\equiv 1 \pmod{p}\\
x^2 - 1 &\equiv 0 \pmod{p}\\
(x + 1) (x - 1) &\equiv 0 \pmod{p}
\end{aligned}

となり、$p$ は素数なので $(x + 1)$ と $(x - 1)$ の少なくとも一方は $p$ で割り切れます。しかし、$x$ は $１$ でも $p-1$ でもないので

\begin{aligned}
(x + 1) \not \equiv 0 \pmod{p}\\
(x - 1) \not \equiv 0 \pmod{p}
\end{aligned}

となります。つまり $(x + 1)$ と $(x - 1)$ は共に $p$ で割り切れません。よって矛盾が生じたので素数 $p$ を法とする非自明な平方根は存在しないことが示されました。

1.6. 奇素数に対するフェルマーの小定理

判定したい自然数が２の場合は素数であることが自明なので、事前に処理したとします。このとき、素数は必ず奇数になるので、奇素数 $p$ に対するフェルマーの小定理を考えます。いま、$p-1$ は偶数なので自然数 $s$ と奇数 $d$ を用いて

p - 1 = 2^s \cdot d

と表せます。よってフェルマーの小定理は

a^{2^s \cdot d} \equiv 1 \pmod{p}

となり、これは

(a^{2^{s-1} \cdot d})^2 \equiv 1 \pmod{p}

と見ることもできます。前節で示した通り、素数 $p$ を法とする $1$ の非自明な平方根は存在しないので、

\begin{aligned}
a^{2^{s-1} \cdot d} &\equiv \;\;\:1 \pmod{p}\\
a^{2^{s-1} \cdot d} &\equiv -1 \pmod{p}
\end{aligned}

のどちらかになります。１つ目であった場合は $s - 1 > 0$ であればさらに

(a^{2^{s-2} \cdot d})^2 \equiv 1 \pmod{p}

と見ることができ、これもやはり

\begin{aligned}
a^{2^{s-2} \cdot d} &\equiv \;\;\:1 \pmod{p}\\
a^{2^{s-2} \cdot d} &\equiv -1 \pmod{p}
\end{aligned}

のどちらかとなります。これを繰り返していくといつかは

\begin{aligned}
a^{d} &\equiv \;\;\:1 \pmod{p}\\
a^{d} &\equiv -1 \pmod{p}
\end{aligned}

となります。
ここまでをまとめると次のようになります。

$p$ が奇素数のとき、自然数 $s$ と奇数 $d$ を用いて

p - 1 = 2^s \cdot d

と書くことができ、以下の $p$ を法とした合同式のうち必ずどれか一つを満たす。

\begin{aligned}
a^{2^{s-1} \cdot d} &\equiv -1\\
a^{2^{s-2} \cdot d} &\equiv -1\\
\cdots\\
a^{2 \cdot d} &\equiv -1\\
a^{d} &\equiv -1\\
a^d &\equiv \;\;\:1
\end{aligned}

この合同式は高々 $\log_2{p}$ 個程度ですので全てチェックしてしまいましょう。

以上より判定したい自然数 $n$ に対するテストはこのようになります。

$n=1$ ならば不合格
$n=2$ ならば合格
$n \geq 3$ が偶数ならば不合格
$n \geq 3$ が奇数ならば $n - 1 = 2^s \cdot d$ とし、$n$ と互いに素な自然数 $a$ に対し $n$ を法として
```
\begin{aligned}
a^{2^{s-1} \cdot d} &\equiv -1\\
a^{2^{s-2} \cdot d} &\equiv -1\\
\cdots\\
a^{2 \cdot d} &\equiv -1\\
a^{d} &\equiv -1\\
a^d &\equiv 1
\end{aligned}
```
のうち一つでも満たせば合格、一つも満たさなければ不合格

このテストを用いた確率的素数判定法はミラー-ラビン素数判定法と呼ばれています。

1.7. テストの実装法

３以上の奇数 $n$ に対する判定を整理します。
確率的素数判定法は複数の底でテストを行い全てに合格した場合のみ素数であると判定するので、テストの過程では不合格の場合のみを検出すれば良いです。
いま $n - 1 = 2^s \cdot d$ とし、テストの底を $a$ とします。このとき $0 \leq r < s$ となる $r$ について $a^{2^rd}$ を計算します。まず、$r=0$ のとき

a^d \equiv 1\; or\; -1 \pmod{n}

であれば底 $a$ のテストは終了です。
そうでない場合、$r=1, 2, \cdots, s-1$ について

a^{2^rd} \not \equiv 1\; or\; -1 \pmod{n}

である限り計算します。また、$r$ についての終了条件 $r < s$ は $2^rd < n - 1$ に対応します。
もし $a^{2^rd} \equiv -1 \pmod{n}$ となった場合は底 $a$ のテストは終了です。
しかし、 $a^{2^rd} \equiv 1 \pmod{n}$ となった場合はどうでしょうか。このとき、$a^{2^{r-1}d}$ は $1$ でも $-1$ でもないことはすでに確認されているので $a^{2^{r-1}d}$ は $1$ の非自明な平方根です。よって $n$ は素数でないので不合格となります。
また、最後まで $a^{2^rd} \not \equiv 1\; or\; -1 \pmod{n}$ であった場合も不合格となります。
以上をまとめると下図のようになります。

1.8. 底の選び方

ミラー-ラビン素数判定法では一般的には底の数を自分で決めて $a < n$ となる自然数 $a$ をランダムに選び底とします。計算速度と判定の精度はトレードオフの関係にあるので必要な精度を確保できる中でできるだけ少ない数の底が望ましいです。そのため効率的に精度を上げることができる底の組み合わせが研究されてきました。
ACLでは底として $a = {2, 7, 61}$ を採用しています。Jaeschke(1993)によると、この底のテストを全て通過する最小の合成数は $4759123141(=48781 \cdot 97561 > 2^{32})$ です。よって $n < 2^{32}$ の範囲(符号なし４バイト整数の範囲)では $100\%$ の精度で判定できます。

1.9. 実装

それでは実装します。なおACLにおいてpow_modを用いている部分は、同等の機能であるPythonの組み込み関数powで代用しています

def is_prime(n):
    # 自明な部分
    if (n <= 1): return False
    if (n == 2 or n == 7 or n == 61): return True
    if (n % 2 == 0): return False

    d = n - 1  # nは奇数
    while (d % 2 == 0): d //= 2  # 奇数dを求める

    for a in (2, 7, 61):
        t = d  # dは他の底でも使うので保持
        y = pow(a, t, n)

        # a^d = 1, -1ならこのループは入らない
        # a^t = 1, -1となるまで繰り返す
        while (t != n - 1 and y != 1 and y != n - 1):
            y = y * y % n
            t <<= 1
        
        # a^d = 1, -1は通過　(t%2 == 0)
        # a^t = -1は通過　(y != n - 1)
        if (y != n - 1 and t % 2 == 0):
            return False
    return True



print(is_prime(17))  # True
print(is_prime(1000000007))  # True
print(is_prime(121))  # False
print(is_prime(561))  # False (561はカーマイケル数のひとつ)

# 底{2, 7, 61}のテストを通過する最小の合成数
print(is_prime(4759123141))  # True

1.10. 計算速度の比較

時間計算量 $O(\sqrt{n})$ の決定的素数判定法と比較してみます。
比較に使用したコードは以下のものです。

import random

# ミラー-ラビン素数判定法（確率的素数判定法）
def pro_is_prime(n):
    if (n <= 1): return False
    if (n == 2 or n == 7 or n == 61): return True
    if (n % 2 == 0): return False
    d = n - 1
    while (d % 2 == 0): d //= 2
    for a in (2, 7, 61):
        t = d
        y = pow(a, t, n)
        while (t != n - 1 and y != 1 and y != n - 1):
            y = y * y % n
            t <<= 1
        if (y != n - 1 and t % 2 == 0):
            return False
    return True

# 決定的素数判定法
def det_is_prime(n):
    if n < 2: return False
    if n == 2: return True
    if n % 2 == 0: return False
    for i in range(3, int(n ** 0.5) + 1):
        if n % i == 0: return False
    return True


def random_1():
    l1, r1 = [3, 1 << 16]
    return random.randrange(l1, r1, 2)

def random_2():
    l2, r2 = [(1 << 16) + 1, 1 << 32]
    return random.randrange(l2, r2, 2)

def random_3():
    l3, r3 = [3, 1 << 32]
    return random.randrange(l3, r3, 2)


def main():
    """
    seed_list = [111, 222, 333, 444, 555, \
                 666, 777, 888, 999, 101010]
    """
    random.seed(111)  # 乱数固定
    loop = 10000  # ループ回数 10^4 or 10^6

    for _ in range(loop):
        n = random_1()
        #n = random_2()
        #n = random_3()

        pro_is_prime(n)  # 確率的
        #det_is_prime(n)  # 決定的
    

if __name__ == "__main__":
    main()

計測方法

AtCoderのコードテスト(Python 3.8.2)を使用しました。
10種類のシード値で3種類の範囲から乱数を生成し、素数判定10000回あたりの実行時間を調べました。
偶数の場合はどちらも最初に処理されるので奇数のみとしました。

random_1
- 乱数の範囲 : $[3, 2^{16}]$ の奇数
- ループ回数 : $10^4, 10^6$ ($10^4$ ではimportなどのオーバーヘッドがメインになる可能性があるため)
random_2
- 乱数の範囲 : $[2^{16}, 2^{32}]$ の奇数
- ループ回数 : $10^4$
random_3
- 乱数の範囲 : $[3, 2^{32}]$ の奇数
- ループ回数 : $10^4$

測定結果

結果は下図のようになりました。数値は素数判定10000回あたりの時間で単位はmsです。

	ミラー-ラビン	決定的素数判定法
random_1(10^4)	63	59
random_1(10^6)	34	30
random_2	100	4060
random_3	99	4113

判定したい数が $10^5$ 程度の小さな数であれば決定的素数判定法の方がわずかに速いようです。
しかし、数が大きくなるとミラー-ラビン素数判定法の優位性が顕著になります。

2. primitive_root

素数 $m$ を法とする原始根のうち最小のものを求めます。

2.1. 位数

最初に、原始根を説明する上で重要な用語である「位数」についてみていきます。定義はこうです。

２以上の自然数 $m$ および $m$ と互いに素な整数 $a$ に対して

a^{n} \equiv 1 \pmod{m}

となる最小の自然数 $n$ を $m$ を法とする $a$ の位数という。

具体例をみてみましょう。
$m = 5, a = 3$ のとき、

\begin{aligned}
3^1  &= 3 & \not \equiv 1 \pmod{5}\\
3^2  &= 9 & \not \equiv 1 \pmod{5}\\
3^3  &= 27 & \not \equiv 1 \pmod{5}\\
3^4  &= 81 & \equiv 1 \pmod{5}
\end{aligned}

なので $5$ を法とする $3$ の位数は $4$ です。
また、$m = 12, a = 5$ のとき、

\begin{aligned}
5^1  &= 2 & \not \equiv 1 \pmod{12}\\
5^2  &= 25 & \equiv 1 \pmod{12}
\end{aligned}

なので $12$ を法とする $5$ の位数は $2$ です。

位数の性質

$m$ を法として $m$ と互いに素な整数 $a$ の位数 $e$ は正整数 $n$ に対して次の性質を持ちます。

a^n \equiv 1 \pmod{m} \Leftrightarrow e\; は \;n\; の約数

(証明)
$\Leftarrow:$
$e$ は $n$ の約数なので自然数 $k$ を用いて $n = ke$ とかける。よって $m$ を法として

\begin{aligned}
a^n &=a^{ke}\\
&= (a^e)^k\\
&\equiv 1
\end{aligned}

となり、成立。
$\Rightarrow:$
$n$ を非負整数 $q, r$ を用いて $n=qe + r ;(0 \leq r < e)$ と表す。このとき $m$ を法として

\begin{aligned}
a^r &\equiv (a^{e})^qa^r\\
&\equiv a^{qe+r}\\
&\equiv 1
\end{aligned}

となる。もし $0 < r < e$ だとすると、$e$ が位数であることと矛盾する（位数は $a^e\equiv 1 \pmod{m}$ を満たす最小の整数）。よって $r = 0$ であり、すなわち $e$ は $n$ の約数となる。

(証明終)

特に $m$ が素数の場合、フェルマーの小定理より位数は $m - 1$ の約数となります。

2.2. 原始根

$m$ が素数の場合を考えます。フェルマーの小定理から $m$ を法として $m$ と互いに素な整数 $a$ の位数は必ず $m - 1$ 以下になります。そこで、位数がちょうど $m - 1$ となる $a$ を特別視し、これを $m$ を法とする原始根と呼びます。

例えば、$m = 7$, $a = 3$ のとき

\begin{aligned}
3^1  &= 3  &\not \equiv 1 \pmod{7}\\
3^2  &= 9  &\not \equiv 1 \pmod{7}\\
3^3  &= 27  &\not \equiv 1 \pmod{7}\\
3^4  &= 81  &\not \equiv 1 \pmod{7}\\
3^5  &= 243  &\not \equiv 1 \pmod{7}\\
3^6  &= 729  &\equiv 1 \pmod{7}\\
\end{aligned}

なので $7$ を法とする $3$ の位数は $6$ です。よって $3$ は $7$ を法とする原始根です。
原始根は１つとは限りません。位数のところで示した例の１つ目から、$3$ は $5$ を法とする原始根であることがわかります。一方、$a = 2$ のとき

\begin{aligned}
2^1  &= 2 & \not \equiv 1 \pmod{5}\\
2^2  &= 4 & \not \equiv 1 \pmod{5}\\
2^3  &= 8 & \not \equiv 1 \pmod{5}\\
2^4  &= 16 & \equiv 1 \pmod{5}
\end{aligned}

となるので $2$ もまた $5$ を法とする原始根です。

なお、$m$ が素数のとき原始根は必ず存在します。
証明は「原始根定理」で調べてみてください。

原始根の性質

素数 $m$ を法とする原始根 $g$ は次の性質を持ちます.

$m$ を法とした $g$ のべき乗の集合

\{g\pmod{m}, g^2\pmod{m}, \cdots , g^{m - 1}\pmod{m}\}

と整数を $m$ で割ったあまりの集合から $0$ を除いた集合

\{1, 2, \cdots , m - 1\}

は一致する。

これは次のように言い換えられます。
$1 \leq k \leq m-1$ の自然数 $k$ について

$g^k \not \equiv 0$
$g^k$ はすべて異なる。

１つ目は $g$ と $m$ が互いに素であることから明らかです。
よって２つ目を証明します。

(証明)
いま、$k < l \leq m - 1$ なる自然数 $k, l$ が存在し、$g^k \equiv g^l\pmod{m}$ であったとします。このとき

\begin{aligned}
g^l - g^k &\equiv 0 \pmod{m}\\
g^k(g^{l-k} - 1) &\equiv 0 \pmod{m}\\
g^{l-k} - 1 &\equiv 0 \pmod{m}\\
g^{l-k} &\equiv 1 \pmod{m}
\end{aligned}

となります。$l - k < m - 1$ なので原始根である $g$ の位数は $m-1$ であることと矛盾します。
したがって、$1 \leq k \leq m-1$ に対し、$g^k$ はすべて異なります。
(証明終)

2.3. 原始根であるかの確認方法1

ある２以上の自然数 $g$ が $m$ を法とする原始根であるかどうかを判定するためには、原始根の定義から $m$ を法として $g, g^2, \cdots, g^{m - 2}$ がすべて $1$ と合同でないことを確かめれば良いです。したがって、最小の原始根を求めるアルゴリズムの実装は以下のようになります。

def naive_primitive_root(m):
    g = 2
    while True:
        for i in range(1, m - 1):
            if pow(g, i, m) == 1:
                g += 1
                break
        else:
            return g


print(naive_primitive_root(7))  # 3
print(naive_primitive_root(23))  # 5
#print(naive_primitive_root(1000000007))  # 非常に時間がかかります

この方法では $g$ の $m - 2$ までのべき乗をすべて調べるので $m$ が大きくなると非常に時間がかかります。

2.4. 原始根であるかの確認方法2

実は原始根の性質を用いることでより簡単に原始根であるかの確認をすることができます。用いる性質は次のものです。

素数 $p_i$ および自然数 $e_i$ を用いて $3$ 以上の素数 $m$ に対して $m-1$ が

m-1 = \prod_{i = 1}^n{p_i^{e_i}}

と素因数分解されたとき、$2$ 以上の自然数 $g$ について次のことが成り立つ。

g\;が\;m\;を法とする原始根
\Leftrightarrow 1 \leq i \leq n\; に対し\; g^{\frac{m-1}{p_i}}\not\equiv 1 \pmod{m}

(証明)
$\Rightarrow:$
$p_i$ は $m-1$ の素因数なので $\frac{m-1}{p_i}$ は整数であり、$1 \leq \frac{m-1}{p_i} < p - 1$ を満たす。$g$ が原始根のとき $g^{p-1}$ が初めて $1$ と合同になるべき乗なので $1 \leq i \leq n\;に対し\; g^{\frac{m-1}{p_i}}\not\equiv 1 \pmod{m}$ となる。

$\Leftarrow:$
この主張の対偶は

gは原始根でない \Rightarrow g^{\frac{m-1}{p_i}}\equiv 1 \pmod{m}となるiが存在する

なのでこれを示す。
いま $g$ は原始根でないので位数 $e$ は $m-1$ 未満である。位数の性質で示した通り、$e$ は $m-1$ の約数なので $m-1$ の素因数と同じ $p_i$ を用いて

e = \prod_{j = 1}^n{p_j^{e_j^{\prime}}}\;\;\;(e_j^{\prime} \leq e_j)

と素因数分解でき、とくに $e_j^{\prime} < e_j$ となる $j$ が少なくとも１つ存在する。このような $j$ のうち１つを $i$ とする。このとき

\begin{aligned}
\frac{m-1}{p_i} &= \frac{1}{p_i}\prod_{j}{p_j^{e_j}}\\
&= \frac{1}{p_i}\left(\prod_{j}{p_j^{e_j - e_j^{\prime}}}\right)\left(\prod_j{p_j^{e_j^{\prime}}}\right)\\
&= p_i^{e_i - e_i^{\prime} - 1}\left(\prod_{j \ne i}{p_j^{e_j - e_j^{\prime}}}\right)\left(\prod_j{p_j^{e_j^{\prime}}}\right)
\end{aligned}

ここで、$e_i - e_i^{\prime} - 1 \geq 0$ かつ $e_j - e_j^{\prime} \geq 0$ なので

\frac{m-1}{p_i} = (自然数) \times e

となる。
したがって、$g$ が原始根でないならば

g^{\frac{m-1}{p_i}} \equiv \left(g^e\right)^{自然数}　\equiv 1 \pmod{m}

となる $i$ が存在する。

(証明終)

以上より

$m-1$ の素因数を求める。
各素因数 $p_i$ について $g^{\frac{m-1}{p_i}}\not\equiv 1 \pmod{m}$ となるか調べる。

という手順で $g$ が原始根であるかを判定できます。

2.5. 実装

それでは実装します。is_primeと同様にpow_modは組み込み関数powで代用します。

# @param m must be prime
def primitive_root(m):
    if m == 2: return 1
    if m == 167772161: return 3
    if m == 469762049: return 3
    if m == 754974721: return 11
    if m == 998244353: return 3

    # m-1の素因数抽出
    divs = [2]
    x = (m - 1) // 2
    while x % 2 == 0: x //= 2
    i = 3
    while i * i <= x:
        if x % i == 0:
            divs.append(i)
            while x % i == 0: x //= i
        i += 2
    if x > 1: divs.append(x)

    # 全ての素因数で1と合同でない最小のgを探す
    g = 2
    while True:
        if all(pow(g, (m - 1) // div, m) != 1 for div in divs): return g
        g += 1


print(primitive_root(7))  # 3
print(primitive_root(23))  # 5
print(primitive_root(1000000007))  # 5

最初に判定しているいくつかの $m$ はconvolutionで使われるmodのようです。

3. おわりに

今回は素数判定法と原始根について見てきました。とくに確率的素数判定法の考え方が面白いと感じました。

internal_mathのなかで今回触れなかったものについてはinternal_math編①で書いていますので、よろしければそちらもご覧ください。

説明の間違いやバグ、アドバイス等ありましたらお知らせください。