● 講義資料
▼ サンプルプログラム
計算機内部の浮動小数演算では,我々の通常の感覚とは異った計算になることがある. この例では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である.
【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よりもわずかに小さくなる.
▼ いくつかの数値の計算
ここでは,今後の授業のイントロダクションとして,いくつかの数値を簡単に（手抜きで）計算するプロ グラムを紹介する. これらの計算方法およびその問題点は,今後の授業の中で解説する.
★ 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 4.142135623730951e-01 1.500000000000000 8.578643762690485e-02 1.250000000000000 1.642135623730951e-01 1.375000000000000 3.921356237309515e-02 1.437500000000000 2.328643762690485e-02 1.406250000000000 7.963562373095145e-03 1.421875000000000 7.661437626904855e-03 1.414062500000000 1.510623730951455e-04 1.417968750000000 3.755187626904855e-03 1.416015625000000 1.802062626904855e-03 1.415039062500000 8.255001269048545e-04 1.414550781250000 3.372188769048545e-04 1.414306640625000 9.307825190485453e-05 1.414184570312500 2.899206059514547e-05 1.414245605468750 3.204309565485453e-05 1.414215087890625 1.525517529854525e-06
【ニュートン法】
/*
* 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 8.578643762690485e-02 1.416666666666667 2.453104293571373e-03 1.414215686274510 2.123901414519125e-06 1.414213562374690 1.594724352571575e-12 1.414213562373095 2.220446049250313e-16 1.414213562373095 2.220446049250313e-16 1.414213562373095 2.220446049250313e-16
★ 自然対数の底の値の計算
【定義をそのまま計算】
/*
* 自然対数の底 $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 2.593742460100002 1.245393683590428e-01 100 2.704813829421529 1.346799903751572e-02 1000 2.716923932235598 1.357896223446620e-03 10000 2.718145926824898 1.359016341466734e-04 100000 2.718268237192295 1.359126675026801e-05 1000000 2.718280469095936 1.359363108743850e-06
10000000 2.718281694132010 1.343270348286296e-07
【テイラー展開を利用する】
/*
* 自然対数の底 $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 2.718281525573192 3.028858528431044e-07 20 2.718281828459046 4.440892098500626e-16 30 2.718281828459046 4.440892098500626e-16
★ 円周率の値の計算
【内接正多角形の長さを利用】
/*
* 円周率 $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 3.000000000000000 1.415926535897931e-01 12 3.105828541230250 3.576411235954335e-02 24 3.132628613281237 8.964040308556243e-03 48 3.139350203046872 2.242450542921048e-03 96 3.141031950890530 5.607026992633379e-04 192 3.141452472285344 1.401813044488165e-04 384 3.141557607911622 3.504567817103066e-05 768 3.141583892148936 8.761440857263381e-06 1536 3.141590463236762 2.190353031394920e-06 3072 3.141592106043048 5.475467448334825e-07 6144 3.141592516588155 1.370016384782957e-07 12288 3.141592618640789 3.494900369105380e-08 24576 3.141592645321216 8.268577378345299e-09 49152 3.141592645321216 8.268577378345299e-09 98304 3.141592645321216 8.268577378345299e-09 196608 3.141592645321216 8.268577378345299e-09
【テイラー展開を利用する】
/*
* 円周率 $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 3.339682539682540 1.980898860927471e-01 100 3.121594652591011 1.999800099878213e-02 1000 3.139592655589785 1.999998000008052e-03 10000 3.141392653591791 1.999999980020206e-04 100000 3.141572653589781 2.000000001167734e-05 1000000 3.141590653589692 2.000000101087807e-06
10000000 3.141592453589780 2.000000134394497e-07
【注意０】 情報メディア教育センターの 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 arc cosine function
acosh inverse hyperbolic cosine function asin arc sine function
asinh hyperbolic sine function atan arc tangent function
atanh inverse hyperbolic tangent function ceil ceiling value function
cos cosine function
cosh hyperbolic cosine function exp exponential function fabs absolute value function floor floor function
log natural logarithm function
pow power function
sin sine function
sinh hyperbolic sine function sqrt square root function
tan tangent function
tanh hyperbolic tangent function
【注意２】 math.hの中では,以下の数学定数の値が浮動小数点定数として定義されている.
M_E The base of natural logarithms (e).
M_LOG2E The base-2 logarithm of e.
M_PI pi, the ratio of the circumference of a circle to its diameter.
M_PI_2 pi/2.
M_PI_4 pi/4.
M_1_PI 1/pi.
M_2_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などはこれに該当す る.)
