● 講義資料

● 講義資料

▼ 講義予定

• 浮動小数点数とは何か（先週の例を元に）

• ニュートン法

● 前回の講義のまとめ

• 二分法を使って√

2の近似値を計算した.

• 二分法で√

2の近似値を計算するためのアルゴリズムは以下の通り.

1. l0= 0,u0= 2,k= 0とおく.

2. l=lk,u=uk とおき,c= (l+r)/2 を計算する.

3. c²<2ならばlk+1=c,rk+1=rk,c²>2 ならばlk+1=lk, rk+1=cとおく.

4. k:=k+ 1とおき, 2に戻る.

• この方法は,kが指定の回数になったときに終了させない限り,無限ループに陥ってしまう.

• このアルゴリズムはf(x) =x²−2に対して,中間値の定理を適用していることに注意が必 要である.

▼ 前回の計算例

• 前回のプログラム１(example_01_1.c)では,繰り返し回数を指定した. √

2を求めるために は, 53回の繰り返し終了時が最も近似の精度が良くなっていたはず.

• 前回のプログラム１(example_01_2.c)のループ継続条件を“fabs(lower - upper) > 3.0E-16”

程度に設定すると,最も近似の精度が良くなっていたはず.

● 講義資料

▼ 前回のプログラム２を書き換えてみる

★ プログラム例１

(example_02_1.c)

/**

* 二分法を用いて Sqrt(2) の値を計算する

*/

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#define EPSILON (3.0E-16) void calc_sqrt(double x) ;

int main(int argc, char **argv) {

if (argc > 1) calc_sqrt(atof(argv[1])) ; return 0 ;

}

void calc_sqrt(double x) {

double lower, upper, square_root ;

lower = 0.0 ; upper = x ;

while(fabs(lower-upper) > EPSILON) { square_root = (lower + upper)/2.0 ;

if (square_root*square_root > x) upper = square_root ; else lower = square_root ;

printf("%e\n", fabs(lower-upper)) ; }

return ; }

• このプログラムは,コマンドライン引数に数値xを与えて, √xの近似値を二分法で計算す る. ただし,このプログラムもxの値によっては停止しない.

• ループ終了条件は,|l−r| ≤ǫであり,ǫの値はEPSILONで与えられている.

– x= 2.0の時には,ǫ= 3.0×10⁻¹⁶でよい.

– x= 20.0 の時には,|l−r| ∼8.9×10⁻¹⁶より小さくならない.

– x= 200.0の時には,|l−r| ∼1.8×10⁻¹⁵より小さくならない.

• したがって,|l−r|を判定基準することは間違いであることがわかる.

★ プログラム例２

(example_02_2.c)

/**

* 二分法を用いて Sqrt(2) の値を計算する

*/

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#define EPSILON (4.0E-16) void calc_sqrt(double x) ;

int main(int argc, char **argv) {

if (argc > 1) calc_sqrt(atof(argv[1])) ; return 0 ;

}

void calc_sqrt(double x) {

double lower, upper, square_root ;

lower = 0.0 ; upper = x ; square_root = x ;

while(fabs(lower-upper) > EPSILON*fabs(square_root)) { square_root = (lower + upper)/2.0 ;

if (square_root*square_root > x) upper = square_root ; else lower = square_root ;

printf("%e %e\n", fabs(lower - upper), fabs((lower - upper)/square_root)) ; }

printf("%.16e %.16e\n", square_root, square_root*square_root) ; return ;

}

• 例１を書き換えて, ループ終了条件として,|l−r| ≤ǫαと設定した. (ǫ= 4.0×10⁻¹⁶, αは その時点での近似値)

• この場合,x= 2.0,x= 20.0,x= 200.0,x= 5.0×10³⁰⁰,x= 5.0×10⁻³⁰⁰など,すべての代 入可能な数値に対してプログラムが停止し,求めたx >1.0の場合には,√

xは「ほぼ」望ま しい近似値を持っていることがわかる.

• したがって,近似の判断基準として|l−r|< ǫαは適切なものであると考えられる.

• ǫ= 4.0×10⁻¹⁶ととる理由は2⁻⁵²<4.0×10⁻¹⁶<2⁻⁵¹であるから.

▼ 注意事項

• 浮動小数点数をprintf関数を表示する場合には,以下の２種類の表示方法がある

– %f: 通常の12.345などと表示する – %e: 12.345を 1.2345E02+と表示する.

– %.Neとすると,小数点以下N桁を表示する

– %M.Neとすると,全体で M桁,小数点以下N桁を表示する(+M+は省略可能) – これ以外にも書式指定があるので,それらはprintf関数のマニュアルを参照

• 一般に,ある値xの近似値xが与えられたとき,

– ǫ=|x−x|を「絶対誤差」と呼び,

– |x−x|=ǫ|x|となるǫを「相対誤差」と呼ぶ.

• 実際の計算の相対誤差は, いくら小さくてもǫ∼10⁻¹⁶程度にしかならないので,%e表示の 時に意味のある桁は,高々15桁である. （もっと長い桁を表示させることは可能であるが, 意味のない数値が並ぶだけである）

▼

IEEE 784

s eN−1 · · · e5 e4 e3 e2 e1 e0 f1 f2 · · · fM−2fM−1

パラメータ 単精度 倍精度 ４倍精度

N 8 11 15

M 24 53 113

ビット長 32 64 128

最大指数 +127 +1023 +16383

最小指数 −126 −1022 −16382

バイアスb 127 1024 16383

ǫM = 2^1−M 2⁻²³ 2⁻⁵² 2⁻¹¹¹

1.19×10⁻⁷ 2.22×10⁻¹⁶ 3.85×10⁻³⁴

ǫM = 2^−M 2⁻²⁴ 2⁻⁵³ 2⁻¹¹²

5.96×10⁻⁸ 1.11×10⁻¹⁶ 1.92×10⁻³⁴

x= (−1)^s(1 +

MX−1

j=1

fj2^−j)×2^e−b, e=

NX−1

k=0

ek2^k, b=

NX−2

k=0

2^k. (1)

いずれの場合にでも,

• 全てのkに対してek= 1ならば「オーバフロー」,

• 全てのkに対してek= 0ならば±0または「非正規数」

を示す.

仮数部２桁,指数部 −1≤e≤1 の２進浮動小数点数で正確に表示可能な値

0=0.00

q

1/8=0.01×2−1

a

2/8=0.10×2−1

a

3/8=0.11×2−1

a

−1/8=−0.01×2−1

a

−2/8=−0.10×2−1

a

−3/8=−0.11×2−1

a

4/8=1.00×2−1

q

−4/8=−1.00×2−1

q

5/8=1.01×2−1

q

−5/8=−1.01×2−1

q

6/8=1.10×2−1

q

−6/8=−1.10×2−1

q

7/8=1.11×2−1

q

−7/8=−1.11×2−1

q

4/4=1.00×20

q

−4/4=−1.00×20

q

5/4=1.01×20

q

−5/4=−1.01×20

q

6/4=1.10×20

q

−6/4=−1.10×20

q

7/4=1.11×20

q

−7/4=−1.11×20

q

4/2=1.00×21

q

−4/2=−1.00×21

q

5/2=1.01×21

q

−5/2=−1.01×21

q

6/2=1.10×21

q

−6/2=1.10×21

q

7/2=1.11×21

q

−7/2=−1.11×21

q

Underflow

subnormal

e=−1

e=−1 e= 0

e= 0

e= 1 e= 1

-

overflow

overflow

この浮動小数点数の体系で表現できる数は,図のなかの黒円で示したものに限られる. ま

た, “subnormal”と示した数は,有効桁が少なくなっていることに注意.

したがって, 1/2 + 1/8 = 5/8と正確に計算することが可能であるが, 1/2 + 1/16は正確 な計算ができないことがわかる.

● 実習内容

（以下で「★」の数は推奨の程度を示します. 多いほど推奨の度合いが大きい. また,「†」は 難易度またはマニアックな程度を表します. 多いほど難しいかマニアックかです.）

1. （★★★）プログラム２を入力し, 種々のxの値に対して, ǫ= 4.0×10⁻¹⁶ を使うと,実際 に√

xの近似値を求めることができることを確かめなさい.

2. （★）プログラム２は, 0< x <1なる入力に対しては正しく√

xの近似値を計算できない.

どのように修正すればよいか？

3. （★★★）ニュートン法を用いて√

2 の近似値を相対誤差4.0×10⁻¹⁶ で求めるプログラム を書きなさい.

4. （★★★）√

2の近似値をニュートン法で求める際に,各繰り返しでの近似値と√

2との誤差

を以下のように出力しなさい. さらに,その結果からgnuplotを用いて, 誤差の振る舞いのグ ラフを作成しなさい.

0 8.578644e-02 1 2.453104e-03 2 2.123901e-06 3 1.594724e-12 4 2.220446e-16 5 2.220446e-16

1e-16 1e-14 1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1

0 1 2 3 4 5

Newton method for sqrt(2)

5. （★★★）関数f(x) = (x²−2)^k (k= 2,3, . . .)にニュートン法を適用して√

2 の近似値を 求め,また,それぞれのkに対して,各繰り返しでの近似値と√

2との誤差の振る舞いをグラ フに表しなさい.

6. （★★）ニュートン法を用いて√

xの近似値を相対誤差4.0×10⁻¹⁶ で求めるプログラムを 書きなさい. さらに, その近似値αx に対して|α²_x−x| を計算し, 求めた近似値が相対誤差 4.0×10⁻¹⁶以内であることを確かめるプログラムを書きなさい.

7. 「大きな」xに対して,ニュートン法のプログラムを適用することを考える.

(a) 繰り返し回数が多くなることを確かめなさい.

(b) （†）繰り返しの初期段階において,近似の様子（誤差の挙動）があまり良くないこと がわかる. その理由を考察しなさい.

(c) そのような挙動を示すxは「大きな」xに限るのかを考察しなさい.

8. （★†）MacOSX の unsigned int 型整数は 32 bits integer である. このことを用いて,

unsigned int型整数として格納可能な自然数nの平方根の整数部分を求めるプログラムを

書きなさい. また,その繰り返し回数何回程度になるか？さらに,求めた値が,実際にnの平 方根の整数部分になっていることを,確かめなさい.