講義資料前回の講義のまとめ

2014 年度・前期・数理解析・計算機数学３・第４回
1
● 講義資料
▼ 講義予定
• 円周率の近似値の計算と浮動小数点演算の誤差
• テイラー級数による円周率の近似値の計算
● 前回の講義のまとめ
▼ 浮動小数点演算
• ここでは, 考える浮動小数点数の体系を 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 の指数部をそのまま使えば良い.
Oct. 22, 2014, Version: 1.0
[email protected]
2014 年度・前期・数理解析・計算機数学３・第４回
2
この手順で「誤差」が生じるのは, y （小さい方の数）の最初の桁移動, 仮数部の和の桁
移動と丸めの部分である.
– 同一の符号をもつ数に対する減算の場合（x ⊖ y の計算）も加算と手順は同じである
が, このままではうまく行かない.
∗ F (10, 3) において, x = 1.00 × 100 , 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ǫ = 21−p である. これは, 「２進数で書いたとき, 最
後の１桁が信用できない」ということをあらわしている. これを, 「最後の１桁が汚染されて
いる」, 「汚染桁は 1 桁」, 「1 ulp の誤差がある」などと表現する.
Oct. 22, 2014, Version: 1.0
[email protected]
2014 年度・前期・数理解析・計算機数学３・第４回
3
▼ 自然対数の底 e の近似値の計算について
• 今回考えるのは,
n
1
an = 1 +
n
の値を計算することによって, 自然対数の底 e の近似値を求める. 数列 {an } は, 単調増加か
つ上に有界である. また, 簡単な計算から, 任意の n ≥ 2 に対して, 2 < an < 3 が成り立つ.
よって, e = limn→∞ an が存在し, 2 < e ≤ 3 が成り立つ.
• 数学的には
an = e −
e
+ O(n−2 )
2n
が成り立つ.
証明は, log(1 + x) のテイラー展開に x = 1/n を代入し, さらにその結果を exp(x) のテイ
ラー展開に適用すればよい.
• 上の評価式は（n に関する −2 次以上の無視すれば
|an − e| ≤
e
2n
と書けるので,
e
– e の近似値を誤差 ǫ 以下で求めたい場合,
< ǫ をみたす最小の自然数 N をとり,
2n
n ≥ N となる an を計算すればよいことがわかる. これは, 数列 {an } が e に収束する
3
e
≤
であるこ
ことの定義である ǫ-N 論法そのものである. この場合, N としては
2ǫ
2ǫ
3
とから, N ≥ ⌊ ⌋ ととれば十分である.
2ǫ
1
< ǫ をみたす最小の自然数 N をとり,
– e の近似値を相対誤差 ǫ 以下で求めたい場合,
2n
1
n ≥ N となる an を計算すればよいことがわかる. この場合, N としては N ≥ ⌊ ⌋ と
2ǫ
とれば十分である.
• しかし, 浮動小数点演算の演算誤差によって, 実際の an の計算値 an は,
|e − an | ∼
e
+ Cn
2n
という挙動を示し, ある程度以上大きな n に対しては, an は e を近似しなくなる.
• 実際に計算結果として得ることができるのは, an の浮動小数点演算の結果であり, (1 + 1/n)
を n 回乗算することによって得られる結果 an は,
|an − an | ≤ 2(n − 1)ǫ|an |
となる. 実際, δ = 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 |
Oct. 22, 2014, Version: 1.0
[email protected]
2014 年度・前期・数理解析・計算機数学３・第４回
4
k
となる. 繰り返し２乗法（高速乗算法）を用いても, (1 + 1/n)2 の段階で, おおよそ 2k δ の
誤差が含まれているので, 同一の結果を得ることになる. （この計算は, (1 + 1/n) そのもの
のが F で正確に表現できていると仮定しているが, 実際には, (1 + 1/n) が F で相対誤差 ǫ
を含んでいると考えるべきである.）
• 以上により, 実際に得られた (1 + 1/n)n の値 an は,
|an − e| ∼
e
+ nδ,
2n
δ = 2n
となることがわかる.
p
√
e
+ nδ ∼ δ となる.
• 上式の右辺は, n = e/(2δ) ∼ δ の時に最小となり, その時の値は, 2n
−52
−8
いま, δ ∼ 2
であるので, およそ 10 程度の精度でしか e の近似値を求めることはでき
ない.
● 講義資料
▼ 円周率の近似値を求める
半角の公式を利用して, π の近似値を求めた例（π と近似値の値との絶対誤差を表示している）
Calculate Pi (1)
1
inner
outer(1)
outer(2)
0.1
0.01
0.001
0.0001
1e-05
1e-06
1e-07
1e-08
1e-09
1
10
100
ℓ2n = 2n
1000
s
1−
10000
100000
1e+06
1e+07
1e+08
1e+09
p
1 − (ℓn /n)2
,
2
2n2 p
1 + (Ln /n)2 − 1 ,
Ln
π
π
ℓn = n sin( ), Ln = n tan( ),
n
n
L2n =
Oct. 22, 2014, Version: 1.0
[email protected]
2014 年度・前期・数理解析・計算機数学３・第４回
5
半角の公式を書き直して π の近似値を求めた例
Calculate Pi (2)
1
inner(3)
outer(3)
inner(4)
outer(4)
0.01
0.0001
1e-06
1e-08
1e-10
1e-12
1e-14
1e-16
1
10
100
1000
10000
100000
1e+06
1e+07
1e+08
1e+09
√
2ℓn
= q
,
p
1 + 1 − (ℓ/n)2
ℓ2n
2Ln
L2n = p
1 + (Ln /n)2 + 1
1
1 1
1
,
=
+
L2n
2 ℓn
Ln
p
ℓ2n = ℓn L2n ,
arctan(x) の近似値をテイラー級数を用いて計算した例
1
x=1
x = 0.999
x = 0.99
x = 0.9
0.01
0.0001
1e-06
1e-08
1e-10
1e-12
1e-14
1e-16
1
10
arctan(x) = x −
Oct. 22, 2014, Version: 1.0
100
1000
10000
x5
(−1)2n x2n+1
x3
+
+ ···+
+ ···
3
5
2n + 1
[email protected]
2014 年度・前期・数理解析・計算機数学３・第４回
6
● 実習内容
（以下で「★」の数は推奨の程度を示します. 多いほど推奨の度合いが大きい. また, 「†」は
難易度またはマニアックな程度を表します. 多いほど難しいかマニアックかです.）
1. （★★★）半角の公式を書き直した式
ℓ2n
√
2ℓn
=q
,
p
1 + 1 − (ℓ/n)2
2Ln
L2n = p
1 + (Ln /n)2 + 1
を用いて π の近似値を計算しなさい. さらに, π とそれぞれの近似値の絶対誤差のグラフを
書きなさい.
2. （★★★）半角の公式とそれから派生する公式
1
1 1
1
,
=
+
L2n
2 ℓn
Ln
p
ℓ2n = ℓn L2n ,
を用いて π の近似値を計算しなさい. さらに, π とそれぞれの近似値の絶対誤差のグラフを
書きなさい.
3. （★★★）ライプニッツの公式で π の近似値を計算しなさい.
4. （★★★）マチンの公式で π の近似値を計算しなさい. また, ライプニッツの公式での計算
時間（計算回数）との比較をしなさい.
5. （★★★）自然対数の底 e の近似値を相対誤差 10−14 以内で計算しなさい.
6. （★★†）x ∈ [0, 1] に対して, ex の近似値を相対誤差 10−14 以内で計算しなさい. プログラ
ムの第一引数に x の値を与えて, ex の近似値を出力するように書きなさい.
7. （★★）sin(1.0) の近似値を相対誤差 10−14 以内で計算しなさい.
8. （★†）x ∈ [−π/2, π/2] に対して, sin(x) の近似値を相対誤差 10−14 以内で計算しなさい.
プログラムの第一引数に x の値を与えて, sin(x) の近似値を出力するように書きなさい.
9. （★†††）log(2.0) の近似値を相対誤差 10−14 以内で計算しなさい. ただし, 「高速」に
計算する方法を考えて計算すること.
Oct. 22, 2014, Version: 1.0
[email protected]

Download Report