● 講義資料
▼ 講義予定
•
円周率の近似値の計算と浮動小数点演算の誤差（続き）•
テイラー級数による円周率・自然対数の底などの近似値の計算● 前回の講義のまとめ
▼ 浮動小数点演算
•
ここでは, 考える浮動小数点数の体系をF
とおく. （特に,β
進p
桁の浮動小数点を考えるときには
F = F (β, p)
と書くこともある.）倍精度浮動小数点数の体系はβ = 2, p = 53
である.
•
実数x ∈ R
を,浮動小数点数として代入すると,β, p
に依存する定数ǫ > 0
が存在して,| x − x | < ǫ | x |
が成り立つ. ここで,
x ∈ F
は,x
をF = F (β, p)
に代入した値である.•
一般に使われる「丸めモード」は,通常の「四捨五入」に対応する「nearest rounding」（一 番近い数への近似）であり, その際には,ǫ = (β/2)β − p = β 1 − p /2
となる. 特にF (2, 53)
に 対しては,ǫ = 2 − 53
となる.•
浮動小数点数x, y ∈ F
に対して,その四則演算を·
とおく. すなわち,x · y
は,x + y, x − y, x × y, x/y
のいずれかである. この時, 実際に行われる演算の結果をx ⊙ y (x ⊕ y, x ⊖ y, x ⊗ y, x ⊘ y
のいずれか)とおき,x · y
とx ⊙ y
の相対誤差を評価することを考える.•
実際に行われる浮動小数点演算の概要は以下の通り.–
乗算（および除算）の場合（x⊗ y
の計算）1. x ∈ F
の仮数部とy ∈ F
の仮数部の積（除算）を計算する.2.
仮数部の積に桁上がりが無い場合,仮数部の積をp
桁に丸める.仮数部の積に桁上がりがある場合,必要な桁移動を行ったうえで
p
桁に丸める.3. x ∈ F
の指数部とy ∈ F
の指数部の和（差）を計算する.4.
得られた仮数部と指数部の結果から浮動小数点数x ⊗ y (x ⊘ y)
を構成する.この手順で「誤差」が生じるのは,仮数部の積の桁移動と丸めの部分である.
–
同一の符号をもつ数に対する加算の場合（x⊕ y
の計算）以下では,x ≥ y
と仮定する.1. y ∈ F
の仮数部をx ∈ F
の仮数部の桁にそろうように桁移動する.2.
桁移動した後の仮数部の和を計算する.3.
仮数部の和に桁上がりが無い場合,仮数部の和をp
桁に丸める.仮数部の和に桁上がりがある場合,必要な桁移動を行ったうえで
p
桁に丸める.4.
得られた仮数部の結果から浮動小数点数x ⊕ y
を構成する. 指数部は（桁上がりが あれば,それを調整して）x
の指数部をそのまま使えば良い.この手順で「誤差」が生じるのは,
y
（小さい方の数）の最初の桁移動,仮数部の和の桁 移動と丸めの部分である.–
同一の符号をもつ数に対する減算の場合（x⊖ y
の計算） も加算と手順は同じである が,このままではうまく行かない.∗ F(10, 3)
において,x = 1.00 × 10 0 , y = 9.99 × 10 − 1
に対して,x ⊖ y
を上の手順で 計算すると,x − y = 0.001, x ⊖ y = 0.01
となり,相対誤差が90
となってしまう.∗
このような問題を解消するために「保護桁」が用いられる. 保護桁とは,浮動小数 点数演算を行う際に仮数部の桁数を（演算の時のみ）余分に与えたものである.•
浮動小数点演算の誤差については,次の事実が成り立つ.保護桁が無くても,加算・乗算・除算に関しては,
| (x · y) − (x ⊙ y) | ≤ 2ǫ | x · y | , x, y ∈ F
が成り立つ. また,保護桁が
1
桁あれば,減算に対しても同じ評価が成り立つ.通常の浮動小数点演算に関しては,保護桁が
1
桁用意されているので,四則演算については, 相対誤差2ǫ
で計算できていると考えて良い.•
加算の場合の上の評価の証明は,以下の通り.いま,
x ≥ y > 0
と仮定する. この時,y
を桁移動すると,β 1 − p
の丸め誤差が発生する. いま, 和の結果に桁あがりがないと仮定すると,和x + y ≥ 1
であり,β − p+1 = 2ǫ
であるので,| (x ⊕ y) − (x + y) | ≤ 2ǫ ≤ 2ǫ | x + y |
が成り立つ.次に桁上がりがあると仮定すると,
y
のシフト時のrounding error 2ǫ
の他に, 桁上がりの桁 移動時に, (1/2)β− 2+p
の誤差が発生する. （1/2は丸めを行うため）また,x + y ≥ β
である ので,| (x ⊕ y) − (x + y) | ≤ (β − p+1 + (1/2)β − p+2 ) = β − p (1 + β/2)/β ≤ β − p (1 + β/2) | x + y |
が成り立つ. ここで,β − p (1 + β/2) ≤ β 1 − p = 2ǫ
を使えば,| (x ⊕ y) − (x + y) | ≤ 2ǫ | x + y |
が成り立つ.• β = 2, nearest rounding
の場合には, 2ǫ= 2 1 − p
である. これは,「２進数で書いたとき, 最 後の１桁が信用できない」ということをあらわしている. これを,「最後の１桁が汚染されて いる」,「汚染桁は1
桁」,「1 ulpの誤差がある」などと表現する.▼ 自然対数の底
e
の近似値の計算について（続き）•
前回,a n =
1 + 1 n
n
の値を計算することによって,自然対数の底
e
の近似値を求めることを考え,数学的にはa n = e − e
2n + O(n − 2 )
が成り立つことを示した.•
実際に計算結果として得ることができるのは,a n
の浮動小数点演算の結果であり, (1 + 1/n) をn
回乗算することによって得られる結果a n
は,| a n − a n | ≤ 2(n − 1)ǫ | a n |
となる. 実際,δ = 2ǫ
とおけば,| ((1 + 1/n) ⊗ (1 + 1/n)) − ((1 + 1/n) × (1 + 1/n)) | < δ | (1 + 1/n) 2 | ,
| (((1 + 1/n) ⊗ (1 + 1/n)) ⊗ (1 + 1/n)) − (((1 + 1/n) × (1 + 1/n)) × (1 + 1/n)) |
< δ | ((1 + 1/n) ⊗ (1 + 1/n)) × (1 + 1/n) | < 2δ | (1 + 1/n) 3 |
となる. 繰り返し２乗法（高速乗算法）を用いても, (1 + 1/n)
2
k の段階で,おおよそ2 k δ
の 誤差が含まれているので, 同一の結果を得ることになる. （この計算は, (1 + 1/n)そのもの のがF
で正確に表現できていると仮定しているが, 実際には, (1 + 1/n)がF
で相対誤差ǫ
を含んでいると考えるべきである.）•
以上により,実際に得られた(1 + 1/n) n
の値a n
は,| a n − e | ∼ e
2n + nδ, δ = 2n
となることがわかる.•
上式の右辺は,n = p
e/(2δ) ∼ √
δ
の時に最小となり,その時の値は,2n e + nδ ∼ δ
となる.いま,
δ ∼ 2 − 52
であるので,およそ10 − 8
程度の精度でしかe
の近似値を求めることはでき ない.▼ 参考文献
•
浮動小数点演算について詳細な解説がある文献で,容易に入手可能なものとしては,Sun Microsystems, “What every computer scientists should know about floating-point arithmetic”, http://dlc.sun.com/pdf/800-7895/800-7895.pdf
がある.
•
多倍長数値演算パッケージとしては,The GNU multiple precision arithmetic library (GNU MP), http://gmplib.org/
が代表的なものである.
● 講義資料
▼ 円周率の近似値を求める
以下の図は,半角の公式を利用して,
π
の近似値を求めた例である.（πと近似値の値との絶対誤 差を表示している）1e-09 1e-08 1e-07 1e-06 1e-05 0.0001 0.001 0.01 0.1 1
1 10 100 1000 10000 100000 1e+06 1e+07 1e+08 1e+09
Calculate Pi (1)
inner outer(1) outer(2)
次の図は,半角の公式を書き直して
π
の近似値を求めた例である.1e-16 1e-14 1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1
1 10 100 1000 10000 100000 1e+06 1e+07 1e+08 1e+09
Calculate Pi (2)
inner(3) outer(3) inner(4) outer(4)
次の図は, arctan(x)の近似値をテイラー級数を用いて計算した例である.
1e-16 1e-14 1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1
1 10 100 1000 10000
x = 1 x = 0.999 x = 0.99 x = 0.9
いろいろな方法で
π
の近似値を計算した際の計算時間1e-05 0.0001 0.001 0.01 0.1 1 10 100 1000 10000 100000
2^10 2^20 2^30
circle length Machine Gauss-Legendre Quadratic Borwein
● 実習内容
（以下で「★」の数は推奨の程度を示します. 多いほど推奨の度合いが大きい. また,「†」は 難易度またはマニアックな程度を表します. 多いほど難しいかマニアックかです.）
1.
（★★★）半角の公式を書き直した式ℓ 2n =
√ 2ℓ n
q 1 + p
1 − (ℓ/n) 2 ,
L 2n = 2L n
p 1 + (L n /n) 2 + 1
を用いて
π
の近似値を計算しなさい. さらに,π
とそれぞれの近似値の絶対誤差のグラフを 書きなさい.2.
（★★★）半角の公式とそれから派生する公式1
L 2n
= 1 2
1 ℓ n
+ 1 L n
, ℓ 2n = p
ℓ n L 2n ,
を用いて
π
の近似値を計算しなさい. さらに,π
とそれぞれの近似値の絶対誤差のグラフを 書きなさい.3.
（★★★）ライプニッツの公式でπ
の近似値を計算しなさい.4.
（★★★）マチンの公式でπ
の近似値を計算しなさい. また, ライプニッツの公式での計算 時間（計算回数）との比較をしなさい.5.
（★★★）自然対数の底e
の近似値を相対誤差10 − 14
以内で計算しなさい.6.
（★★†）x∈ [0, 1]
に対して,e x
の近似値を相対誤差10 − 14
以内で計算しなさい. プログラ ムの第一引数にx
の値を与えて,e x
の近似値を出力するように書きなさい.7.
（★★）sin(1.0)の近似値を相対誤差10 − 14
以内で計算しなさい.8.
（★†）x∈ [ − π/2, π/2]
に対して, sin(x)の近似値を相対誤差10 − 14
以内で計算しなさい.プログラムの第一引数に
x
の値を与えて, sin(x)の近似値を出力するように書きなさい.9.
（★†††）log(2.0) の近似値を相対誤差10 − 14
以内で計算しなさい. ただし, 「高速」に 計算する方法を考えて計算すること.10.
（★††）いろいろな方法でπ
の近似値を計算しなさい. たとえば,π 2
6 = X
n=1
1 n 2 , π 4
90 = X
n=1
1 n 4 , π
6 = arcsin(1/2),
などがある. その場合,マチンの公式と比較してどちらが高速に計算できるかを確かめなさ
