回加えて も

● 講義資料

▼ サンプルプログラム

計算機内部の浮動小数演算では, 通常とは違った計算になることがある. この例では

を

回加えて も

にならないことがわかる.

【C (1)】

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

* $Id: ex_01_01.c,v 1.2 2006-04-10 13:37:00+09 naito Exp $

*

浮動小数点演算の例（倍精度浮動小数点）

* 0.1

を

回加えても

にはならない

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

#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

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

•

「単精度浮動小数点型」は

である.

•

「浮動小数点定数」

などの表記は

型定数であり,

型定数として

を書

くときには

など書く.

【C (2)】

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

* $Id: ex_01_01_1.c,v 1.3 2006-04-10 13:37:07+09 naito Exp $

*

浮動小数点演算の例（単精度浮動小数点）

* 0.1

を

回加えても

にはならない

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

#include <stdio.h>

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

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.0) 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

注意事項

•

倍精度浮動小数点演算では

を

回加算しても

よりも小さな値となる.

•

単精度浮動小数点演算では

を

回加算したら

よりも大きな値となる.

【Java】

//

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

import java.io.* ;

/**

*

浮動小数点演算の例：

を

回加えても

にはならない

* @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

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

•

「単精度浮動小数点型」は

である.

•

「浮動小数点定数」

などの表記は

型定数であり,

型定数として

を書 くときには

など書く.

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

y = 1.0000001

y does not equals to 1.0

【Fortran】

C $Id: ex_01_01.f,v 1.1 2006-04-10 21:17:19+09 naito Exp $ C

浮動小数点演算の例：

を

回加えても

にはならない

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

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

•

「倍精度浮動小数点型」は

である.

•

「浮動小数点定数」

などの表記は

型定数であり,

型定数とし て

を書くときには

など書く.

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

y = 0.9999999999999999 y does not equals to 1.0

【OCaml】

(* ex_01_01.ml *)

(* $Id: ex_01_01.ml,v 1.3 2006-04-05 18:57:54+09 naito Exp $ *) (*

浮動小数点演算の例：

を

回加えても

にはならない

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

における

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

•

「単精度浮動小数点型」は

には存在しない.

【Pascal】

{

$Id: ex_01_01.p,v 1.3 2006-04-05 17:34:26+09 naito Exp $

浮動小数点演算の例

0.1

を

回加えても

にはならない

program ex_01_01 (output);

var

x, y : real ; n : integer ; begin

x := 0.1 ; y := 0.0 ;

for n := 1 to 10 do y := y + x ; writeln (’y = ’, y) ;

if (y = 1.0) then writeln(’y equals to 1.0’) else writeln(’y does not equals to 1.0’);

end.

コンパイル・実行方法

gpc ex_01_01.p -o ex_01_01 ./ex_01_01

結果

y = 1.0000000e+00 y does not equals to 1.0

注意事項

• Pascal

では

は処理系定義の「浮動小数点型」である.

• Gnu Pascal

では

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

【PostScript】

% $Id: ex_01_01.ps,v 1.4 2006-04-10 13:50:43+09 naito Exp $

%

浮動小数点演算の例

1 1 10 {

pop X add } for

==

実行方法

gs -dNODISPLAY -dBATCH -q ex_01_01.ps

結果

1.00000012

注意事項

• PostScript

における浮動小数点数は「単精度浮動小数点型」である.

【Mathematica】

(* $Id: ex_01_01.math,v 1.2 2006-04-10 21:21:51+09 naito Exp $ *) (* 0.1

を

回加えても

とはならない

x := 0.1 ; y := 0.0 ;

For[i=0,i<10,i++,y=y+x];

Print[N[y,100]];

If[y == 1.0,

Print["Y equals to 1.0"],

Print["Y does not equals to 1.0"]

];

実行方法

math -batchinput -batchoutput < ex_01_01.math

結果

0.9999999999999999

"Y equals to 1.0"

注意事項

• Mathematica

は「無限精度」の演算を行なっているように思えるが, 数値演算に関しては, 倍精度浮

動小数点演算となっている.

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

•

どんな言語でも

を

回加算しても

と等しくはならない.

•

この計算は, 「言語」や「実行する機器」には依存せず, 浮動小数点演算の精度のみに依存した結果 となる. （正確には「丸めモード」には依存している.）

•

単精度浮動小数点演算で

を

回加算すると

よりもわずかに大きくなり, 倍精度浮動小数

点演算で

を

回加算すると

よりもわずかに小さくなる.

▼ 計算例

次のページの表はいくつかの方法を利用して

n→∞(1 + 1/n)ⁿ

を計算しようとしたものである.

なお, よく知られている通り

n→∞(1 + 1/n)ⁿ=e

であり,

の小数点以下５０桁までの真の値は,

である.

n ^{方法１ 方法２ 方法３ 方法４}

1 2.0000000000 2.0000000000 2.000000000000000000000000 2.000000000000000000000000 5 2.4883201122 2.4883203506 2.488319999999999865281097 2.488319999999999421191887 10 2.5937430859 2.5937428474 2.593742460100001867573383 2.593742460100002311662593 50 2.6915805340 2.6915857792 2.691588029073604726448821 2.691588029073607835073290 100 2.7048130035 2.7048108578 2.704813829421528925678331 2.704813829421528481589121 500 2.7155508995 2.7155337334 2.715568520651697959067405 2.715568520651728157133675 1000 2.7170419693 2.7170493603 2.716923932235520311451182 2.716923932235593586170808 5000 2.7180233002 2.7184610367 2.718010050101671737365905 2.718010050101554941903714 10000 2.7180233002 2.7185957432 2.718145926824356184425824 2.718145926824925506792852 50000 2.7215235233 2.7219350338 2.718254646133538621199932 2.718254646126732065880560 100000 2.7215235233 2.7219107151 2.718268237197528414128556 2.718268237192297487325732 500000 2.7544298172 2.7532434464 2.718279110219865568609521 2.718279110260366060458637 1000000 2.5946049690 2.5898523331 2.718280469156427514576535 2.718280469095753382191560 5000000 3.2929797173 3.2012779713 2.718281555159559825796123 2.718281555200129151472765 1.0e+07 3.2929797173 2.8841857910 2.718281693980372448748994 2.718281694132081760528763 5.0e+07 1.0000000000 1.0000000000 2.718281816570668141253009 2.718281814934938811489928 1.0e+08 2.718281786395797539057639 2.718281798347357725020856 5.0e+08 2.718281729032110849431092 2.718281748862504176855737 1.0e+09 2.718282030814509475646901 2.718282052011560256943312 5.0e+09 2.718282043475248155317558 2.718282053098873163321514 1.0e+10 2.718282043475248155317558 2.718282053234787554174545 5.0e+10 2.718282030116025538291069 2.718282053343518800403444 1.0e+11 2.718282030116025538291069 2.718282053357110150670906 5.0e+11 2.718221668654040357182566 2.718221696051625002610308 1.0e+12 2.718523469568279615771189 2.718523496037237752176452 5.0e+12 2.719127181674036553005180 2.719127196535721324721635 1.0e+13 2.716110017431422107847538 2.716110034086900881789006 5.0e+13 2.716110008485271620770618 2.716110034087009239556210 1.0e+14 2.716110008485271620770618 2.716110034087022562232505 5.0e+14 2.716110004857417603574277 2.716110034087033664462751 1.0e+15 3.035035181396292358613209 3.035035206549261843633758 5.0e+15 3.035035181414545313316466 3.035035206549263619990597 1.0e+16 1.000000000000000000000000 1.000000000000000000000000

てなわけで, プログラムを組んで調べてみると

n→∞lim(1 + 1/n)ⁿ = 1

であることがわかる. あれぇ？

方法１ 演算精度：単精度浮動小数点演算

乗算方法：高速乗算法

方法２ 演算精度：単精度浮動小数点演算

乗算方法：n 回の積

方法３ 演算精度：倍精度浮動小数点演算

乗算方法：高速乗算法

方法４ 演算精度：倍精度浮動小数点演算

乗算方法：pow 関数

各方法について, (1 + 1

の値および

をプロットしたもの

1 1.5 2 2.5 3 3.5

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

"method1"

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

方法１

1 1.5 2 2.5 3 3.5

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

"method2"

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

方法２

1 1.5 2 2.5 3 3.5

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

"method3"

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

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

方法３

1.5 2 2.5 3 3.5

"method4"

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

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

● 実習内容

• 1.0

を

回加算したら, その結果は

と等しくなるかを調べなさい. また, その理由はなぜかを 考察しなさい.

• lim

n→∞(1 + 1/n)ⁿ

を計算するために,

を上の表のように動かして, (1 + 1

の値を出力するプログ ラムを書き, 得られたデータから前ページのようなグラフを

を用いて書きなさい. 乗算の方 法としては「高速乗算法」を用いなさい. また, このような結果となる理由を考察しなさい.

なお, 大学院生が利用する

上では

は

バイトであるが, 学部生が利用する

上では

バイトである. Linux 上では

が

バイト整数である.

• lim

n→∞n^1/n

や

x→0sin(x)/x

は, 上記のような方法で正しく計算できるかを考えなさい. また, その理 由を考察しなさい.

•

電子メールで「今日の講義の感想や意見」を送ってください.

【注意１】

において

などの「数学関数ライブラリ」を用いる場合には,

gcc -Wall -ansi -pedantic foo.c -lm

とする. （-lm がついていることに注意.）

【注意２】 この授業の実習・課題やレポートの回答では, 「IEEE754 浮動小数点演算標準」に準拠して

いる言語であれば, 何を使っても構わないこととする. （上記の

はこれに該当する.）