実習メモ

２００８年度・前期・数理解析・計算機数学２・第１回
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 である.
• 「浮動小数点定数」 0.1, 1.0E-1 などの表記は double 型定数であり, float 型定数として 0.1 を書
くときには 0.1F, 1.0E-1F など書く.
ex01.tex,v 1.5 2008-04-14 18:00:24+09 naito Exp
２００８年度・前期・数理解析・計算機数学２・第１回
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>
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 よりも大きな値となる.
▼ サンプルプログラムのまとめ
• 0.1 を 10 回加算しても 1.0 と等しくはならない.
• この計算は, 「言語」や「実行する機器」には依存せず, 浮動小数点演算の精度のみに依存した結果
となる. （正確には「丸めモード」には依存している.）
• 単精度浮動小数点演算で 0.1 を 10 回加算すると 1.0 よりもわずかに大きくなり, 倍精度浮動小数
点演算で 0.1 を 10 回加算すると 1.0 よりもわずかに小さくなる.
ex01.tex,v 1.5 2008-04-14 18:00:24+09 naito Exp
２００８年度・前期・数理解析・計算機数学２・第１回
3
▼ いくつかの数値の計算
ここでは, 今後の授業のイントロダクションとして, いくつかの数値を簡単に（手抜きで）計算するプロ
グラムを紹介する. これらの計算方法およびその問題点は, 今後の授業の中で解説する.
★ 2 の平方根の値の計算
【区間縮小法】
/*
* SQRT[2] の値を近似的に求める.
方法１: 区間縮小法
*/
#include <stdio.h>
#include <math.h>
#define N
(20)
int main(int argc, char **argv)
{
double left, right, sqrt2 ;
int i ;
left
= 0.0 ;
right = 2.0 ;
for(i=0;i<N;i++) {
sqrt2 = (left + right)/2.0 ;
printf("%.15f\t%.15e\n", sqrt2, fabs(sqrt2-M_SQRT2)) ;
if (sqrt2*sqrt2 <= 2.0) left = sqrt2 ;
else right = sqrt2 ;
}
return 0 ;
}
1.000000000000000
1.500000000000000
1.250000000000000
1.375000000000000
1.437500000000000
1.406250000000000
1.421875000000000
1.414062500000000
1.417968750000000
1.416015625000000
1.415039062500000
1.414550781250000
1.414306640625000
1.414184570312500
1.414245605468750
1.414215087890625
1.414199829101562
1.414207458496094
1.414211273193359
1.414213180541992
4.142135623730951e-01
8.578643762690485e-02
1.642135623730951e-01
3.921356237309515e-02
2.328643762690485e-02
7.963562373095145e-03
7.661437626904855e-03
1.510623730951455e-04
3.755187626904855e-03
1.802062626904855e-03
8.255001269048545e-04
3.372188769048545e-04
9.307825190485453e-05
2.899206059514547e-05
3.204309565485453e-05
1.525517529854525e-06
1.373327153264547e-05
6.103877001395475e-06
2.289179735770475e-06
3.818311029579746e-07
ex01.tex,v 1.5 2008-04-14 18:00:24+09 naito Exp
２００８年度・前期・数理解析・計算機数学２・第１回
4
【ニュートン法】
/*
* SQRT[2] の値を近似的に求める. 方法２: ニュートン法
*/
#include <stdio.h>
#include <math.h>
#define N
(8)
int main(int argc, char **argv)
{
double sqrt2 ;
int i ;
sqrt2 = 2.0 ;
for(i=0;i<N;i++) {
printf("%.15f\t%.15e\n", sqrt2, fabs(sqrt2-M_SQRT2)) ;
sqrt2 = sqrt2/2.0 + 1.0/sqrt2 ;
}
return 0 ;
}
2.000000000000000
5.857864376269049e-01
1.500000000000000
1.416666666666667
8.578643762690485e-02
2.453104293571373e-03
1.414215686274510
1.414213562374690
2.123901414519125e-06
1.594724352571575e-12
1.414213562373095
2.220446049250313e-16
1.414213562373095
1.414213562373095
2.220446049250313e-16
2.220446049250313e-16
ex01.tex,v 1.5 2008-04-14 18:00:24+09 naito Exp
２００８年度・前期・数理解析・計算機数学２・第１回
5
★ 自然対数の底の値の計算
【定義をそのまま計算】
/*
* 自然対数の底 $e$ の値を近似的に求める. 方法１の１: (1+1/n)^n の値を計算する. ひたす
らかけ算する.
*/
#include <stdio.h>
#include <math.h>
double calc_e(unsigned long n) ;
unsigned long n[] = {
10UL, 100UL, 1000UL, 10000UL, 100000UL, 1000000UL, 10000000UL,
100000000UL, 1000000000UL, 10000000000UL, 100000000000UL
} ;
int main(int argc, char **argv)
{
int i, m ;
double e ;
m = sizeof(n)/sizeof(unsigned long) ;
for(i=0;i<m;i++) {
e = calc_e(n[i]) ;
printf("%lu\t%.15f\t%.15e\n", n[i], e, fabs(e-M_E)) ;
}
return 0 ;
}
double calc_e(unsigned long n) {
double x = 1.0 ;
double y ;
unsigned long int i ;
y = 1.0 + 1.0/n ;
for(i=0UL;i<n;i++) x *= y ;
return x ;
}
10
100
2.593742460100002
2.704813829421529
1.245393683590428e-01
1.346799903751572e-02
2.716923932235598
2.718145926824898
1.357896223446620e-03
1.359016341466734e-04
100000 2.718268237192295
1000000 2.718280469095936
1.359126675026801e-05
1.359363108743850e-06
1000
10000
10000000
2.718281694132010
1.343270348286296e-07
100000000
1000000000
2.718281798346360
2.718282052011900
3.011268523422928e-08
2.235528544503040e-07
10000000000
100000000000
2.718282053233160
2.718282053194587
2.247741144323356e-07
2.247355421758357e-07
ex01.tex,v 1.5 2008-04-14 18:00:24+09 naito Exp
２００８年度・前期・数理解析・計算機数学２・第１回
6
【テイラー展開を利用する】
/*
* 自然対数の底 $e$ の値を近似的に求める. 方法２の１: exp(x) のテイラー展開を計算する.
*/
#include <stdio.h>
#include <math.h>
unsigned long n[] = {
5UL, 10UL, 20UL, 30UL,
} ;
double calc_e(unsigned long n) ;
int main(int argc, char **argv)
{
int i, m ;
double e ;
m = sizeof(n)/sizeof(unsigned long) ;
for(i=0;i<m;i++) {
e = calc_e(n[i]) ;
printf("%lu\t%.15f\t%.15e\n", n[i], e, fabs(e-M_E)) ;
}
return 0 ;
}
double calc_e(unsigned long n) {
double x = 1.0 ;
double y = 1.0 ;
unsigned long int i ;
for(i=1UL;i<n;i++) {
y *= 1.0/(double)i ;
x += y ;
}
return x ;
}
5
2.708333333333333
9.948495125712054e-03
10
20
2.718281525573192
2.718281828459046
3.028858528431044e-07
4.440892098500626e-16
30
2.718281828459046
4.440892098500626e-16
ex01.tex,v 1.5 2008-04-14 18:00:24+09 naito Exp
２００８年度・前期・数理解析・計算機数学２・第１回
7
★ 円周率の値の計算
【内接正多角形の長さを利用】
/*
* 円周率 $pi$ の値を近似的に求める. 方法３の１: 多角形近似. 単なる倍角の公式.
*/
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#define N
(20)
unsigned long pi_montecalro(unsigned long n) ;
int main(int argc, char **argv)
{
int i ;
unsigned long n ;
unsigned int p ;
double pi ;
double inner ;
inner = 3.0 ;
for(i=0;i<N;i++) {
n = (1UL<<i)*6 ;
printf("%lu\t%.15f\t%.15e\n", (1UL<<i)*6, inner, fabs(inner-M_PI)) ;
inner = 2*n*sqrt((1.0 - sqrt(1.0-(inner/n)*(inner/n)))/2.0) ;
}
return 0 ;
}
6
12
24
48
96
192
384
768
1536
3072
6144
12288
24576
49152
98304
196608
393216
786432
1572864
3145728
3.000000000000000
3.105828541230250
3.132628613281237
3.139350203046872
3.141031950890530
3.141452472285344
3.141557607911622
3.141583892148936
3.141590463236762
3.141592106043048
3.141592516588155
3.141592618640789
3.141592645321216
3.141592645321216
3.141592645321216
3.141592645321216
3.141593669849427
3.141592303811738
3.141608696224804
3.141586839655041
1.415926535897931e-01
3.576411235954335e-02
8.964040308556243e-03
2.242450542921048e-03
5.607026992633379e-04
1.401813044488165e-04
3.504567817103066e-05
8.761440857263381e-06
2.190353031394920e-06
5.475467448334825e-07
1.370016384782957e-07
3.494900369105380e-08
8.268577378345299e-09
8.268577378345299e-09
8.268577378345299e-09
8.268577378345299e-09
1.016259633779271e-06
3.497780554084784e-07
1.604263501064906e-05
5.813934751852656e-06
ex01.tex,v 1.5 2008-04-14 18:00:24+09 naito Exp
２００８年度・前期・数理解析・計算機数学２・第１回
8
【テイラー展開を利用する】
/*
* 円周率 $pi$ の値を近似的に求める. 方法１の１: arctan(x) のテイラー展開を x = 1 で求
める.
*/
#include <stdio.h>
#include <math.h>
unsigned long n[] = {
10UL, 100UL, 1000UL, 10000UL, 100000UL, 1000000UL, 10000000UL
} ;
double atan_taylor(unsigned long n, double x) ;
int main(int argc, char **argv)
{
int i, m ;
double pi ;
m = sizeof(n)/sizeof(unsigned long) ;
for(i=0;i<m;i++) {
pi = atan_taylor(n[i], 1.0)*4.0 ;
printf("%lu\t%.15f\t%.15e\n", n[i], pi, fabs(pi-M_PI)) ;
}
return 0 ;
}
double atan_taylor(unsigned long n, double x) {
double v = 0.0 ;
double y ;
unsigned long int i ;
int sign = 1 ;
y = x ;
for(i=1UL;i<n;i+=2) {
if (sign) v += (1.0/(double)i) * y ;
else v -= (1.0/(double)i) * y ;
y *= x*x ;
sign ^= 1 ;
}
return v ;
}
10
100
3.339682539682540
3.121594652591011
1.980898860927471e-01
1.999800099878213e-02
1000
10000
3.139592655589785
3.141392653591791
1.999998000008052e-03
1.999999980020206e-04
100000
3.141572653589781
2.000000001167734e-05
1000000 3.141590653589692
2.000000101087807e-06
3.141592453589780
2.000000134394497e-07
10000000
ex01.tex,v 1.5 2008-04-14 18:00:24+09 naito Exp
２００８年度・前期・数理解析・計算機数学２・第１回
9
【注意０】 情報メディア教育センターの MacOSX では, C コンパイラ利用時には -m64 オプションをつ
けると, long, unsigned long が 8 バイトとなります.
gcc -m64 foo.c
とする. （-m64 がついていることに注意.）
【注意１】 C において pow などの「数学関数ライブラリ」を用いる場合には,
gcc foo.c -lm
または
gcc -m64 foo.c -lm
とする. （-lm がついていることに注意.）代表的な数学関数は以下の通り.
acos
acosh
arc cosine function
inverse hyperbolic cosine function
asin
asinh
arc sine function
hyperbolic sine function
atan
arc tangent function
atanh
ceil
inverse hyperbolic tangent function
ceiling value function
cos
cosh
cosine function
hyperbolic cosine function
exp
fabs
exponential function
absolute value function
floor
log
floor function
natural logarithm function
pow
sin
power function
sine function
sinh
sqrt
hyperbolic sine function
square root function
tan
tanh
tangent function
hyperbolic tangent function
【注意２】 math.h の中では, 以下の数学定数の値が浮動小数点定数として定義されている.
M_E
M_LOG2E
The base of natural logarithms (e).
The base-2 logarithm of e.
M_LOG10E
M_LN2
The base-10 logarithm of e.
The natural logarithm of 2.
M_LN10
The natural logarithm of 10.
ex01.tex,v 1.5 2008-04-14 18:00:24+09 naito Exp
２００８年度・前期・数理解析・計算機数学２・第１回
10
M_PI
pi, the ratio of the circumference of a circle to its diameter.
M_PI_2
M_PI_4
pi/2.
pi/4.
M_1_PI
M_2_PI
1/pi.
2/pi.
M_2_SQRTPI 2 over the square root of pi.
M_SQRT2
The positive square root of 2.
M_SQRT1_2
The positive square root of 1/2.
【注意３】 この授業の実習・課題やレポートの回答では, 「IEEE754 浮動小数点演算標準」に準拠してい
る言語であれば, 何を使っても構わないこととする. (C, Java, Fortran, OCaml などはこれに該当す
る.)
ex01.tex,v 1.5 2008-04-14 18:00:24+09 naito Exp