#10 多項式の因数分解 その2 Cantor-Zassenhaus の方法

投稿日:2018年1月8日,最終更新:2025年4月5日

今回は、有限体上の因数分解アルゴリズムの一つである、 Cantor-Zassenhaus の方法を取り扱う。

この記事では、特に断らない限り、多項式の係数は有限体 \mathbb{F}_q で考える。 ただし、時々、その拡大体 \mathbb{F}_{q^d} が登場するので注意してほしい。

Cantor-Zassenhaus の方法

Cantor-Zassenhaus の方法は、無平方モニック多項式を対象にし、2つのステップからなる。 (無平方とは限らない多項式を因数分解する際は、先に無平方分解をしておく必要がある)

最初のステップは次数別因数分解 (DDF; distinct-degree factorization) である。 入力は無平方モニック多項式 f\in\mathbb{F}_q[x] とすると、出力は、 f=f_1f_2\cdots f_k を満たす多項式の列 f_1,f_2,\ldots,f_k である。 ただし、各 f_if_i=f_{i1}\cdots f_{il_i}, \deg f_{ij}=i, f_{ij} は既約モニック多項式、 となる。 なお、各 f_i について、 \deg f_i=i であればこの時点で f_i が既約であることがわかる。

次のステップは等次数因数分解 (EDF; equal-degree factorization) である。 このアルゴリズムは、次数別因数分解で得られた各 f_i の自明でない因数 f_{i}^* を見つけ出す。 このアルゴリズムは確率的である。 つまり、ある確率で f_i の自明でない因数を見つけ出すが、ある確率で失敗する(自明な因数しか見つけられない)。 このアルゴリズムを再帰的に適用すれば、完全な因数分解 f_i=f_{i1}\cdots f_{il_i} が得られる。

(なお、 D. G. Cantor と H. Zassenhaus の提案したアルゴリズムは EDF の方であり、 DDF の方はそれ以前から知られていたようである。 しかし、 Cantor-Zassenhaus のアルゴリズムを適用するにはどのみち DDF が必要なため、ここでは DDF もアルゴリズムの一部として紹介する。 また、 EDF を「等次数因数分解」とするのはこの連載独自の訳である。)

次数別因数分解

次数別因数分解は、フェルマーの小定理に立脚する。

定理x^{q^n}-x は、その次数が n の約数であるような既約なモニック多項式全ての積である。 つまり、以下の3つが成り立つ:

  1. x^{q^n}-x は無平方多項式である。
  2. d\mid n なる自然数 d について、次数 d の既約多項式 gx^{q^n}-x の因数である。
  3. x^{q^n}-x の既約因数 g の次数を d とすると、 dn の約数である。

証明. 1について: f(x)=x^{q^n}-x について f'(x)=-1 なので、 \gcd(f,f')=1 である。 よって、 f は無平方である。

2について: \mathbb{F}_q[x]/(g)\cong\mathbb{F}_{q^d} である。 x\in\mathbb{F}_q[x] の、 \mathbb{F}_{q^d} における像を \alpha とおく。 \begin{array}{ccccc} \mathbb{F}_q[x] & \longrightarrow & \mathbb{F}_q[x]/(g) & \stackrel{\sim}{\longrightarrow} & \mathbb{F}_{q^d} \\ x & \longmapsto & \overline{x} & \longmapsto & \alpha \end{array} フェルマーの小定理より、 \alpha^{q^d}=\alpha である。 d\mid n より、 \alpha^{q^n}=\alpha^{q^{d+d+\cdots+d}}=\alpha^{q^d q^d\cdots q^d}=\Bigl(\cdots\bigl((\alpha^{q^d})^{q^d}\bigr)^{\cdots}\Bigr)^{q^d}=\alpha である。 よって \alpha^{q^n}-\alpha=0 であり、これは \mathbb{F}_q[x]/(g) で見ると \overline{x}^{q^n}-\overline{x}=0 で、これは \mathbb{F}_q[x] で見ると x^{q^n}-x\in(g) を意味する。 よって、 x^{q^n}-xg で割り切れる。

3について: まず、 \mathbb{F}_{q^n} は多項式 x^{q^n}-x の最小分解体である。 つまり、 \mathbb{F}_{q^n}\mathbb{F}_{q}x^{q^n}-x の全ての根を添加したものである。 このことは、次の3つからわかる:\mathbb{F}_{q^n} の元は全て多項式 x^{q^n}-x の根であり、 x^{q^n}-x は無平方であり、 x^{q^n}-x の根の個数は最大でも q^n 個である。

次に、\mathbb{F}_{q^d}\mathbb{F}_{q^n} の部分体であることを言う。 同型 \mathbb{F}_q[x]/(g)\cong\mathbb{F}_{q^d} を使い、 \mathbb{F}_{q^d} の元 \alpha を、次のように取る: \begin{array}{ccccc} \mathbb{F}_q[x] & \longrightarrow & \mathbb{F}_q[x]/(g) & \stackrel{\sim}{\longrightarrow} & \mathbb{F}_{q^d} \\ x & \longmapsto & \overline{x} & \longmapsto & \alpha \end{array} gx^{q^n}-x の因数なので、 \alphax^{q^n}-x の根である。 従って、 \mathbb{F}_q(\alpha)=\mathbb{F}_{q^d}\mathbb{F}_{q^n} の部分体である。

\mathbb{F}_{q^n}\mathbb{F}_{q^d} 上の線形空間とみなした時の次元を m:=\dim_{\mathbb{F}_{q^d}}\mathbb{F}_{q^n} とすると、 元の個数について q^n=\bigl(q^d\bigr)^m が成り立つので、 n=md である。 つまり、 dn の約数である。

例:n=1 の場合。 1 の約数は 1 だけだから、 x^q-x は次数1の既約なモニック多項式(つまり、任意の1次式 x-a)全ての積である: x^q-x=\prod_{a\in\mathbb{F}_q} (x-a) これは、フェルマーの小定理の帰結である。

例:n=2 の場合。 x^{q^2}-x は 次数1の既約なモニック多項式全てと、次数2の既約なモニック多項式全ての積である: x^{q^2}-x=\left(\prod_{a\in\mathbb{F}_q} (x-a)\right)\left(\prod_{a,b\in\mathbb{F}_q,x^2+ax+b\text{ は既約}} (x^2+ax+b)\right)

例:q=3 の場合 x^3-x=x(x-1)(x-2) \mathbb{F}_3 係数の2次の既約なモニック多項式は x^2+1, x^2+x+2, x^2+2x+2 の3つであるので、 x^9-x=x(x-1)(x-2)(x^2+1)(x^2+x+2)(x^2+2x+2) となる(演習:この等式が正しいことを、右辺を展開することにより確かめよ。係数は mod 3 で考える)。

この定理に基づいて、次数別因数分解のアルゴリズムを構成できる:

  1. g_0:=f とおく。
  2. k=1,2,\ldots に関して、順番に i.〜iii. を実行する
    1. k 次の既約因数を全て掃き出すために、次を計算する: f_k:=\gcd(g_{k-1},x^{q^k}-x)f に含まれる k 次の既約因数全ての積となる。
    2. g_{k-1}f_k で割ったものが g_k となる: g_k:=g_{k-1}/f_k g_k は、 k 次の既約因数が全て掃き出され、 k+1 次以上の既約因数のみからなる(f は無平方と仮定していることに注意せよ)。
    3. もしも g_k が 1 となれば、既約因数を掃き出し終えたということなので、アルゴリズムを停止する。
  3. f=f_1f_2\cdots f_k であり、 f_i の既約因数は全て i 次である。

x^{q^k}-x の計算であるが、どうせ f (の因数)との GCD にしか使わないので、計算の際に \mod f で計算して良い。 具体的には、 u_{k-1}= x^{q^{k-1}}\bmod{f} に対して \mod fk 乗したものを u_k:=u_{k-1}^k\bmod{f} とし、 f_k の計算の際に f_k:=\gcd(g_{k-1},u_k-x) とする。

終了条件であるが、 g_k の次数 d2(k+1) 未満であれば f_{k+1}=\cdots=f_{d-1}=1,f_d=g_k としてアルゴリズムを打ち切って良い。

等次数因数分解

ここで紹介するアルゴリズムは、標数が奇数である必要がある。

入力となる多項式を h とする(これは先のアルゴリズムで得られた f_i のことである)。 h は無平方なモニック多項式であり、 d 次の既約多項式 h_jl 個の積として h=h_1\cdots h_l と書ける(つまり \deg h=dl である)。

成功した場合の出力は、 h の自明でない因数 h^* である。 この場合、 h^*h/h^* に再帰的に等次数因数分解を適用すれば、いつか \deg h=d となって完全な因数分解を得る。

数学的な説明の前に、アルゴリズムを説明してしまおう。

  1. u\in\mathbb{F}_q[x] を、 u\ne 0\deg u<\deg h となるようにランダムに選ぶ。(このような選び方は q^{dl}-1 通りある)
  2. \gcd(h,u)\ne 1 となった場合は \gcd(h,u)h の自明でない因数なので、 h^*=\gcd(h,u) と置いて、アルゴリズムを終了する(成功)。
  3. そうでない場合、つまり \gcd(h,u)=1 の場合。 mod hu^{(q^d-1)/2} を計算する。 それを使って \gcd(h,u^{(q^d-1)/2}-1)\ne 1,h となれば h^*=\gcd(h,u^{(q^d-1)/2}-1) と置いてアルゴリズムは成功、そうでなければアルゴリズムは失敗である。

では数学的な説明に入ろう。

定理. 上記のアルゴリズムは期待通りに動作し、 1/2 以上の確率で成功する。

証明h=h_1\cdots h_l の状況では、中国剰余定理より、環の同型として \mathbb{F}_q[x]/(h) \cong\mathbb{F}_q[x]/(h_1)\times\cdots\times\mathbb{F}_q[x]/(h_l) \cong\mathbb{F}_{q^d}\times\cdots\times\mathbb{F}_{q^d} が成り立つ(前半が中国剰余定理、後半は直積の各成分について有限体の性質を使った)。 この対応を \varphi\colon\mathbb{F}_q[x]/(h)\stackrel{\sim}{\longrightarrow}\mathbb{F}_{q^d}\times\cdots\times\mathbb{F}_{q^d} と書き、さらにその第 j 射影を \varphi_j\colon\mathbb{F}_q[x]/(h)\to\mathbb{F}_{q^d} と書くことにする。 (もちろん、我々は個々の h_j を知らないので、この \varphi\varphi_j は実際の計算アルゴリズムには利用できない。しかし、アルゴリズムが成功する確率の計算には利用できる)

ここで、 \mathbb{F}_q[x]dl 次未満の多項式を一つ選ぶことは、環 \mathbb{F}_q[x]/(h) の元を一つ選ぶことに相当する。 そこで、 u\in\mathbb{F}_q[x]u\ne 0, \deg u<dl となるように取ったとする。 \begin{array}{ccccc} \mathbb{F}_q[x]&\to&\mathbb{F}_q[x]/(h)&\stackrel{\sim}{\longrightarrow}&\mathbb{F}_{q^d}\times\cdots\times\mathbb{F}_{q^d} \\ u&\mapsto&\overline{u}&\mapsto&\varphi(\overline{u})=(\varphi_1(\overline{u}),\ldots,\varphi_l(\overline{u})) \end{array} を考える。 この第 j 成分 \varphi_j(\overline{u})u_j とおく。

ここでもし、ある j について u_j=0 となり、別の j' について u_{j'}\ne0 となったならば、 u から h の自明でない因数を構成できる。 具体的には、 uh_j で割り切れ、 h_{j'} で割り切れない。 つまり、 \gcd(h,u)h_j を因数として含み、 h_{j'} を含まない h の因数である。

しかしそのようなラッキーはなかなか起きない。 そこで、全ての j について u_j\ne 0 となった場合を考える。 これは \gcd(h,u)=1 の場合に相当する。

u_j\in\mathbb{F}_{q^d} はフェルマーの小定理 u_j^{q^d}=u_j を満たす。 さらに、 u_j\ne 0 であれば、 u_j^{q^d-1}=1 である。 q が奇数であるという仮定より、その平方根 u_j^{(q^d-1)/2} を考えることができ、これは \pm1 である: u_j^{(q^d-1)/2}=\pm 1

\mathbb{F}_{q^d} の 0 でない元であって、 \frac{q^d-1}{2} 乗が 1 となるようなものはちょうど半分、つまり \frac{q^d-1}{2} 個存在し、 \frac{q^d-1}{2} 乗が -1 となるようなものは残りの半分、つまり \frac{q^d-1}{2} 個存在する(後述の補題)。

よって、 u をランダムに選べば、 u_j^{(q^d-1)/2} が 1 となる確率と -1 となる確率はちょうど半々(1/2 ずつ)である(u をランダムに選んだ時に u_j がランダムに選ばれるのというのは、中国剰余定理の帰結である)。

特に、 \varphi(\overline{u})^{(q^d-1)/2}=(u_1^{(q^d-1)/2},\ldots,u_l^{(q^d-1)/2})(\pm1,\ldots,\pm1) となり、しかも成分の符号が全て一致((1,\ldots,1) または (-1,\ldots,-1) となる)する確率は \frac{2}{2^{l}}=2^{1-l} である。 \begin{array}{ccc} \mathbb{F}_q[x]/(h)&\stackrel{\sim}{\longrightarrow}&\mathbb{F}_{q^d}\times\cdots\times\mathbb{F}_{q^d} \\ \overline{u}^{(q^d-1)/2}&\mapsto&\varphi(\overline{u})^{(q^d-1)/2}=(\pm1,\ldots,\pm1) \\ \overline{u}^{(q^d-1)/2}-1&\mapsto&\varphi(\overline{u})^{(q^d-1)/2}-1=(0\text{ or }{-2},\ldots,0\text{ or }{-2}) \end{array} さらに、 \varphi(\overline{u})^{(q^d-1)/2} に 1 を足すと、そこそこの確率で「ある j について第 j 成分が 0 となり、別の j' について第 j' 成分が非 0 となる」 (0\text{ or }2,\ldots,0\text{ or }2) を得る。 こうなればアルゴリズムは成功で、 \gcd(h,u^{(q^d-1)/2}+1) を計算すれば h の自明でない因数が得られる。

アルゴリズムが失敗、つまり (0,\ldots,0)(2,\ldots,2) が得られた場合は、 \gcd(h,u^{(q^d-1)/2}+1) はそれぞれ h または 1 となる。

補題q は奇素数のべきとする。 \mathbb{F}_{q}^\times の元のうち、 (q-1)/2 乗が 1 となるものと -1 になるものはちょうど半分ずつ((q-1)/2 個ずつ)存在する。

この補題は、次の定理の系として得られる:

補題q は素数べき、 dq-1 の約数とする。 \mathbb{F}_{q}^\times の元のうち、 \frac{q-1}{d} 乗が 1 となるものはちょうど \frac{q-1}{d} 個存在する。

証明. 乗法群 \mathbb{F}_q^\times 上の群の自己準同型 f, g\begin{array}{cccc} &\mathbb{F}_q^\times&\to&\mathbb{F}_q^\times \\ f\colon&x&\mapsto&x^d \\ g\colon&x&\mapsto&x^{(q-1)/d} \end{array} と定める。

\operatorname{Ker}f=\{x\in\mathbb{F}_q^\times\mid x^d=1\}, \operatorname{Ker}g=\{x\in\mathbb{F}_q^\times\mid x^{(q-1)/d}=1\} より、元の個数について \#(\operatorname{Ker}f)\le d, \#(\operatorname{Ker}g)\le \frac{q-1}{d} が成り立つ。

f\circ g=g\circ f=1 であるので、 \operatorname{Im}g\subset\operatorname{Ker}f, \operatorname{Im}f\subset\operatorname{Ker}g である。

準同型定理より、 \mathbb{F}_q^\times/\operatorname{Im}g\cong \operatorname{Ker}g であり、 元の個数について q-1=\#(\mathbb{F}_q^\times)=\#(\operatorname{Im}g)\cdot \#(\operatorname{Ker}g)\le\#(\operatorname{Ker}f)\cdot\#(\operatorname{Ker}g)\le d\cdot\frac{q-1}{d}=q-1 が成り立つ。 よって \#(\operatorname{Ker}f)=d, \#(\operatorname{Ker}g)=\frac{q-1}{d} である。

この \#(\operatorname{Ker}g)=\frac{q-1}{d} が示したいことであった。

実装

乱数の生成には random パッケージの System.Random モジュールを使うことにする。

今回のソースファイルは src/Numeric/AlgebraicReal/Factor/CantorZassenhaus.hs とする:

module Numeric.AlgebraicReal.Factor.CantorZassenhaus 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 GHC.TypeLits (KnownNat)
import Numeric.AlgebraicReal.UniPoly

まず、よくある powMod アルゴリズムを実装する:

powMod :: (Eq k, Fractional k) => UniPoly k -> Integer -> UniPoly k -> UniPoly k
powMod _ 0 _ = 1
powMod a n m = let a' = a `modP` m in loop a' (n-1) a'
  where loop a 0 acc = acc
        loop a n acc = case n `quotRem` 2 of
          (n',0) -> loop ((a * a) `modP` m) n' acc
          (n',1) -> loop ((a * a) `modP` m) n' ((acc * a) `modP` m)

【1月13日 修正】powMod の実装が不適切だったので修正した。

最初は次数別因数分解である。 こちらのアルゴリズムは乱数を使わずに書ける。

-- Input: nonzero, monic, squarefree
distinctDegreeFactorization :: (FiniteField k, Eq k) => UniPoly k -> [(Int,UniPoly k)]
distinctDegreeFactorization f = loop 1 ind f
  where loop k u g | degree' g == 0 = []
                   | degree' g < 2*k = [(degree' g,g)]
                   | f' == 1 = loop (k+1) u' g'
                   | otherwise = (k,f') : loop (k+1) u' g'
          where u' = powMod u q f
                f' = toMonic $ gcdP g (u' - ind)
                g' = g `divP` f'
        q = order (leadingCoefficient f)

Haskell は純粋なので、 Haskell の流儀で乱数を扱うには乱数生成器を持って回る形になる。 まず、「n 次未満の多項式をランダムに生成する」関数を書く。 係数の型は Random クラスのインスタンスである必要がある。

randomPolyOfDegreeLessThan :: (Eq k, Num k, Random k, RandomGen g) => Int -> g -> (UniPoly k, g)
randomPolyOfDegreeLessThan n g = case takeR n [] g of
                                   (xs,g) -> (fromCoeff (V.fromList xs), g)
  where takeR 0 acc g = (acc,g)
        takeR n acc g = case random g of
                          (a,g) -> takeR (n-1) (a:acc) g

等次数因数分解の、「自明でない因数を一つ見つけ出す」操作を書く。 引数 h の既約因数の次数は全て d であると仮定する。 返り値は Maybe 型にくるみ、自明な因数しか見つからなかった場合は Nothing を返す。

-- Input: nonzero, monic, squarefree, equal-degree, reducible
equalDegreeFactorizationOne :: (FiniteField k, Eq k, Random k, RandomGen g) => Int -> UniPoly k -> g -> (Maybe (UniPoly k), g)
equalDegreeFactorizationOne d h g
  = let (u,g') = randomPolyOfDegreeLessThan (degree' h) g
        m = do
          guard (u /= 0)
          let v = toMonic $ gcdP h u
          if v /= 1
            then return v
            else do
            let w = powMod u ((q^d-1) `div` 2) h
                h' = toMonic $ gcdP h (w-1)
            guard (h' /= 1 && h' /= h)
            return h'
    in (m,g')
  where q = order (leadingCoefficient h)

等次数因数分解によって、与えられた多項式を完全に因数分解する関数を書く。 状態を持って回るのがさすがに面倒になってくるので、中で State モナドを使う。 試すときの利便性を考え、グローバルな StdGen を使う版 (equalDegreeFactorizationIO) も用意しておく。

-- Input: nonzero, monic, squarefree, equal-degree
equalDegreeFactorization :: (FiniteField k, Eq k, Random k, RandomGen g) => Int -> UniPoly k -> g -> ([UniPoly k], g)
equalDegreeFactorization d h = runState (loop h [])
  where loop h acc | degree' h == 0 = return acc
                   | degree' h == d = return (h:acc)
                   | otherwise = do
                       m <- state (equalDegreeFactorizationOne d h)
                       case m of
                         Nothing -> loop h acc
                         Just h' -> do -- h' is a nontrivial factor of h
                           acc' <- loop h' acc
                           loop (h `divP` h') acc'

equalDegreeFactorizationIO :: (FiniteField k, Eq k, Random k) => Int -> UniPoly k -> IO [UniPoly k]
equalDegreeFactorizationIO d h = getStdRandom (equalDegreeFactorization d h)

次数別因数分解と等次数因数分解を組み合わせ、与えられた無平方モニック多項式を完全に因数分解する関数を書く。 これも State モナドを使っている。 試すときの利便性を考え、グローバルな StdGen を使う版 (factorCZIO) も用意しておく。

-- Input: nonzero, monic, squarefree
factorCZ :: (FiniteField k, Eq k, Random k, RandomGen g) => UniPoly k -> g -> ([UniPoly k], g)
factorCZ f = runState
             $ fmap concat
             $ forM (distinctDegreeFactorization f)
             $ \(d,h) -> state (equalDegreeFactorization d h)

factorCZIO :: (FiniteField k, Eq k, Random k) => UniPoly k -> IO [UniPoly k]
factorCZIO f = getStdRandom (factorCZ f)

finite-field パッケージの PrimeField 型は Random クラスのインスタンスにはなっていないようである。 そこで、 PrimeField 型に対する Random クラスのインスタンスを用意してやる(本当は、いち利用者にすぎない我々のパッケージでこういうインスタンスを定義するのは行儀が良くない)。

instance KnownNat p => Random (PrimeField p) where
  randomR (lo,hi) g = case randomR (toInteger lo, toInteger hi) g of
                        (a,g') -> (fromInteger a, g')
  random = randomR (minBound, maxBound)

algebraic-real.cabal の exposed-modules に今回の Numeric.AlgebraicReal.Factor.CantorZassenhaus モジュールを追加し、

  exposed-modules:     Numeric.AlgebraicReal.UniPoly
                     , ...
                     , Numeric.AlgebraicReal.Factor.CantorZassenhaus

build-depends に random と mtl を追加しておく:

  build-depends:       base >= 4.7 && < 5
                     , vector
                     , finite-field
                     , random
                     , mtl

さあ試してみよう。 真っ先に :set prompt "> " して、プロンプトを読みやすくする。 次に :m + System.Random Data.FiniteField して、 System.RandomData.FiniteField モジュールが見えるようにする。 その次に、 :set -XDataKinds で DataKinds 拡張を有効にする。

$ stack repl
...
...> :set prompt "> "
> :m + System.Random Data.FiniteField
> :set -XDataKinds

グローバルな乱数生成器を初期化してやる。

> setStdGen (mkStdGen 325235)

さあ始めよう。

> let x = ind
> factorCZIO ((x^7-1) :: UniPoly (PrimeField 37))
[UniPoly [36,1],UniPoly [36,8,9,1],UniPoly [36,28,29,1]]

x^7-1\mathbb{F}_{37} で因数分解したところ、 x^7-1=(x+36)(x^3+9x^2+8x+36)(x^3+29x^2+28x+36) という結果が得られた。 先に与えた乱数のシードが違えば、同じ次数の因数の順番が違うかもしれない(つまり、 [UniPoly [36,1],UniPoly [36,28,29,1],UniPoly [36,8,9,1]] が出力されるかもしれない)。

考える体を変えれば、異なる因数分解が出てくる:

> factorCZIO ((x^7-1) :: UniPoly (PrimeField 11))
[UniPoly [10,1],UniPoly [10,4,5,1],UniPoly [10,6,7,1]]

一方、同じ多項式 x^7-1\mathbb{F}_7 で分解しようとすると、結果が返ってこない:

> factorCZIO ((x^7-1) :: UniPoly (PrimeField 7))

これは、「引数が無平方多項式である」という仮定を満たさないためである(この場合 x^7-1=(x-1)^7 となる)。