#7 整数係数
これまで代数的実数の定義多項式には有理数係数多項式（Haskellの型で言うと UniPoly Rational）を使っていたが、今後は整数係数多項式（Haskellの型で言うと UniPoly Integer）を使うことにする。
また、そのついでに、整域やユークリッド整域を表す型クラスを用意する。
可換環
この連載では、単に環と言ったら単位元を持つ可換環を指す。
整域 (integral domain) とは、 \(1\ne 0\) であって、非ゼロな元どうしの積が非ゼロとなる環のことである。 つまり、環 \(R\) が整域であるとは、任意の元 \(a,b\in R\) に対して \[ab=0 \Rightarrow a=0\text{ or }b=0\] を満たすことである。
整域でない環の例としては、合成数 \(n\) についての \(\mathbb{Z}/n\mathbb{Z}\) がある。 例えば、 \(\mathbb{Z}/6\mathbb{Z}\) において、 \(2\ne 0\), \(3\ne 0\) であるが \(2\cdot 3=0\) である。
別の整域でない環の例としては、 \(n\ge 2\) について \(n\) 次の正方行列のなす環 \(M(n,K)\) がある。
以後、もっぱら整域を扱う。
一意分解整域 (UFD; unique factorization domain) とは、素因数分解の類似である素元分解ができる整域である。
単項イデアル整域 (PID; principal ideal domain) とは、全てのイデアルが単項生成であるような整域である。
単項イデアル整域であれば一意分解整域である。
ユークリッド整域 (Euclidean domain) とは、ユークリッド除算ができる整域である。
ユークリッド整域であれば単項イデアル整域である。
体であればユークリッド整域である。
ややマイナーな概念ではあるが、GCD 整域 (GCD domain) を、任意の2元 \(a,b\) に対して最大公約元 \(\gcd(a,b)\) が存在する整域とする（もちろん、最大公約元は一意とは限らない）。
一意分解整域であれば GCD 整域である。
（TODO：図）
例 整数 \(\mathbb{Z}\) はユークリッド整域である。
例 多項式環 \(R[x]\) は、
- \(R\) が体の時はユークリッド整域、
- \(R\) がUFDの時はUFD,
- \(R\) がGCD整域の時はGCD整域、
- \(R\) が整域の時は整域
である。
例 ガウス整数のなす環 \(\mathbb{Z}[\sqrt{-1}]\) や、アイゼンシュタイン整数のなす環 \(\mathbb{Z}[\omega]\) （ただし \(\omega^2+\omega+1=0\)）はユークリッド整域である。
（TODO：GCD整域であるがUFDでない例）
演算
我々が抽象化して扱いたい演算について確認する。
除算
\(a\) が \(b\) (\(b\ne 0\)) の倍数であることが分かっている場合に \(a=bc\) となる \(c\) を計算する操作（除算）があると良い。
このような除算が行える環のクラスとしては、体やユークリッド整域があるが、実際にはもっと広いクラスでこのような操作を行える。 例えば、整数係数多項式環 \(\mathbb{Z}[x]\) で \(2x+2\) を \(x+1\) で割ると商が \(2\) となる。
一方で、一般の可換環ではこのような除算は行えない（well-defined ではない）。 例えば、 \(\mathbb{Z}/6\mathbb{Z}\) においては \[3=3\cdot 1=3\cdot 3=3\cdot 5\] なので、3 を 3 で「割った」商として 1, 3, 5 の3個が取れ、一意性が成り立たない。
環が整域であれば商の一意性が成り立つので、このような除算は整域に対して定義するのが妥当であろう。
ユークリッド除算
GCD の計算に使うユークリッドの互除法では、ユークリッド除算を使う。
しかし、多項式の GCD の場合は、簡約剰余列や部分終結式剰余列など、ユークリッド除算に頼らないアルゴリズムが利用できる（#4 を参照）。
GCD
GCD があると、多項式の内容と原始部分を計算できる。
TODO: アルゴリズム
Haskell での型クラス
Haskell 標準の Prelude で定義されている型クラスの階層は雑なので、 divMod や gcd などの関数を多項式について使えないという問題がある。 もちろん、自分で型クラスを定義すれば良いのだが、じゃあどういう階層を定義するんだという話になる。
参考として、 algebraパッケージ では、 Numeric.Domain 以下のモジュールとして Numeric.Domain.Integral, Numeric.Domain.Euclidean などが定義されている。
我々の用途では、抽象化したい演算として、整域の除算と GCD、それにユークリッド除算があげられる。 よって、それらの演算に対応する型クラスを定義する。
数学的な環のクラスに対応する Haskell の型クラスは、次のようにする：
- 環:
Num a（Prelude のものをそのまま使う） - 整域:
(Num a) => IntegralDomain a- 除算を定義する
- GCD 整域:
(IntegralDomain a) => GCDDomain a- GCD を定義する
- ユークリッド整域:
(GCDDomain a) => EuclideanDomain a- ユークリッド除算を定義する
- 体:
(Fractional a, GCDDomain a)（2つの型クラスを組み合わせる）- 本当は
GCDDomainあたりをFractionalのスーパークラスとしたいが、 Prelude の型クラスにそういう風に手を加えることはできない。
- 本当は
これらのクラス階層を定義するために、 src/Numeric/AlgebraicReal/Class.hs を作成しよう：
このモジュールで定義する型クラスについて、有理数 Ratio a と整数 Integer に対してのインスタンスを定義することにする。 多項式環 UniPoly a に関しては、 UniPoly.hs でインスタンスを定義する。
整域
除算を定義する。
-- | 整域
--
-- \(1 \ne 0\) となる可換環で、非自明な零因子を持たないものを整域という。
class (Num a) => IntegralDomain a where
-- | 除算
--
-- @a@ が @b@ の倍数であるとき、 @divide a b@ は @a = b * c@ となる @c@ を返す。
-- @a@ が @b@ の倍数でない時の挙動は規定しない。
divide :: a -> a -> a
infixl 7 `divide`
instance (Integral a) => IntegralDomain (Ratio a) where
divide = (/)
instance IntegralDomain Integer where
divide = div
{-
With UndecidableInstances:
instance (Fractional a) => IntegralDomain a where
divide = (/)
-}GCD 整域
gcd という名前はすでに Prelude で使われているので、衝突を避けるためにここでは gcdD という関数名を使う（Prelude との名前の衝突上等、という立場もあるだろうが、ここではそうしない）。
GCDDomain クラスのメソッドは本来は gcdD :: a -> a -> a だけで良いが、多項式の内容 (content) を計算するための関数も用意しておく。
内容を計算する関数の引数の型は UniPoly a となるべきだが、このモジュール (Class.hs) からは UniPoly 型が見えないので、 Vector a で代用する。
このクラスで内容を計算する意図としては、体の場合に contentV 関数の定義をいじることによって、整数係数多項式の原始多項式と体係数多項式のモニック多項式を統一的に扱えるようにするというものがある。 まあ、 Prelude の showList 並みにせこいトリックである。
なお、体の場合の GCD は、 0 同士の場合を除き、常に 1 とする。
class (IntegralDomain a) => GCDDomain a where
gcdD :: a -> a -> a
contentV :: V.Vector a -> a
contentV = V.foldr gcdD 0
-- | 体の場合の gcdD 関数のデフォルト実装
fieldGcd :: (Eq a, Fractional a) => a -> a -> a
fieldGcd 0 0 = 0
fieldGcd _ _ = 1
-- | 体の場合の contentV 関数のデフォルト実装
fieldContentV :: (Eq a, Fractional a) => V.Vector a -> a
fieldContentV xs | V.null xs = 0
| otherwise = V.last xs
instance (Integral a) => GCDDomain (Ratio a) where
gcdD = fieldGcd
contentV = fieldContentV
instance GCDDomain Integer where
gcdD = gcd
contentV xs | V.null xs = 0
| V.last xs < 0 = negate (gcdV 0 xs) -- 内容の符号に、最高次の係数の符号を反映させる
| otherwise = gcdV 0 xs -- 短絡評価を考えなければ foldr gcd 0 でも良い
where
-- foldl/foldr と gcd の組み合わせでは GCD が 1 になっても残りの部分が評価される。
-- 列の途中で GCD が 1 になれば全体の GCD は 1 で確定なので、そういう短絡評価する。
gcdV 1 _ = 1
gcdV a v | V.null v = a
| otherwise = gcdV (gcdD (V.last v) a) (V.init v)
{-
With UndecidableInstances:
instance (Eq a, Fractional a) => GCDDomain a where
gcdD = fieldGcd
contentV = fieldContentV
-}【12月14日 更新】 fieldGcd 関数と fieldContentV 関数を追加した。
ユークリッド整域
ユークリッド整域を表すクラスも一応用意しておく：
-- ユークリッド整域
class (Num a) => EuclideanDomain a where
-- a, b に対して、 (q,r) = divModD a b は a = q * b + r を満たす
divModD :: a -> a -> (a, a)
modD :: a -> a -> a
modD x y = snd (divModD x y)
infixl 7 `modD`
instance EuclideanDomain Integer where
divModD = divMod
{- あってもなくても良い：
instance (Integral a) => EuclideanDomain (Ratio a) where
divModD x y = (x / y, 0)
With UndecidableInstances:
instance (Fractional a) => EuclideanDomain a where
divModD x y = (x / y, 0)
-}ただ、今後はもっぱら整数係数多項式（これはユークリッド整域ではない）を扱いたいので、今更ユークリッド整域を追加しても実はそんなに嬉しくない。
一変数多項式環についての実装
ここからは UniPoly.hs を編集していく。 まずは先頭付近に import Numeric.AlgebraicReal.Class を追加する。
定数倍する操作は scaleP という関数ですでに提供しているが、今回 IntegralDomain という型クラスを導入したので、定数で割る操作も書けるようになった：
unscaleP :: (Eq a, IntegralDomain a) => a -> UniPoly a -> UniPoly a
unscaleP a f = mapCoeff (`divide` a) f以前は UniPoly Integer に対して content_int, contentAndPrimitivePart_int, primitivePart_int などの関数を定義したが、我々は GCD 整域を表す型クラスを定義したので、これらの関数をより一般化できる：
-- | 多項式の内容を計算する
content :: (GCDDomain a) => UniPoly a -> a
content (UniPoly xs) = contentV xs
-- | 多項式の内容と原始部分を計算する
contentAndPrimitivePart :: (Eq a, GCDDomain a) => UniPoly a -> (a, UniPoly a)
contentAndPrimitivePart f
| c == 1 = (c, f)
| otherwise = (c, unscaleP c f)
where c = content f
-- | 多項式の原始部分を計算する
primitivePart :: (Eq a, GCDDomain a) => UniPoly a -> UniPoly a
primitivePart = snd . contentAndPrimitivePartあとは、ファイルの最後の方に型クラスのインスタンスを書いていく。 まずは整域である。
instance (Eq a, IntegralDomain a) => IntegralDomain (UniPoly a) where
divide f g
| g == 0 = error "divide: divide by zero"
| degree f < degree g = zeroP -- f should be zero
| otherwise = loop zeroP f
where
l = degree' f - degree' g + 1
b = leadingCoefficient g
-- invariant: f == q * g + r
loop q r | degree r < degree g = q -- r should be zero
| otherwise = loop (q + q') (r - q' * g)
where q' = unscaleP b (UniPoly (V.drop (degree' g) (coeff r)))【注】この記事を公開した時点での上記 divide 関数にはバグがあった。 現在のバージョンでは修正済みなので、バグ持ちのコードをコピペしたという人は修正しておくこと。
多項式に関しては、係数が体の場合のユークリッド除算 divModP と、一般の場合の擬除算 pseudoDivModP をすでに実装しているが、係数が整域の場合の普通の除算を実装するのは今回が初めてである。
次は GCD 整域である。#4 でいくつか多項式剰余列を見たが、ここでは部分終結式剰余列を使うことにする。
gcd_subresultantPRS :: (Eq a, IntegralDomain a) => UniPoly a -> UniPoly a -> UniPoly a
gcd_subresultantPRS f 0 = f
gcd_subresultantPRS f g
| degree f < degree g = gcd_subresultantPRS g f
| otherwise = case pseudoModP f g of
0 -> g
rem -> let !d = degree' f - degree' g
!s = (-1)^(d + 1) * rem
in loop d (-1) g s
where
loop !_ _ f 0 = f
loop !d psi f g = case pseudoModP f g of
0 -> g
rem -> let !d' = degree' f - degree' g
!c = leadingCoefficient f
!psi' | d == 0 = psi
| d > 0 = ((-c)^d) `divide` (psi^(d-1))
!beta = -c * psi' ^ d'
!s = unscaleP beta rem
in loop d' psi' g s
instance (Eq a, GCDDomain a) => GCDDomain (UniPoly a) where
gcdD x y = let (xc,xp) = contentAndPrimitivePart x
(yc,yp) = contentAndPrimitivePart y
in scaleP (gcdD xc yc) $ primitivePart (gcd_subresultantPRS xp yp)
contentV xs | V.null xs = 0
| otherwise = scaleP (contentV (V.map fst ys)) $ gcdV 0 (V.map snd ys)
where
ys = V.map contentAndPrimitivePart xs
gcdV 1 _ = 1
gcdV a v | V.null v = a
| otherwise = gcdV (primitivePart $ gcd_subresultantPRS (V.last v) a) (V.init v)【12月14日 更新】 gcd_subresultantPRS 関数が、不適当に degree f >= degree g を仮定していたのを修正。
【2018年2月4日 更新】 gcd_subresultantPRS 関数が degree f == degree g の場合にエラーを吐くことがあったので、 psi' の定義を修正。
内容 content の実装でトリックを使っているので、有理数係数の多項式について gcdD を計算した場合はモニック多項式が得られる。
ユークリッド整域も実装しておく：
型の変更
これまでに定義した多項式についての関数のうちいくつかは、今回定義した型クラスを使って一般化できる。
UniPoly.hs の content_int, primitivePart_int, contentAndPrimitivePart_int の各関数はすでに一般化した。
Resultant.hs の resultant_int, resultant_poly, primitivePRS_int, reducedPRS_int, subresultantPRS_int も、一般化できる（読者への課題とする）。
AlgReal.hs の negativePRS も一般化できるが、スツルム列では符号が重要になるので、その辺を善処してやる必要がある。 その準備として、多項式剰余列の計算の際に係数も計算する変種を用意しておく。
-- Resultant.hs に追加
-- subresultantPRS' f g = (b,r) : subresultantPRS' g r
-- where lc(g)^l * f = q * g + b * r, l = deg f - deg g + 1
subresultantPRS' :: (Eq a, IntegralDomain a) => UniPoly a -> UniPoly a -> [(a,UniPoly a)]
subresultantPRS' _ 0 = []
subresultantPRS' f g = case pseudoModP f g of
0 -> []
rem -> let !d = degree' f - degree' g
!s = (-1)^(d + 1) * rem
in ((-1)^(d + 1), s) : loop d (-1) g s
where
loop !_ _ _ 0 = []
loop d psi f g = case pseudoModP f g of
0 -> []
rem -> let !d' = degree' f - degree' g
!c = leadingCoefficient f
!psi' = (-c)^d `divide` psi^(d-1)
!beta = -c * psi' ^ d'
!s = unscaleP beta rem
in (beta,s) : loop d' psi' g sこれを利用すると、係数が整域の場合の負係数多項式剰余列を次のように実装できる：
-- AlgReal.hs の negativePRS を置き換え
-- | Negative polynomial remainder sequence using subresultant PRS
--
-- Assumption: 'degree f > degree g'
negativePRS :: (Ord a, IntegralDomain a) => UniPoly a -> UniPoly a -> [UniPoly a]
negativePRS f g = f : g : loop 1 f 1 g (subresultantPRS' f g)
where
loop !_ _ !_ _ [] = []
loop !s f !t g ((b,x):xs)
-- b * (lc g)^(degree f - degree g + 1) * s > 0
| sign b * (sign $ leadingCoefficient g)^(degree' f - degree' g + 1) * s > 0 = -x : loop t g (-1) x xs
| otherwise = x : loop t g 1 x xsUniPoly.hs の squareFree :: (Eq a, Fractional a) => UniPoly a -> UniPoly a は、次のように一般化できる：
squareFree :: (Eq a, GCDDomain a) => UniPoly a -> UniPoly a
squareFree f = f `divide` gcdD f (diffP f)整数係数へ
型クラスについての準備が整ったので、代数的実数のコードも書き換えていこう。 まず、
data AlgReal = FromRat !Rational
| AlgReal !(UniPoly Rational) !Int !Rational !Rational
deriving (Show)を
data AlgReal = FromRat !Rational
| AlgReal !(UniPoly Integer) !Int !Rational !Rational
deriving (Show)に、
definingPolynomial :: AlgReal -> UniPoly Rational
definingPolynomial (FromRat x) = ind - constP x
definingPolynomial (AlgReal p _ _ _) = pは
definingPolynomial :: AlgReal -> UniPoly Integer
definingPolynomial (FromRat x) = fromCoeff (V.fromList [- numerator x, denominator x])
definingPolynomial (AlgReal p _ _ _) = pに書き換える。
当然、このままでは型の不整合が大量に発生するので、頑張って不整合を潰していく。
整数係数多項式と有理数
これまでの多項式に有理数を代入するという操作を何箇所かで行なっている。 有理数係数多項式に有理数を代入する操作については valueAt :: (Num a) => a -> UniPoly a -> a という関数を利用できたが、整数係数多項式と有理数では型が違うのでそのまま利用できない。
そこで、場当たり的ではあるが「整数係数多項式に有理数を代入する」次のような関数を追加する：
-- UniPoly.hs に追加
valueAtZQ :: Rational -> UniPoly Integer -> Rational
valueAtZQ t (UniPoly xs) = V.foldr' (\a b -> fromInteger a + t * b) 0 xs多項式 \[a_n X^n+\cdots+a_1 X+a_0\] に \(u/v\) を代入した時の分子と分母を返す関数があると便利なので、ついでに追加しておく：
-- UniPoly.hs に追加
-- homogeneousValueAt x y (a_n X^n + ... + a_1 X + a_0)
-- = (a_n x^n + a_{n-1} x^{n-1} y + ... + a_1 x y^{n-1} + a_0 y^n, y^n)
homogeneousValueAt :: (Eq a, Num a) => a -> a -> UniPoly a -> (a, a)
homogeneousValueAt num den f@(UniPoly coeff)
| f == 0 = (0, 1)
| otherwise = (V.foldr' (\x y -> x + num * y) 0 $ V.zipWith (*) coeff denseq, V.head denseq)
where
-- numseq = V.iterateN (V.length coeff) (* num) 1
denseq = V.reverse (V.iterateN (V.length coeff) (* den) 1)【2018年2月4日 追加】 また、多項式 \[f(X)=a_n X^n+\cdots+a_1 X+a_0\] に別の多項式 \(g(x)/v\) を代入した時の分子と分母を返す関数があると便利なので、追加しておく：
-- homogeneousCompP (a_n X^n + ... + a_1 X + a_0) g d
-- = (a_n g^n + a_{n-1} g^{n-1} d + ... + a_1 g d^{n-1} + a_0 d^n, d^n)
homogeneousCompP :: (Eq a, Num a) => UniPoly a -> UniPoly a -> a -> (UniPoly a, a)
homogeneousCompP f@(UniPoly coeff) g den
| f == 0 = (0, 1)
| otherwise = (V.foldr' (\x y -> constP x + g * y) 0 $ V.zipWith (*) coeff denseq, V.head denseq)
where
denseq = V.reverse (V.iterateN (V.length coeff) (* den) 1)「整数係数多項式に有理数を代入した時の符号」も valueAtZQ を使って実装できる：
-- AlgReal.hs に追加
signAtZQ :: Rational -> UniPoly Integer -> Int
signAtZQ x p = sign (valueAtZQ x p)
isZeroAtZQ :: Rational -> UniPoly Integer -> Bool
isZeroAtZQ x p = valueAtZQ x p == 0サフィックス ZQ は、「整数」 \(\mathbb{Z}\) 「有理数」 \(\mathbb{Q}\) のつもりであり、深い意味はない。
ちなみに、 signAtZQ は整数演算のみで実装できる（演習問題）。
signAtZQ と同様に、
signAtZQX :: ExtReal Rational -> UniPoly Integer -> IntvarianceAtZQ :: Rational -> [UniPoly Integer] -> IntvarianceAtZQX :: ExtReal Rational -> [UniPoly Integer] -> IntcountRealRootsBetweenZQ :: Rational -> Rational -> UniPoly Integer -> IntcountRealRootsBetweenZQX :: ExtReal Rational -> ExtReal Rational -> UniPoly Integer -> Int
も実装する（コード略）。
以下の関数については、 UniPoly Rational を UniPoly Integer に変え、 ZQ 版の関数を使うようにする：
intervalsWithSign :: UniPoly Integer -> Int -> Rational -> Rational -> [Interval]
mkAlgReal :: UniPoly Integer -> Rational -> Rational -> AlgReal
mkAlgRealWithInterval :: UniPoly Integer -> Interval -> AlgReal
realRoots :: UniPoly Integer -> [AlgReal]
realRootsBetween :: UniPoly Integer -> ExtReal Rational -> ExtReal Rational -> [AlgReal]
mkAlgRealWithCReal :: UniPoly Integer -> CReal -> AlgReal多項式の実根を求める操作に関しては、有理数係数に関するものも引き続き用意しておく：
realRootsBetweenQ :: UniPoly Rational -> ExtReal Rational -> ExtReal Rational -> [AlgReal]
realRootsBetweenQ f lb ub
| f == 0 = error "realRoots: zero"
| degree' f == 0 = []
| otherwise = let commonDenominator = V.foldl' (\a b -> lcm a (denominator b)) 1 (coeff f)
fz = primitivePart $ mapCoeff (\x -> numerator x * (commonDenominator `div` denominator x)) f
in realRootsBetween fz lb ubその他のコード変更
根の限界については、 UniPoly Integer -> Rational の関数に変える：
rootBound :: UniPoly Integer -> Rational
rootBound f | f == 0 = error "rootBound: polynomial is zero"
| otherwise = 1 + (V.maximum $ V.map (abs . (% lc)) $ V.init $ coeff f)
where lc = leadingCoefficient fNum クラスのインスタンスは、次のように書き換える：
instance Num AlgReal where
...
FromRat k + AlgReal f s a b = mkAlgReal (fst $ homogeneousCompP f (definingPolynomial (FromRat k)) (denominator k)) (a + k) (b + k)
...
FromRat k - AlgReal f s a b = mkAlgReal (fst $ homogeneousCompP f (- definingPolynomial (FromRat k)) (denominator k)) (k - b) (k - a)
AlgReal f s a b - FromRat k = mkAlgReal (fst $ homogeneousCompP f (definingPolynomial (FromRat (-k))) (denominator k)) (a - k) (b - k)
...
FromRat k * AlgReal f s a b
| k == 0 = 0
| k > 0 = AlgReal f_x_k s (a * k) (b * k)
| k < 0 = AlgReal f_x_k (-s) (b * k) (a * k)
where f_x_k = fst $ homogeneousValueAt (scaleP (denominator k) ind) (fromInteger $ numerator k) (mapCoeff fromInteger f) -- f(x/k)
...【12月17日 修正】FromRat と AlgReal の和がゼロになる場合を考慮し忘れていた件と、 FromRat _ - AlgReal _ _ _ _ の実装が間違っていたのを修正。
【2018年2月4日 修正】FromRat と AlgReal の足し算と引き算の実装が間違っていたので、 homogeneousCompP を使うように修正。
sqrtQ / nthRootQ では realRootsBetween の代わりに realRootsBetweenQ を使うようにする。
powIntA と valueAtA では、定義多項式の型が変わったことに追従するため、 definingPolynomial x を mapCoeff fromInteger $ definingPolynomial x :: UniPoly Rational のような感じで有理数係数に直す。
【12月14日 追記】 IntegralDomain, GCDDomain の各クラスに対する AlgReal のインスタンスを追加する：