#8 多項式の因数分解 その1 Kronecker の方法と無平方分解

投稿日：2017年12月16日，最終更新：2025年4月5日

今回から数回に分けて、多項式の因数分解を取り扱う。 多項式の因数分解は、代数的数の最小多項式を求めるために必要となる。 #6 までで見たように、因数分解しなくても代数的実数についての演算は一通り実装できるが、その場合、演算するにつれて多項式の次数が不必要に増大する。 例えば、 (1+\sqrt{2})-\sqrt{2} は有理数 1 だが、前回までの実装ではこれは4次の多項式の根として表されることになる。

前回は代数的数を表すのに使う多項式を有理数係数から整数係数にしたので、因数分解したい多項式は整数係数である。

今回は、 Kronecker の方法と、無平方分解を紹介する。 Kronecker の方法は結局採用しないので、無駄に文章を読みたくない方は「無平方分解」のところまで読み飛ばして欲しい。

Kronecker の方法

Kronecker の方法は、整数係数多項式を因数分解するための最も古くから知られたアルゴリズムである。 ただし、効率が悪すぎるので現代的な計算機代数の教科書には載っていない。 英語版の Wikipedia でも参照してもらいたい。 …で済ませたいところだが、一応ここでも説明しておく。

因数分解したい整数係数多項式を f とする。 整数係数多項式 g が f の自明でない因数であったとする。 すると、整数 k を固定した時、整数 g(k) は f(k) の約数である。 f(k) が 0 でなければ、 g(k) としてありうる値の種類は有限個である。 この事実から、 g としてありうる多項式を有限個に絞り込むことができる。 （本当は f の因数分解が整数係数の範囲と有理数係数の範囲で一致することを確認しなければならないが、ここでは省略する）

アルゴリズム的に記そう： f(k) が 0 にならないような k を m 個選び、それを k_1,\ldots,k_m とする。 各 k_i について f(k_i) の約数を全て列挙し、それぞれ1つずつ選ぶ。 つまり、 d_i\mid f(k_i) となるような整数の組 (d_1,\ldots,d_m) を1つ選ぶ。 この組 (d_1,\ldots,d_m) に対して、 g(k_i)=d_i となるような m-1 次以下の多項式 g がただ1つ存在する。 この g が f を割り切れば f の因数が見つかったことになるし、どのような (d_1,\ldots,d_m) を選んでも g が f を割り切らなければ、 f は m-1 次以下の因数を持たないと結論づけられる。

具体例で実践してみよう。

例 f(x)=x^3+1

f が既約でなければ、次数が \frac{\deg f}{2} 以下の因数が存在するはずである。 今回は \frac{\deg f}{2}=\frac{3}{2} なので、つまり、次数が 1 の因数があるかどうかを判断し、それを求めたい。

値が 0 でないような x をいくつか挙げてみると、 x=0,1 がある。 f(0)=1, \quad f(1)=2 であり、 f(0) の約数は \{1,{-1}\} の2個、 f(1) の約数は \{1,{-1},2,{-2}\} の4個である。 そこで、 g(0),g(1) の値としてありうる組み合わせ8個と、それを実際に与える g(x) を列挙する：

g(0) g(1) g(x) f(x) の因数？
1 1 1 自明
1 -1 -2x+1 No
1 2 x+1 Yes
1 -2 -3x+1 No
-1 1 2x-1 No
-1 -1 -1 自明
-1 2 3x-1 No
-1 -2 -x-1 Yes

非自明な因数が欲しいので、次数が 0 の多項式（定数）は適格ではない。 定数でない g で、 f を実際に割り切るのは x+1 と -x-1 である。 というわけで、 f の自明でない因数として \pm(x+1) が見つかった。

f を x+1 で割ると f_2(x)=x^2-x+1 が得られるので、 f_2 について同様に次数 1 の因数を探す。 今度は、ありうる場合を全て試しても f_2 の因数は見つからないので、 f_2 は既約だとわかる。

多項式補間

以下、標数 0 の体およびその部分環を係数とする多項式について考える。

多項式が通る点 (x_i,y_i) （ただし i\ne j なら x_i\ne x_j）を n+1 個指定すると、次数 n 以下の多項式がただ一つ決まる。 というのは、多項式 f を f(x)=a_nx^n+\cdots+a_1x+a_0 と書くと、 \{a_i\} に関しての連立方程式 \begin{aligned} a_nx_0^n+\cdots+a_1x_0+a_0&=y_0, \\ a_nx_1^n+\cdots+a_1x_1+a_0&=y_1, \\ \vdots& \\ a_nx_n^n+\cdots+a_1x_n+a_0&=y_n \end{aligned} の解が一意的に存在するからである。 （ちなみに、この連立方程式の係数行列はファンデルモンドの行列 (Vandermonde matrix) と呼ばれている。）

ただ、上記の連立方程式の解を直接求めるのは大変である。 そのため、ラグランジュ補間をはじめとする、いくつかの多項式補間の方法が考案されている。 このような、指定された点を通る多項式のことを補間多項式と呼ぶ。

実装例

Kronecker の方法を実装してみよう。 まず、整数に関するラグランジュ補間を適当にでっち上げる（ラグランジュ補間についての説明はここではしない）。

{-# LANGUAGE BangPatterns #-}
module Numeric.AlgebraicReal.Factor.Kronecker where
import Numeric.AlgebraicReal.Class
import Numeric.AlgebraicReal.UniPoly
import qualified Data.Vector as V
import Data.List (inits)
import Data.Maybe (listToMaybe,maybeToList)
import Data.Ratio
import Control.Monad (sequence,guard)

-- | ラグランジュ補間
lagrangePoly :: (Eq a, Fractional a) => [(a,a)] -> UniPoly a
lagrangePoly [] = 0
lagrangePoly s = loop [] s 0
  where loop ss [] m = m
        loop ss (t@(x,y):ts) m = loop (t:ss) ts $ m + scaleP y (product [scaleP (recip (x - x')) (ind - constP x') | x' <- map fst ss ++ map fst ts])

-- | 有理数の分母が 1 の場合に整数 ('Integer') として返す
ratToInt :: Rational -> Maybe Integer
ratToInt x | denominator x == 1 = Just (numerator x)
           | otherwise = Nothing

-- | 有理数係数多項式の係数の分母がすべて 1 の場合に、整数係数多項式として返す
ratPolyToIntPoly :: UniPoly Rational -> Maybe (UniPoly Integer)
ratPolyToIntPoly f = fromCoeff <$> V.mapM ratToInt (coeff f)

-- | 整数係数の範囲でラグランジュ補間多項式を計算する
lagrangePolyI :: [(Integer,Integer)] -> Maybe (UniPoly Integer)
lagrangePolyI [] = Nothing
lagrangePolyI ts = ratPolyToIntPoly $ lagrangePoly $ map (\(x,y) -> (fromInteger x, fromInteger y)) ts

整数の約数を列挙する関数も用意する：

-- | 整数 n の全ての約数をリストとして返す
intDivisors :: Integer -> [Integer]
intDivisors 0 = error "intDivisors: zero"
intDivisors n | n < 0 = intDivisors (-n)
intDivisors n = loop 1 []
  where loop !i xs | i*i > n = xs
                   | i*i == n = i : (-i) : xs
                   | (q,0) <- divMod n i = i : (-i) : loop (i+1) (q : (-q) : xs)
                   | otherwise = loop (i+1) xs

いよいよ Kronecker の方法の本体である。 まず、「d 次以上の自明でない因数を一つ見つける」関数 oneFactor を書き、それを利用して完全な因数分解 factor を実装する。

「約数の組み合わせを全て試す」処理を書くのには、リストモナドが利用できる（Haskell を使っていて良かったと思える瞬間である）。 コードの詳しい説明は、コメントを参照して欲しい。

-- | Kronecker の方法で、整数係数多項式の自明でない因数を探す
oneFactor :: Int                      -- ^ 探すべき因数の次数の下限
          -> UniPoly Integer          -- ^ 因数を探す多項式
          -> Maybe (UniPoly Integer)  -- ^ 非自明な因数。見つからなかったら 'Nothing'
oneFactor !d 0 = Nothing
oneFactor !d f
  | n == 0 = Nothing
  | n == 1 = Nothing -- already irreducible
  | otherwise = listToMaybe $ do
      -- リストモナドを使い、非決定的計算っぽく書く

      -- nz は f(k) /= 0 となるような (k,f(k)) のリストである
      let nz = [(k,v) | k <- [0..], let v = valueAt k f, v /= 0]

      -- nz' は [[(k0,d) | d は f(k0) の約数],[(k1,d) | d は f(k1) の約数],..] の形のリストである
      let nz' = [map ((,) x) (intDivisors y) | (x,y) <- nz]

      -- inits によって、 nz' の先頭のいくつかを取り出したリストが得られる
      -- 非自明な因数が存在するなら n を f の次数として n/2 次以下のはずである
      -- d-1 次以下の因数は全て掃き出したという前提なので、補間によって d 次以上の多項式を作りたい
      -- 補間によって m 次 (d <= m <= n/2) の多項式を作りたい
      (m,ev) <- zip [d..] $ drop (d + 1) $ take (n `div` 2 + 2) $ inits nz'
      -- ここで length ev == m+1, d <= m <= n/2 が成り立つ

      -- 各 f(ki) について、その約数の組み合わせを全て試す
      ev' <- sequence ev
      -- ここでも length ev' == m+1 が成り立つ
      -- ev' は [(k0,d0),(k1,d1),...,(km,dm) | di は f(ki) の約数] の形のリストである

      -- ラグランジュ補間して多項式を得る
      g <- maybeToList (lagrangePolyI ev')

      -- 欲しいのは m 次の多項式である
      guard (degree g == Just m)

      -- g が f を割り切ることをチェックする
      guard (f `pseudoModP` g == 0)

      -- 最高次の係数を正にして返す
      if leadingCoefficient g < 0
        then return (-g)
        else return g

        where n = degree' f

-- | Kronecker の方法で、整数係数多項式を因数分解する
factor :: UniPoly Integer -> [UniPoly Integer]
factor 0 = error "factor: zero"
factor f = loop 1 f
  where loop !d f
          | degree f == Just 0 = [] -- constant
          | otherwise = case oneFactor d f of
                          Just g -> g : loop (degree' g) (f `divide` g)
                          Nothing -> [f]

以上のコードを src/Numeric/AlgebraicReal/Factor/Kronecker.hs に保存し、 .cabal ファイルに Numeric.AlgebraicReal.Factor.Kronecker を追加しよう。 それができたら早速、 stack repl で因数分解を試してみよう。 例によって :set prompt "> " を実行して長ったらしいプロンプトを省略する。

まずは簡単なところからいってみよう。

> let x = ind
> factor $ x^2 - 1
[UniPoly [-1,1],UniPoly [1,1]]
> factor $ x^3 - 2*x + 1
[UniPoly [-1,1],UniPoly [-1,1,1]]

x^2-1 の因数分解として (x-1)(x+1) が、 x^3-2x+1 の因数分解として (x-1)(x^2+x-1) が得られた。

3以上の n について、 x^n-1 の因数分解をやってみよう。

> factor $ x^3 - 1
[UniPoly [-1,1],UniPoly [1,1,1]]
> factor $ x^4 - 1
[UniPoly [-1,1],UniPoly [1,1],UniPoly [1,0,1]]
> factor $ x^5 - 1
[UniPoly [-1,1],UniPoly [1,1,1,1,1]]
> factor $ x^6 - 1
[UniPoly [-1,1],UniPoly [1,1],UniPoly [1,-1,1],UniPoly [1,1,1]]
> factor $ x^7 - 1
[UniPoly [-1,1],UniPoly [1,1,1,1,1,1,1]]
> factor $ x^8 - 1
[UniPoly [-1,1],UniPoly [1,1],UniPoly [1,0,1],UniPoly [1,0,0,0,1]]
> factor $ x^9 - 1
[UniPoly [-1,1],UniPoly [1,1,1],UniPoly [1,0,0,1,0,0,1]]
> factor $ x^10 - 1
[UniPoly [-1,1],UniPoly [1,1],UniPoly [1,-1,1,-1,1],UniPoly [1,1,1,1,1]]
> factor $ x^11 - 1
[UniPoly [-1,1],UniPoly [1,1,1,1,1,1,1,1,1,1,1]]

筆者の環境では、 x^{11}-1 の因数分解に40秒以上かかった。 「俺のハイスペックマシンでは一瞬だったぜ」という方は、 n を 13 や 17 など、より大きな素数にして試して欲しい。 （ちなみに、 ghci で式の評価にかかる時間を測定したいなら :set +s とすれば良いようだ。）

無平方分解

これまでに使ったアルゴリズムの中には、多項式が無平方であることを要求するものがあった。 今後登場するアルゴリズムの中にも、多項式が無平方であることを要求するものがある。 これまでは無平方成分（#1 で紹介）を使って多項式を無平方にしていたが、実は、多項式の微分と GCD の組み合わせを上手く利用すれば、因数の重複度ごとに因数分解することができる。

例を使って説明しよう。 f(x)=x^6+7x^5+20x^4+31x^3+29x^2+16x+4 という多項式があった時に、以前紹介した無平方成分の計算（squareFree 関数）を使うと x^4+4x^3+6x^2+5x+2 と無平方成分が得られる。 これは f の因数なので、 f(x)=(x^4+4x^3+6x^2+5x+2)g(x) という因数分解が得られることになる。

しかし、実はアルゴリズムを少し改変すれば (x^2+x+1)(x^2+3x+2)^2 という風に無平方多項式の積として因数分解できるのである。

【2018年9月2日 修正】無平方成分の計算が間違っていたので修正した。

なお、 Kronecker の方法と違い、無平方分解は、完全な（既約多項式への）因数分解を行う訳ではない。 実際、今の例では、 x^2+3x+2 はさらに (x+1)(x+2) と因数分解できる。

一般には、多項式 f が互いに素な無平方多項式 f_i の積として f=f_1f_2^2\cdots f_m^m と書ける時、無平方成分 f/\gcd(f,f') は指数を全て1に変えたもの f_1f_2\cdots f_m であった。 一方、無平方分解のアルゴリズムでは、各 f_i を個別に求めることができる。

素朴な無平方分解のアルゴリズム

このアルゴリズムの入力は（標数 0 の GCD 整域を係数とする）多項式であり、出力はその無平方分解、つまり「f=f_1 f_2^2\cdots f_m^m となる無平方多項式 f_i のリスト」である。

このアルゴリズムでは、 f の次数に関する帰納法（再帰）を使う。

まず、定数の場合は既に因数分解が完了している。

そうでない場合、 \gcd(f,f') の次数は \deg f よりも真に小さいので、 \gcd(f,f') に対して無平方分解を行うことができる。 無平方分解の結果は \gcd(f,f')=g_1 g_2^2\cdots g_{m-1}^{m-1} の形をしている。

分解したい多項式 f を f=f_1 f_2^2 \cdots f_m^{m} と書くと \gcd(f,f')=f_2 f_3^2\cdots f_m^{m-1} なので、 1 以上の i について f_{i+1}=g_i である。 あとは f_1 を求めたいが、これは \gcd(f,f')=f_2 f_3^2\cdots f_m^{m-1} の指数に1ずつ加えたもの f_2^2 \cdots f_m^{m} で f を割れば良い。

実装例は次のようになる (src/Numeric/AlgebraicReal/Factor/SquareFree.hs):

module Numeric.AlgebraicReal.Factor.SquareFree where
import Numeric.AlgebraicReal.Class
import Numeric.AlgebraicReal.UniPoly
import qualified Data.Vector as V

squareFreeFactorization :: (Eq a, GCDDomain a) => UniPoly a -> [(UniPoly a,Int)]
squareFreeFactorization 0 = [(0,1)]
squareFreeFactorization f = mf (primitivePart f)
  where
    mf f | degree f == Just 0 = [] -- constant
         | degree f == degree p = u -- t = 1 の場合
         | otherwise = (t,1) : u
      where r = mf (sqPart f)
            u = map (\(g,i) -> (g,i+1)) r
            p = sqPart f * product (map fst r) -- = product (map (uncurry (^)) u)
            t = f `divide` p
    sqPart f = primitivePart $ gcdD f (diffP f)

補助関数や変数の名前には特に意味はない。

発展：Yun のアルゴリズム

無平方分解のアルゴリズムの改良版として、 Yun のアルゴリズムを紹介しておく。

入力は f とする。

  1. v_1:=f/\gcd(f,f'), w_1:=f'/\gcd(f,f') を計算する。（v_1 は f の無平方成分となる）
  2. i=1,2,3,\ldots について、 v_i が 1 になるまで、以下を繰り返す：
    1. h_i:=\gcd(v_i,w_i-v_i') を計算する。
    2. v_{i+1}:=\frac{v_i}{h_i}, w_{i+1}:=\frac{w_i-v_i'}{h_i} を計算する。
  3. f=h_1 h_2^2 \cdots h_k^k である。

証明だとか、実行時間がどうなるかについては、 Modern Computer Algebra 等の教科書を参照してほしい（マルナゲ！）。 もちろん、自分で証明を試みても良いだろう。

Haskell での実装例は次のようになる：

yun :: (Eq a, GCDDomain a) => UniPoly a -> [(UniPoly a,Int)]
yun 0 = [(0,1)]
yun f = let f' = diffP f
            u = gcdD f f'
        in loop 1 (f `divide` u) (f' `divide` u)
  where loop !i v w | degree' v == 0 = []
                    | h == 1 = loop (i+1) v s
                    | otherwise = (h,i) : loop (i+1) (v `divide` h) (s `divide` h)
          where s = w - diffP v
                h = gcdD v s

REPL で試してみる。

> :set +s
> let x = ind
(0.00 secs, 344,424 bytes)
> let f = (x^10-1)*(x^5-1)*(x^14-1)^4*(x^29-1)^7 :: UniPoly Integer
(0.02 secs, 344,640 bytes)
> squareFreeFactorization f
[(UniPoly [1,-1,1,-1,1],1),(UniPoly [1,1,1,1,1],2),(UniPoly [1,0,1,0,1,0,1,0,1,0,1,0,1],4),(UniPoly [1,1],5),(UniPoly [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1],7),(UniPoly [-1,1],13)]
(4.56 secs, 2,090,679,664 bytes)
> yun f
[(UniPoly [1,-1,1,-1,1],1),(UniPoly [1,1,1,1,1],2),(UniPoly [1,0,1,0,1,0,1,0,1,0,1,0,1],4),(UniPoly [1,1],5),(UniPoly [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1],7),(UniPoly [-1,1],13)]
(2.81 secs, 1,251,632,416 bytes)

確かに、素朴なアルゴリズムよりも Yun のアルゴリズムの方が早いようである。 （Haskell は遅延評価なので最初の squareFreeFactorization f には f 自体の評価にかかる時間も含まれているが、それは無視できる。）

正標数の場合

上記のアルゴリズムでは f が定数でなければ \gcd(f,f') の次数が f の次数よりも小さくなるということを仮定しているが、正標数ではそれは成り立たない。 どういうことかというと、標数 p>0 の場合 f が定数でないのに f'=0 となる場合がある（従って \gcd(f,f')=f となる）のだ。

例えば、標数 p の体で f(x)=x^{2p}+2x^p+3 を考えると、 f は定数でないにも関わらず f'=0 となる。 つまり、係数が 0 でない項の次数がすべて p の倍数だった場合に微分が消える。

ただ、係数体が有限体の場合は多項式 f\in\mathbb{F}_q[x] が f'=0 を満たせばある多項式 g\in\mathbb{F}_q[x] について f=g^p と書ける。 例えば、係数体が有限素体 \mathbb{F}_p であれば f(x)=(x^2+2x+3)^p となるし、有限体 \mathbb{F}_q (q=p^e) であれば f(x)=(x^2+2^{q/p}x+3^{q/p})^p となる。

そのような事情で、正標数の場合は、無平方成分や無平方分解のアルゴリズムの見直しが必要になる。

無平方分解したい多項式を f=f_1\cdots f_r^r, 各 f_i は無平方で、互いに素、 f の次数は n とする。 この時、 f'=\sum_{i=1}^r i f_i' f/f_i=\sum_{1\le i\le r,p\nmid i} i f_i' f/f_i である。 f の因数 f_i^i に注目すると \gcd(f_i^i,f')=\begin{cases} f_i^{i-1} & (p\nmid i) \\ f_i^i & (p \mid i) \end{cases} となるので \gcd(f,f')=\left(\prod_{1\le i\le r,p\nmid i}f_i^{i-1}\right) \left(\prod_{1\le i\le r,p\mid i}f_i^{i}\right) となる。 以後これを u:=\gcd(f,f') とおく。 （ただし、縦棒 \mid は「割り切る」という記号であり、 p\mid i は「p は i を割り切る」つまり「i は p の倍数である」という意味である）

標数 0 の場合の無平方成分は f/\gcd(f,f') で計算できたが、正標数の場合は f/\gcd(f,f')=\prod_{1\le i\le r,p\nmid i}f_i となる。 つまり、因数のうち、指数が p の倍数でないものだけが取り出されている。 以後これを v:=f/\gcd(f,f') とおく。

指数が p の倍数のものを取り出すにはどうするか。 v が p\nmid i についての f_i から成っていることに注目すると \gcd(u,v^n)=\prod_{1\le i\le r,p\nmid i}f_i^{i-1} を得る。 よって u をこれで割れば p\mid i についての f_i たちを取り出すことができる： u/\gcd(u,v^n)=\prod_{1\le i\le r,p\mid i}f_i^{i}

この u/\gcd(u,v^n) の微分は 0 となる。 よって、係数体が有限体であれば w^p=u/\gcd(u,v^n) となる多項式 w を取れる。

あとは、無平方成分を計算したいなら w について再帰的にアルゴリズムを適用して v と合わせれば良い。 無平方分解を（素朴に）計算したいなら、 \gcd(u,v^n) と w に再帰的にアルゴリズムを適用すれば良い。

（TODO: 「完全体」等の用語を使って説明する）

次回予告：モジュラーな方法

今後紹介する現代的な因数分解の手法は、モジュラーな方法、つまり係数を \mod n で考える手法である。 しかし我々が因数分解したい多項式は、あくまで整数係数の多項式である。 n を法とする因数分解が分かっても、それがそのまま整数係数の因数分解とはならない。

だが、もしも、因数分解したい多項式の係数の上限 N があらかじめ分かっていれば、 2N < n となる法 n で因数分解を行い、係数を \pm N に収まるように整数に戻してやる、ということができる。

モジュラーな方法の基礎になるのは、有限体 \mathbb{F}_q 係数での因数分解である。 というわけで、これから有限体をガンガン使って行くことになるが、筆者の感覚では、数学科以外の読者は有限体にそれほど馴染みがない。 そのため、このままでは今後の議論についてこられるのが数学科だけになってしまう可能性が高い。

そこで、次回は有限体の入門的な内容を取り扱い、数学科以外の読者でも有限体の基本的な性質を理解してもらえるようにする。 （要するに「予備知識なしでも読める」という謳い文句を使えるようにするためのアリバイ作りである）

それが済んだら、何回かに分けて有限体上の多項式の因数分解のアルゴリズムを解説する。

有限体での因数分解が実装できたら、それを利用して整数での因数分解を実装する。