● 講義資料
▼ 講義予定
•
円周率の近似値の計算と浮動小数点演算の誤差•
テイラー級数による円周率の近似値の計算● 前回の講義のまとめ
▼ 浮動小数点演算
•
ここでは, 考える浮動小数点数の体系を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
桁に丸める.この手順で「誤差」が生じるのは,
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 }
は,単調増加か つ上に有界である. また,簡単な計算から,任意のn ≥ 2
に対して, 2< a n < 3
が成り立つ.よって,
e = lim n→∞ a n
が存在し, 2< e ≤ 3
が成り立つ.•
数学的にはa n = e − e
2n + O(n −2 )
が成り立つ.証明は, log(1 +
x)
のテイラー展開にx = 1/n
を代入し, さらにその結果をexp(x)
のテイ ラー展開に適用すればよい.•
上の評価式は（nに関する− 2
次以上の無視すれば| a n − e | ≤ e 2n
と書けるので,– e
の近似値を誤差ǫ
以下で求めたい場合,e
2n < ǫ
をみたす最小の自然数N
をとり,n ≥ N
となるa n
を計算すればよいことがわかる. これは,数列{ a n }
がe
に収束する ことの定義であるǫ-N
論法そのものである. この場合,N
としてはe
2ǫ ≤ 3
2ǫ
であることから,
N ≥ ⌊ 3
2ǫ ⌋
ととれば十分である.– e
の近似値を相対誤差ǫ
以下で求めたい場合,1
2n < ǫ
をみたす最小の自然数N
をとり,n ≥ N
となるa n
を計算すればよいことがわかる. この場合,N
としてはN ≥ ⌊ 1
2ǫ ⌋
と とれば十分である.•
しかし,浮動小数点演算の演算誤差によって,実際のa n
の計算値a n
は,| e − a n | ∼ e 2n + Cn
という挙動を示し,ある程度以上大きな
n
に対しては,a n
はe
を近似しなくなる.•
実際に計算結果として得ることができるのは,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)
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
の近似値を求めることはでき ない.● 講義資料
▼ 円周率の近似値を求める
半角の公式を利用して,
π
の近似値を求めた例（πと近似値の値との絶対誤差を表示している）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)
ℓ 2n = 2n s
1 − p
1 − (ℓ n /n) 2
2 ,
L 2n = 2n 2 L n
p 1 + (L n /n) 2 − 1 , ℓ n = n sin( π
n ), L n = n tan( π
n ),
半角の公式を書き直して
π
の近似値を求めた例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)
ℓ 2n =
√ 2ℓ n
q 1 + p
1 − (ℓ/n) 2 ,
L 2n = 2L n
p 1 + (L n /n) 2 + 1 1
L 2n
= 1 2
1 ℓ n
+ 1 L n
, ℓ 2n = p
ℓ n L 2n ,
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
arctan(x) = x − x 3 3 + x 5
5 + · · · + ( − 1) 2n x 2n+1
2n + 1 + · · ·
● 実習内容
（以下で「★」の数は推奨の程度を示します. 多いほど推奨の度合いが大きい. また,「†」は 難易度またはマニアックな程度を表します. 多いほど難しいかマニアックかです.）
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
以内で計算しなさい.プログラムの第一引数に
