#11 多項式の因数分解 その3 Berlekamp の方法
前回紹介した Cantor-Zassenhaus のアルゴリズムとは別の方法として、 Berlekamp の方法も紹介しておこう。 Berlekamp のアルゴリズムは、有限体上の多項式を多項式時間で因数分解するアルゴリズムの草分けであり、取り上げることに歴史的意義もあるだろう（と言いつつ筆者は原論文を確認したわけではないのだが……）。
Cantor-Zassenhaus の EDF アルゴリズムを適用するには事前に次数別因数分解を行っておく必要があったが、 Berlekamp の方法では無平方多項式を直接因数分解できる。 Berlekamp の方法もまた、確率的なアルゴリズムである。
前回の Cantor-Zassenhaus のアルゴリズムは割と可換環論的な話で済んだが、 Berlekamp のアルゴリズムでは（有限体上の）線形代数を使う。 とはいえ、線形代数は理系なら大学1年で扱うはずなので、特に問題にはならないだろう。
Berlekamp の方法
因数分解したい無平方モニック多項式を f とおく。 f は既約なモニック多項式 f_1,\ldots,f_l を使って f=f_1\cdots f_l と因数分解されたとする。 それぞれの f_i の次数を d_i:=\deg f_i とおく。 この状況では、中国剰余定理より、環の同型として \mathbb{F}_q[x]/(f) \cong\mathbb{F}_q[x]/(f_1)\times\cdots\times\mathbb{F}_q[x]/(f_l) \cong\mathbb{F}_{q^{d_1}}\times\cdots\times\mathbb{F}_{q^{d_l}} が成り立つ。 この対応を \chi\colon\mathbb{F}_q[x]/(f)\stackrel{\sim}{\longrightarrow}\mathbb{F}_{q^{d_1}}\times\cdots\times\mathbb{F}_{q^{d_l}} と書くことにする。
前回の EDF アルゴリズムでも似たような同型を考えたが、前回は各既約因子 f_i の次数は一定だった。 そのため、前回の EDF アルゴリズムは \mathbb{F}_q[x]/(f) の元を一つ選んで \frac{q^d-1}{2} 乗すれば、右辺の各 \mathbb{F}_{q^d} 成分が \pm1 になった。 だが、今回は各成分の次数が異なるため、 \frac{q^d-1}{2} 乗すれば良いというものではない。
それでも、どの \mathbb{F}_{q^{d_i}} も部分体として \mathbb{F}_q を含む。 よって、 \chi で \mathbb{F}_q\times\cdots\times\mathbb{F}_q に対応するような \mathbb{F}_q[x]/(f) の部分空間 \mathcal{B} を決定し、 \mathcal{B} から取った元を \frac{q-1}{2} 乗すれば、 \chi の行き先において \pm1 になる。 この \mathcal{B} は Berlekamp 部分代数 (Berlekamp subalgebra) と呼ばれる。 「部分空間」と書いているのは \mathbb{F}_q-代数としてである。
\begin{array}{rccc} \chi\colon & \mathbb{F}_q[x]/(f) & \stackrel{\sim}{\longrightarrow} & \mathbb{F}_{q^{d_1}}\times\cdots\times\mathbb{F}_{q^{d_l}} \\ & \cup & & \cup \\ & \mathcal{B} & \stackrel{\sim}{\longrightarrow} & \mathbb{F}_q\times\cdots\times\mathbb{F}_q \\ & a^{(q-1)/2} & \longmapsto & (\pm1,\ldots,\pm1) \end{array}
ではどうやって \mathcal{B} からランダムな元を取ってくるかという話になる。 これは \mathcal{B} の \mathbb{F}_q-線形空間としての基底 b_1,\ldots,b_l (b_i\in\mathbb{F}_q[x]/(f)) を決定してしまえば簡単で、それにランダムな係数をかけた線形結合 c_1b_1+\cdots+c_lb_l (c_i\in\mathbb{F}_q, 各 c_i はランダム) は \mathcal{B} のランダムな元となる。
というわけで残る問題は \mathcal{B} の基底 b_1,\ldots,b_l (b_i\in\mathbb{F}_q[x]/(f)) を決定する、ということになる。 なお、 \mathcal{B} が1次元であれば、 f の既約性が言える。
ここで、写像 \begin{array}{cccc} \beta\colon&\mathbb{F}_q[x]/(f)&\longrightarrow&\mathbb{F}_q[x]/(f) \\ & a & \longmapsto & a^q-a \end{array} を考えると、 \beta は \mathbb{F}_q-線形写像である（演習：これを示せ）。 さらに、 \beta の核 \operatorname{Ker}\beta は \mathcal{B} と一致する。
\mathbb{F}_q[x]/(f) の基底は、 n:=\deg f として \overline{x^{n-1}},\overline{x^{n-2}},\ldots,\overline{x},\overline{1} と取れる。 よって、線形写像 \beta\colon\mathbb{F}_q[x]/(f)\rightarrow\mathbb{F}_q[x]/(f) をこの基底について行列表示するということもできる。 \beta の表現行列について掃き出し法（ガウスの消去法）を実行してやれば、 \operatorname{Ker}\beta の基底を計算できる。 特に、 \operatorname{Ker}\beta の次元は、準同型定理より、 l=n-\operatorname{rank}\beta である。
アルゴリズムとしてまとめると、次のようになる：
- \beta の表現行列 M を計算する。 つまり、 \beta(x^j)=x^{qj}-x^j\equiv \sum_{i=0}^{n-1} m_{ij}x^i\bmod f \quad (0\le j<n) とおいて m_{ij}\in\mathbb{F}_q を求める。 この時の x^{qj} の計算は、先に x^q\bmod f を計算し、それを j 乗すると良い。
- 行列 M にガウスの消去法を適用し、 \operatorname{Ker}\beta の次元 l と基底 b_1,\ldots,b_l を計算する。 l=1 であれば、 f はすでに既約なのでアルゴリズムを終了する。
- c_1,\ldots,c_l (c_i\in\mathbb{F}_q) をランダムに選び、 a:=c_1b_1+\cdots+c_lb_l\in\mathbb{F}_q[x]/(f) とおく。
- a^{(q-1)/2}\in\mathbb{F}_q[x]/(f) を計算する。
- f^*:=\gcd(f,a^{(q-1)/2}-1) を計算する。
- 運が良ければ、 f^*\ne 1,f^*\ne f となり、 f^* は f の自明でない因数を与える。（成功）
- f^*=1 または f^*=f の場合は、アルゴリズムは失敗である。
定理． 上記のアルゴリズムは期待通りに動作し、 1/2 以上の確率で成功する。
証明． 正当性については、すでに書いた通りである。
成功する確率は、つまり（一様）ランダムに選んだ a\in\mathcal{B} に対して \chi(a^{(q-1)/2}-1)\in\mathbb{F}_q\times\cdots\times\mathbb{F}_q の各成分がどうなるかである。 前回と同様の議論により、 (0,\ldots,0) が出てくる確率は 1/2^l, (-2,\ldots,-2) が出てくる確率は 1/2^l である。 よって失敗する確率は 2/2^l=1/2^{l-1} で、成功する確率は 1-1/2^{l-1}\ge 1/2 である。
計算例
簡単な例で手計算してみよう。 q=3, f=(x+2)(x^2+x+2)=x^3+x+1\in\mathbb{F}_3[x] とする。 \beta は \begin{array}{cccc} \beta\colon&\mathbb{F}_3[x]/(f) & \longrightarrow & \mathbb{F}_3[x]/(f) \\ & a & \longmapsto & a^3-a \\ & \overline{1} & \longmapsto & \overline{0} \\ & \overline{x} & \longmapsto & \overline{x^3-x}=\overline{x+2} \\ & \overline{x^2} & \longmapsto & \overline{x^6-x^2}=\overline{2x+1} \end{array} なので、 \mathbb{F}_3[x]/(f) の基底 \overline{1},\overline{x},\overline{x^2} に関する表現行列 M は \begin{pmatrix} 0 & 2 & 1 \\ 0 & 1 & 2 \\ 0 & 0 & 0 \end{pmatrix} となる。 これの下に単位行列をくっつけて（拡大係数行列）、 \begin{pmatrix} 0 & 2 & 1 \\ 0 & 1 & 2 \\ 0 & 0 & 0 \\ %\hline 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} 列基本変形を行うと \begin{pmatrix} 2 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \\ %\hline 0 & 1 & 0 \\ 1 & 0 & 1 \\ 0 & 0 & 1 \end{pmatrix} を得る。 これの右2列の下半分 \begin{pmatrix}1\\0\\0\end{pmatrix}, \begin{pmatrix}0\\1\\1\end{pmatrix} が Berlekamp 部分代数 \mathcal{B} の基底である。 わかりやすいように x とかを使って書くと \begin{aligned} b_1&=1\cdot\overline{1}+0\cdot\overline{x}+0\cdot\overline{x^2}=\overline{1}, \\ b_2&=0\cdot\overline{1}+1\cdot\overline{x}+1\cdot\overline{x^2}=\overline{x^2+x} \end{aligned} となる。 この時点で、 f は2個の既約因数からなるということがわかる。
あとはこいつらのランダムな線形結合を計算する。 例えば c_1=2,c_2=0 として a=c_1b_1+c_2b_2=2\cdot\overline{1} としよう。 a^{(q-1)/2}-1=a^1-1=0 なので、 f^*=\gcd(f,a^{(q-1)/2}-1)=\gcd(f,0)=f である。 これは f の自明な因数なので、アルゴリズムは失敗である。
気を取り直して、別の c_i で試してみよう。 今度は c_1=2,c_2=1 として a=c_1b_1+c_2b_2=\overline{x^2+x+2} としよう。 a^{(q-1)/2}-1=a^1-1=\overline{x^2+x+1} なので、 f^*=\gcd(f,a^{(q-1)/2}-1)=\gcd(f,x^2+x+1)=x+2 である。 これは f の非自明な因数なので、アルゴリズムは成功である。 あとは x+2 と f/(x+2)=x^2+x+2 に再帰的にアルゴリズムを適用すれば、 f の完全な因数分解を得る。
実装
今回の実装は src/Numeric/AlgebraicReal/Factor/Berlekamp.hs に書き込む：
module Numeric.AlgebraicReal.Factor.Berlekamp where
import Prelude hiding (toInteger)
import Data.FiniteField
import System.Random
import qualified Data.Vector as V
import Control.Monad
import Control.Monad.State
import Data.Array
import Data.List
import Numeric.AlgebraicReal.UniPoly
import Numeric.AlgebraicReal.Factor.CantorZassenhaus (powMod)
type Mat a = Array (Int,Int) a前回実装した powMod 関数を使うので、前回のモジュールを import する（powMod は UniPoly モジュールに書いた方が良かったかもしれない）。
行列は Data.Array の2次元配列で表すことにする。 今回の用途では、行列の足し算や掛け算は特に必要ない。
\beta の表現行列を計算する関数 betaMat を作る（名前が雑だ）：
-- beta : x \mapsto x^q-x の表現行列を計算する。
betaMat :: (FiniteField k, Eq k) => UniPoly k -> Mat k
betaMat f
= let q = order (leadingCoefficient f)
n = degree' f
xq = powMod ind q f -- x^q mod f
ys = iterate (\y -> (y * xq) `modP` f) 1
in array ((0,0),(n-1,n-1))
[ ((i,j), mij - d)
| (j,xqj) <- zip [0..n-1] ys
, (i,mij) <- zip [0..n-1] (V.toList (coeff xqj) ++ repeat 0)
, let d | i == j = 1
| i /= j = 0
]変数 ys は 1,x^q \bmod f,x^{2q} \bmod f,x^{3q} \bmod f,\ldots の形のリストとなる。 それの係数を使って行列を作る。
列基本変形を行う関数 elim を適当に実装する：
-- 列基本変形を行う
-- 非 0 な列の index と、変形後の行列を返す
elim :: (Fractional k, Eq k) => Int -> Mat k -> ([Int],Mat k)
elim imax m = loop [] i0 m
where
b@((i0,j0),(_,jn)) = bounds m
loop seen i m
| i > imax = (seen,m)
| otherwise = case [k | k <- [j0..jn] \\ seen, m!(i,k) /= 0] of
[] -> loop seen (i+1) m
k:_ -> loop (k:seen) (i+1) $ array b
[ ((i',j),v)
| (i',j) <- indices m
, let v | j /= k = m!(i',j) - m!(i',k)*m!(i,j)/m!(i,k)
| j == k = m!(i',k)
]Berlekamp 部分代数の基底を計算する関数 berlekampBasis を書く：
-- Input: nonzero, monic, squarefree
berlekampBasis :: (FiniteField k, Eq k) => UniPoly k -> [UniPoly k]
berlekampBasis f
= let n = degree' f
m = betaMat f
-- m' は m の拡大係数行列（m の下に単位行列をくっつけたもの）
m' = array ((0,0),(2*n-1,n-1)) (assocs m ++ [ ((i+n,j),v)
| (i,j) <- indices m
, let v | i == j = 1
| i /= j = 0
])
(ix,l) = elim (n-1) m'
in [fromCoeff (V.fromList [l!(n+i,j) | i<-[0..n-1]]) | j <- [0..n-1] \\ ix]Berlekamp のアルゴリズムにより、非自明な因数を見つけ出す関数 berlekampOne を書く：
-- Input: nonzero, monic, squarefree
-- Returns Nothing if f is already irreducible
berlekampOne :: (FiniteField k, Eq k, Random k, RandomGen g) => UniPoly k -> g -> (Maybe (UniPoly k), g)
berlekampOne f = let bs = berlekampBasis f
r = length bs
q = order (leadingCoefficient f)
-- challenge :: State g (UniPoly k)
-- f の非自明な因数を見つけて返す
challenge = do
cs <- sequence $ replicate r $ state random -- F_q の元をランダムに r 個選ぶ
let a = sum (zipWith scaleP cs bs) -- 線形結合を作る
w = powMod a ((q-1) `div` 2) f -- a^((q-1)/2) mod f の計算
f' = toMonic $ gcdP f (w-1)
if f' /= 1 && f' /= f
then return f' -- f' は f の非自明な因数である
else challenge -- f' は非自明な因数ではなかった；係数を選ぶところからやり直し
in runState (if r == 1
then return Nothing -- f is already irreducible
else fmap Just challenge)
berlekampOneIO :: (FiniteField k, Eq k, Random k) => UniPoly k -> IO (Maybe (UniPoly k))
berlekampOneIO = getStdRandom . berlekampOneこの関数は、入力 f がすでに既約だったら Nothing を返し、そうでない場合は非自明な因数が見つかるまでアルゴリズムを繰り返してそれを返す。
前回と同様、グローバルな StdGen を利用する版 (berlekampOneIO) も用意しておく。
これを使って、与えられた無平方モニック多項式を完全に因数分解する関数 factorBerlekamp を書く：
-- Input: nonzero, monic, squarefree
factorBerlekamp :: (FiniteField k, Eq k, Random k, RandomGen g) => UniPoly k -> g -> ([UniPoly k], g)
factorBerlekamp f = runState (loop f [])
where loop f acc | degree' f == 0 = return acc
| otherwise = do
mf' <- state (berlekampOne f)
case mf' of
Nothing -> return (f:acc) -- f is already irreducible
Just f' -> do
-- f' is a nontrivial factor of f
acc' <- loop f' acc
loop (f `divP` f') acc'
factorBerlekampIO :: (FiniteField k, Eq k, Random k) => UniPoly k -> IO [UniPoly k]
factorBerlekampIO = getStdRandom . factorBerlekampここでも、グローバルな StdGen を利用する版 (factorBerlekampIO) も用意しておく。
algebraic-real.cabal には今回の Numeric.AlgebraicReal.Factor.Berlekamp モジュールを追加し、依存パッケージとして array を追加する：
exposed-modules: Numeric.AlgebraicReal.UniPoly
, ...
, Numeric.AlgebraicReal.Factor.CantorZassenhaus
, Numeric.AlgebraicReal.Factor.Berlekamp
build-depends: base >= 4.7 && < 5
, vector
, finite-field
, random
, mtl
, array実装を使う
今回書いた実装を REPL で試してみよう。 stack repl を起動し、前回同様のコマンドを打ち込む：
$ stack repl
...
...> :set prompt "> "
> :m + System.Random Data.FiniteField
> :set -XDataKinds
> setStdGen (mkStdGen 4539814)さあ始めよう。
> let x = ind
> factorBerlekampIO (x^3+x+1 :: UniPoly (PrimeField 3))
[UniPoly [2,1,1],UniPoly [2,1]]さっきの手計算でも使った x^3+x+1\in\mathbb{F}_3[x] を因数分解させたところ、 x^3+x+1=(x+2)(x^2+x+2) という分解が得られた。 良さそうである。
次に、前回もやった x^7-1\in\mathbb{F}_{37}[x] を因数分解してみよう。
> factorBerlekampIO (x^7-1 :: UniPoly (PrimeField 37))
[UniPoly [36,28,29,1],UniPoly [36,8,9,1],UniPoly [36,1]]因数分解として x^7-1=(x^3+29x^2+28x+36)(x^3+9x^2+8x+36)(x+36) という結果が得られた。 前回のアルゴリズムでは先に次数別因数分解を行なっているため、因数の次数が低い順に並んで出てきたが、今回はそんなことはない。
もちろん、乱数を使っているので、因数の順番は試行ごとに変わる。
> factorBerlekampIO (x^7-1 :: UniPoly (PrimeField 37))
[UniPoly [36,28,29,1],UniPoly [36,8,9,1],UniPoly [36,1]]
> factorBerlekampIO (x^7-1 :: UniPoly (PrimeField 37))
[UniPoly [36,1],UniPoly [36,8,9,1],UniPoly [36,28,29,1]]2回目は1回目と同じ順番で出てきたが、3回目は順番が変わった。
他の例も試してみる。
> factorBerlekampIO (x^7-1 :: UniPoly (PrimeField 11))
[UniPoly [10,4,5,1],UniPoly [10,1],UniPoly [10,6,7,1]]
> factorBerlekampIO (x^20-1 :: UniPoly (PrimeField 11))
[UniPoly [1,0,1],UniPoly [2,1],UniPoly [5,0,1],UniPoly [10,1],UniPoly [5,1],UniPoly [4,1],UniPoly [4,0,1],UniPoly [7,1],UniPoly [3,0,1],UniPoly [3,1],UniPoly [8,1],UniPoly [1,1],UniPoly [6,1],UniPoly [9,1],UniPoly [9,0,1]]
> factorBerlekampIO (x^50-1 :: UniPoly (PrimeField 11))
[UniPoly [5,1],UniPoly [4,0,0,0,0,1],UniPoly [7,0,0,0,0,1],UniPoly [8,1],UniPoly [8,0,0,0,0,1],UniPoly [7,1],UniPoly [2,1],UniPoly [10,1],UniPoly [3,1],UniPoly [2,0,0,0,0,1],UniPoly [6,1],UniPoly [9,0,0,0,0,1],UniPoly [5,0,0,0,0,1],UniPoly [1,1],UniPoly [9,1],UniPoly [4,1],UniPoly [6,0,0,0,0,1],UniPoly [3,0,0,0,0,1]]
> factorBerlekampIO (x^100-1 :: UniPoly (PrimeField 11))
[UniPoly [7,1],UniPoly [7,0,0,0,0,1],UniPoly [8,1],UniPoly [8,0,0,0,0,1],UniPoly [1,1],UniPoly [6,1],UniPoly [9,0,0,0,0,1],UniPoly [4,0,1],UniPoly [4,1],UniPoly [4,0,0,0,0,0,0,0,0,0,1],UniPoly [2,0,0,0,0,1],UniPoly [10,1],UniPoly [1,0,1],UniPoly [9,0,0,0,0,0,0,0,0,0,1],UniPoly [5,0,0,0,0,0,0,0,0,0,1],UniPoly [5,0,1],UniPoly [4,0,0,0,0,1],UniPoly [3,0,0,0,0,0,0,0,0,0,1],UniPoly [6,0,0,0,0,1],UniPoly [5,1],UniPoly [5,0,0,0,0,1],UniPoly [3,1],UniPoly [3,0,1],UniPoly [3,0,0,0,0,1],UniPoly [2,1],UniPoly [9,1],UniPoly [9,0,1]]x^{50}-1 や x^{100}-1 といった、比較的次数の高い多項式も、実用的な時間で因数分解できた。 係数環が違うといはいえ、 #8 で扱った Kronecker の方法では考えられないことである。
前回と今回扱った有限体上の多項式の因数分解アルゴリズムを利用し、次回は、いよいよ整数係数多項式の因数分解を行う。