● 講義資料
▼ 講義予定
•
円周率の近似値の計算と浮動小数点演算の誤差（続き）•
テイラー級数による円周率の近似値の計算● 前回の講義のまとめ
▼ 浮動小数点演算
•
ここでは, 考える浮動小数点数の体系を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
以内で計算しなさい.プログラムの第一引数に
