#6 代数的数の演算 まとめ
#3 代数的数の演算と終結式 以降、代数的数の演算についての話題を引っ張ってきたが、ここらでまとめたいと思う。
#3 の最後で取り上げた課題としては、終結式の計算と、演算結果の分離区間についてのものがあったが、前者は #4 擬除算と多項式剰余列 で一応実装できた。 後者の方針の一つとして、前回（#5 区間演算と計算可能実数）は区間列の扱いを掘り下げた。 今回は、もう一つの方針として、実根の距離の最小値の評価を紹介する。
判別式と根の分離
実根の距離の評価に入る前に、判別式について論じておく。
判別式
二次方程式の判別式はみなさんご存知だろう。 二次方程式 \[ax^2+bx+c=0\] の判別式は \(D=b^2-4ac\) で、 \(D=0\) なら重根を持つ、というやつだ（実係数の場合は、 \(D\) の値によって実根の個数がどうのこうのということも言える）。 この方程式の（重複を込めて）2つの根を \(\alpha_1,\alpha_2\) と書けば、 \(D=a(\alpha_1-\alpha_2)^2\) となる。
三次方程式 \(x^3+px+q=0\) にも判別式と呼ばれる量があり、 \(D=0\) であればこの三次方程式は重根を持つ： \[D=-4p^3-27q^2\]
一般に、 \(n\) 次方程式 \(f(x)=0\) の根を重複を込めて \(\alpha_1,\ldots,\alpha_n\) と書いたとき、判別式 (discriminant) を \[D(f)=\operatorname{lc}(f)^{2n-2}\prod_{i=1}^n \prod_{j=i+1}^n (\alpha_i-\alpha_j)^2\] と定義する。 基本的な性質として、多項式 \(f\) が重根を持つ時、かつその時に限り \(D(f)=0\) となる。
演習問題：\(n=3\) の場合の判別式を、係数を使って表せ（ヒント：根と係数の関係を使う）。 特に、 \(f(x)=x^3+px+q\) の判別式が \(-4p^3-27q^2\) となることを確かめよ。
さて、多項式 \(f\) が重根を持つということは、 \(f\) とその微分 \(f'\) の GCD が非自明となることと同値であった。 さらに、終結式の性質より、それは \(\operatorname{res}(f,f')=0\) と同値だ。
この事実、および、判別式が「根に関する積」で書けることから、判別式と終結式に密接な関係があると察した読者も多いだろう。 実際、終結式を使って判別式を表すことができる。定理． \(n\) 次多項式 \(f(x)=ax^n+\cdots+a_0\) について、 \(f\) の判別式 \(D(f)\) は次のように書ける： \[D(f)=\frac{({-1})^{n(n-1)/2}}{a}\operatorname{res}(f,f').\]
証明． \[f(x)=a(x-\alpha_1)\cdots(x-\alpha_n)\] とおく。 この時、それぞれの \(\alpha_i\) について \[f'(\alpha_i)=a(\alpha_i-\alpha_1)\dots(\alpha_i-\alpha_{i-1})(\alpha_i-\alpha_{i+1})\dots(\alpha_i-\alpha_n)\] が成り立つ。 なぜなら、 \(f(x)\) は \(a(x-\alpha_1)\dots(x-\alpha_{i-1})(x-\alpha_{i+1})\dots(x-\alpha_n)\) と \(x-\alpha_i\) の積として書けるので、微分のライプニッツ則より \[\begin{aligned} f'(x)=&\bigl(a(x-\alpha_1)\cdots(x-\alpha_{i-1})(x-\alpha_{i+1})\cdots(x-\alpha_n)\bigr)'(x-\alpha_i) \\ &\quad+a(x-\alpha_1)\cdots(x-\alpha_{i-1})(x-\alpha_{i+1})\cdots(x-\alpha_n)\bigl(x-\alpha_i\bigr)' \\ =&(\cdots\cdots)(x-\alpha_i)+a(x-\alpha_1)\cdots(x-\alpha_{i-1})(x-\alpha_{i+1})\dots(x-\alpha_n) \end{aligned}\] となり、それに \(x=\alpha_i\) を代入すれば良い。
終結式と根の関係より、 \[\operatorname{res}(f,f')=a^{n-1}\prod_{i=1}^n f'(\alpha_i)\] が成り立つ。 これを使うと \[\begin{aligned} D(f)&=a^{2n-2}\prod_{i=1}^n \prod_{j=i+1}^n (\alpha_i-\alpha_j)^2 \\ &=a^{2n-2}({-1})^{n(n-1)/2}\prod_{i=1}^n \prod_{j=i, i\ne j}^n (\alpha_i-\alpha_j) \\ &=a^{2n-2}({-1})^{n(n-1)/2}\prod_{i=1}^n f'(\alpha_i)/a \\ &=a^{n-2}({-1})^{n(n-1)/2}\prod_{i=1}^n f'(\alpha_i) \\ &=\frac{({-1})^{n(n-1)/2}}{a}\operatorname{res}(f,f') \end{aligned}\] を得る。
\(D(f)\) が終結式を最高次の係数 \(a\) で割っているのは、 \(\operatorname{res}(f,f')\) が \(a\) で割り切れるからである。 というのは、シルベスター行列の行列式で書いたときに一番上の行に \(a\) が単独で現れるからだ。 つまり、次が言える：
系． \(R\)-係数多項式 \(f\in R[x]\) の判別式 \(D(f)\) は \(R\) の元となる。
「\(R\) の元となる」よりは「\(R\) の元として定義できる」と言った方が正確かもしれない。
\(R=\mathbb{Z}\) とすれば次が得られる：
系． 整数係数多項式の判別式は整数である。 特に、 \(D(f)\ne0\) （つまり、 \(f\) が無平方）の場合、 \(\left\lvert D(f)\right\rvert\ge 1\) が成り立つ。
根の距離の最小値
\(f\) の根の距離の最小値（英語では root separation と呼ぶ）を \[\operatorname{sep}(f)=\min_{i\ne j} \left\lvert \alpha_i-\alpha_j\right\rvert\] とおく。 最小値を実現する \(i\) と \(j\) を \(i=1\), \(j=2\) とおいても一般性を失わない。
\(f\) の実根を含む区間の幅を \(\operatorname{sep}(f)\) 未満に取れば、実根を確実に分離することができる。 そのため、 \(\operatorname{sep}(f)\) を下から評価したい。
根の限界（根の絶対値の上界；#2 で取り扱った）を \(M(f)\) で表すと、2つの根の距離を \[\left\lvert \alpha_i-\alpha_j\right\rvert\le 2M(f)\] と評価できる。 それを踏まえて \[D(f)=a^{2n-2}\prod_{i=1}^n \prod_{j=i+1}^n (\alpha_i-\alpha_j)^2\] の絶対値を取ると \[\begin{aligned} \left\lvert D(f)\right\rvert&=\left\lvert a\right\rvert^{2n-2}\prod_{i=1}^n \prod_{j=i+1}^n \left\lvert \alpha_i-\alpha_j\right\rvert^2 \\ &=\left\lvert a\right\rvert^{2n-2}\cdot\left\lvert \alpha_1-\alpha_2\right\rvert^2\cdot\prod_{i=1}^n \prod_{j=i+1,(i,j)\ne(1,2)}^n \left\lvert \alpha_i-\alpha_j\right\rvert^2 \\ &\le\left\lvert a\right\rvert^{2n-2}\cdot\left\lvert \alpha_1-\alpha_2\right\rvert^2\cdot\prod_{i=1}^n \prod_{j=i+1,(i,j)\ne(1,2)}^n (2M(f))^2 \\ {}&=\left\lvert a\right\rvert^{2n-2}\cdot\left\lvert \alpha_1-\alpha_2\right\rvert^2\cdot(2M(f))^{n(n-1)-2} \end{aligned}\] より、 \[\operatorname{sep}(f)=\left\lvert \alpha_1-\alpha_2\right\rvert \ge\sqrt{\frac{\left\lvert D(f)\right\rvert}{\left\lvert a\right\rvert^{2n-2}(2M(f))^{n^2-n-2}}} =\frac{\sqrt{\left\lvert D(f)\right\rvert}}{\left\lvert a\right\rvert^{n-1}(2M(f))^{n(n-1)/2-1}} \] という評価を得る。 右辺は \(f\) の係数のみから計算できる。
さらに、 \(f\) が整数係数で無平方の場合、 \(\left\lvert D(f)\right\rvert\ge 1\) なので、判別式の計算を省いて \[\operatorname{sep}(f)\ge\frac{1}{\left\lvert a\right\rvert^{n-1}(2M(f))^{n(n-1)/2-1}}\] という評価ができる。
こういう「代数的数同士の距離」というトピックは、超越数論にも関わってくるようである。 今回はこれ以上は掘り下げない。
代数的実数の演算の実装
四則演算
まずは和から考えよう。
定義多項式 \(f\) を持つ \(\alpha\) と、定義多項式 \(g\) を持つ \(\beta\) を考える。 #3 の結果より、和 \(\alpha+\beta\) を根に持つ多項式として \(h(x)=\operatorname{res}_y(f(x-y),g(y))\) を取ることができる。 必要であれば \(h\) を \(\gcd(h,h')\) で割ることによって無平方とできる。 あとは、その分離区間を求めたい。 分離区間を求める方法として、3通り考えられる。
方法1・実根の数え上げを使う： \(\alpha\) は計算可能実数、 \(\beta\) も計算可能実数なので、 \(\alpha+\beta\) も計算可能実数である。 \(\alpha+\beta\) に収束する縮小区間列 \(\{[a_n,b_n]\}\) を取る。 \(h\) のスツルム列を計算すれば、それぞれの区間 \([a_n,b_n]\) に何個実根が存在するかがわかる。 \(n\) を十分大きく取れば、区間 \([a_n,b_n]\) に存在する実根を1個にできる（つまり、 \([a_n,b_n]\) が分離区間になる）。
方法2・根の距離を使う方法1： 方法1と同様に、 \(\alpha+\beta\) に収束する縮小区間列 \(\{[a_n,b_n]\}\) を取る。 さっきの節の結果を使うと、 \(\operatorname{sep}(h)\) の下からの評価 \(m\) が計算できる。 区間 \([a_n,b_n]\) の幅が \(m\) 未満になるように \(n\) を十分大きく取れば、 \([a_n,b_n]\) が分離区間になる。
方法3・根の距離を使う方法2： \(\alpha\) の分離区間 \((a,b)\) を幅 \(m/2\) 未満で計算、 \(\beta\) の分離区間 \((c,d)\) を幅 \(m/2\) 未満で計算する。 すると、 \(\alpha+\beta\) の分離区間は \((a+c,b+d)\) と取れる。 （この方法は和や差の場合は良いが、積の場合は幅を事前に見積もるのが難しく、使いづらい。）
ここでは、方法1を採用することにする。 この方法の欠点は、スツルム列の計算が必要になることである。
なお、単項マイナス negate と逆数 recip に関しては、分離区間を計算し直す必要はない。
AlgReal 型に対する各種演算の実装は AlgReal.hs に追加で書き込んでいくことにする。
多項式と分離区間 \((a,b]\) から AlgReal 型の代数的実数を構築する関数 mkAlgReal はすでに用意しているが、追加で補助関数を用意しておく：
mkAlgRealWithInterval :: UniPoly Rational -> Interval -> AlgReal
mkAlgRealWithInterval f (Iv a b) = mkAlgReal f a b
-- | 与えられた無平方多項式と、その実根を表す計算可能実数から、代数的実数を構築する
mkAlgRealWithCReal :: UniPoly Rational -> CReal -> AlgReal
mkAlgRealWithCReal f (CReal xs) = mkAlgRealWithInterval f (head $ dropWhile (\(Iv a b) -> countRealRootsBetween a b f >= 2) xs)では、 AlgReal 型に対する Num クラスのインスタンスを作っていこう。 先にコードを提示し、そのあとに解説する。
instance Num AlgReal where
negate (FromRat x) = FromRat (negate x)
negate (AlgReal f s a b) = AlgReal (compP f (-ind)) (-s) (-b) (-a)まずは単項マイナスである。 単項マイナスに関しては、分離区間は単純に区間として符号を逆転させるだけで良い：\((a,b)\mapsto({-b},{-a})\)
定義多項式は \(f(x)\) を \(f({-x})\) に変える。コードで書くと compP f (-ind) であるが、この実装に関しては改善の余地がある。
前回、 AlgReal データ構築子の定義に「多項式の符号の変化」を表すフィールド（\(f'(\alpha)\) の符号と同じ）を追加したが、それも逆になる（符号を反転させる）。
-- Num クラス 続き（和）
FromRat x + FromRat y = FromRat (x + y)
FromRat k + AlgReal f s a b = mkAlgReal (compP f (ind - constP k)) (a + k) (b + k)
x@(AlgReal _ _ _ _) + y@(FromRat _) = y + x
x + y = mkAlgRealWithCReal (squareFree $ resultant_poly f_x_y g) (algRealToCReal x + algRealToCReal y)
where f = mapCoeff constP $ definingPolynomial x
f_x_y = compP f (constP ind - ind) -- f(x-y)
g = mapCoeff constP $ definingPolynomial y【12月17日 修正】FromRat と AlgReal の和がゼロになる場合を考慮し忘れていたので修正。
和の実装を見ていこう。 有理数同士の場合は Rational 型の和を使う。
有理数 \(k\) と代数的実数 \(\alpha\) を足す場合は、定義多項式は \(f(x)\) から \(f(x-k)\) になる（\(f(\alpha)=0\) ならば \(f((\alpha+k)-k)=0\) なので）。 \(f'((\alpha+k)-k)=f'(\alpha)\) なので、符号の変化はそのままである。 分離区間は、 \(k\) だけ平行移動する：\((a,b)\mapsto(a+k,b+k)\)
代数的実数同士の和が肝心で、定義多項式は \(\operatorname{res}_y(f(x-y),g(y))\) （を無平方にしたもの）となる。 この際、2変数多項式環 \(\mathbb{Q}[x][y]\) が登場するわけだが、2つの不定元は、内側の \(x\) が constP ind :: UniPoly (UniPoly Rational)、外側の \(y\) が ind :: UniPoly (UniPoly Rational) と書ける。
終結式の計算には、前々回の演習問題とした resultant_poly 関数を使う。 分離区間の計算には、一旦計算可能実数を経由する。
-- Num クラス 続き（差）
FromRat x - FromRat y = FromRat (x - y)
FromRat k - AlgReal f s a b = mkAlgReal (compP f (constP k - ind)) (k - b) (k - a)
AlgReal f s a b - FromRat k = mkAlgReal (compP f (ind + constP k)) (a - k) (b - k)
x - y = mkAlgRealWithCReal (squareFree $ resultant_poly f_x_y g) (algRealToCReal x - algRealToCReal y)
where f = mapCoeff constP $ definingPolynomial x
f_x_y = compP f (constP ind + ind) -- f(y+x)
g = mapCoeff constP $ definingPolynomial y【12月17日 修正】FromRat と AlgReal の差がゼロになる場合を考慮し忘れていたので修正。
差も、和とほとんど同じである。
-- Num クラス 続き（積）
FromRat x * FromRat y = FromRat (x * y)
FromRat k * AlgReal f s a b
| k == 0 = 0
| k > 0 = AlgReal (compP f (scaleP (recip k) ind)) s (a * k) (b * k)
| k < 0 = AlgReal (compP f (scaleP (recip k) ind)) (-s) (b * k) (a * k)
x@(AlgReal _ _ _ _) * y@(FromRat _) = y * x
x * y = mkAlgRealWithCReal (squareFree $ resultant_poly y_f_x_y g) (algRealToCReal x * algRealToCReal y)
where f = definingPolynomial x
y_f_x_y = UniPoly $ V.reverse $ V.imap (\i a -> constP a * ind^i) $ coeff f -- y^n f(x/y)
g = mapCoeff constP $ definingPolynomial y積も似たような感じで実装するが、使う終結式は \(\operatorname{res}_y(y^nf(x/y),g(y))\) と、やや複雑になる。 一見すると分数（有理式）が出てくるように見えるが、 \[f(x)=a_nx^n+\cdots+a_1x+a_0\] と書けば \[y^nf(x/y)=a_nx^n+\cdots+a_ix^iy^{n-i}+\cdots+a_1xy^{n-1}+a_0y^n\] なので、有理式にはならない。 計算には、注意深く多項式の係数を操作すれば良い。
-- Num クラスの残り
abs x | x >= 0 = x
| otherwise = negate x
signum x | x > 0 = 1
| x == 0 = 0
| x < 0 = -1
fromInteger n = FromRat (fromInteger n)Num クラスの残りの演算は適当に実装する。
instance Fractional AlgReal where
recip (FromRat x) = FromRat (recip x)
recip (AlgReal f s a b) = AlgReal (UniPoly $ V.reverse $ coeff f) s' (recip b) (recip a)
where s' | even (degree' f) || 0 < a = -s
| otherwise = s
fromRational = FromRatFractional クラスも実装する。 fromRational は FromRat データ構築子をそのまま使えるので、あとは逆数関数 recip の実装である（除算に関してはデフォルト実装を使えば良い）。
逆数に関する定義多項式は、係数を反転させる。 区間も、区間として逆数を取れば良い（AlgReal データ構築子における区間は、ゼロを含まないようなものになっていることに注意）。
符号の変化はどうだろうか。 \(\alpha\) の定義多項式が次数 \(n\) の \(f\) で \(\beta=1/\alpha\) の定義多項式を \(h(x)=x^nf(1/x)\) とすると、 \[h'(x)=nx^{n-1}f(1/x)-x^{n-2}f'(1/x)\] より \[h'(\beta)=n\beta^{n-1}f(\alpha)-\alpha^{n-2}f'(\alpha)=-\alpha^{n-2}f'(\alpha)\] となる。つまり、
- \(n\) が偶数または \(\alpha\) が正の場合は、 \(h'(\beta)\) の符号は \(f'(\alpha)\) と逆
- \(n\) が奇数かつ \(\alpha\) が負の場合は、 \(h'(\beta)\) の符号は \(f'(\alpha)\) と同じ
である（s' の計算）。
平方根、べき根
平方根や冪根（\(n\) 乗根）を計算する関数があると便利である。 これらは「多項式の実根を求める操作」の一種ではあるが、使用頻度が高いと考えられるので、用意しておく価値はあるだろう。
多項式の根は一般に複数あるが、 \(\sqrt[n]{\alpha}\) の符号については
- \(\alpha\) が非負の場合は、 \(\sqrt[n]{\alpha}\) も非負
- \(n\) が奇数で \(\alpha\) が負の場合は \(\sqrt[n]{\alpha}\) は負
- \(n\) が偶数で \(\alpha\) が負の場合は、エラー
とするのが普通だろう。
有理数の平方根
有理数 \(a\) に対して、多項式 \(x^2-a\) の実根のうち、正のもの（区間 \((0,+\infty)\) にあるもの）を返す。
sqrtQ :: Rational -> AlgReal
sqrtQ a | a > 0 = case realRootsBetween (ind^2 - constP a) (Finite 0) PositiveInfinity of
[sqrt_a] -> sqrt_a
_ -> error "sqrt: none or multiple roots"
| a == 0 = 0
| otherwise = error "sqrt: negative"有理数のべき根
\(x^n-a\) の実根のうち、 \(x\) と符号が同じもの（区間 \((-\infty,0]\) または \((0,+\infty)\) にあるもの）を返す。
nthRootQ :: Int -> Rational -> AlgReal
nthRootQ !n !a
| n == 0 = error "0th root"
| n < 0 = nthRootQ (-n) (recip a)
| a > 0 = case realRootsBetween (ind^n - constP a) (Finite 0) PositiveInfinity of
[b] -> b
l -> error ("nthRoot: none or multiple roots " ++ show l)
| a == 0 = 0
| odd n = case realRootsBetween (ind^n - constP a) NegativeInfinity (Finite 0) of
[b] -> b
l -> error ("nthRoot: none or multiple roots " ++ show l)
| otherwise = error "nthRoot: negative"代数的実数の平方根
定義多項式 \(f\) に \(x^2\) を代入したもの \(f(x^2)\) は、 \(\pm\sqrt{\alpha}\) を根に持つ多項式となる。 \(f(x^2)\) の実根の中から、2乗したものが \(\alpha\) の分離区間に入るものを返せば良い。
sqrtA :: AlgReal -> AlgReal
sqrtA (FromRat x) = sqrtQ x
sqrtA x = case filter (\y -> FromRat a < y^2 && y^2 <= FromRat b) $ realRootsBetween (compP f (ind^2)) (Finite 0) PositiveInfinity of
[sqrtx] -> sqrtx
r -> error $ "sqrt: none or multiple roots" ++ show r
where f = definingPolynomial x
(a,b) = isolatingInterval x代数的実数のべき根
定義多項式 \(f\) に \(x^n\) を代入したもの \(f(x^n)\) は、 \(\sqrt[n]{\alpha}\) を根に持つ多項式となる。 \(f(x^n)\) の実根の中から、 \(n\) 乗したものが \(\alpha\) の分離区間に入るものを返せば良い。
nthRootA :: Int -> AlgReal -> AlgReal
nthRootA !n (FromRat x) = nthRootQ n x
nthRootA !n x
| n == 0 = error "0th root"
| n < 0 = nthRootA (-n) (recip x)
-- now n must be positive
| x == 0 = 0
| x > 0 = case filter (\x -> FromRat a < x^n && x^n <= FromRat b) $ realRootsBetween (compP f (ind^n)) (Finite 0) PositiveInfinity of
[rx] -> rx
_ -> error "nthRoot: none or multiple roots"
-- x must be negative
| odd n = case filter (\x -> FromRat a < x^n && x^n <= FromRat b) $ realRootsBetween (compP f (ind^n)) NegativeInfinity (Finite 0) of
[rx] -> rx
_ -> error "nthRoot: none or multiple roots"
| otherwise = error "nthRoot: negative"
where f = definingPolynomial x
(a,b) = isolatingInterval xべき乗、有理数係数多項式への代入
かけ算を実装済みなので、それを使った冪乗（\(n\) 乗）の実装ももちろんできる。 やってみよう。
例えば、 \(x^3-4x+2\) の根 \(\alpha\) の10乗を計算するとしたら、皆さんはどうするだろうか？ 指数10を二進法で表して \[\alpha^{10}=\left(\alpha\cdot\left(\alpha^2\right)^2\right)^2\] により4回のかけ算で済ませる！……と答えた方は不正解である。 かけ算を9回行う、と答えた方ももちろん不正解である。
たった4回のかけ算、と思われるかもしれないが、現在の素朴な（因数分解をしない）実装では、結果として得られる代数的数の定義多項式の次数は \((3\cdot(3^2)^2)^2=3^{10}=59049\) となってしまう。 このように次数が大きいと、計算には非常に時間がかかる。 （演算に必要な時間や空間が一定ならばこのような繰り返し乗算による方法はうまくいくが、代数的数には通用しない）
賢い読者は、関係式 \(\alpha^3-4\alpha+2=0\) を利用することを思いついたかもしれない。 この関係式を使って \(\alpha\) の冪乗を計算すると、 \[\begin{aligned} \alpha^3&=4\alpha-2, \\ \alpha^5&=(4\alpha-2)\alpha^2=-2\alpha^2+16\alpha-8, \\ \alpha^{10}&=(\alpha^5)^2=(-2\alpha^2+16\alpha-8)^2=8(38\alpha^2-65\alpha+24) \end{aligned}\] となる。 手計算でやってみるとなかなか大変だが、やっていることは有理数係数の多項式除算である。
あとは \(8(38\alpha^2-65\alpha+24)=8((38\alpha-65)\alpha+24)\) に対して代数的実数のかけ算等を使えば良い。 この例では、代数的実数どうしの演算は1回しか使わない。
「関係式 \(\alpha^3-4\alpha+2=0\) を使って \(\alpha^{10}\) を計算する」というのは、言い方を変えれば「\(\alpha\) についての多項式 \(\alpha^{10}\) を \(\alpha^3-4\alpha+2\) で割った余りを計算する」ということである。 多項式除算はすでに実装済みで、何回も使っている。
Haskell で実装すると次のようになるだろう：
powIntA :: AlgReal -> Int -> AlgReal
powIntA _ 0 = 1
powIntA x n | n < 0 = recip $ powIntA x (-n)
powIntA (FromRat x) n = FromRat (x^n)
powIntA x n = case (ind^n) `modP` f of
g -> valueAt x (mapCoeff FromRat g)
where f = definingPolynomial xcase 式以下が、多項式除算を使う部分である。
同じ理屈で、「次数の高い有理数係数多項式に代数的実数を代入する操作」を多項式除算を使って実装できる。
演習問題：有理数係数多項式に代数的実数を代入したものを計算する関数 valueAtA :: AlgReal -> UniPoly Rational -> AlgReal を実装せよ。
有理数乗
冪根と（整数による）冪乗を実装したので、有理数乗 \(\alpha^{\frac{p}{q}}\) も実装できる。
powRatA :: AlgReal -> Rational -> AlgReal
powRatA x y = nthRootA (fromInteger $ denominator y) (powIntA x (fromInteger $ numerator y))ちなみに、有理数乗 \(\alpha^{\frac{p}{q}}\) の定義多項式は、終結式を使って \(\operatorname{res}_y(f(y),x^q-y^p)\) と計算できる。
遊んでみよう
実装が済んだら、早速試してみよう。 $ stack repl で REPL を起動する。
$ stack repl
algebraic-real-0.1.0.0: initial-build-steps (lib)
Configuring GHCi with the following packages: algebraic-real
GHCi, version 8.0.2: http://www.haskell.org/ghc/ :? for help
[1 of 5] Compiling Numeric.AlgebraicReal.UniPoly ( .../src/Numeric/AlgebraicReal/UniPoly.hs, interpreted )
...
[5 of 5] Compiling Numeric.AlgebraicReal.AlgReal ( .../src/Numeric/AlgebraicReal/AlgReal.hs, interpreted )
Loaded GHCi configuration from /private/var/folders/0c/.../T/ghci41849/ghci-script
*Numeric.AlgebraicReal.UniPoly Numeric.AlgebraicReal.AlgReal Numeric.AlgebraicReal.CReal Numeric.AlgebraicReal.Interval Numeric.AlgebraicReal.Resultant Numeric.AlgebraicReal.UniPoly> :set prompt "> "
> sqrtQ 2
AlgReal (UniPoly [(-2) % 1,0 % 1,1 % 1]) 1 (3 % 4) (3 % 2)
> 2 :: AlgReal
FromRat (2 % 1)
> (sqrtQ 2)^2
AlgReal (UniPoly [(-4) % 1,0 % 1,1 % 1]) 1 (9 % 16) (9 % 4)
> (sqrtQ 2)^2 == 2
True\(\sqrt{2}\) を2乗すると2となることが確認できた（ただし、 2 と (sqrtQ 2)^2 では表示させた時の見た目は異なる）。
> sqrtQ 2 + sqrtQ 8
AlgReal (UniPoly [36 % 1,0 % 1,(-20) % 1,0 % 1,1 % 1]) 1 (3 % 1) (6 % 1)
> sqrtQ 2
AlgReal (UniPoly [(-2) % 1,0 % 1,1 % 1]) 1 (3 % 4) (3 % 2)
> 2 * sqrtQ 2
AlgReal (UniPoly [(-2) % 1,0 % 1,1 % 4]) 1 (3 % 2) (3 % 1)
> sqrtQ 8
AlgReal (UniPoly [(-8) % 1,0 % 1,1 % 1]) 1 (9 % 4) (9 % 2)
> 2 * sqrtQ 2 == sqrtQ 8
True等式 \(2\sqrt{2}=\sqrt{8}\) が確認できた。
係数膨張
もっと複雑な代数的数も扱ってみよう。
> let x = ind
> let [a0,a1,a2,a3] = realRoots $ (x-4)*(x-3)*(x-2)*(x+3)+1
> a0
AlgReal (UniPoly [(-71) % 1,54 % 1,(-1) % 1,(-6) % 1,1 % 1]) (-1) ((-9) % 2) ((-9) % 4)
> a1
AlgReal (UniPoly [(-71) % 1,54 % 1,(-1) % 1,(-6) % 1,1 % 1]) 1 (9 % 8) (9 % 4)
> a2
AlgReal (UniPoly [(-71) % 1,54 % 1,(-1) % 1,(-6) % 1,1 % 1]) (-1) (9 % 4) (27 % 8)
> a3
AlgReal (UniPoly [(-71) % 1,54 % 1,(-1) % 1,(-6) % 1,1 % 1]) 1 (27 % 8) (9 % 2)多項式 \((x-4)(x-3)(x-2)(x+3)+1\) から適当に代数的数をでっち上げる。
根どうしの和を計算してみよう。
…
筆者の環境ではしばらく時間がかかり、画面が数字で埋め尽くされた。
これは、すでに何回か触れた、係数膨張である。 この膨大な数字が虚仮威しに過ぎないことは、モニック多項式に変換する（定数倍する）ことでわかる：
> toMonic (definingPolynomial it)
UniPoly [11778048 % 5,(-97967232) % 5,191816064 % 5,(-23103648) % 1,(-57623376) % 5,130774488 % 5,(-17815540) % 1,32400936 % 5,(-6375679) % 5,399798 % 5,111231 % 5,(-31524) % 5,3631 % 5,(-42) % 1,1 % 1]これまでに載せた通りの実装では、終結式の計算と、無平方成分の計算で係数膨張が起こる。 そこで、それぞれの計算後に、モニック多項式への変換 toMonic を挟む。
演習問題 適宜モニック多項式に変換することにより、係数膨張を軽減せよ。 具体的には、 UniPoly.hs の squareFree 関数と、さっき実装した各演算の resultant_poly の呼び出しの部分を編集する。
因数分解へ
これで代数的実数についての演算を一通り実装できた。 だが、大きな課題として、定義多項式の因数分解をしていないために多項式の次数が不必要に増大してしまうという欠点がある。
例えば、4の平方根 \(\sqrt{4}\) を計算させると有理数 FromRat 2 ではなく2次の代数的数 AlgReal (UniPoly [(-4) % 1,0 % 1,1 % 1]) 1 (5 % 4) (5 % 2) が得られる（もちろん、有理数 2 と比較すれば「等しい」と判定される）。 また、 \(\sqrt{2}+\sqrt{8}\) の次数は2（\(3\sqrt{2}\) なので）となるべきであるが、今の実装で計算すると次数が4となってしまう。
これらの問題を解決するには、終結式の計算で得られた多項式を因数分解し、演算結果の代数的実数に対する最小多項式を得ることが必要である。 （ちなみに、有理数かどうかの判定に関しては、因数分解しなくても判別する方法はある。）
というわけで早速次回から多項式の因数分解に取り掛かりたいところだが、その前に、多項式の係数を整数係数に変えておきたい。 今の実装では、代数的実数の定義に使う多項式は有理数係数としているが、実際のところ、そこは整数係数で十分である。 それに、整数係数の方が何かと効率的であるし、有理数係数多項式の因数分解においても、結局は整数係数多項式の因数分解に帰着させることになる。
そこで、次回は多項式の係数を整数係数化し、その次から多項式の因数分解に取り掛かることにする。