IEEE 784

● 講義資料

▼ 講義予定

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

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

● 講義資料

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

★ プログラム例１

(example_02_1.c)

/**

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

*/

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#define EPSILON (1.0E-15) double calc_sqrt(double x) ;

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

if (argc > 1) printf("%.15e\n", calc_sqrt(atof(argv[1]))) ; return 0 ;

}

double 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)) ; }

return square_root ; }

• このプログラムは,コマンドライン引数に数値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) double calc_sqrt(double x) ;

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

if (argc > 1) printf("%.15e\n", calc_sqrt(atof(argv[1]))) ; return 0 ;

}

double 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)) ; }

return square_root ; }

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

• この場合,x= 2.0,x= 20.0,x= 200.0,x= 5.0×10³⁰⁰,x= 5.0×10⁻³⁰⁰など,すべての代入可能

√

▼ 注意事項

• 浮動小数点数をprintf関数を表示する場合には,以下の２種類の表示方法がある – %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と近似値の値との絶対誤差を表 示している）

����������

����������

����������

����������

����������

����������

����������

����������

����������

��� ���� ����� ������ ������� �����^� �����^� �����^� �����^� �����^��

★ プログラム例

(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 = 10UL ; double e ;

for(int 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 ;

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

}

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

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

▼

IEEE 784

s eN−1 · · · e₅ e₄ e₃ e₂ e₁ e₀ 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 +

M∑−1

j=1

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

N∑−1

k=0

ek2^k, b=

N∑−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

● 実習内容

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

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

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

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

x−α|/√

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

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

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

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