原始根のアルゴリズム

1 原始根とは
2 原始根の存在定理
- 2.1 例題
3 原始根を求めるアルゴリズム
4 原始根を求めないアルゴリズム
- 4.1 Tonelli-Shanksのアルゴリズム
- 4.2 部分乗法群
5 円分体
- 5.1 FFT
6 演習問題
7 参考文献

原始根とは

$\bmod p$ において、$g^0,g^1,\ldots,g^{p-2}$ が相異なるような $g\neq 0\pmod p$ を原始根と呼ぶ。つまり、$g^0,g^1,\ldots,g^{p-2}$ を並び替えると $1,2,\ldots,p-1$ になるような $g$ である。

例：例えば、mod 5 では
\begin{align}
2^0&=1 \\
2^1&=2 \\
2^2&=4 \\
2^3&=3
\end{align}であるから、mod 5 において $2$ は原始根である。また、
\begin{align}
3^0&=1 \\
3^1&=3 \\
3^2&=4 \\
3^3&=2
\end{align}であるから、mod 5 において $3$ も原始根である。

原始根 $g$ が存在すると、対数のように $g$ の指数部分に着目することで掛け算を足し算に帰着することができる。そして、次節で示すが、必ず $g$ は存在する。

原始根 $g$ の存在を認めると、各 $a \in (\mathbb{Z}/p\mathbb{Z})^\times$ に対して、$a=g^i$ となる唯一の $0 \leq i < p-1$ が存在する。この $i$ を「$\bmod p$ における、$g$ を底とする $a$ の指数 $\mathrm{ind_g}(a)$」と呼ぶ。

原始根の存在定理

原始根の存在定理として知られている次の定理を示す（難しくはないが、道具として使う分には読み飛ばしても問題ない）。

定理：$\bmod p$ において、原始根は必ず存在する。

証明：原始根が存在しない、つまり位数 $p-1$ の元が存在しないと仮定して矛盾を導く。

$d$ を $p-1$ の約数とする。Fermatの小定理より $x^{p-1}=1$ はちょうど $p-1$ 個の解を持つ。$\mathbb{Z}/p\mathbb{Z}$ は体だから $n$ 次方程式は高々 $n$ 個の解を持つ。因数分解 $x^{p-1}-1=(x^d-1)\sum_{k=0}^{(p-1)/d-1} x^{kd}$ で左辺はちょうど $p-1$ 個の根を持つから、$x^d=1$ はちょうど $d$ 個の解を持つ（原始根 $g$ の存在を認めれば、解は $x=g^{(p-1)/d},g^{2(p-1)/d},\ldots,g^{(d-1)(p-1)/d}$ である）。

$x^{p^{s-1}}=1$ の解は $x^{p^s}=1$ の解でもあるから、前者に含まれない後者の $p^s-p^{s-1}$ 個の解は位数 $p^s$ である。よって $p-1=\prod p_i^{e_i}$ と素因数分解できるとき、位数 $p_i^{e_i}$ の元が存在する（原始根 $g$ の存在を認めるならば、位数 $p_i^{e_i}$ の解は $x=g^{n(p-1)/p_i^{e_i}}$ s.t. $\gcd(n,p_i)=1$ である）。

元 $a, b$ の位数 $s,t$ が互いに素ならば、$ab$ の位数が $st$ になることを示す。$(ab)^{st} = (a^s)^t(b^t)^s=1$ であるから、$ab$ の位数は $st$ を割り切る。逆に $(ab)^m=1$ ならば $a^{m}$ と $b^{-m}$ がそれぞれ生成する部分群は同一で、その位数は $s, t$ の約数である。$s,t$ は互いに素であったから、$a^m=b^m=1$ となる。つまり $(st)|m$ を得る。以上より、$ab$ の位数は $st$ である。

位数 ${p_i}^{e_i},{p_j}^{e_j}$ の２つの元から位数 ${p_i}^{e_i}{p_j}^{e_j}$ の元が生成できることが分かった。この操作を全ての素冪の因数について行うことで位数 $p-1$ の元、すなわち原始根を構成することができる。

これで安心して原始根が使える。

例題

原始根が具体的にどれか分からなくても、その存在を仮定するだけで多くの性質を導くことができる。$\bmod p$ の問題に対して、取り敢えず原始根を考えるというのは（計算量的に原始根を求めるのが難しくとも）有用な手だと思う。

問題１：$1^k+2^k+\dots+(p-1)^k \pmod p$ を計算せよ。

解法：原始根 $g$ により、$1,2,\ldots,p-1$ は並び替えると $g^0,g^1,\ldots,g^{p-2}$ になる。すると、べき乗和は等比級数に変換できる：
$$1^k+2^k+\dots+(p-1)^k=g^{0}+g^{k}+\dots+g^{(p-2)k}$$
等比級数の和の公式により次を得る：
$$g^{0}+g^{k}+\dots+g^{(p-2)k}=
\begin{cases}
p-1 & ( k = 0 \pmod {p-1} ) \\
0 & (\text{otherwise})
\end{cases}\\
$$

問題２：$(p-1)!= -1 \pmod p$ を示せ。（出典：Wilsonの定理）

証明：$p = 2$ のとき、$(p-1)! = -1 \pmod 2$ より成立。以下、$p$ を奇素数とする。原始根 $g$ により $$(p-1)! = \prod_{i=0}^{p-1} g^{i} = g^{(p-1)p/2} \pmod p$$ である。$g^{p-1}=1$ より $g^{(p-1)/2} = \pm 1$ であるが、$g^{(p-1)/2} = 1$ だと $g$ の位数が $p-1$ であることに反する。従って、$g^{(p-1)/2} = -1$ である。$p$ は奇素数であったから、$g^{(p-1)p/2} = -1$ である。

問題３：$g^a, g^b$ の２元が生成する乗法群を求めよ。

解法：$g^a,g^b$ の生成する乗法群の元の指数全体は $a\mathbb{Z}+b\mathbb{Z}+(p-1)\mathbb{Z}=\gcd(a,b,p-1)\mathbb{Z}$ である。よって、$g^{\gcd(p-1,a,b)}$ が生成する乗法群に等しい。

また、この事実から位数 $a,b$ の２元が生成する乗法群の位数は $\mathrm{lcm}(a,b)$ であることが導かれる。

類題：CodeForces360D

原始根を求めるアルゴリズム

前節では原始根の存在を示した。次は原始根を求めるアルゴリズムを構築しよう。原始根は乱択により、効率的に求めることが出来る。

問題４：$\bmod p$ における原始根を一つ求めるアルゴリズムを構築せよ。

原始根が出るまでランダムに整数 $a$ を選び続ける。定義から、 $a$ が原始根であることと $a$ の位数が $p-1$ であることは同値である。Fermatの小定理（あるいはLagrangeの定理）より $(\mathbb{Z}/p\mathbb{Z})^\times$ の任意の元の位数は $p-1$ の約数である。従って、各 $p-1$ の素因数 $p_i$ について $a^{(p-1)/p_i} \neq 1 \pmod p$ であることを確かめれば十分。この部分の素直な計算量は素因数分解を前計算しておけば $O((\log p)^2)$ であるが、工夫すると $O(\log(p) \log(\log p))$ にできる（参考文献[1]を参照せよ）。

試行回数を評価しよう。原始根は $\phi(p-1)$ 個あるから、原始根を引くまでの試行回数の期待値は $\frac{p-1}{\phi(p-1)}$ である。$\frac{n}{\phi(n)} = O(\log(\log{n}))$ であるから、試行回数は expect $O(\log(\log{p}))$ 回である。$O(n/\phi(n))$ の評価はこちらの記事「素数が関する計算量」を参照せよ。

計算量は $p-1$ の素因数分解がボトルネックになる。素因数分解の計算量は $\sqrt{p}$ まで試し割りする $O(p^{1/2})$ や ρ法とミラー・ラビン素数判定法を組み合わせた expect $O(p^{1/4})$ [3]がある。

原始根が求まれば、その原始根を底とする指数が次のように求められる：

問題５：$a \neq 0 \pmod p$ とする。原始根 $g$ に対して、$a=g^x \pmod p$ なる $x$ を求めるアルゴリズムを構築せよ。

この問題は離散対数問題そのものであるから、Baby-Step Giant-Step により $O(\sqrt{p})$ で求まる。

原始根を求めないアルゴリズム

原始根を求めるには $O(p^{1/4})$ 、さらにその原始根を底とする指数を求めるには一回当たり $O(p^{1/2})$ かかる。原始根を実際に求めたり、さらにその指数を実際に求めたりするには、この計算量は不可避である。従って、アルゴリズムを高速にするには、「原始根を求めない」ということが重要になる。

どのように原始根を迂回するかは、臨機応変に考える必要がある（が、その意識を持つだけで変わると思う）。実例を２つ紹介する。

Tonelli-Shanksのアルゴリズム

次の問題を原始根と指数を陽に求めることで解いてみよう。

問題６：$p$ を奇素数として、$a\neq0\pmod p$ とする。$x^2=a \pmod p$ を満たす $x$ を求めよ。

解法（改善前）：原始根を $g$ として、 $x=g^y$ だとする。すると、$x^{2y}=a$ は離散対数問題そのものだから、$O(p^{1/2})$ で解ける。

原始根と指数を実際に求める方針では、これ以上の計算量改善は不可能である。しかし、原始根を陽に求めず、その情報を部分的に引き出せば、計算量を $O(\log(p)^2)$ にできる。このアルゴリズムはTonelli-Shanksのアルゴリズムと呼ばれている。

Tonelli-Shanksのアルゴリズムの準備として次の問題を考えよう。この問題は、後述する誤差項の指数を 0 にする操作と対応している。

問題７：$w = 1 \pmod 2$ なる未知数 $w$ と $z$ が与えられる。$$w + zx = 0 \pmod {2^Q}$$ となる $zx$ を求めたい。許された操作は次の３つである。

bitShift(u, v)
- $u2^v \pmod {2^Q}$ を返す。
isZero(u)
- $u = 0 \pmod {2^Q}$ ならば True、さもなくば False を返す。
add(u, v)
- $u + v \pmod {2^Q}$ を返す。

各関数の計算量は bitShift が $O(Q)$、それ以外が $O(1)$ である。$zx$ を $O(Q^2)$ で求めるアルゴリズムを構築せよ。

解法：$w+zx$ が２進法表記でいくつ末尾に $0$ が連なっているかは、$w+zx$ をいくつ左に bit-shift すると $0$ になるかで分かる。つまり、末尾に $k$ 個 $0$ が連なるならば、$(w+zx)2^m = 0 \pmod {2^Q}$ となる最小の $m$ が $Q-k$ に等しい。このとき、$z$ を k-bit だけ左に shift した数を $w+zx$ に加えると、末尾に連なる 0 の個数を少なくとも１つ増やすことができる。つまり、$zx\leftarrow zx+z2^k$ とする。この操作を繰り返すと $O(Q^2)$ で $w+zx=0\bmod {2^Q}$ となる。Python によるコードは次のようになる。

zx = 0
for i in range(Q):
   if not isZero(bitShift(add(w, zx), Q - 1 - i)):
	zx = add(zx, bitShift(z, i)

この問題の関数は、原始根や指数を知らずとも可能な、指数の世界の計算と対応している。次の通りである。

$q$ を奇数として $p-1=q2^Q$ だとする。中国剰余定理より $(p-1)\mathbb{Z}\cong q\mathbb{Z}\times 2^Q\mathbb{Z}$ と分解できる。簡単のため、$\bmod q$ では指数が $0$ になっているとして、$q\mathbb{Z}$ をひとまず無視する。

$$\mathrm{ind}(u)=0 \pmod {2^Q} \Leftrightarrow u=1 \pmod p$$であるから、isZero(ind(u)) は $O(1)$ で可能である。$$\mathrm{ind}(x)=2^v\mathrm{ind}(u) \pmod {2^Q} \Leftrightarrow x=u^{2^v} \pmod p$$であるから、繰り返し二乗法によって bitShift(ind(u), v) は $O(Q)$ で可能である。$$\mathrm{ind}(x)=\mathrm{ind}(u)+\mathrm{ind}(v) \pmod {2^Q} \Leftrightarrow x=uv \pmod p$$ であるから、add(ind(u), ind(v)) は $O(1)$ で可能である。

以上、原始根や指数を知らずとも、$\bmod q$ で指数が $0$ ならば、これらの操作が指数の世界で可能なことを示した。

問題７を踏まえ、Tonelli-Shanksのアルゴリズムの方針は次の通りである。

$x^2=a \pmod p$ の解の第一近似として $\mathrm{ind}(y^2a^{-1})=0 \pmod q$ なる $y$ を取る。
$\mathrm{ind}(b)=0 \pmod q, \mathrm{ind}(b) = 1 \pmod 2$ なる $b$ を取る。
問題７のようにして、$y$ を $b$ のべき乗を乗算した値で更新することで、$\mathrm{ind}(y^2a^{-1})=0 \pmod {2^Q}$ とする。このとき $y^2=a \pmod {p}$ となって解が得られる。

$x$ の近似 $y$ を求めよう。$a^{(q+1)/2} \pmod p$ を２乗すると、その指数は $(q+1)\mathrm{ind}(a) = \mathrm{ind}(a)\pmod q$ となって、$\bmod q$ で $a$ と等しい。よって、$x$ の近似 $y$ として $a^{(q+1)/2}$ を選び、誤差項 $y^2a^{-1}$ を修正していく。

次に $b$ を求めよう。指数が $\bmod (p-1)$ で奇数である数 $b’$ をランダムに取る。奇数は $(p-2)/2$ 個あるので確率 $(p-2)/(2(p-1))$ で見つかる。指数の偶奇は $\mathrm{ind}(x) = 0\pmod 2 \Leftrightarrow x^{(p-1)/2}=1 \pmod p$ で分かる（問題１０）。$q$ 乗した $b’^q$ の指数は、 $\bmod q$ では $0$ に、$\bmod 2^Q$ では奇数になる。よって、$b=b’^Q$ とすればよい。

以上、$y$ と $b$ が $O(\log p)$ で得られた。誤差項を修正する部分も併せると、全体で $O(\log(p)^2)$ になる。

部分乗法群

ある元 $a=g^b$ の生成する部分乗法群の元の指数全体は $b\mathbb{Z}+(p-1)\mathbb{Z}=\gcd(b,p-1)\mathbb{Z}$ となる。従って、この部分乗法群を求めるには、$a$ の原始根 $g$ を底とする指数を求めればよい。指数を経由せずに部分乗法群を求められないだろうか。部分乗法群の位数 $d$ が与えられると、その部分乗法群は一意に定まる。なぜなら、$x^d=1\pmod p$ は高々 $d$ 個の解しか持たないからである。従って、$p-1$ の素因数分解が前計算されている場合、$a$ の位数を計算することで、指数を求めずに済む。位数の計算は原始根判定と同様にすれば、$O(\log(\log(p))^2\log(p))$ である。

次の例題のように、生成元の個数が多い場合、原始根の計算を迂回することで計算量を落とすことができる。

問題８：$a_1,a_2,\ldots,a_n \in (\mathbb{Z}/p\mathbb{Z})^\times$ から生成される乗法群に $b$ が含まれるか判定せよ。つまりある整数 $m_1,m_2,\ldots,m_n$ が存在して、$b=\prod a_i^{m_i} \pmod p$ とできるか判定せよ。

$2 \leq p < 2^{31}$
$1 \leq n < 10^{5}$

（出典：RitsCamp2019Day2 L）

解法：原始根 $g$ を用いて、$a_i=g^{c_i}, b=g^{d}$ と書けるとする。すると、$a_1,a_2,\ldots,a_n$ が生成する乗法群の元の指数全体は $(p-1)\mathbb{Z}+\sum c_i\mathbb{Z}=\gcd(p-1,c_1,\ldots,c_n)\mathbb{Z}$ と書ける。よって $d$ が $\gcd(p-1,a_1,\ldots,a_n)$ の倍数であるか判定すればよい。ところで $a_i$ の生成する乗法群の位数 $q_i$ から、その乗法群が一意に定まる（なぜなら $x^{q_i}=1$ の解は高々 $q_i$ 個）。位数 $q_i$ ならば、その指数全体は $\frac{p-1}{q_i}\mathbb{Z}/(p-1)\mathbb{Z}=c_i\mathbb{Z}/(p-1)\mathbb{Z}$ である。$c_i$ を $(p-1)/q_i$ で置き換えてもよい。$q_i$ は原始根判定と同様にすれば $O(\log(\log(p))^2\log(p))$ で求まるから、全体では $O(\sqrt{p}+n\log(\log(p))^2\log(p))$ である。

円分体

$\mathbb{F}_p$ での計算を円分体 $\mathbb{F}_p(\zeta_n)$ に持ち上げて行った後に $\mathbb{F}_p$ に戻すという操作を、$\zeta_n=g^{(p-1)/n}$ で代用することができる。$n \mid (p-1)$ でなければならない。そのような素数 $p$ はどの程度の割合で見つかるだろうか。初項 $1$、公差 $n$ の等差数列の $x$ 以下の項のうち、素数であるものの個数は算術級数の素数定理より
$$\Theta\left(\frac{x}{\phi(n)\log(x)} \right)$$
である。ただし $\phi$ はオイラーのトーティエント関数。$\phi(n)\leq n$ だから $x$ を第 $k$ 項目とすると素数の個数は $\Theta\left(\frac{kn}{n\log(kn)}\right)$である。従って、平均 $\Theta(\log(n))$ 項目まで探索すれば素数が見つかる。

FFT

素数 $p=2^{23}\times119+1$ における $\mathbb{Z}/p\mathbb{Z}$ で FFT を計算する際に、原始根 $g$ が活用できる。FFT では 2 冪乗根の回転子が必要であった。そこで、$g^{(p-1)/119}$ を $2^{23}$ 乗根の回転子として活用することができる[2]。$fg$ の巡回畳み込みの $k$ 次の項は $\mathbb{F}_p(\zeta_k)$ での離散フーリエ変換を利用することで線形の空間計算量で求まる。すなわち、$\frac{1}{k}\sum_i f(\zeta^i)g(\zeta^i)$ を求めればよい（yuki3024）。

演習問題

※出典の問題の一部分の要素だけを取り出していたり、表現を変えたりしています。具体的な数値を課している場合、$10^8$ step 以内に収まるように解いて下さい。

問題９：数列 $f_1,f_2,\ldots$ が
$$f_i = \prod_{j=1}^k f_{i-j}^{b_j} \pmod p$$
の関係を満たしている。$f_1=f_2=\ldots=f_{k-1}=1,f_n=m\neq0$ と $b_1,b_2,\ldots,b_k$ が分かっている。そのような数列が存在するか判定し、存在するならば、$f_k$ としてあり得る値を求めよ。

$1 \leq k \leq 100$
$k < n \leq 10^{9}$
$p=998244353$

（出典：CodeForces 1106F）

解法：$f_i\neq0$より原始根 $g$ を用いて $f_i=g^{a_i}$ とできる。このとき、漸化式は
$$a_i = \sum_{j=1}^k b_ja_{i-j} \pmod {p-1}$$
となる。$a_1=a_2=\ldots=a_{k-1}=0$ であるから、繰り返し二乗法によって $a_n=Aa_{k} \pmod {p-1}$ となるような $A$ が求められる。あとは $a_k$ をこの式から求めればよい。

問題１０：$a\neq 0\pmod p$ とする。次の３つは同値であることを示せ。
(i) $x^k = a \pmod p$ を満たす $x$ が存在する。
(ii) $g$ を原始根とすると $a \in \{g^{1\cdot\gcd(p-1,k)},\ldots,g^{\frac{p-1}{\gcd(p-1,k)}\cdot\gcd(p-1,k)}\}$ である。
(iii) $a^{(p-1)/\gcd(p-1,k)} = 1 \pmod p$ である。

証明：(i) $\Rightarrow$ (ii) を示す。$g$ の位数は $p-1$ であるから $x = g^n$ とすると、$x^k = g^{\{kn\pmod{p-1}\}}$ である。ユークリッドの互除法より $k\mathbb{Z} + (p-1)\mathbb{Z} = \gcd(k,p-1)\mathbb{Z}$ であるから、(i) $\Rightarrow$ (ii) である。

(ii) $\Rightarrow$ (iii) はFermatの小定理より成立する。

(iii) $\Rightarrow$ (i) を示す。$a=g^n$ とすると $a^{(p-1)/\gcd(p-1,k)}=g^{n(p-1)/\gcd(p-1,k)}=1$ より $n(p-1)/\gcd(p-1,k)=0 \pmod {p-1}$ である。これは $n=0 \pmod {\gcd(p-1,k)}$ を導く。

問題１１：素数 $p$ と２つの数列 $(a_i)_{1 \leq i < p}$ と $(b_i)_{1 \leq i < p}$ が与えられる。$c_k = \sum_{i \times j = k \bmod p}a_ib_j$ によって定められる数列 $(c_i)$ を $o(p^2)$ で計算するアルゴリズムを構築せよ。（出典：yukicoder : No.931）

解法：原始根を一つ求めて $g$ と置く。$i = g^l, j=g^m$ とすると $k = g^{l + m}$ であるから、$g$ を不定元に置き換えた多項式の掛け算で数列 $c$ が計算できる。多項式の掛け算は FFT を用いれば $O(p\log(p))$ である。

問題１２：$$a_1x_1^k+a_2x_2^k+\dots+a_nx_n^k=b \pmod p$$を満たす $(x_1,\ldots,x_n)\in(\mathbb{Z}/p\mathbb{Z})^n$ の個数を求めるアルゴリズムを構築せよ

$1 \leq n \leq 500$
$2 < p < 10^{5}$

（出典：yukicoder : No.1025）

解法：$k\mathbb{Z}+(p-1)\mathbb{Z}=\gcd(k,p-1)\mathbb{Z}$ より $k$ を $\gcd(k,p-1)$ で置き換えても解の個数は変わらない。$x_i \neq 0$ ならば変数変換 $x_i=g^{km_i}$ としても一般性を失わない。$g^k$ の位数が $(p-1)/k$ だから $0 \leq m_i < (p-1)/k$ の範囲の遷移だけ真面目に計算すればよい。$$\begin{align}
&a_1x_1^{k}+\ldots+a_ix_i^{k}=b’ \\ \Leftrightarrow&a_1(x_1g^N)^k+\ldots+a_i(x_ig^N)^k=b’g^{kN}
\end{align}
$$ より、任意の整数 $N$ について、$i$ 項目までの和が $b$ になるものと $bg^{kN}$ になるものの個数は等しい。つまり「dp[i][j] = i 項目までの和が j になる場合の数」という dp をする際に $j=0,g^0,g^1,\ldots,g^{k-1}$ の部分だけ更新すればよい。以上より、全体で $O(n(1+(p-1)/k)(k+1))=O(np)$ となる。

参考文献

熨斗袋氏の Tweet：その１
- 原始根判定に使う冪乗の計算の高速化。位数に使う場合は二分探索すればよい。
熨斗袋氏の Tweet：その２
- NTT-friendly な素数の個数の見積り。
素因数分解 $O(n^{1/4})$ by Kiri8128氏
けんちょんの競プロ精進記録：AOJ 3062 Product (RUPC 2019 day2-L)
- 丁寧に説明されている。僕の記事が分からなかった人でもこっちなら、理解できるはず。
数学の女王 ―歴史から見た数論入門
- 原始根の存在定理の証明を参考にした。