#5 区間演算と計算可能実数
この連載では、代数的実数を分離するのに区間を使っている。 今回は、区間について掘り下げよう。
区間演算
#0 イントロダクション では、計算機上で実数を正確に取り扱うための方法として、計算可能実数を紹介した。 それと関連するが少し別の方向性として、区間演算 (interval arithmetic) というものがある。
計算機上で浮動小数点数を使って実数を近似することの問題点は、誤差が発生することだった。 しかし、誤差の大きさを正確にコントロールできているのであれば、ある意味で「正確に実数を取り扱えている」と言える。 例えば、ある実数を計算した結果が 3.14\pm 0.01 の範囲にある、ということが確実に分かるのなら、「その実数は3よりも大きい」という性質は確実に言える。 ただし「その実数は3.14より大きい」という性質に関しては「成り立つかどうかわからない」。 （単なる近似計算の場合は、誤差の大きさが不明なので、これらの性質が成り立つかどうかを保証することはできない）
数学的には、このような誤差のコントロールは区間を使って実現できる（有効数字などという雑な概念は使わない！）。 つまり 3.14\pm 0.01 ではなく [3.13,3.15] と表す。 計算機上では、区間の端点としては有理数や浮動小数点数を使う。 この文脈では、浮動小数点数は「実数の近似」ではなく、「実数の離散部分集合」として取り扱うことになる（特に、浮動小数点数同士の等号比較が意味を持つ）。
区間に対して四則演算等を定義することによって、既存の数式やアルゴリズムを機械的に変換して区間に対応させることができる（区間演算）。
まずは、有理数の区間を使った区間演算を実装しよう。
四則演算
区間に対しての四則演算は次のように定義できる：
\begin{aligned} {[a,b]}+[c,d]&=[a+c,b+d], \\ {[a,b]}-[c,d]&=[a-d,b-c], \\ {[a,b]}\cdot[c,d]&=[\min\{ac,ad,bc,bd\},\max\{ac,ad,bc,bd\}] \end{aligned}
除算は、符号が確定している（区間がゼロを含まない）場合のみ行える： \frac{1}{[a,b]}=\begin{cases} [1/b,1/a] & (0<a\text{ または } b<0) \\ \text{ゼロ除算} & (\text{その他}) \end{cases}
いずれも、 [a,b]\star[c,d]\supseteq\{x\star y\in\mathbb{R}\mid x\in[a,b],y\in[c,d]\} となるように定めている。
Haskellで実装するならば、次のようになる（src/Numeric/AlgebraicReal/Interval.hs に保存）：
module Numeric.AlgebraicReal.Interval where
-- Iv a b represents the interval [a,b]
data Interval = Iv !Rational !Rational deriving (Show)
width :: Interval -> Rational
width (Iv a b) = b - a
instance Num Interval where
negate (Iv a b) = Iv (-b) (-a)
Iv a b + Iv c d = Iv (a + c) (b + d)
Iv a b - Iv c d = Iv (a - d) (b - c)
Iv a b * Iv c d = Iv (minimum [a*c,a*d,b*c,b*d]) (maximum [a*c,a*d,b*c,b*d])
abs x@(Iv a b) | 0 <= a = x
| b <= 0 = -x
| otherwise = Iv 0 (max (-a) b)
signum (Iv a b) | 0 < a = 1 -- 正の場合
| b < 0 = -1 -- 負の場合
| a == 0 && b == 0 = 0 -- 0 の場合
| 0 <= a = Iv 0 1 -- 非負の場合
| b <= 0 = Iv (-1) 0 -- 非正の場合
| otherwise = Iv (-1) 1 -- 符号不明の場合
fromInteger n = Iv (fromInteger n) (fromInteger n)
instance Fractional Interval where
recip (Iv a b) | 0 < a || b < 0 = Iv (recip b) (recip a)
| otherwise = error "divide by zero"
fromRational x = Iv x x比較
通常の実数は全順序が定まっており、2つの実数に対して大小等しいのいずれかが確定するが、区間に対してはそうはいかない。 ある実数 \alpha を計算した時に区間 [3.14,3.15] が得られ、別の実数 \beta を計算した時に区間 [3.13,3.15] が得られたとしても、 \alpha=\beta が成り立つかはわからないし、 \alpha<\beta や \alpha>\beta が成り立つかもわからない。
しかし、次は言える：
- b<c ならば [a,b]<[c,d],
- b\le c ならば [a,b]\le[c,d],
- a>d ならば [a,b]>[c,d],
- a\ge d ならば [a,b]\ge[c,d],
- a=b=c=d ならば [a,b]=[c,d] （つまり、区間の幅が 0 の場合は正確な比較ができる）
それ以外の場合は順序が定まらない。
Haskell で実装するなら、 Maybe Ordering を返すような比較関数となるだろう：
maybeCompareInterval :: Interval -> Interval -> Maybe Ordering
maybeCompareInterval (Iv a b) (Iv c d)
| b < c = Just LT
| d < a = Just GT
| a == d && b == c = Just EQ
-- b <= c = LE
-- d <= a = GE
| otherwise = NothingHaskell の Ordering 型は EQ (equal; =), LT (less than; <), GT (greater than; >) の3つの値しか取れず、 \le や \ge を表現できないので、それらの場合は実装していない。
計算例
関数 f(x)=x^3+2x^2-2 に区間 [3.13,3.15] を代入して計算してみる：
> let f x = x^3 + 2 * x^2 - 2
> f (Iv 3.13 3.15)
Iv (48258097 % 1000000) (392807 % 8000)計算結果として [48258097/1000000, 392807/8000]=[48.258097, 49.100875] という区間が得られた。
なお、数学的に同じ関数でも、計算方法によって結果が変わる場合がある。 例えば、次の g1 と g2 はいずれも g(x)=x^3-2x^2-2 という関数を表すが、得られる区間は後者の方が狭い：
> let g1 x = x^3 - 2 * x^2 - 2
> let g2 x = (x - 2) * x^2 - 2
> g1 (Iv 3.13 3.15)
Iv (8819297 % 1000000) (386483 % 40000)
> g2 (Iv 3.13 3.15)
Iv (9070497 % 1000000) (75287 % 8000)もっと極端な例を挙げると、 h(x)=x-x という関数は数学的には常に0だが、区間演算では（引数の幅が0でなければ）結果は有限の幅を持った区間となる：
> let h x = x - x
> h (Iv 3.13 3.15)
Iv ((-1) % 50) (1 % 50)区間に対する関数値の計算 f([a,b]) において、区間幅の増大を抑えるために平均値形式 (mean value form) と呼ばれる方法が知られている： f(c)+f'([a,b])([a,b]-c), \quad c=\frac{a+b}{2}
浮動小数点数の場合
有理数の場合は、区間の端点の演算を正確に行える。 しかし、浮動小数点数を使って計算する場合は、そもそも区間の端点を正確に表せない。
例えば、区間 [a,b]=[1,1.5] と [c,d]=[2^{-53}+2^{-55},2^{-51}+2^{-53}] の和を考えよう： [a,b]+[c,d]=[a+c,b+d]=[1+2^{-53}+2^{-55},1.5+2^{-51}+2^{-53}] 有理数の場合はこれで良いが、精度53ビットの2進浮動小数点数を使う場合は a+c=1+2^{-53}+2^{-55} や b+d=1.5+2^{-51}+2^{-53} を正確に表現できないので困る。 具体的には、（最近接丸めによって） a+c\leadsto 1+2^{-52}, b+d\leadsto 1.5+2^{-51} となり、下端が不正に大きく、上端が不正に小さくなってしまう。
![区間 [c,d]](../images/05-fig-1.svg)
例えば Haskell では次のように確かめられる：
Prelude> let a = 1 :: Double
Prelude> let b = 1.5 :: Double
Prelude> let c = 2^^(-53) + 2^^(-55) :: Double
Prelude> let d = 2^^(-51) + 2^^(-53) :: Double
Prelude> a + c == 1 + 2^^(-52)
True
Prelude> b + d == 1.5 + 2^^(-51)
True演習問題 普段使いのプログラミング言語で同じことを試せ。つまり、 a=1, b=1.5, c=2^{-53}+2^{-55}, d=2^{-51}+2^{-53} についての a+c と b+d を浮動小数点数を用いて計算せよ。 上記の Haskell での実行例と異なる結果が得られた場合は、なぜそうなるかを考察せよ。
しかし、浮動小数点数の丸めモードを制御できれば、下端を計算する時に「なるべく小さくなるように丸める」、上端を計算する時に「なるべく大きくなるように丸める」、ということができる： \begin{aligned} a+c&\stackrel{\text{down}}{\leadsto} 1, \\ b+d&\stackrel{\text{up}}{\leadsto} 1.5+2^{-51}+2^{-52}. \end{aligned}
式と図の中で \stackrel{\text{down}}{\leadsto} と \stackrel{\text{up}}{\leadsto} を使って表したものは方向付き丸め (directed rounding) で、前者は負の無限大方向の丸め、後者は正の無限大方向の丸めと呼ばれる。
これを利用すれば、浮動小数点数を端点とする区間演算を（有理数の場合よりは幅が広がってしまうものの）正しく行うことができる。
これらの丸めモードについては、浮動小数点数の規格である IEEE754 で定められている。 IEEE754 では、丸めモードとして
- 最近接（偶数への）丸め (roundTiesToEven)
- 負の無限大方向の丸め (roundTowardNegative)
- 正の無限大方向の丸め (roundTowardPositive)
の他に、
- ゼロ方向の丸め (roundTowardZero) と
- 最近接（無限大への）丸め (roundTiesToAway; IEEE754-2008 より)
が定められており、以下の各演算について正しく丸められることが要請されている：
- 四則演算および FMA （融合積和; Fused Multiply-Add）
- 平方根 (
sqrt) - 剰余（C言語の
remainder関数） - 整数との変換
- 浮動小数点数の異なる形式間の変換
これ以外の数学関数、例えば exp や sin は値が正しく丸められるとは限らない。
浮動小数点数の丸めモードを制御するには
プログラミング言語において浮動小数点数の丸めモードを制御する方法であるが、言語によって利用できる場合と利用できない場合がある。
C言語では <fenv.h> で提供される fegetround / fesetround 関数を使って丸め方向を取得、設定できる（Visual C++ の古いバージョンではこれらの関数が正しく実装されていないので注意すること。最新の VS2017 ではどうやら修正されてしまったらしく、ディスりの材料にできない）。
Fortran の場合は IEEE_ARITHMETIC モジュールの IEEE_GET_ROUNDING_MODE, IEEE_SET_ROUNDING_MODE を使うようだ（筆者は Fortran を書かないので、実際に試したわけではない）。
最近の言語では、 Julia でも丸め方向を制御できる（setrounding 関数）。
残念ながら、この連載で主に使っている Haskell では、丸め方向を制御できない。 浮動小数点数の丸めモードを変更することはグローバルな状態を触ることであり、純粋な計算とは相性が悪い。 （もちろん、四則演算と丸め方向の指定をセットにした関数をC言語で実装して、HaskellからFFIでそれを呼び出すということはできる。）
課題 読者の普段使いの言語が浮動小数点数の丸めモードの変更に対応しているかどうか調べよ。
なお、言語組み込みの（ハードウェアで実装されることの多い）浮動小数点数にこだわらなければ、多倍長浮動小数点数のポピュラーな実装である GNU MPFR が IEEE754 準拠の丸めモードの実装を謳っている。
この連載では、当面は浮動小数点数による区間演算は実装しないことにする。
計算可能実数
計算可能実数については #0 イントロダクション でふわっとした紹介をした。 ここではもう少し精密な形で説明する。
まず、計算可能関数であるが、ざっくりと、自然数から自然数への数学的な関数のうちプログラミング言語を使って実装できるようなもの、と思ってもらって良いだろう。 より厳密な定義が知りたければ、ゲーデルの不完全性定理を扱っているようなまともな論理学の教科書なら必ず計算可能関数（または帰納的関数、再帰的関数）の定義が書かれているはずなので、そういう本を参照せよ（本当に入門的な本だと、ここで書いたようなざっくりとした定義で済ませている場合もあるが…）。
計算可能性は自然数から自然数への関数について定義されるが、それを拡張して整数や有理数およびその組についての関数（\mathbb{Z}\to\mathbb{Q} や \mathbb{Z}\to\mathbb{Q}\times\mathbb{Q} など）についても計算可能性を定義できる。 また、整数の列や有理数列は、自然数から整数や有理数への関数だと思えるので、計算可能性が定義できる。
実数 \alpha が計算可能実数 (computable real number, または単に computable number) であるとは、ある計算可能な有理数列 \{a_n\} と計算可能関数 e\colon\mathbb{N}\to\mathbb{N} の組が存在して、任意の自然数 n_0, n に対して e(n_0)\le n ならば \left\lvert\alpha-a_n\right\rvert<2^{-n_0} が成り立つこととする。 要するに、 \varepsilon-N 論法の \varepsilon を 2^{-n_0} で、「数列」を計算可能な数列で、「存在」を計算可能な関数で置き換えたものである。
別の言い方をすれば、 \alpha が計算可能実数であるとは、計算可能な有理数の区間列 \{[a_n,b_n]\} であって \begin{gathered} a_n\le a_{n+1}\le b_{n+1}\le b_n \\ b_n-a_n < 2^{-n} \\ a_n\to\alpha, b_n\to\alpha \end{gathered} を満たすもの、と定義することもできる（TODO: 定義の同値性？）。
（余談：計算可能実数の定義に関して、「その実数の2進小数表示（10進小数でも良い）を出力するプログラムが存在する」の方が直感的だと思われるかもしれないが、それでは小数表示の一意性が成り立たないこと (1=0.999\cdots) による面倒が発生する。）
ただ、計算機上で実装する立場としては、そもそも「実数」なんてものはないので、この定義をそのまま使うわけにはいかない。 代わりに、「適切な条件を満たす組 (\{a_n\},e)」なり「適切な条件を満たす区間列 \{[a_n,b_n]\}」を適当な同値関係で割ったものを計算可能実数と呼ぶことにする（有理数列を使って実数を構成する方法を思い出せ）。
定義のやり方はいくつか考えられる（コーシー列、区間列、デデキント切断など）が、ここでは、区間列を使って計算可能実数を定義する。
実装例 (src/Numeric/AlgebraicReal/CReal.hs):
module Numeric.AlgebraicReal.CReal where
import Numeric.AlgebraicReal.Interval
import Data.Maybe
newtype CReal = CReal [Interval]演算は、列のそれぞれの要素について区間演算を行えば良い。
instance Num CReal where
negate (CReal xs) = CReal (map negate xs)
CReal xs + CReal ys = CReal (zipWith (+) xs ys)
CReal xs - CReal ys = CReal (zipWith (-) xs ys)
CReal xs * CReal ys = CReal (zipWith (*) xs ys)
abs (CReal xs) = CReal (map abs xs)
signum (CReal xs) = error "signum on CReal is not computable"
fromInteger n = CReal $ repeat (fromInteger n)さっき書いた計算可能実数の定義では n 番目の区間の幅が 2^{-n} 以下になることを要請していたが、この実装例ではその辺を横着している。 また、計算可能実数の符号は一般には計算可能でない（ゼロ判定ができない）ので、 signum 関数は実装しない。
除算についてだが、区間演算の場合は分母の区間にゼロが含まれるとゼロ除算となるのだった。 計算可能実数の場合は、符号が確定するまで区間の幅を狭めて精度を上げるという操作ができる。 ゼロ除算の場合は無限ループに陥るが、そもそもゼロ除算はやってはいけない操作なのでその結果何が起ころうとも問題ないだろう。
instance Fractional CReal where
recip (CReal xs) = CReal $ map recip (dropWhile (\(Iv a b) -> a <= 0 && 0 <= b) xs)
fromRational x = CReal $ repeat (fromRational x)【2018年8月31日 修正】 recip の実装が間違っていたのを修正した。
計算可能実数の比較
区間演算においては順序が定まらない場合があったが、計算可能実数の場合は、区間を狭めて近似の精度を上げていけば順序を定めることができる。 ただし、2つの計算可能実数が等しい場合には、一般には無限ループに陥ってしまう。 なので、計算可能実数の実装する比較関数には名前に unsafe をつけることにする。
unsafeCompareCReal :: CReal -> CReal -> Ordering
unsafeCompareCReal (CReal xs) (CReal ys) = head $ catMaybes (zipWith maybeCompareInterval xs ys)より安全な（確実に停止する）操作としては、有理数に対しての「のりしろ付きの」比較ができる（TODO: 用語）。 例えば、計算可能実数 x について「0<x か x\le 0 か」は判断できないが、「0<x か x<0.01 か」は、どちらか成り立つ方を判断できる。 近似する区間の幅が 0.01 未満になるようにすれば少なくとも一方が成り立つからだ。 もちろん両方成り立つ場合は結果が一意的に定まらず、計算可能実数 x と x' が実数として同じでも、一方は「0<x」を返しもう一方は「x'<0.01」が帰ってくるかもしれない。 つまり、この比較操作は実数上の演算としては well-defined ではない。
TODO: 実装
浮動小数点数への変換
実数を浮動小数点数に丸める方法として IEEE754 で何種類か定められているのは先に書いたが、計算可能実数についてそれらの演算を（無限ループせずに）正しく実装することはできない。 これもゼロ判定（等号比較）が不可能なことに由来する。 もしも正しい丸め演算を実装できたら、それを利用してゼロ判定ができてしまう。
次善の策として、さっき書いた「のりしろ付きの」比較を応用して、精度 (0.5+\varepsilon)ulp の最近接（とは限らない）丸めと、精度 (1+\varepsilon)ulp の（最良とは限らない）方向付き丸めを実装できる（\varepsilon は正）。
注：ulp とは unit in the last place のことで、浮動小数点数の仮数の一番下の桁の大きさを表す。
TODO: 実装
代数的実数と計算可能実数
代数的実数の区間縮小：補足
前の記事 (#2) では分離区間を求めた後もスツルム列を使って近似値を計算していた。 しかし、一旦分離区間が分かってしまえば、定義多項式がその区間内に重複度1の根をただ1個もつということ、つまり区間内で1回だけ定義多項式の符号が変化するということが言えるので、スツルム列を使わなくても、定義多項式の符号を見るだけで実根との大小が分かってしまう。
AlgReal.hs の冒頭に次を追加し：
import Numeric.AlgebraicReal.Interval
import Numeric.AlgebraicReal.CReal次を追加する：
-- s: この代数的実数における f' の符号（正なら区間 (a,b) において f は負から正に変わり、負なら区間 (a,b) において f は正から負に変わる）
intervalsWithSign :: UniPoly Rational -> Int -> Rational -> Rational -> [Interval]
intervalsWithSign !f !s !a !b = Iv a b : ivs a b
where
ivs !a !b | signAt c f == 0 = repeat (Iv c c)
| s * signAt c f < 0 = Iv c b : ivs c b
| s * signAt c f > 0 = Iv a c : ivs a c
where c = (a + b) / 2#2 で書いた AlgReal 型の定義を、この「符号の変化」を保持するように書き換えておく：
-- | 代数的実数を表す型
data AlgReal = FromRat !Rational
| AlgReal !(UniPoly Rational) !Int !Rational !Rational
deriving (Show)
-- | 定義多項式
definingPolynomial :: AlgReal -> UniPoly Rational
definingPolynomial (FromRat x) = ind - constP x
definingPolynomial (AlgReal p _ _ _) = p
-- | 実根の分離区間
isolatingInterval :: AlgReal -> (Rational, Rational)
isolatingInterval (FromRat x) = (x - 1, x + 1)
isolatingInterval (AlgReal _ _ x y) = (x, y)
-- | 近似する区間の列
intervals :: AlgReal -> [Interval]
intervals (FromRat x) = repeat (Iv x x)
intervals (AlgReal f s a b) = intervalsWithSign f s a b
-- | 与えられた定義多項式と、分離区間 (a,b] から、代数的実数を構築する。
mkAlgReal :: UniPoly Rational -> Rational -> Rational -> AlgReal
mkAlgReal f a b
-- 0 の場合は FromRat 0 を使う
| a < 0 && b >= 0 && valueAt 0 f == 0 = FromRat 0
-- 区間が空の場合はエラー
| b <= a = error "mkAlgReal: empty range"
-- 区間の端点で多項式の値が 0 でないようにする
| valueAt b' f' == 0 = FromRat b'
-- 定数項の 0 を取り除き、また、区間の符号が確定したものをデータ構築子として使う
| otherwise = AlgReal f' s a' b'
where nonZeroPart xs | V.head xs == 0 = nonZeroPart (V.tail xs)
| otherwise = xs
f' = UniPoly $ nonZeroPart (coeff f)
s | signAt b f' > 0 = 1
| signAt b f' < 0 = -1
| otherwise = signAt b (diffP f')
Just (Iv a' b') = find (\(Iv x y) -> 0 < x || y < 0) (intervalsWithSign f' s a b)ついでに、 intervals 関数の型と実装を書き換えた。
代数的実数を計算可能実数とみなす
さっき書いた intervals :: AlgReal -> [Interval] 関数がまさに、代数的実数を計算可能実数としてみなすことに相当している。
algRealToCReal :: AlgReal -> CReal
algRealToCReal x = CReal (intervals x)代数的実数どうしの比較演算（compare 関数）は、等しくないことを確認した後は計算可能実数としての比較を利用できる：
instance Ord AlgReal where
compare (FromRat x) (FromRat y) = ...
compare (FromRat x) y = ...
compare x (FromRat y) = ...
-- 代数的実数同士の比較
compare x y
| b <= a' = LT -- 区間が重なっていない場合1（y の方が大きい）
| b' <= a = GT -- 区間が重なっていない場合2（x の方が大きい）
| countRealRootsBetween a'' b'' g == 1 = EQ -- 等しいかどうか？
-- x と y が等しくないことが確定した場合、計算可能実数として比較する
| otherwise = unsafeCompareCReal (algRealToCReal x) (algRealToCReal y)
where f = definingPolynomial x -- x の定義多項式
(a,b) = isolatingInterval x -- x の区間
f' = definingPolynomial y -- y の定義多項式
(a',b') = isolatingInterval y -- y の区間
g = gcdP f f'
a'' = max a a' -- x の区間と y の区間の共通部分の、下限
b'' = min b b' -- 同、上限代数的実数を浮動小数点数に変換する
TODO: 実装