● 講義資料
▼ サンプルプログラム
計算機内部の浮動小数演算では,通常とは違った計算になることがある. この例では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である.
• 「浮動小数点定数」0.1,1.0E-1などの表記はdouble型定数であり,float型定数として 0.1を書
【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>
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.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よりも大きな値となる.
【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
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 浮動小数点演算の例： 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
【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には存在しない.
▼ サンプルプログラムのまとめ
• どんな言語でも0.1を10回加算しても1.0と等しくはならない.
• この計算は, 「言語」や「実行する機器」には依存せず,浮動小数点演算の精度のみに依存した結果 となる. （正確には「丸めモード」には依存している.）
• 単精度浮動小数点演算で 0.1を10 回加算すると1.0よりもわずかに大きくなり, 倍精度浮動小数 点演算で 0.1を 10回加算すると1.0よりもわずかに小さくなる.
▼ 計算例
次のページの表はいくつかの方法を利用して 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 てなわけで,プログラムを組んで調べてみると
n→∞lim(1 + 1/n)n = 1
であることがわかる. あれぇ？
方法１ 演算精度：単精度浮動小数点演算(float),乗算方法：高速乗算法 方法２ 演算精度：倍精度浮動小数点演算(double),乗算方法：高速乗算法 方法３ 演算精度：倍精度浮動小数点演算(double),乗算方法：pow関数
各方法について,|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
単精度 倍精度
★ 方法２のプログラム(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 main(int argc, char **argv) {
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 ; }
★ 方法２のプログラム(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 main(int argc, char **argv) {
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 ; }
● 実習内容
• 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→∞n1/n や lim
x→0sin(x)/x は, 上記のような方法で正しく計算できるかを考えなさい. また, その理 由を考察しなさい.
【注意１】 Cにおいて powなどの「数学関数ライブラリ」を用いる場合には,
gcc foo.c -lm
とする. （-lmがついていることに注意.）
【注意２】 この授業の実習・課題やレポートの回答では,「IEEE754浮動小数点演算標準」に準拠してい る言語であれば,何を使っても構わないこととする. （上記のC, Java, Fortran, OCamlなどはこれ に該当する.）
