高信頼・高精度・高可搬な数値計算法の研究

(1)

高信頼・高精度・高可搬な数値計算法の研究

Studies on Numerical Algorithms Retaining

High Reliability, High Accuracy and High Portability

2011 年 2 月

早稲田大学大学院基幹理工学研究科数学応用数理専攻数値解析研究

山中脩也

Naoya Yamanaka

(2)

(3)

3

序論

1.1 本論文の概要

本論文は高信頼・高精度・高可搬な数値計算法に関する研究成果をまとめたものである．本論文において高信頼な数値計算法とは，数値計算結果に対して数学的に厳密な誤差上限が与えられる計算法を指す．また，高精度な数値計算法とは，環境によらず精度のよい結果を返す計算法を指す．高可搬な数値計算法とは様々な環境で所望の値を得る計算法を指す．本論文ではこれらの三つの特徴を持つ数値計算法について述べ，特にそれらの特徴を持つ区間演算法について詳細に述べる．またこれらに関連する内容として，高精度な内積計算法の並列化法と，高信頼かつ高速な数値積分法について説明する．

本論文は，現在の電子計算機による数値計算の主流である，IEEE 754-1985 標準規格の浮動小数点数を基礎としている．IEEE 754-1985 標準規格は浮動小数点規格の一つであり，「電子計算機の中における数の表現」や「四則演算と平方根の演算結果精度」，「丸め」，「例外」などを定めている．数値計算では実数を直接扱えない代わりに浮動小数点数を扱っているため多くの計算において誤差が伴い，電子計算機の誕生当初から問題点として認識されてきた．その結果，計算法の安定性等について多くの知見が得られ，多くの優れた計算法は誤差の影響がなるべく小さくなるように設計されている．しかし実際の計算において丸め誤差の影響を完全に把握することは難しく，既存の計算法の大半は計算結果の精度について確実な保証を与えることはできない．

これらの現状を踏まえて，本論文の前半では，著者が提案した丸めの変更を用いない区間演算法について述べる．浮動小数点数を利用した計算機では全ての実数を表すことはできないため，入力された値や計算結果は丸められることになる．そのため計算機の中で数学的な計算結果を正しく保持するために，区間の概念が自然に導入される．区間は数の集合を表し，その区間内の全ての数を表す．数学的な計算結果より大きい(または小さい) 浮動小数点数の中で最も計算結果に近い値を得るのに，広く丸めの変更を用いた方法が

(8)

利用されてきた．丸めとは IEEE 754-1985 標準規格で規定される浮動小数点演算のパラメータの一つで，数学的な計算結果に対して，どの浮動小数点数を返すかを定めている．

IEEE 754-1985 標準規格では，最近点への丸め，上向き丸め，下向き丸め，ゼロへの丸め

の4つが定義されており，浮動小数点数を用いた区間の四則演算では，上向き丸めと下向き丸めを利用してそれぞれ計算を行なえば，数学的な結果を含む最も区間幅の小さな区間が得られることが知られていた．一方で，丸めの変更はそれなりに実行時間がかかるだけでなく，利用しているプラットフォームやコンパイラによって利用命令が大きく異なり，

計算環境によっては最近点への丸めしか無い物も多く存在する．これらの問題点に対し著者は，現在の数値計算環境で最も利用されている最近点への丸めのみを利用し，丸めを利用した区間演算と同じ精度を達成する，区間演算法を提案した．提案法は無誤差変換を用いて計算誤差を把握し，その値を判定することで，数学的な結果がどの区間に含まれるかを判定する．提案法を利用することで丸めの変更が不要となり，また，最近点への丸めに基づく既存の高精度計算ルーチンが直接利用可能となったことで，高精度な区間演算の構築が容易となった．

本論文の後半では，区間演算に関連する二つの事柄について述べる．

まず，著者が提案した高精度な内積計算法の並列化法について述べる．これは，荻田・

Rump・大石が提案した高精度内積計算法に対し，分散型並列計算機上においても，一台の計算機で計算した結果と同じ，または，より高精度な結果を得る並列化法である．荻田・Rump・大石による高精度内積計算法は，無誤差変換を用いた高精度な内積計算法で，

倍精度浮動小数点数のみを利用し，それより長い仮数部を持つ浮動小数点数を利用し計算を行なったときと同等の結果を返す特徴を持つ．一方でこの計算法は，誤差無しで計算を行なう性質上，データ依存性が非常に強く，分散型並列計算機における並列計算は難しいと考えられていた．この問題点に対し著者は，一台の計算機での内積計算結果を複数の倍精度浮動小数点数の形で保持することで，データの情報を出来る限り失うこと無く，また，各計算機間の通信情報量を大幅に増やすこと無く，一台の計算機で計算した結果と同じ，または，より高精度な結果を並列計算機で得ることに成功した．また申請者は，提案法の数学的解析によって誤差上限を示し，数値実験においてもその特徴が表れていることを示した．

次に著者が提案した，高信頼かつ高速な数値積分法について述べる．これは，高橋・森が提案した二重指数関数型積分公式を利用した，一変数関数に対する高信頼かつ高速な自動積分法である．一般に，数値積分の実行時間は，必要となる分点数の関数値計算時間が多くを占める．近似計算において利用者が指定した要求精度を満たす分点数を事前に把握することは容易ではなく，また，多くの数値積分法では得られた結果が要求精度を満足するか否かの判定を繰り返しているため，実際に必要な分点数より多くの関数値計算を行なうことが多く，実行時間に直接的な影響を与えていた．この問題点に対し著者が提案した

(9)

1.2 本論文の構成 9 計算法は，数値積分を行なう過程で発生する全ての誤差を事前に把握するよう試み，要求精度を満たす分点数を予め求める．これにより，二重指数関数型積分公式を用いた数値計算法において，高信頼で，従来の計算法に比べて同程度か，時にそれより高速な計算法を作成することに成功した．

本論文を通して，オーバフローとアンダーフローは起きないものと仮定する．また，計算法を表現するのに，MATLAB プログラムで用いる簡易的な表記を利用する．なお，浮動小数点加算，減算，乗算は浮動小数点演算回数 (ﬂop) を用いて数えられる．

1.2 本論文の構成

第1章では，本研究が行われた背景として高信頼・高精度・高可搬な数値計算法に関する事柄が記述され，本論文の概要と構成について述べる．

第2章では，高信頼・高精度・高可搬な区間演算法について説明する．まず，最近点への丸めを用いた区間演算の方針として，二つの手法があることを説明する．続いて必要となるいくつかの関数について具体的な計算手順が示し，それぞれの方針における四則演算の詳細な計算法について述べる．特に乗算に関しては，区間行列積の計算を考慮し，従来広く行なわれてきた区間行列積の計算法に比べ，より高精度な結果を返す計算法についても述べる．高精度区間行列積計算法の詳細な計算手順が示し，数値実験により従来の計算法に比べ，高精度な結果が得られていることを示す．

第3章では，高精度内積計算法の並列化法について議論する．まず荻田・Rump・大石によって提案された高精度な内積計算法について，計算手順とその誤差解析を述べる．続いて，それらの計算法に対応する並列化法について詳細な計算手順とその誤差解析が証明と共に述べる．最後に数値実験において，その特徴が表れていることや，提案された並列化法の並列化効率について述べる．

第4章では，高橋・森が提案した二重指数関数型積分公式を利用した，高信頼かつ高速な自動積分法について述べる．はじめに，二重指数関数型積分公式についての基礎的な事実を述べ，続いて，岡山・松尾・杉原による誤差解析定理を説明する．次に，この誤差解析を改良し，二重指数関数型積分公式のより適切な分点数と区間幅を計算する定理を示し，その有効性について説明する．また，高速化を達成するための手法として，二重指数関数型積分公式に対する丸め誤差の事前誤差計算法とその有効性を述べる．最後に提案手法を詳細に述べ，数値実験によって高信頼な計算結果が，従来の計算法に比べて同程度か，時にそれより高速に計算が行なえることを示す．

第5章では，本論文の結論として，まとめと今後の展望を述べる．

(10)

1.3 数値計算法の基礎

本論文は現在の電子計算機による数値計算の主流である，IEEE 754-1985標準規格の浮動小数点数を基礎としている. 本節では数値計算法の基礎として，IEEE 754-1985 標準規格について述べる．

1.3.1 浮動小数点数

IEEE 754-1985標準規格に基づく浮動小数点数システムは，浮動小数点数の集合とその

上の演算によって定義される．IEEE 754-1985 標準規格によって規定される浮動小数点数としては，正規化数（Normalized Numbers），零（Zeros），非正規化数（Denormalized Numbers），浮動小数点例外（無限（Inﬁnities），NaNs）が用意されている．

正規化数

0 か1 である di に対し a=±

( 1 2⁰ + d1

2¹ + d2

2² + d3

2³ +· · ·+ dN−1

2^N⁻¹ )

×2^e, (1.1)

と書ける数を正規化数という．e_min を負の整数，e_max を正の整数として，e は e_min 5 e 5emax となる整数である．

m=± ( 1

2⁰ + d₁ 2¹ + d₂

2² + d₃

2³ +· · ·+ d_N₋₁ 2^N⁻¹

)

(1.2) を符号付き仮数（signed mantissa）といい，e を指数（exponent）という．指数 e も二進数で表される．通常，単精度（4 byte = 32 bit），倍精度（8 byte = 64 bit），拡張倍精度

（16 byte = 128 bit）の浮動小数点システムがあるが，それらは次のような浮動小数点シ

ステムである．

単精度 N = 24, −1265e5127 倍精度 N = 53, −10225e51023 拡張倍精度 N = 64, −163825e516383

(1.3)

本論文では，特に断ることのない限り，浮動小数点数といえば，IEEE 754-1985標準規格に基づく倍精度浮動小数点のことを指すこととする．

正規化数の絶対値の最大値は x_max =

( 1 2⁰ + 1

2¹ + 1 2² + 1

2³ +· · ·+ 1 2^N⁻¹

)

×2^e^max (1.4) であり，その最小値は

xmin = 2^e^min (1.5)

(11)

1.3 数値計算法の基礎 11 である．倍精度ではそれぞれ

( 1 2⁰ + 1

2¹ + 1 2² + 1

2³ +· · ·+ 1 2⁵²

)

×2¹⁰²³= 1.79769· · · ×10³⁰⁸ (1.6) 2⁻¹⁰²²= 2.22507· · · ×10⁻³⁰⁸ (1.7) である．|x|> xmax のときにオーバ―フロー（overﬂow）が生じたという．

倍精度浮動小数点数においては，仮数部が53ビットであるから

2⁻⁵³ = 1.110223· · · ×10⁻¹⁶ (1.8) より倍精度浮動小数点数は十進法で約16桁の精度がある．

また，本論文において倍精度浮動小数点数の集合を F ^とする．

零

零は規格化され +

( 0 2⁰ + 0

2¹ + 0 2² + 0

2³ +· · ·+ 0 2^N⁻¹

)

×2^e^min (1.9)

と表される．

非正規化数

IEEE 754-1985 標準規格では，浮動小数点数は指数部が emin となったとき，仮数

部の最初の桁が 1 より小さい数を表すために，デフォルトで最初の桁を 1とすることをやめ，ここが 0 となる数を置くことを許す規格となっている．これを，非正規化数

（denormalized number）という．非正規化数の範囲に数が入ることを漸化アンダーフロー

（gradual underﬂow）という．

非正規化数の正で最も小さな数は

2⁻¹⁰⁷⁴= 4.940656· · · ×10⁻³²⁴ (1.10) これ以下の数になると，アンダーフロー（underﬂow）が生じたという．

本論文における非正規化数の集合を U とする．

浮動小数点例外

この他に，IEEE 754-1985 標準規格では次のような特別な数が用意されている．

• Inf (Inﬁnity) はオーバーフローの結果や零で割った結果として得られる．

• NaN（Not a Number）は √

−5,∞

∞, +∞, + (−∞) などの不当な演算の結果として得られる．

• ±0 はアンダーフローから ±∞ での割り算の結果として得られる．

(12)

1.3.2 丸めと四則演算

丸め

IEEE 754-1985 標準規格では，浮動小数点数の演算は丸めを用いて定義されている．c

を実数とする．

• ^{上向きの丸め（}round upward） 4:R→F c 以上の浮動小数点数の中で最も小さい数に丸める．

• 下向きの丸め（round downward） 5:R→F c 以下の浮動小数点数の中で最も大きい数に丸める．

• 最近点への丸め（round to nearest） or fl(· · ·) :R→F

cに最も近い浮動小数点数に丸める．もしも，このような点が2点ある場合は，仮数部の最終ビットが偶数である浮動小数点に丸める（これを偶数への丸めという）．

• 切り捨て（round toward 0）

絶対値が c以下の浮動小数点数の中で最も c に近いものに丸める．

四則演算

IEEE 754-1985標準規格において浮動小数点数の四則演算の結果は，数学的に正しい実

数としての四則演算の結果を指定された丸めを行なって得られた数に一致するよう定められている．ここで，数学的な四則演算を ◦ ∈ {+,−,×, /} とし，最近点への丸めを用いた四則演算を ∈ {⊕, ,⊗,} ^{とすれば任意の} x, y∈F ^{について，}

xy=(x◦y) (1.11)

が成立する．また，平方根も同様の性質を満たす．すなわち，数値計算において，演算結果の精度が保証されている演算規則は四則演算と平方根のみである．注意すべきことは，

指数関数や三角関数などはこのような規格をみたすようには規格されていないことである．これは，これら初等関数の値を精度保証付きで求めるためには別途工夫が必要であることを意味する．

また，十進数を二進数に変換すること，また，その逆の変換をすることは，しばしば数値計算で必要となる．十進数整数は二進整数へ誤差無しで変換できること，および十進小数を二進小数へ変換すると誤差が生じることに注意する．

(13)

1.3 数値計算法の基礎 13

1.3.3 浮動小数点数とその演算の特徴

eps・eta

epsは相対丸め誤差定数のことであり，本論文では1.0 から次に小さな浮動小数点数までの距離を示す^*1．また，etaはアンダーフロー定数であり，正の非正規化数の中で最も小さな数のことである．IEEE 754-1985標準規格浮動小数点数の場合，

eps= 2⁻⁵³ (1.12)

eta= 2⁻¹⁰⁷⁴ (1.13)

である．よって，1

2eps⁻¹eta は正の規格化浮動小数点数の中で最も小さな浮動小数点数であるので，f ∈F に対して，

f ∈U ⇐⇒ 05|f|5 1

2eps⁻¹eta (1.14)

が成立する．

a, b∈F に対して浮動小数点加減算は，

|ε|5eps (1.15)

を満たす ε を用いて

fl(a±b) = (a±b) (1 +ε) (1.16)

と書ける．アンダーフロー近くの加減算には

|a±b|5eps⁻¹eta =⇒ fl(a±b) =a±b (1.17)

fl(a±b) = 0 ⇐⇒ a =−b (1.18)

が成立するため，式(1.16) には etaが必要ないことに注意する．

ufp・ulp

ufp とは ”unit in the ﬁrst place” のことであり，

06=r ∈R =⇒ ufp(r) := 2^b^log²^|^r^|c (1.19) で定義される．なお，ufp(0) := 0 とする．この記号を用いると，正規化数 f のビットがどの領域に入っているのかを簡単に表すことができるようになり，f のビットの幅は，先

*1epsの定義は幾つか存在することに注意する．広く知られた異なる定義として，1.0から次に大きな浮動小数点数までの距離がある．[1]

(14)

頭ビットが ufp(f) であり，最後のビットは 2eps·ufp(f) と表せる．なお，最後のビットを ”unit in the last place (ulp)”と呼ぶ．すなわち，

ulp(f) = 2eps·ufp(f) (1.20) が成立する．本論文における解析では，ある浮動小数点数がどの位置にあるかを見ることがあり，次式が有用である；

05r∈R =⇒ ufp(r)5|r|<2ufp(r) (1.21) r, r⁰ ∈R, ufp(r)5|r⁰| =⇒ ufp(r)5ufp(r⁰). (1.22)

pred・succ

浮動小数点数 f より小さい浮動小数点の中で最も大きな浮動小数点数と，浮動小数点数 f より大きい浮動小数点の中で最も小さな浮動小数点数を，次のように定義する：

pred(r) := max{f ∈F : f < r} (1.23) succ(r) := min{f ∈F : r < f}. (1.24) ufp を用いると，pred とsucc は次のように書ける：

f /∈U, |f| 6=ufp(f) =⇒ pred(f) =f−2eps·ufp(f), (1.25) succ(f) =f + 2eps·ufp(f). (1.26) f /∈U, |f|=ufp(f) =⇒ pred(f) = (1−eps)f, (1.27) succ(f) = (1 + 2eps)f. (1.28) f /∈U, |f|=− 6=ufp(f) =⇒ pred(f) = (1 + 2eps)f, (1.29) succ(f) = (1−eps)f. (1.30) f ∈U =⇒ pred(f) =f−eta, (1.31)

succ(f) =f +eta. (1.32)

浮動小数点数 f に対して

pred(f)5f−eps·ufp(f)5f +eps·ufp(f)5succ(f) (1.33) が成立し，f /∈U に対して

f −2eps·ufp(f)5pred(f)5succ(f)5f+ 2eps·ufp(f) (1.34) が成立する．なおここで，非正規化数の集合 U ^{は最も小さな正規化数} ±1

2eps⁻¹eta を含むことに注意する．

次に succ(x) とpred(x) を計算する手法をそれぞれ Algorithm 1，Algorithm 2 に示す：

なおここでabssucc と abspred は非負の入力に対して，Algorithm 3 ，Algorithm 4 で定義される：

(15)

1.3 数値計算法の基礎 15 Algorithm 1 次に大きい浮動小数点数を求める計算法

x∈F に対して，次に大きい浮動小数点数 y∈F を求める計算法．

function y =succ(x) if x=0

y=abssucc(x);

else

y=−abspred(−x);

end end

Algorithm 2 次に小さい浮動小数点数を求める計算法

x∈F に対して，次に小さい浮動小数点数 y∈F を求める計算法．

function y =pred(x) if x=0

y=abspred(x);

else

y=−abssucc(−x);

end end

Algorithm 3 絶対値が次に大きい浮動小数点数を求める計算法

x∈F に対して，絶対値が次に大きい浮動小数点数 y∈F を求める計算法．

function y =abssucc(x) if x=2⁻¹⁰²²

y= x 1−2⁻⁵³; else

y=x+ 2⁻¹⁰⁷⁴; end

end

(16)

Algorithm 4 絶対値が次に小さい浮動小数点数を求める計算法

x ∈F に対して，絶対値が次に小さい浮動小数点数y ∈F を求める計算法．

function y =abspred(x) if x=2⁻¹⁰²²

y =(

1−2⁻⁵³) x;

else

y =x−2⁻¹⁰⁷⁴; end

end

1.4 高信頼数値計算法の基礎

本節では，高信頼な数値計算法の基礎について述べる．高信頼な数値計算法とは，数値計算結果に対して数学的に厳密な誤差上限が与えられる計算手法を指し，広く「精度保証付き数値計算」とも呼ばれる．本節では特に，高信頼な数値計算法において重要となる区間解析について述べる．

浮動小数点数システム上での高信頼な数値計算の原理をまず簡単に説明する．ここで問題は連続数学の問題で，解が実数で与えられるものを考える．浮動小数点数を用いた数値計算では，真の解そのものを計算することは一般には不可能である．そこで，真の値の下限を与える浮動小数点数と上限を与える浮動小数点数を計算する．もし，下限と上限の数値が小数点以下15 桁まで一致すれば，真の解の小数点以下15 桁まで計算できたことになる．すなわち，精度保証付き数値計算の基本的なアイデアは，浮動小数点数を両端とする区間の中に連続数学の問題の解を包み込もうというものである．

精度保証付き数値計算の原理は，実数値で与えられる真の値の上限と下限を浮動小数点演算により計算しようとするものである．このような考え方は，区間解析（interval

analysis）と呼ばれる数学的理論と密接な関係がある．区間解析のアイディアは日本の須

永教授が1950年代に提案したものである．その後，1960年代に米国のムーアが学位論文のテーマとし，専門書を出版するに至って，多くの研究者の研究対象となった．

区間解析では，区間を数の拡張と考える．そして，その四則演算が定義される．これを区間演算という．区間演算を浮動小数点システムの上で展開するものを機械区間演算という．

(17)

1.4 高信頼数値計算法の基礎 17

1.4.1 上限下限型の区間演算

上限下限型の区間解析

まず，実数上での演算を仮定する区間演算の概念を説明する．本節においては四則演算は実数の四則演算とする．したがって，有理数演算を実装して，扱うすべての有理数の場合を除けば，本節の議論は数学的仮想の上のもので，電子計算機上では実現できないことに注意する．

区間解析においては，区間とは x∈R に対して，

[x, x] ={x|x 5x 5x} (1.35) と表される閉区間であるとする．ただし，x 5x∈R で，それぞれの区間の下端，上端と呼ぶことにする．以下，記号の節約のため

[x] = [x, x] (1.36)

と区間を表すことにし，区間 [x] と単に書けば，それは [x] = [x, x] を表すものとする．

x=x となる区間[x] は点区間（point interval）という．点区間は実数であるので，点区間 [x] の表す実数を x と書くことにする．

区間 [x] について以下を定義する：

d([x]) =x−x, (1.37)

r([x]) = x−x

2 , (1.38)

m([x]) = x+x

2 . (1.39)

d([x]), r([x]), m([x]) をそれぞれ区間 [x]の直径，半径，中心という．

区間 [x] の最小絶対値 mig と最大絶対値 magをそれぞれ次で定義する．

mig([x]) =h[x]i= min{|x| |x∈[x]} (1.40) mag([x]) =|[x]|= max{|x| |x∈[x]}= max{−x, x} (1.41) 二つの区間 [x],[y] が与えられたとき，その二つの区間の四則演算を次で定義する．

[x]◦[y] ={x◦y|x ∈[x], y∈[y]} (1.42) ただし，◦ ∈ {+,−,×, /} である．これを区間演算という．区間演算に対して次の Algorithm 5 が成立する．

Algorithm 5から区間演算が有限回の実数演算で実行できることがわかる．乗算と除算

は，場合分けにより少ない手間で計算する規則も簡単に導ける．Algorithm 6，Algorithm 7 に示す．

(18)

Algorithm 5 上端下端型区間同士の実四則演算上端下端型区間 [x] = [x, x], [y] =[

y, y]

の四則演算計算法．

[x] + [y] = [x+y, x+y] (1.43)

[x]−[y] = [x−y, x−y] (1.44)

[x]×[y] = [min{xy, xy, xy, xy}, max{xy, xy, xy, xy}] (1.45) [x]

[y] = [x]× [1

y, 1 y ]

(0∈/ [y]) (1.46)

Algorithm 6 上端下端型実区間乗算上端下端型区間 [x] = [x, x], [y] =[

y, y]

の実区間乗算計算法．

[x]<0 [x]30 [x]>0 [y]<0 [xy, xy] [xy, xy] [xy, xy]

[y]30 [xy, xy] [min(xy, xy), max(xy, xy)] [xy, xy]

[y]>0 [xy, xy] [xy, xy] [xy, xy]

Algorithm 7 上端下端型実区間除算上端下端型区間 [x] = [x, x], [y] =[

y, y]

の実区間除算計算法．

[x]<0 [x]30 [x]>0 [y]<0 [

x/y, x/y]

[x/y, x/y] [

x/y, x/y] [y]>0 [

x/y, x/y] [

x/y, x/y]

区間演算については，包含関係における単調性

[x]⊆[x⁰],[y]⊆[y⁰] =⇒[x]◦[y]⊆[x⁰]◦[y⁰], ◦ ∈ {+,−,×, /} (1.47) が成立する．

また，加算と乗算に関し，交換則と結合則が成立する．

[x]◦[y] = [y]◦[x], ◦ ∈ {+,×} (1.48) [x]◦([y]◦[z]) = ([x]◦[y])◦[z]. ◦ ∈ {+,×} (1.49) しかし，加算と乗算の逆元は存在しない．すなわち −[x] = [−x, −x] であるが，

0 = [0]⊆[x]−[x] = [x−x, x−x] (1.50)

1 = [1]⊆[x]/[x] (1.51)

(19)

1.4 高信頼数値計算法の基礎 19 となることもわかる．上式で等号は [x] が点区間のときのみに成立する．

また，分配則も区間演算に対しては成立しない．そのかわり次の劣分配則が成立する．

[x]×([y] + [z])⊆[x]×[y] + [x]×[z] (1.52) この式で等号は，例えば区間 [y] と [z] が同じ符号を持つ時に成立する．

上限下限型の機械区間演算

前節で説明した区間演算においては，区間の両端の数は実数であるとし，厳密な実数演算に基づいて区間演算が定義された．しかし，浮動小数点数システム上で数値計算するときには，厳密な区間演算は実行できない．そこで，数学的な正解が含まれる区間が演算の結果得られることを念頭に，区間演算を与えられた浮動小数点数システム F ^上に展開する方法を示す．IF ={[x]∈ IR |x, x∈F} とする．IF の要素を機械区間という．

機械区間は区間の定義で上端と下端の数が浮動小数点数で与えられるものをいう．区間 [a, b]∈IR ^を区間 [c, d]∈IF ^{に丸める丸め} ♦:IR →IF ^を

[x]⊂♦[x] (1.53)

および，すべての [x],[y]∈IR について次を満たすように定義する．

♦[x] = [x] (1.54)

[x]⊂[y] =⇒♦[x]⊂♦[y] (1.55)

♦(−[x]) =−♦[x] (1.56)

IEEE 754-1985 標準規格においては

♦[x] = [5(x), 4(x)] (1.57)

とすると，上記条件は満たされる．上向き丸めと下向き丸めを利用した，上端下端型の機械区間演算の四則演算は Algorithm 8–Algorithm 11 のように書ける．

(20)

Algorithm 8 上限下限型機械区間加算上限下限型区間 x = [x, x] と y = [

y, y]

の丸めを用いた機械区間加算計算法．このとき，上端下端型区間 z は次式を満たす：

z = [z, z]⊃[x, x] +[ y, y]

. (1.58)

function z =infsup-roundupdown-plus(x, y) z =5(

x+y)

; z =4(x+y) ; end

Algorithm 9 上限下限型機械区間減算上限下限型区間 x = [x, x] と y = [

y, y]

の丸めを用いた機械区間減算計算法．このとき，上端下端型区間 z は次式を満たす：

z = [z, z]⊃[x, x]−[ y, y]

. (1.59)

function z =infsup-roundupdown-minus(x, y) z =5(x−y) ;

z =4( x−y)

; end

(21)

Algorithm 10 上限下限型機械区間乗算上限下限型区間 x= [x, x] とy=[

y, y]

の丸めを用いた区間乗算計算法．このとき，上端下端型区間 z は次式を満たす：

z = [z, z]⊃[x, x]×[ y, y]

. (1.60)

function z =infsup-roundupdown-times(x, y) if x >0

if y >0 z =5(

xy)

; z =4(xy) ; elseif y <0

z =5( xy)

; z =4(xy) ; else

z =5( xy)

; z =4(xy) ; end

elseif x <0 if y >0

z =5(xy) ; z =4( xy)

; elseif y <0

z =5(xy) ; z =4( xy)

; else

z =5(xy) ; z =4( xy)

; end

else

if y >0 z =5(

xy)

; z =4( xy)

; elseif y <0

z =5(xy) ; z =4(xy) ; else

z = min(

5(xy),5( xy))

; z = max( 4(

xy)

,4(xy))

; end

end end

(22)

Algorithm 11 上限下限型機械区間除算上限下限型区間 x = [x, x] とy =[

y, y]

の丸めを用いた区間除算計算法．このとき，上端下端型区間 z は次式を満たす：

z = [z, z]⊃[x, x]/[ y, y]

. (1.61)

function z =infsup-roundupdown-divide(x, y) if x >0

if y >0 z =5(

x/y)

; z =4(x/y) ; elseif y <0

z =5( x/y)

; z =4(x/y) ; else

Error;

end elseif x <0 if y >0

z =5(x/y) ; z =4( x/y)

; elseif y <0

z =5(x/y) ; z =4( x/y)

; else

Error;

end else

if y >0 z =5(

x/y)

; z =4( x/y)

; elseif y <0

z =5(x/y) ; z =4(x/y) ; else

Error;

end end end

(23)

1.4.2 中心半径型の区間演算

中心半径型の区間解析

上端下端型区間演算における掛け算は，場合分けも多く複雑である．これをシンプルに計算するために区間の中心と半径による表示を考える．中心半径型区間演算の四則演算としては，以下の計算則（Algorithm 12）が知られている．

Algorithm 12 中心半径型区間同士の実四則演算

中心半径型区間 < α >=< α_c, α_r >, < β >=< β_c, β_r > の四則演算計算法．

< α >+< β >=< αc +βc, αr+βr> (1.62)

< α >−< β >=< αc −βc, αr+βr> (1.63)

< α >×< β >⊂< α_cβ_c, |α_c|β_r+|β_c|α_r+α_rβ_r > (1.64)

< α >

< β > =< α >× 1

< β > = 1

β_c²−β_r² (< α >×< β >) (0∈/< β >) (1.65)

中心半径型の機械区間演算

中心半径型区間の四則演算は，上向き丸めと下向き丸めを利用して Algorithm 13–

Algorithm 16 のように書ける．

Algorithm 13 中心半径型機械区間加算

中心半径型区間 x =< x_c, x_r > と y =< y_c, y_r > の丸めを用いた区間加算計算法．このとき，中心半径型区間 z は次式を満たす：

z =< z_c, z_r >⊃< x_c, x_r >+< y_c, y_r > . (1.66) function z =midrad-roundupdown-plus(x, y)

z_c =5(x_c+y_c) e_c =4(x_c+y_c−z_c) zr =4(xr+yr+ec) end

(24)

Algorithm 14 中心半径型機械区間減算

中心半径型区間 x =< xc, xr > と y =< yc, yr > の丸めを用いた区間減算計算法．このとき，中心半径型区間 z は次式を満たす：

z =< zc, zr >⊃< xc, xr>−< yc, yr > . (1.67) function z =midrad-roundupdown-minus(x, y)

zc =5(xc−yc) e_c =4(x_c−y_c−z_c) zr =4(xr+yr+ec) end

Algorithm 15 中心半径型機械区間乗算

中心半径型区間 x =< xc, xr > と y =< yc, yr > の丸めを用いた区間乗算計算法．このとき，中心半径型区間 z は次式を満たす：

z =< zc, zr >⊃< xc, xr>×< yc, yr > . (1.68) function z =midrad-roundupdown-times(x, y)

zc =5(xcyc) ec =4(xcyc−zc)

z_r =4(|x_c|y_r+x_r|y_c|+x_ry_r+e_c) end

Algorithm 16 中心半径型機械区間除算

中心半径型区間 x =< xc, xr > と y =< yc, yr > の丸めを用いた区間除算計算法．このとき，中心半径型区間 z は次式を満たす：

z =< z_c, z_r>⊃< x_c, x_r > / < y_c, y_r > . (1.69) function z =midrad-roundupdown-divide(x, y)

zc =5

( xcyc

(x_c −x_r) (x_c+x_r) ) t =4

( 1

(x_c−x_r) (x_c+x_r) )

e_c =4(tx_cy_c−z_c)

zr =4(t(|xc|yr+xr|yc|+xryr) +ec) end

(25)

1.4.3 丸め誤差の事前誤差解析

区間演算を用いることで関数値 f の計算で発生する丸め誤差を事前に見積もることができる．与えられた一変数関数 f(x) に対して，a 5x5bを満たす全ての x に対して，

max

a5x5b|fl(f(x))−f(x)|5ε, (1.70) を満たす定数 ε を計算する．この手法は “オペレータオーバーロード” という概念に基づいている．関数 f(x) は再帰的に

(1) 四則演算 {+,−,×, /},

(2) sin,cos,exp,log などの初等関数, (3) それらの合成関数，

で書き表されていると仮定する．このとき，データ構造として (I, ε) を考える．I が区間，ε は正の定数である．このデータ構造に対し，四則演算はAlgorithm 17のように定義される：

Algorithm 17 (I, ε) に対する四則演算データ構造 (I, ε) に対する四則演算計算法．

(1) (加算)

I =I1+I2 (1.71)

(I1, ε1) + (I2, ε2) = (I, 4(ε1+ε2+eps·mag(I))), (1.72) (2) (減算)

I =I₁−I₂ (1.73)

(I1, ε1)−(I2, ε2) = (I, 4(ε1+ε2+eps·mag(I))), (1.74) (3) (乗算)

I =I1×I2 (1.75)

(I₁, ε₁)·(I₂, ε₂) = (I, 4(ε₁mag(I₂) +ε₂mag(I₁) +eps·mag(I))), (1.76) (4) (除算)

I⁰ =I₁/I₂, I⁰⁰ = 1/I₂, I⁰⁰⁰ =I⁰×I⁰⁰ (1.77) (I1, ε1)/(I2, ε2) = (I⁰, 4(ε1mag(I⁰⁰) +ε2mag(I⁰⁰⁰) +eps·mag(I⁰))). (1.78)

(26)

f(x) を連続的に微分可能な関数とする．ここで，関数 f(x) を次の精度で計算できるようなライブラリを持っていると仮定する^*2：

|fl(f(x))−f(x)|5 1

2ulp((f(x))). (1.79)

I_f を，f(x) の x に区間 I を入力し区間四則演算を利用して計算した f の計算区間とする．また同様に，I_f⁰ を，f⁰(x) の x に区間 I を入力し区間四則演算を利用して計算した f⁰ の計算区間とする．このとき f((I, ε)) は次のように定義される．

I_f =f(I), I_f⁰ =f⁰(I) (1.80) f((I, ε)) = (I_f, 4(

ε·mag( I_f⁰)

+eps·mag(I_f))

). (1.81)

ここで，C++, MATLAB などのオペレータオーバーロードがサポートされている言

語を使用していると仮定する．このとき，与えられた関数 f に対して，(J, ε) を使って Algorithm 18 のようにボトムアップで計算行なっていく事で，

a5x5bmax |fl(f(x))−f(x)|5ε (1.82) を満たす定数 ε が計算される．

Algorithm 18 関数の丸め誤差の事前計算法

f(ξ) の計算を行なったときの丸め誤差の事前計算法 (a5ξ 5b, ξ ∈F)．手順 1 区間 I = [a, b] とする.

手順 2 x= (I, eps·max (|a|,|b|)) とする.

手順 3 (Iy, εy) =f(x) をオペレータオーバーロードを使って計算する．

手順 4 ε_y を出力する.

*2初等関数に対して式 (1.79)の特徴を満たすライブラリとしてCRlibmライブラリ[2] がある．

(27)

1.5 高精度数値計算法の基礎 27

1.5 高精度数値計算法の基礎

1.5.1 演算精度と結果精度

精度という言葉は，数値計算の中でも誤解を招きやすいものなので，ここで定義を明確にしたい．数値計算における精度には，下記の2つの意味がある．

• 演算精度（precision）

• ^{結果精度（}accuracy）

演算精度は，1回の浮動小数点演算の相対精度のことであり，すなわち，定数である．一方，結果精度は，真値に対する数値計算結果の精度であるため問題に依存し，演算精度からだけではわからない．数値計算の分野で高精度計算というと，多くの場合，演算精度が高い，という意味で使われる．したがって，高精度計算と言っても，必ずしも結果精度が高いとは限らない．

1.5.2 準備

p= (p₁, . . . , p_n)^T ∈Fⁿ とする．このとき，次式が成立する [3]：

˜ s :=fl

( _n

∑

i=1

pi

)

=⇒

s˜−

∑n i=1

pi

5γn−1

∑n i=1

|pi|. (1.83)

このとき，γn は

γ_n:= n·eps

1−n·eps (1.84)

であり，n·eps<1 を暗に仮定している．

s と S をそれぞれ次のように定義する：

s:=

∑n i=1

pi, (1.85)

S :=

∑n i=1

|p_i|. (1.86)

また，ベクトル p の総和に関する条件数は次のように定義される：

cond ( _n

∑

i=1

pi

) := S

|s|, s6= 0. (1.87)

(28)

同様に，ベクトルx= (x₁, . . . , x_n)^T ∈Fⁿ^，y = (y₁, . . . , y_n)^T ∈Fⁿ ^{に対して内積} x^Ty の条件数は次のように定義される：

cond( x^Ty)

:= 2|x^T||y|

|x^Ty| , x^Ty6= 0. (1.88)

1.5.3 加乗算に関する無誤差変換

1969 年，Knuth は2つの浮動小数点数 a, b の加算は，誤差なく計算できることを示し

た（Algorithm 19）[4]．アルゴリズム19 は，下記の事実を示している：

Algorithm 19 加算の無誤差変換加算に関する無誤差変換．このとき

a+b=x+y (1.89)

|y|5 1

2ulp(x) (1.90)

が成立する．

function [a, b] =TwoSum(x, y) a =x⊕y

c=a x

b= (x (a c))⊕(y c) end

• x は a+b の最も良い近似である．

• x の誤差 y =a+b−x は浮動小数点数である．

• y は上記の手法を用いて浮動小数点演算のみで計算できる．

1971 年に，乗算に対しても同様な性質が成り立つことを Dekker [5]が示した．これを説明するためにまず，浮動小数点数を無誤差で二つの浮動小数点数に分割する計算法を Algorithm 20 に示す．ここで，x_h と x_t は，それぞれx の仮数部上位半分の 26 ビットと下位半分の 26 ビットに対応するものである．ただし，必ずしもビットパターンが一致するわけではない．この式は，x を xh と xt の和に誤差なしで変換できることを意味している．

Algorithm 20 を用いて Dekker は次のことを示した．まず，手法 TwoProduct を Algorithm 21 のように定義する．

Algorithm 19やAlgorithm 21によって，2つの浮動小数点数a, bの加算や乗算は，いずれも，その近似xとその誤差y よいう2つの浮動小数点数の和x+y の形に誤差なく変

(29)

1.5 高精度数値計算法の基礎 29 Algorithm 20 無誤差の分割

浮動小数点数 x を無誤差で二つに分割する計算法．このとき

x=xh+xt (1.91)

|x_h|=|x_t| (1.92)

が成立する．

function [xh, xt] =Split(x) c= (2²⁷+ 1)⊗x x_h =c (c x) xt =x xh

end

Algorithm 21 乗算の無誤差変換乗算に関する無誤差変換．このとき

a×b=x+y (1.93)

|y|5 1

2ulp(x) (1.94)

が成立する．

function [a, b] =TwoProduct(x, y) a =x⊗y

[x₁, x₂] =Split(x) [y1, y2] =Split(y)

b=x2⊗y2 (((a x1⊗y1) x2⊗y1) x1⊗y2) end

換できることを分かった．このような変換を，無誤差変換（Error-Free Transformations）と呼ぶ．

1.5.4 ベクトルの総和の高精度な計算法

計算誤差を捨てない計算法 VecSum

TwoSumは，ベクトルに対する誤差無しの加算に応用できる（Algorithm 22）：

(30)

Algorithm 22 ベクトルの総和に関する無誤差変換ベクトルの総和に関する誤差無し変換．

function q =VecSum(p) q1 =p1

for i= 2 :n

[qi, qi−1] =TwoSum(pi, qi−1) end

end

図 1.1は n= 5 の時のVecSum のイメージ図である．

qn は明らかに

q_n=

∑n i=1

p_i, i.e. q_n =fl ( _n

∑

i=1

p_i )

(1.95)

である．そして，q₁, q₂, . . . , q_n₋₁ は fl(∑n

i=1p_i) の残差である．VecSum はベクトル p= (p1, p2, . . . , pn)^T ∈Fⁿ (1.96) を新たなベクトル

q = (q1, q2, . . . , qn)^T ∈Fⁿ (1.97) へ次式を満たすように変換する．すなわち

∑n i=1

pi =

∑n i=1

qi (1.98)

が成立する．

図1.1 n= 5 の時のVecSum のイメージ図

(31)

1.5 高精度数値計算法の基礎 31 式(1.83)によって，VecSum は次式を満たす [6, (4.4)]：

n−1

∑

i=1

|qi|5γn−1

∑n i=1

|pi|. (1.99)

擬似的に K 倍精度で総和を計算する SumK

荻田・Rump・大石は文献[6] の中で，VecSumを再帰的に用いて総和を計算する SumK を提案した（Algorithm 23）．

Algorithm 23 K 倍精度の総和計算

(K −1) 回のベクトルに関する誤差無し変換を用いた，ベクトル p ∈Fⁿ に対する K 倍精度の総和計算．

function res=SumK(p, K) for i= 1 :K−1

p=VecSum(p) end

res =fl ( _n

∑

i=1

pi

)

end

図1.2 はn= 5，K = 4 の時の SumK のイメージ図を表している．

図1.2 n= 5，K= 4 の時のSumK のイメージ図

(32)

このときSumK による結果 res の誤差上限は以下のように表される [6]：定理 1 (荻田ら[6])

res をSumK で得られた結果とするとき，res は次式を満たす．

|res−s|5(eps+ 3γ_n²₋₁)|s|+γ_2(n^K ₋₁₎S. (1.100) さらに，もし s6= 0 ならば，次式を満たす．

res−s s

5eps+ 3γ_n²₋₁+γ_2(n−1)^K ·cond ( _n

∑

i=1

pi

)

. (1.101)

この定理は次式

res−s s

.eps+O(eps^K)·cond ( _n

∑

i=1

p_i )

, (1.102)

が成立する事を述べており，SumKによって内部で K 倍精度演算を行ない（この事実は，

右辺の二項目と一致する），その値を倍精度演算に丸めたかのような結果 res を返すこと

（この事実は，右辺の一項目と一致している）を意味している．すなわち res は，もし cond

( _n

∑

i=1

pi

)

< γ_2(n⁻^K₋₁₎eps≈(2neps)⁻^Keps (1.103) ならば，倍精度においてほとんど最大の精度を達成できていることになる．

擬似的な L 倍精度で総和を計算し，L 個で和を表わす計算法 SumL

次に，L個の浮動小数点数の和によって結果を表す総和計算法 SumL（Algorithm 24）を示す．

Algorithm 24 L 倍精度総和計算

L 個の浮動小数点数の和によって表された L 倍精度総和計算．

function q =SumL(p, L) for k = 1 :L−1

p(1 :n−k+ 1) =VecSum(p(1 :n−k+ 1)) qk =pn−k+1

end q_L=fl

(_n₋_L+1

∑

i=1

p_i )

end

図 1.3は n= 5，L= 3 のときの SumL のイメージ図を表している．

SumL に関して次の定理を示す：

(33)

1.5 高精度数値計算法の基礎 33

図1.3 n= 5，L= 3 のときのSumLのイメージ図

定理 2

SumL は次の不等式を満たす．

∑L k=1

q_k−

∑n i=1

p_i

5γ_n−1^L

∑n i=1

|p_i| (1.104)

∑L k=1

|qk|5(1 +γ_2(n₋₁₎)

∑n i=1

|pi|. (1.105)

式(1.104) から SumL の結果ベクトル q は s=∑n

i=1p_i とするとき，

∑L

k=1q_k−s s

.O(eps^L)·cond ( _n

∑

i=1

pi

)

, (1.106)

を満たすことがわかる．これは，s が総和 ∑L

k=1qk をL 倍精度を使って計算した近似値になっていることを示している．

1.5.5 ベクトルの内積の高精度な計算法

計算誤差を捨てない内積計算法 VecProduct

続いて，内積計算について述べる．x, y ∈Fⁿ とする．内積 x^Ty を x^Ty =

∑2n i=1

t_i (1.107)

(34)

Algorithm 25 ベクトルの内積に関する無誤差変換二つのベクトル x, y ∈ Fⁿ ^をx^Ty =∑2n

i=1t_i を満たす新たなベクトル t ∈F²ⁿ ^に変換する計算手法．

function t =VecProduct(x, y) r= 0

for i= 1 :n

[h_i, t_i] =TwoProduct(x_i, y_i) [r, tn+i−1] =TwoSum(r, hi) end

t_2n =r end

となるように変換した次の計算手法 VecProduct （Algorithm 25）を示す．

図 1.4は n= 4 の時のVecProduct のイメージ図を表している．

VecProduct は次の定理を満たす．

定理 3

VecProduct は次の不等式を満たす：

2n∑−1 i=1

|t_i|5γ_n

∑n i=1

|x_iy_i| (1.108)

|t_2n|5(1 +γ_n)

∑n i=1

|x_iy_i|. (1.109)

!

" ! " " "

# ! # # #

$&% $(' $() $+*

図1.4 n= 4の時の VecProductのイメージ図

高信頼・高精度・高可搬な数値計算法の研究