● 講義資料
▼ 講義予定
• 円周率の近似値の計算と浮動小数点演算の誤差
• テイラー級数による円周率の近似値の計算
● 前回の講義のまとめ
▼ 浮動小数点演算
• ここでは, 考える浮動小数点数の体系を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×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の誤差がある」などと表現する.
▼ 自然対数の底
e
の近似値の計算について（続き）• 前回,
an=
1 + 1 n
n
の値を計算することによって,自然対数の底eの近似値を求めることを考え,数学的には an =e− e
2n+O(n−2) が成り立つことを示した.
• 実際に計算結果として得ることができるのは,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|
となる. 繰り返し２乗法（高速乗算法）を用いても, (1 + 1/n)2k の段階で,おおよそ2kδの 誤差が含まれているので, 同一の結果を得ることになる. （この計算は, (1 + 1/n)そのもの のがF で正確に表現できていると仮定しているが, 実際には, (1 + 1/n)が F で相対誤差ǫ を含んでいると考えるべきである.）
• 以上により,実際に得られた(1 + 1/n)n の値an は,
|an−e| ∼ e
2n+nδ, δ= 2n となることがわかる.
• 上式の右辺は, n=p
e/(2δ)∼√
δ の時に最小となり,その時の値は, 2ne +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)
次の図は,半角の公式を書き直してπの近似値を求めた例である.
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-14 1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1
x = 1 x = 0.999 x = 0.99 x = 0.9
● 実習内容
（以下で「★」の数は推奨の程度を示します. 多いほど推奨の度合いが大きい. また,「†」は 難易度またはマニアックな程度を表します. 多いほど難しいかマニアックかです.）
1. （★★★）半角の公式を書き直した式
ℓ2n=
√2ℓn
q 1 +p
1−(ℓ/n)2 ,
L2n= 2Ln
p1 + (Ln/n)2+ 1
を用いてπ の近似値を計算しなさい. さらに,πとそれぞれの近似値の絶対誤差のグラフを 書きなさい.
2. （★★★）半角の公式とそれから派生する公式 1
L2n
=1 2
1 ℓn
+ 1 Ln
, ℓ2n =p
ℓnL2n,
を用いてπ の近似値を計算しなさい. さらに,πとそれぞれの近似値の絶対誤差のグラフを 書きなさい.
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 以内で計算しなさい. ただし, 「高速」に 計算する方法を考えて計算すること.
● 繰り返し２乗法のプログラム例
int power(int a, unsigned int n) {
int p = 1 ; while(n > 0) {
if (n%2 == 1) p = p*a ; a = a*a ;
n = n/2 ; }
return p ; }
int power(int a, unsigned int n) {
int p = 1 ; while(n) {
if (n&1) p *= a ; a *= a ;
n >>= 1 ; }
return p ; }
int power(int a, unsigned int n) {
int p = 1 ; for(;n;n>>=1) {
if (n&1) p *= a ; a *= a ;
}
return p ; }
