ex01.dvi

(1)

● 講義資料

▼ サンプルプログラム

計算機内部の浮動小数演算では, 通常とは違った計算になることがある. この例では 0.1 を 10 回加えても 1.0 にならないことがわかる. 【C (1)】 /*******************************

* $Id: ex_01_01.c,v 1.2 2006-04-10 13:37:00+09 naito Exp $ * 浮動小数点演算の例（倍精度浮動小数点）

* 0.1 を 10 回加えても 1.0 にはならない *******************************/ #include <stdio.h>

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

double x = 0.1, y = 0.0 ;

int n ;

for(n=0;n<10;n++) y += x ; printf("y = %f\n", y) ;

if (y == 1.0) printf ("y equals to 1.0\n") ; else printf ("y does not equals to 1.0\n") ; printf("y = %.16f\n", y) ; return 0 ; } コンパイル・実行方法 gcc ex_01_01.c -o ex_01_01 ./ex_01_01 結果 y = 1.000000

y does not equals to 1.0 y = 0.9999999999999999 注意事項

• double は「倍精度浮動小数点型」である. • 「単精度浮動小数点型」は float である.

(2)

【C (2)】

/*******************************

* $Id: ex_01_01_1.c,v 1.3 2006-04-10 13:37:07+09 naito Exp $ * 浮動小数点演算の例（単精度浮動小数点）

* 0.1 を 10 回加えても 1.0 にはならない *******************************/ #include <stdio.h>

float x = 0.1F, y = 0.0F ;

int n ;

for(n=0;n<10;n++) y += x ; printf("y = %f\n", y) ;

if (y == 1.0F) printf ("y equals to 1.0\n") ; else printf ("y does not equals to 1.0\n") ; printf("y = %.8f\n", y) ; return 0 ; } 結果 y = 1.000000

y does not equals to 1.0 y = 1.00000012

注意事項

• 倍精度浮動小数点演算では 0.1 を 10 回加算しても 1.0 よりも小さな値となる. • 単精度浮動小数点演算では 0.1 を 10 回加算したら 1.0 よりも大きな値となる.

(3)

【Java】

//

// $Id: ex_01_01.java,v 1.2 2006-04-05 16:38:59+09 naito Exp $ //

import java.io.* ;

/**

* 浮動小数点演算の例： 0.1 を 10 回加えても 1.0 にはならない * @author Hisashi Naito

*/

public class ex_01_01 {

public static void main (String args[]) throws IOException { double x = 0.1, y = 0.0 ;

int n ;

for(n=0;n<10;n++) y += x ; System.out.println("y = " + y) ;

if (y == 1.0) System.out.println("y equals to 1.0") ; else System.out.println("y does not equals to 1.0") ; } } コンパイル・実行方法 javac ex_01_01.java java ex_01_01 結果 y = 0.9999999999999999 y does not equals to 1.0 注意事項

• double は「倍精度浮動小数点型」である. • 「単精度浮動小数点型」は float である.

• 「浮動小数点定数」 0.1, 1.0E-1 などの表記は double 型定数であり, float 型定数として 0.1 を書 くときには 0.1F, 1.0E-1F など書く.

結果（単精度浮動小数点演算）

y = 1.0000001

(4)

【Fortran】

C $Id: ex_01_01.f,v 1.1 2006-04-10 21:17:19+09 naito Exp $ C 浮動小数点演算の例： 0.1 を 10 回加えても 1.0 にはならない PROGRAM EX_01_01 REAL X, Y INTEGER N X = 0.1 ; Y = 0.0 ; DO N=1,10,1 Y = Y + X ENDO WRITE(*,601) ’Y = ’, Y IF (Y .EQ. 1.0) THEN

PRINT *, ’Y EQUALS TO 1.0’ ELSE

PRINT *, ’Y DOES NOT EQUALS TO 1.0’ ENDIF 601 FORMAT(A,F10.8) STOP END コンパイル・実行方法 g77 ex_01_01.f ./ex_01_01 結果 y = 1.00000012

y does not equals to 1.0 注意事項

• real は「単精度浮動小数点型」である.

• 「倍精度浮動小数点型」は double precision である.

• 「浮動小数点定数」 0.1, 1.0E-1 などの表記は real 型定数であり, double precision 型定数とし て 0.1 を書くときには 0.1D0, 1.0D-1 など書く.

結果（倍精度浮動小数点演算）

y = 0.9999999999999999 y does not equals to 1.0

(5)

【OCaml】

(* ex_01_01.ml *)

(* $Id: ex_01_01.ml,v 1.3 2006-04-05 18:57:54+09 naito Exp $ *) (* 浮動小数点演算の例： 0.1 を 10 回加えても 1.0 にはならない *)

let ex_01_01 () =

let x = 0.1 and y = ref 0.0 in

for n = 1 to 10 do y := x +. !y done ; Printf.printf "%f\n" !y ;

if (!y == 1.0) then Printf.printf "y equals to 1.0\n" else Printf.printf "y does not equals to 1.0\n"; Printf.printf "%.16f\n" !y ; ;; ex_01_01() ;; コンパイル・実行方法 ocamlc ex_01_01.ml

ocamlc ex_01_01.cmo -o ./ex_01_01 ./ex_01_01

結果

1.000000

y does not equals to 1.0 0.9999999999999999 注意事項

• OCaml における float は「倍精度浮動小数点型」である. • 「単精度浮動小数点型」は OCaml には存在しない.

(6)

▼ サンプルプログラムのまとめ

• どんな言語でも 0.1 を 10 回加算しても 1.0 と等しくはならない. • この計算は, 「言語」や「実行する機器」には依存せず, 浮動小数点演算の精度のみに依存した結果 となる. （正確には「丸めモード」には依存している.） • 単精度浮動小数点演算で 0.1 を 10 回加算すると 1.0 よりもわずかに大きくなり, 倍精度浮動小数 点演算で 0.1 を 10 回加算すると 1.0 よりもわずかに小さくなる.

(7)

▼ 計算例

次のページの表はいくつかの方法を利用して lim n→∞(1 + 1/n) n _{を計算しようとしたものである.} なお, よく知られている通り lim n→∞(1 + 1/n) n_{= e であり, e の小数点以下５０桁までの真の値は,} 2.71828182845904523536028747135266249775724709369995 である. n 方法１方法２方法３ 1 2.0000000000 2.0000000000000000 2.0000000000000000 5 2.4883201122 2.4883199999999999 2.4883199999999994 10 2.5937430859 2.5937424601000019 2.5937424601000023 50 2.6915805340 2.6915880290736047 2.6915880290736078 100 2.7048130035 2.7048138294215289 2.7048138294215285 500 2.7155508995 2.7155685206516980 2.7155685206517282 1000 2.7170419693 2.7169239322355203 2.7169239322355936 5000 2.7180233002 2.7180100501016717 2.7180100501015549 10000 2.7180233002 2.7181459268243562 2.7181459268249255 50000 2.7215235233 2.7182546461335386 2.7182546461267321 100000 2.7215235233 2.7182682371975284 2.7182682371922975 500000 2.7544298172 2.7182791102198656 2.7182791102603661 1000000 2.5946049690 2.7182804691564275 2.7182804690957534 5000000 3.2929797173 2.7182815551595598 2.7182815552001292 1.0e+07 3.2929797173 2.7182816939803724 2.7182816941320818 5.0e+07 1.0000000000 2.7182818165706681 2.7182818149349388 1.0e+08 2.7182817863957975 2.7182817983473577 5.0e+08 2.7182817290321108 2.7182817488625042 1.0e+09 2.7182820308145095 2.7182820520115603 5.0e+09 2.7182820434752482 2.7182820530988732 1.0e+10 2.7182820434752482 2.7182820532347876 5.0e+10 2.7182820301160255 2.7182820533435188 1.0e+11 2.7182820301160255 2.7182820533571102 5.0e+11 2.7182216686540404 2.7182216960516250 1.0e+12 2.7185234695682796 2.7185234960372378 5.0e+12 2.7191271816740366 2.7191271965357213 1.0e+13 2.7161100174314221 2.7161100340869009 5.0e+13 2.7161100084852716 2.7161100340870092 1.0e+14 2.7161100084852716 2.7161100340870226 5.0e+14 2.7161100048574176 2.7161100340870337 1.0e+15 3.0350351813962924 3.0350352065492618 5.0e+15 3.0350351814145453 3.0350352065492636 1.0e+16 1.0000000000000000 1.0000000000000000 てなわけで, プログラムを組んで調べてみると lim n→∞(1 + 1/n) n _{= 1} であることがわかる. あれぇ？方法１演算精度：単精度浮動小数点演算 (float), 乗算方法：高速乗算法方法２演算精度：倍精度浮動小数点演算 (double), 乗算方法：高速乗算法方法３演算精度：倍精度浮動小数点演算 (double), 乗算方法：pow 関数

(8)

各方法について,|e − (1 + 1/n)n| をプロットしたもの 1e-09 1e-08 1e-07 1e-06 1e-05 0.0001 0.001 0.01 0.1 1 10

1 100 10000 1e+06 1e+08 1e+10 1e+12 1e+14 1e+16 Method 2 1e-09 1e-08 1e-07 1e-06 1e-05 0.0001 0.001 0.01 0.1 1 10

1 100 10000 1e+06 1e+08 1e+10 1e+12 1e+14 1e+16 Method 2 方法１方法２ 1e-09 1e-08 1e-07 1e-06 1e-05 0.0001 0.001 0.01 0.1 1 10

1 100 10000 1e+06 1e+08 1e+10 1e+12 1e+14 1e+16 Method 3 0 1 2 3 4 5

1 100 10000 1e+06 1e+08 1e+10 1e+12 1e+14 1e+16

方法３単精度及び倍精度の方法について, 誤差の評価直線を追加したもの 1e-09 1e-08 1e-07 1e-06 1e-05 0.0001 0.001 0.01 0.1 1 10

1 100 10000 1e+06 1e+08 1e+10 1e+12 1e+14 1e+16

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

1 100 10000 1e+06 1e+08 1e+10 1e+12 1e+14 1e+16

(9)

★ 方法２のプログラム(Solaris 用) /* sample_01_2.c * (1+1/n)^n を計算する * 倍精度浮動小数点演算 */ #include <stdio.h>

double power(double, unsigned long) ;

unsigned long print_value[] = { 1, 5, 10, 50, 100, 500, 1000, 5000, 10000, 50000, 100000, 500000, 1000000, 5000000, 10000000, 50000000, 100000000, 500000000, 1000000000, 5000000000, 10000000000, 50000000000, 100000000000, 500000000000, 1000000000000, 5000000000000, 10000000000000, 50000000000000, 100000000000000, 500000000000000, 1000000000000000, 5000000000000000, 10000000000000000, 50000000000000000 } ;

int i ;

int limit = sizeof(print_value)/sizeof(print_value[0]) ;

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

fprintf(stdout, "%-16lu\t%18.16f\n", print_value[i],

power((1.0+1.0/(double)print_value[i]), print_value[i])) ; }

return 0 ; }

double power(double x, unsigned long n) { double x0 = 1.0, x1 = x ; while(n) { if (n&1) x0 *= x1 ; x1 *= x1 ; n >>= 1 ; } return x0 ; }

(10)

★ 方法２のプログラム(Linux 用) /* sample_01_2.c * (1+1/n)^n を計算する * 倍精度浮動小数点演算 */ #include <stdio.h>

double power(double, unsigned long long) ;

unsigned long long print_value[] = { 1, 5, 10, 50, 100, 500, 1000, 5000, 10000, 50000, 100000, 500000, 1000000, 5000000, 10000000, 50000000, 100000000, 500000000, 1000000000, 5000000000, 10000000000, 50000000000, 100000000000, 500000000000, 1000000000000, 5000000000000, 10000000000000, 50000000000000, 100000000000000, 500000000000000, 1000000000000000, 5000000000000000, 10000000000000000, 50000000000000000 } ;

int i ;

int limit = sizeof(print_value)/sizeof(print_value[0]) ;

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

fprintf(stdout, "%-16lu\t%18.16f\n", print_value[i],

power((1.0+1.0/(double)print_value[i]), print_value[i])) ; }

return 0 ; }

double power(double x, unsigned long long n) { double x0 = 1.0, x1 = x ; while(n) { if (n&1) x0 *= x1 ; x1 *= x1 ; n >>= 1 ; } return x0 ; }

(11)

● 実習内容

• 1.0 を 10 回加算したら, その結果は 10.0 と等しくなるかを調べなさい. また, その理由はなぜかを 考察しなさい. • lim n→∞(1 + 1/n) n _{を計算するために,}_{n を上の表のように動かして, (1 + 1/n)}n _{の値を出力するプログ} ラムを書き, 得られたデータから前ページのようなグラフを gnuplot を用いて書きなさい. 乗算の方法としては「高速乗算法」を用いなさい. また, このような結果となる理由を考察しなさい.

なお, 大学院生が利用する Solaris 上では unsigned long は 8 バイトであるが, 学部生が利用する Linux 上では 4 バイトである. Linux 上では unsigned long long が 8 バイト整数である. • lim n→∞n 1/n _{や lim} x→0sin(x)/x は, 上記のような方法で正しく計算できるかを考えなさい. また, その理 由を考察しなさい. 【注意１】 C において pow などの「数学関数ライブラリ」を用いる場合には, gcc foo.c -lm とする. （-lm がついていることに注意.）【注意２】この授業の実習・課題やレポートの回答では, 「IEEE754 浮動小数点演算標準」に準拠している言語であれば, 何を使っても構わないこととする. （上記の C, Java, Fortran, OCaml などはこれに該当する.）