● 講義資料

● 講義資料

▼ 講義予定

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

• 自然対数の底 eの近似値の計算

● 前回の講義のまとめ

• 二分法を使って√

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に対して,中間値の定理を適 用していると考えてもよい. すなわち, 連続関数f : [a, b]−→Rが,区間[a, b]上で単調増加であり, f(a)<0,f(b)>0 を満たしているとき,以下の二分法で,区間[a, b]に一意的に存在するf(x) = 0 の解x∞に収束する数列{xk} をつくることができる.

1. l0=a,r0=b,k= 0とおく.

2. xk= (lk+rk)/2 を計算する.

3. f(xk)<0ならばlk+1=xk,rk+1=rk,f(xk)>0ならばlk+1 =lk,rk+1 =xk とおく.

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

この時,すべてのkについてf(lk)<0,f(rk)>0 であるので,x∞∈[lk, rk]が成り立つので,

|x^∞−xk| ≤1

2|rk−lk|= 2^−(k+1)|b−a| が成り立つ. よって,xk は x∞ に収束する.

• この二分法において,数学的にはf(xk) = 0となる可能性がある. しかし,コンピュータによるプロ グラムを考える場合には,f(xk)が演算誤差の範囲内で 0に近いときにプログラムを停止して, その ときの xk を x∞ の近似値と考える. 逆に言えば, x∈[a, b] において, f(x)の値を演算誤差が少な い状況で計算可能な関数 f にしか,この二分法を適用することはできない.

▼ 前回の計算例

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

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

★ プログラム例１

(example_02_1.c)

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

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#define EPSILON (1.0E-15) 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 = (lower+upper)/2.0 ; 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)) ; }

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

}

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

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

• x= 2.0の時には,プログラムは終了するが,x= 200.0の時には,プログラムは終了しない. すなわ

ち,プログラム終了条件として,|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 ; double x0 = x ;

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("%.12e %.12e\n", square_root,

fabs(square_root*square_root - x0)/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|< ǫαは適切なものであると考えられる.

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

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

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

• 一般に,ある値xの近似値xが与えられたとき, – ǫ=|x−x|を「絶対誤差」と呼び,

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

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

▼ 自然対数の底の近似値を求める

以下の図は, (1 + 1/n)ⁿを計算して,eの近似値を求めた例である.（eと近似値の値との絶対誤差を表 示している）

1e-08 1e-07 1e-06 1e-05 0.0001 0.001 0.01 0.1 1

10 100 1000 10000 100000 1e+06 1e+07 1e+08 1e+09 1e+10

★ プログラム例

(example_02_3.c)

/* double を用いて (1.0 + 1.0/n)^n を計算して, e との誤差を調べる */

#include <stdio.h>

#include <math.h>

#define MAX (10)

double power(double k, unsigned long n) ;

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

unsigned long n ; int i ;

double e ; n = 10 ;

for(i=0;i<MAX;i++) {

e = power(1.0+1.0/n, n) ;

printf("%lu %.14e\n", n, fabs(e - M_E)) ; n *= 10 ;

}

return 0 ; }

/* マジメに n 回の乗算を行う */

double power(double k, unsigned long n) {

double x = 1.0 ; unsigned long i ;

for(i=0;i<n;i++) x *= k ; return x ;

}

• このプログラムでは,n= 10からn= 10¹⁰まで,nを10倍しながら, (1 + 1/n)ⁿ を計算している.

その際 nとしてunsigned long型変数を使う必要がある. int型変数では,n= 10¹⁰までは計算 できない.

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

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 +

M−1X

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は正確 な計算ができないことがわかる.

二分法による√

2の計算経過. 初期区間[0,2], c1 = 1.0

1 00111111 11110000 00000000 00000000 00000000 00000000 00000000 00000000 2 00111111 11111000 00000000 00000000 00000000 00000000 00000000 00000000 3 00111111 11110100 00000000 00000000 00000000 00000000 00000000 00000000 4 00111111 11110110 00000000 00000000 00000000 00000000 00000000 00000000 5 00111111 11110111 00000000 00000000 00000000 00000000 00000000 00000000 6 00111111 11110110 10000000 00000000 00000000 00000000 00000000 00000000 7 00111111 11110110 11000000 00000000 00000000 00000000 00000000 00000000 8 00111111 11110110 10100000 00000000 00000000 00000000 00000000 00000000 9 00111111 11110110 10110000 00000000 00000000 00000000 00000000 00000000 10 00111111 11110110 10101000 00000000 00000000 00000000 00000000 00000000 11 00111111 11110110 10100100 00000000 00000000 00000000 00000000 00000000 12 00111111 11110110 10100010 00000000 00000000 00000000 00000000 00000000 13 00111111 11110110 10100001 00000000 00000000 00000000 00000000 00000000 14 00111111 11110110 10100000 10000000 00000000 00000000 00000000 00000000 15 00111111 11110110 10100000 11000000 00000000 00000000 00000000 00000000 16 00111111 11110110 10100000 10100000 00000000 00000000 00000000 00000000 17 00111111 11110110 10100000 10010000 00000000 00000000 00000000 00000000 18 00111111 11110110 10100000 10011000 00000000 00000000 00000000 00000000 19 00111111 11110110 10100000 10011100 00000000 00000000 00000000 00000000 20 00111111 11110110 10100000 10011110 00000000 00000000 00000000 00000000 21 00111111 11110110 10100000 10011111 00000000 00000000 00000000 00000000 22 00111111 11110110 10100000 10011110 10000000 00000000 00000000 00000000 23 00111111 11110110 10100000 10011110 01000000 00000000 00000000 00000000 24 00111111 11110110 10100000 10011110 01100000 00000000 00000000 00000000 25 00111111 11110110 10100000 10011110 01110000 00000000 00000000 00000000 26 00111111 11110110 10100000 10011110 01101000 00000000 00000000 00000000 27 00111111 11110110 10100000 10011110 01100100 00000000 00000000 00000000 28 00111111 11110110 10100000 10011110 01100110 00000000 00000000 00000000 29 00111111 11110110 10100000 10011110 01100111 00000000 00000000 00000000 30 00111111 11110110 10100000 10011110 01100110 10000000 00000000 00000000 31 00111111 11110110 10100000 10011110 01100110 01000000 00000000 00000000 32 00111111 11110110 10100000 10011110 01100110 01100000 00000000 00000000 33 00111111 11110110 10100000 10011110 01100110 01110000 00000000 00000000 34 00111111 11110110 10100000 10011110 01100110 01111000 00000000 00000000 35 00111111 11110110 10100000 10011110 01100110 01111100 00000000 00000000 36 00111111 11110110 10100000 10011110 01100110 01111110 00000000 00000000 37 00111111 11110110 10100000 10011110 01100110 01111111 00000000 00000000 38 00111111 11110110 10100000 10011110 01100110 01111111 10000000 00000000 39 00111111 11110110 10100000 10011110 01100110 01111111 01000000 00000000 40 00111111 11110110 10100000 10011110 01100110 01111111 00100000 00000000 41 00111111 11110110 10100000 10011110 01100110 01111111 00110000 00000000 42 00111111 11110110 10100000 10011110 01100110 01111111 00111000 00000000 43 00111111 11110110 10100000 10011110 01100110 01111111 00111100 00000000 44 00111111 11110110 10100000 10011110 01100110 01111111 00111010 00000000 45 00111111 11110110 10100000 10011110 01100110 01111111 00111011 00000000 46 00111111 11110110 10100000 10011110 01100110 01111111 00111011 10000000 47 00111111 11110110 10100000 10011110 01100110 01111111 00111011 11000000 48 00111111 11110110 10100000 10011110 01100110 01111111 00111011 11100000 49 00111111 11110110 10100000 10011110 01100110 01111111 00111011 11010000 50 00111111 11110110 10100000 10011110 01100110 01111111 00111011 11001000 51 00111111 11110110 10100000 10011110 01100110 01111111 00111011 11001100 52 00111111 11110110 10100000 10011110 01100110 01111111 00111011 11001110 53 00111111 11110110 10100000 10011110 01100110 01111111 00111011 11001101 54 00111111 11110110 10100000 10011110 01100110 01111111 00111011 11001100 55 00111111 11110110 10100000 10011110 01100110 01111111 00111011 11001100 56 00111111 11110110 10100000 10011110 01100110 01111111 00111011 11001100 57 00111111 11110110 10100000 10011110 01100110 01111111 00111011 11001100

またはマニアックな程度を表します. 多いほど難しいかマニアックかです.）

1. （★★★）プログラム２を入力し,種々のxの値に対して,ǫ= 4.0×10⁻¹⁶を使うと,実際に√ xの 近似値を求めることができることを確かめなさい.

2. （★★★）プログラム２の最後の行では, 入力値 x, 近似値αに対して,|x−α²|/α を求めている.

この値と, 実際の近似値の相対誤差 |√

x−α|/√

xとの関係を求めなさい. また,近似値を相対誤差 ǫで求めるためには,プログラム中のEPSILONの値をどのように決めればよいかを考えなさい.

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

xの近似値を計算できない. どのよ うに修正すればよいか？

4. （★★）MacOSXのunsigned int型整数は32 bits integerである. このことを用いて,unsigned int型整数として格納可能な自然数nの平方根の整数部分を求めるプログラムを書きなさい. また, その繰り返し回数が何回程度になるか？さらに,求めた値が,実際にnの平方根の整数部分になって いることを,確かめなさい.

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

unsigned int 型整数の積と商を,乗算・除算を用いずに計算するプログラムを書きなさい. また,

unsigned int型変数の場合に,同様のプログラムを書きなさい.

6. （★★★）example_02_3.cを入力して,|(1 + 1/n)ⁿ−e|の値がnに対してどのように変化するか を調べなさい.

7. （★★★）example_02_3.cでは(1 + 1/n)ⁿ を計算するために, (1 + 1/n)をn−1回乗算してい る. この部分を高速に計算できるように改良しなさい. この場合,nはもう少し大きなところまで計 算できるはずなので, n= 10^k の形（または,n= 2^k の形）でnとしてunsigned long型変数を 使ったときの限界の近くまで計算しなさい.

8. （★★★）example_02_3.cでは(1 + 1/n)ⁿ を計算するために, (1 + 1/n)をn−1回乗算してい る. この部分をpow関数を利用して同様の計算を行いなさい. この場合も前問と同程度のnまで計 算しなさい.