Skip to content →

Cipolla のアルゴリズム

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

フロベニウス写像

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

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

3つ例を挙げる:

  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}\) であるから拡大体は元の体に比べて大きくなる。
  3. \((\mathbb{Z}/13\mathbb{Z})(\sqrt{2},\sqrt{3}) = \{a+b\sqrt{2}+c\sqrt{3}+d\sqrt{6} \mid a,b,c,d\in \mathbb{Z}/13\mathbb{Z}\}\)
    • \(\sqrt{2},\sqrt{3} \not\in \mathbb{Z}/13\mathbb{Z}\) であるから拡大体は元の体に比べて大きくなる。四則演算が閉じるまで体を拡大するので、\(\sqrt{6}\) が新たに含まれる。

位数(要素数)\(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\) で定義する。

定理1:フロベニウス写像 \(\phi\) は自己同型である。すなわち、\(\phi\) は全単射で
\(\phi(x+y)=\phi(x)+\phi(y)\)
\(\phi(xy)=\phi(x)\phi(y)\)
が成り立つ。
証明:\(\phi(xy)=(xy)^p=x^py^p=\phi(x)\phi(y)\) である。\(0 < i < p\) のとき \({p \choose i} =0 \pmod p\) であるから標数 \(p\) と二項展開により \(\phi(x+y)=(x+y)^p=x^p+y^p\) である。以上より \(\phi\) は体の準同型である。

\(x\) の位数は \(p^N-1\) の約数であるから
\(\phi(x^{p^{N-1}}) = x^{p^N} = x\)
となり \(\phi\) は全単射である。

定理2:\(a \in \mathbb{F}_p \Leftrightarrow \phi(a)=a\) である。

\(\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\) は一通りに定まる。

証明:フェルマーの小定理より \(a\in \mathbb{F}_p\) ならば \(a^p=a\)。これは \(x^p=x\) の \(p\) 個の解である。\(p\) 次よりほかに解は存在せず、\(\phi\) は全単射であったから、\(a \in \mathbb{F}_p \Leftrightarrow \phi(a)=a\) となる。
定理3:任意の \(r \in \mathbb{F}_{p^N}\) に対して \(r\) を根に持つ唯一の \(\mathbb{F}_p\) 係数のモニック最小多項式 \(P\) が存在する。更に、\(\phi^d(r)=r\) となる最小の \(d\) を用いて \(P=\prod_{i=0}^{d-1}(x-\phi^i(r))\) と書ける。
証明:\(Q(x):=\sum a_i x^i\) とする。\(r\) が \(Q\) の根ならば定理1、2から \(\phi(Q(r))=Q(\phi(r))\) である。つまり \(r\) が根ならば \(\phi(r)\) も根である。

\(d\) の定義から集合 \(\{r,\phi(r),\ldots,\phi^{d-1}(r)\}\) は \(\phi\) を作用させても変化しない。また、\(\phi^i(r) \neq \phi^j(r)\) \((1\leq i < j < d)\) である。なぜなら、各 \(v\) から \(\phi(v)\) に有向辺を張ったグラフは \(\phi\) の全単射性から有向サイクルの集合であり、\(\phi^i(r) = \phi^j(r)\) ならば \(d\) の最小性に反するからである。従って、 \(P:=\prod_{i=0}^{d-1}(x-\phi^i(r)) = \sum c_ix^i\) に対して \(\sum c_ix^i=\sum \phi(c_i)x^i\) つまり全ての \(i\) に対して \(c_i=\phi(c_i)\) が成り立つ。定理2より \(P\) は \(\mathbb{F}_p\) 係数の多項式である。\(d\) の最小性と因数分解の一意性、\(\phi^i(r)\neq\phi^j(r)\) から \(P(x)\) は唯一のモニック最小多項式である。

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\) の2解 \(x=\alpha,\beta\) を取る。解と係数の関係より \(\alpha\beta=a\) である。\(\alpha,\beta \not\in \mathbb{F}_p\) とする。このとき \(x^2-2b+a\) は \(\mathbb{F}_p\) において既約多項式であるから定理3より \(\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\) であるから、連続で失敗する確率は指数関数的に減衰する。

定理4:\(a \neq 0\) が平方剰余だとする。\(b\in\mathbb{F}_p\) をランダムに取ったとき、\(b^2-a\) が非平方剰余になる確率は \(\frac{p-1}{2p}\) に等しい。
証明:\(\mathbb{F}_p\) における方程式 \(f=0\) の解の個数を\(N(f)\) 、ルジャンドル記号を \(\chi\) と表すことにする。ランダムに \(b \in \mathbb{F}_p\) を取ったとき \(b^2-a\) が非平方剰余になる確率は \(1-N(x^2-y^2-a)/2p\) に等しい。

\(a \neq 0\) のとき \(N(x^2-a) = 1+\chi(a)\) であるから、代入して
\begin{align}
&N(x^2-y^2-a)\\
=&(1+\chi(a))+(1+\chi(-a))+\sum_{A \in \mathbb{F}_p\setminus\{0,a\}} N(x^2-A)N(y^2+a-A)\\
=&3+\chi(-1)+\sum_{A \in \mathbb{F}_p\setminus\{0,a\}} (1+\chi(A))(1+\chi(A-a))\\
=&3+\chi(-1)+\sum_{A \in \mathbb{F}_p\setminus\{0,a\}} (1+\chi(A)+\chi(A-a)+\chi(A)\chi(A-a))\\
\end{align}
を得る。

\(\chi(a)=-1, \sum_{A \in \mathbb{F}_p} \chi(A)=0\) を代入して
\(\sum_{A \in \mathbb{F}_p\setminus\{0,a\}} (1+\chi(A)+\chi(A-a)) = p-3-\chi(-1)\)
を得る。

\(\chi(u)\chi(v)=\chi(uv),\chi(u)=\chi(u^{-1})\) が成り立つことを用いて
\begin{align}
&\sum_{A \in \mathbb{F}_p\setminus\{0,a\}}χ(A)χ(A-a)\\
=&\sum_{A \in \mathbb{F}_p\setminus\{0,a\}}χ(a^2)χ(a^{-1}A)χ(a^{-1}A-1)\\
=&\sum_{A’ \in \mathbb{F}_p\setminus\{0,1\}}χ(A’)χ(A’-1)\\
\end{align}
を得る。この式はルジャンドル記号に関するヤコビ和の \(\chi(-1)\) 倍であり、具体的に計算できる。\(χ(A’)χ(A’-1)=χ(A'(A’-1)^{-1})\) と変形しておく。\(x/(x-1)=y,x\neq0,1\) を \(x\) について解くと \(x=-y/(1-y),y\neq0,1\) である。つまり \(x\) が \(\mathbb{F}_p\setminus\{0,1\}\) を動くと \(y\) は \(\mathbb{F}_p\setminus\{0,1\}\) を動く。従って
\begin{align}
&\sum_{A’ \in \mathbb{F}_p\setminus\{0,1\}}χ(A'(1-A’)^{-1})\\
=&\sum_{A \in \mathbb{F}_p\setminus\{0,1\}}χ(A)\\
=&-1\\
\end{align}
を得る。
以上より \(1-N(x^2-y^2-a)/2p=(p-1)/2p\) である。

コード(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)}\) であるから冪指数は整数。

5次以上の方程式には解の公式が存在しないから、解の公式を使わずとも根が求まる既約多項式を用意する必要がある。そこで、\(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で呟いていたら、教えてもらった。

Published in アルゴリズム

Comments

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

CAPTCHA