\(p\) を奇素数として \(\mathbb{F}_p\) における平方根を \(O(\log p)\) で求める Cipolla のアルゴリズムを解説する。フロベニウス写像との関連を説明して、平方根以外の場合についても拡張する。

フロベニウス写像

まず Cipolla のアルゴリズムにおいて核となるアイディアであるフロベニウス写像について説明する。

体 \(k\) に要素 \(\alpha\) を追加して新たな体 \(K\) を作ることを添加といい、\(K=k(\alpha)\) と書く。\(k(\alpha)\) を \(k\) の拡大体と呼ぶ。

2つ例を挙げる：

1. 複素数 \(\mathbb{C}\) は実数 \(\mathbb{R}\) に \(\sqrt{-1}\) を添加してできた \(\mathbb{R}\) の拡大体。
2. \((\mathbb{Z}/13\mathbb{Z})(\sqrt{2}) = \{a+b\sqrt{2} \mid a,b\in \mathbb{Z}/13\mathbb{Z}\}\)
   • \(\sqrt{2} \not\in \mathbb{Z}/13\mathbb{Z}\) であるから拡大体は元の体に比べて大きくなる。

位数（要素数）\(q\) の有限体は同型を除き、ただ一通りに決まることが知られており、その有限体を \(\mathbb{F}_q\) と書く。このとき、位数 \(q\) は必ず素冪であり、逆に任意の素冪の有限体は必ず存在する。例えば、(2), (3) はそれぞれ \(\mathbb{F}_{13^2}, \mathbb{F}_{13^4}\) に同型である。\(\underbrace{1+\cdots+1}_{n \text{ 個}} = 0\) となる最小の正整数 \(n\) を標数という（存在しないならば 0 と定める）。有限体 \(\mathbb{F}_{p^N}\) の標数は \(p\) に等しいことが知られている。

有限体 \(\mathbb{F}_{p}\) の拡大体 \(\mathbb{F}_{p^N}\) に対して、フロベニウス写像 \(\phi : \mathbb{F}_{p^N} \to \mathbb{F}_{p^N}\) を \(\phi(r)=r^p\) で定義する。

\(\mathbb{F}_{p^N}\) は \(\mathbb{F}_p\) の拡大体として取っているので、この定理の \(\mathbb{F}_p\) は \(\mathbb{F}_{p^N}\) の部分体である。\(\mathbb{F}_{p^N}\) の乗法的単位元 \(1\) の生成する部分体（素体）と考えてもよい。部分体は元の体と同じ乗法的単位元を持つから、逆に \(\mathbb{F}_{p^N}\) が与えられたとき、それを拡大体とし得る部分体 \(\mathbb{F}_p\) は一通りに定まる。

Cipolla のアルゴリズム

(1) 平方根の場合

\(x^2=a \pmod p\) の解 \(x\) を求めたい。解が存在することと、\(a\) が平方剰余であることは同値である。これはオイラーの規準「\(a\) が平方剰余 \(\Leftrightarrow a^{(p-1)/2}=1 \pmod p\)」を使えば判定できる。以下、\(a\) は平方剰余であって解が存在すると仮定する。

必要ならば体 \(\mathbb{F}_p\) を拡大して \(\mathbb{F}_p\) 係数方程式 \(x^2-2b+a=0\) の２解 \(x=\alpha,\beta\) を取る。解と係数の関係より \(\alpha\beta=a\) である。\(\alpha,\beta \not\in \mathbb{F}_p\) とする。このとき \(x^2-2b+a\) は \(\mathbb{F}_p\) において既約多項式であるから定理３より \(\beta=\phi(\alpha)\) である。従って、\(\alpha^{(p+1)/2}\) が \(a\) の平方根となり、累乗法で \(O(\log p)\) で求まる。

\(x^2-2b+a=0\Leftrightarrow x=b \pm \sqrt{b^2-a} \) が \(\mathbb{F}_p\) に含まれないから \(b^2-a\) は非平方剰余である。このような \(b\) はどのぐらいの確率で取れるだろうか。次の定理に示す通り、その成功率は約 \(1/2\) であるから、連続で失敗する確率は指数関数的に減衰する。

コード（Java）は次のようになる。よすぽジャッジで Verify できる：

// a + b √base (MOD p) の形の数を扱うクラス
static class P {
	static long base;
	static long p;
	long a, b;

	public P(long a, long b) {
		this.a = a % p;
		this.b = b % p;
	}

	// 返り値 ： t との積
	P multiply(P t) {
		return new P(a * t.a + b * t.b % p * base, a * t.b + b * t.a);
	}

	// 返り値 ： n 乗
	P pow(long n) {
		if (n == 0)
			return new P(1, 0);
		P ret = multiply(this).pow(n / 2);
		if (n % 2 == 1)
			ret = ret.multiply(this);
		return ret;
	}
}

// 返り値: a^n mod p
long pow(long a, long n, long p) {
	return n != 0 ? pow(a * a % p, n / 2, p) * (n % 2 == 1 ? a : 1) % p : 1;
}

// 返り値： x^2 = a (mod p) の解 x
// そのような x が存在しなければ - 1 を返す。
long cipolla(long a, long p) {
	if (p == 2) return a;
	if (a == 0) return 0;
	if (pow(a, (p - 1) / 2, p) != 1) return -1;
	long b = 0;
	while (pow((b * b + p - a) % p, (p - 1) / 2, p) == 1)
		++b;
	P.base = (b * b + p - a) % p;
	P.p = p;
	P t = new P(b, 1).pow((p + 1) / 2);
	return t.a;
}

(2) 平方根以外の冪根の場合

\(x^n=a\pmod p\) の解 \(x\) を求めたい。解が存在することと、a が n 乗剰余であることは同値である。これはオイラーの規準「a が n 乗剰余 \(\Leftrightarrow\) \(a^{(p−1)/\gcd(p-1,n)}=1 \pmod p\)」を使えば判定できる。以下、a は n 乗剰余であって解が存在すると仮定する。

\(g := \gcd(n,p-1)\) とする。\(n/g\) と \((p-1)/g\) は互いに素であるから、逆元 \((n/g)^{-1} \pmod {(p-1)/g}\) が存在する。両辺を \((n/g)^{-1}_{\pmod{(p-1)/g}}\) 乗することで
\(x^n=a \pmod p \Leftrightarrow x^{\gcd(n,p-1)}=a^{(n/g)^{-1}_{\pmod{(p-1)/g}}} \pmod p\) を得る。従って \(n\) は \(p-1\) の約数と仮定しても一般性を失わない。

平方根の場合を拡張して、次のような方針を取る：

1. 定数項 \((-1)^na\) の \(n\) 次モニック既約多項式 \(f\) を見つける。
2. \(f\) の根 \(\theta\) を \(1\) つ見つける。
3. フロベニウス写像より \(\theta,\theta^p,\ldots,\theta^{p^{n-1}}\) は \(f\) の相異なる根。
4. 解と係数の関係より \(\theta^{(1+p+\ldots+p^{n-1})/n}\) が答え。\(1+p+\ldots+p^{n-1}=0 \pmod {(p-1)}\) であるから冪指数は整数。

５次以上の方程式には解の公式が存在しないから、解の公式を使わずとも根が求まる既約多項式を用意する必要がある。そこで、\(f = (x+b)^n-b^n+a\) とする。根は原始 \(n\) 乗根 \(\zeta\) を用いて \(x=-b+\zeta^i(b^n-a)^{1/n}\) \((i=0,1,\ldots,n-1)\)となる。既約多項式になるには \((b^n-a)^{1/n}\) が \(\mathbb{F}_p\) に含まれなければよい。つまり \(b^n-a\) が非 n 乗剰余になる \(b\) を取ってくる。

\(\mathbb{F}_p((b^n-a)^{1/n})\) で \(\theta^{(1+p+\ldots+p^{n-1})/n}\) の計算量は累乗法と FFT を使えば \(O(n^2\log(n)\log(p))\) である。工夫すればもう少しだけ速くなる。\(\theta^{(1+p+\ldots+p^{n-1})/n} = \theta(\prod_{i=0}^{n-2}\theta^{(1+p+\ldots+p^i)})^{(p-1)/n}\) とする。\((b^n-a)^{(p-1)/n}\) を原始 \(n\) 乗根 \(\zeta\) とすれば \(\theta^{p^i}\) は \(O(1)\) で求まる。従って \(\theta^{1+p+\ldots+p^i}\) を順番に掛け合わせてから \((p-1)/n\) 乗すれば \(O(n^2\log n+n\log n\log p)\) になる。

\(b^n-a\) が \(n\) 乗非剰余になる確率は原始根での表示を考えることでおよそ \(1-1/n\) になるだろうという推測される。実際、\(a≠0\) のとき \(N(x^n−a)=\sum_{i=0}^{n-1} a^{(p-1)i/d}\) であることを用いて、平方非剰余の場合と同様に頑張って計算すると \(1-1/n-O(n/\sqrt{p})\) が成功確率として得られる（らしい。理解してない）[2]。従って、\(n\) が \(\sqrt{p}\) より十分小さければよい。

参考文献

1. Paul Garrett氏のフロベニウス写像の授業資料
2. Williams, Kenneth S., and Kenneth Hardy. “A refinement of HC Williams’ 𝑞th root algorithm.” mathematics of computation 61.203 (1993): 475-483.
   • \(b^n-a\) が \(n\) 乗非剰余になる確率の計算が載っている……が理解できていない。
3. Harasawa, Ryuichi, Yutaka Sueyoshi, and Aichi Kudo. “Root computation in finite fields.” IEICE Transactions on Fundamentals of Electronics, Communications and Computer Sciences 96.6 (2013): 1081-1087.
4. sugarknri さん1, 2
5. Cipollaのアルゴリズムって他の冪根に拡張できるんかなーってTwitterで呟いていたら、教えてもらった。