数値計算のねらいは，できるだけ精確・高速に解を得ることである．誤差 (精度) と収束性 (安定性,速度) が数値計算のキモとなる．前回に説明した収束判定条件による誤差は打ち切り誤差 (truncation error) と呼ばれる．ここでは，誤差のもう一つの代表例である，計算機に特有の丸め誤差 (roundoff error) について見ておこう．
計算機は一般に無限精度の計算をおこなっているわけではない．CPU で足し算をおこなう以上，一般的な計算においては CPUが扱う一番効率のいい数の大きさが存在する．これが,32bit の CPU では 1 ワード,4byte(4x8bits) である．1ワードで表現できる最大整数は,符号に 1bit 必要なので,2^(31)-1 となる．実数は以下のような仮数部と指数を取る浮動小数点数で表わされる．
| 記号 | $s \times f \times B^{e-E}$ |
|---|---|
| $s$ | sign bit(符号ビット:正負の区別を表す) |
| $f$ | fraction portion of the number(仮数部) |
| $e$ | biased exponent(指数部) |
| $B$ | base(基底) で通常は 2 |
| $E$ | bias(下駄)と呼ばれる |
| real(単精度) | $s=1, e=8, f=23, E=127$ |
| double precision(倍精度) | $s=1, e=11, f=52, E=1023$ |
$E$は指数が負となる小数点以下の数を表現するためのもの．演算結果は実際の値から浮動小数点数に変換するための操作「丸め (round-off)」が常に行われる．それに伴って現れる誤差を丸め誤差と呼ぶ．
1 ワードの整数の最大値とその2進数表示．
2**(4*8-1)-1 # => 2147483647
この整数を2進数で表示するように変換するには，bin(n)を用いて，
bin(2**(4*8-1)-1) # => 0b1111111111111111111111111111111
となり，31個の1が並んでいることが分かる．1 ワードの整数の最大桁は，$n$ の長さを戻すコマンドlen(n)を使って，
len(str(2**(4*8-1)-1)) # => 10
となり，たかだか10桁程度であることが分かる．一方，64bit の場合の整数の最大桁．
len(str(2**(8*8-1)-1)) # => 19
である．
python では多倍長計算するので，通常のプログラミング言語で起こるintの最大数あたりでの奇妙な振る舞いは示さない．
2147483647+100 # => 2147483747
単精度の浮動小数点数は,仮数部2進数23bit,2倍長実数で52bitである．この有効桁数は以下の通り．
len(str(2**(23))) # => 7
len(str(2**(52))) # => 16
pythonでは普通に整数を扱う場合には，1ワードによる制約はなく，大きな数を扱える．
しかし，numpyなどでarrayとして作る場合には最大数に対する制約が現れる． 下の例では，
のがわかるでしょう．
import numpy
tmp = numpy.int(2**(8*8-1)-1)
print(tmp) # => 9223372036854775807
print(tmp+100) # => 9223372036854775907
x = numpy.array([0, 0])
print(x.dtype) # => int64
y = numpy.append(x, tmp+1)
print(y[2]) # => 9.22337203685e+18
print(y[2]+10)
print(y.dtype)
x[0] = tmp # OK
print(x[0])
x[1] = tmp+1
9223372036854775807 9223372036854775907 int64 9.22337203685e+18 9.22337203685e+18 float64 9223372036854775807
--------------------------------------------------------------------------- OverflowError Traceback (most recent call last) <ipython-input-1-153177f82482> in <module> 14 x[0] = tmp # OK 15 print(x[0]) ---> 16 x[1] = tmp+1 OverflowError: Python int too large to convert to C long
「丸め」にともなって誤差が生じる. CやFortran等の通常のプログラミング言語では「丸める」仕様なのでプログラマーが気をつけなければならない．以下のようなC programでは，予期したのとは違う結果となる．
//プログラムリスト : 実数のケタ落ち
#include <stdio.h>
int main(void){
float a,b,c;
double x,y,z;
a=1.23456789;
printf(" a= %17.10f\n",a);
b=100.0;
c=a+b;
printf("%20.10f %20.10f %20.10f\n",a,b,c);
x=(float)1.23456789;
y=(double)100;
z=x+y;
printf("%20.12e %20.12e %20.12e\n",x,y,z);
x=(double)1.23456789;
y=(double)100;
z=x+y;
printf("%20.12e %20.12e %20.12e\n",x,y,z);
return 0;
}
//プログラムリスト : 丸め誤差
#include <stdio.h>
int main(void){
float x=77777,y=7,y1,z,z1;
y1=1/y;
z=x/y;
z1=x*y1;
printf("%10.2f %10.2f\n",z,z1);
if (z!=z1){
printf("z is not equal to z1.\n");
}
printf("Surprising?? \n\n\n\n\n%10.5f %10.5f\n",z,z1);
return 0;
}
これを避けるには，EPSILONという小さな数字を定義しておいて，値の差の絶対値を求めるfabsを使って
| \hspace{100mm} |
|---|
とすべき．このときは数学関数であるfabsを使っているので，
maple
> gcc -lm test.c
とmath libraryを明示的に呼ぶのを忘れないように．
x=77777 # change this digits
y=7
y1=1.0/y
z=x/y
z1=x*y1
print("%30.20f %30.20f" % (z, z1))
if (z != z1):
print("no, different")
else:
print("yes, identical")
11111.00000000000000000000 11111.00000000000000000000 yes, identical
e = 10**(-4)
x=77777 # change this digits
y=7
y1=1.0/y
z=x/y
z1=x*y1
print("%30.20f %30.20f" % (z, z1))
if ( abs(z-z1) < e ):
print("yes, identical")
else:
print("no, different")
11111.00000000000000000000 11111.00000000000000000000 yes, identical
上の例では，浮動小数点数で計算した場合に小さい数の差を区別することができなくなるということを示している．つまり，小さい数を足したときにその計算機がその差を認識できなくなる限界ということ．これは，昔はCPU固有の精度で，今でも機械精度(Machine epsilon)と呼ばれる．今は，言語の実装によって違ってくるが，以下のようにして求めることができる．
e=1.0
w=1.0+e
while(w>1.0):
print('%-15.10e %-15.10e %-15.10e' % (e,w,w-1.0))
e = e/2.0
w = 1.0+e
1.0000000000e+00 2.0000000000e+00 1.0000000000e+00 5.0000000000e-01 1.5000000000e+00 5.0000000000e-01 2.5000000000e-01 1.2500000000e+00 2.5000000000e-01 1.2500000000e-01 1.1250000000e+00 1.2500000000e-01 6.2500000000e-02 1.0625000000e+00 6.2500000000e-02 3.1250000000e-02 1.0312500000e+00 3.1250000000e-02 1.5625000000e-02 1.0156250000e+00 1.5625000000e-02 7.8125000000e-03 1.0078125000e+00 7.8125000000e-03 3.9062500000e-03 1.0039062500e+00 3.9062500000e-03 1.9531250000e-03 1.0019531250e+00 1.9531250000e-03 9.7656250000e-04 1.0009765625e+00 9.7656250000e-04 4.8828125000e-04 1.0004882812e+00 4.8828125000e-04 2.4414062500e-04 1.0002441406e+00 2.4414062500e-04 1.2207031250e-04 1.0001220703e+00 1.2207031250e-04 6.1035156250e-05 1.0000610352e+00 6.1035156250e-05 3.0517578125e-05 1.0000305176e+00 3.0517578125e-05 1.5258789062e-05 1.0000152588e+00 1.5258789062e-05 7.6293945312e-06 1.0000076294e+00 7.6293945312e-06 3.8146972656e-06 1.0000038147e+00 3.8146972656e-06 1.9073486328e-06 1.0000019073e+00 1.9073486328e-06 9.5367431641e-07 1.0000009537e+00 9.5367431641e-07 4.7683715820e-07 1.0000004768e+00 4.7683715820e-07 2.3841857910e-07 1.0000002384e+00 2.3841857910e-07 1.1920928955e-07 1.0000001192e+00 1.1920928955e-07 5.9604644775e-08 1.0000000596e+00 5.9604644775e-08 2.9802322388e-08 1.0000000298e+00 2.9802322388e-08 1.4901161194e-08 1.0000000149e+00 1.4901161194e-08 7.4505805969e-09 1.0000000075e+00 7.4505805969e-09 3.7252902985e-09 1.0000000037e+00 3.7252902985e-09 1.8626451492e-09 1.0000000019e+00 1.8626451492e-09 9.3132257462e-10 1.0000000009e+00 9.3132257462e-10 4.6566128731e-10 1.0000000005e+00 4.6566128731e-10 2.3283064365e-10 1.0000000002e+00 2.3283064365e-10 1.1641532183e-10 1.0000000001e+00 1.1641532183e-10 5.8207660913e-11 1.0000000001e+00 5.8207660913e-11 2.9103830457e-11 1.0000000000e+00 2.9103830457e-11 1.4551915228e-11 1.0000000000e+00 1.4551915228e-11 7.2759576142e-12 1.0000000000e+00 7.2759576142e-12 3.6379788071e-12 1.0000000000e+00 3.6379788071e-12 1.8189894035e-12 1.0000000000e+00 1.8189894035e-12 9.0949470177e-13 1.0000000000e+00 9.0949470177e-13 4.5474735089e-13 1.0000000000e+00 4.5474735089e-13 2.2737367544e-13 1.0000000000e+00 2.2737367544e-13 1.1368683772e-13 1.0000000000e+00 1.1368683772e-13 5.6843418861e-14 1.0000000000e+00 5.6843418861e-14 2.8421709430e-14 1.0000000000e+00 2.8421709430e-14 1.4210854715e-14 1.0000000000e+00 1.4210854715e-14 7.1054273576e-15 1.0000000000e+00 7.1054273576e-15 3.5527136788e-15 1.0000000000e+00 3.5527136788e-15 1.7763568394e-15 1.0000000000e+00 1.7763568394e-15 8.8817841970e-16 1.0000000000e+00 8.8817841970e-16 4.4408920985e-16 1.0000000000e+00 4.4408920985e-16 2.2204460493e-16 1.0000000000e+00 2.2204460493e-16
from decimal import *
getcontext().prec=80
e=Decimal(1.0)
w=Decimal(1.0)+e
while(w>1.0):
print('%-15.10e %-15.10e' % (e,w))
e = e/Decimal(2.0)
w = Decimal(1.0)+e
1.0000000000e+00 2.0000000000e+00 5.0000000000e-01 1.5000000000e+00 2.5000000000e-01 1.2500000000e+00 1.2500000000e-01 1.1250000000e+00 6.2500000000e-02 1.0625000000e+00 3.1250000000e-02 1.0312500000e+00 1.5625000000e-02 1.0156250000e+00 7.8125000000e-03 1.0078125000e+00 3.9062500000e-03 1.0039062500e+00 1.9531250000e-03 1.0019531250e+00 9.7656250000e-04 1.0009765625e+00 4.8828125000e-04 1.0004882812e+00 2.4414062500e-04 1.0002441406e+00 1.2207031250e-04 1.0001220703e+00 6.1035156250e-05 1.0000610352e+00 3.0517578125e-05 1.0000305176e+00 1.5258789062e-05 1.0000152588e+00 7.6293945312e-06 1.0000076294e+00 3.8146972656e-06 1.0000038147e+00 1.9073486328e-06 1.0000019073e+00 9.5367431641e-07 1.0000009537e+00 4.7683715820e-07 1.0000004768e+00 2.3841857910e-07 1.0000002384e+00 1.1920928955e-07 1.0000001192e+00 5.9604644775e-08 1.0000000596e+00 2.9802322388e-08 1.0000000298e+00 1.4901161194e-08 1.0000000149e+00 7.4505805969e-09 1.0000000075e+00 3.7252902985e-09 1.0000000037e+00 1.8626451492e-09 1.0000000019e+00 9.3132257462e-10 1.0000000009e+00 4.6566128731e-10 1.0000000005e+00 2.3283064365e-10 1.0000000002e+00 1.1641532183e-10 1.0000000001e+00 5.8207660913e-11 1.0000000001e+00 2.9103830457e-11 1.0000000000e+00 1.4551915228e-11 1.0000000000e+00 7.2759576142e-12 1.0000000000e+00 3.6379788071e-12 1.0000000000e+00 1.8189894035e-12 1.0000000000e+00 9.0949470177e-13 1.0000000000e+00 4.5474735089e-13 1.0000000000e+00 2.2737367544e-13 1.0000000000e+00 1.1368683772e-13 1.0000000000e+00 5.6843418861e-14 1.0000000000e+00 2.8421709430e-14 1.0000000000e+00 1.4210854715e-14 1.0000000000e+00 7.1054273576e-15 1.0000000000e+00 3.5527136788e-15 1.0000000000e+00 1.7763568394e-15 1.0000000000e+00 8.8817841970e-16 1.0000000000e+00 4.4408920985e-16 1.0000000000e+00 2.2204460493e-16 1.0000000000e+00 1.1102230246e-16 1.0000000000e+00 5.5511151231e-17 1.0000000000e+00 2.7755575616e-17 1.0000000000e+00 1.3877787808e-17 1.0000000000e+00 6.9388939039e-18 1.0000000000e+00 3.4694469520e-18 1.0000000000e+00 1.7347234760e-18 1.0000000000e+00 8.6736173799e-19 1.0000000000e+00 4.3368086899e-19 1.0000000000e+00 2.1684043450e-19 1.0000000000e+00 1.0842021725e-19 1.0000000000e+00 5.4210108624e-20 1.0000000000e+00 2.7105054312e-20 1.0000000000e+00 1.3552527156e-20 1.0000000000e+00 6.7762635780e-21 1.0000000000e+00 3.3881317890e-21 1.0000000000e+00 1.6940658945e-21 1.0000000000e+00 8.4703294725e-22 1.0000000000e+00 4.2351647363e-22 1.0000000000e+00 2.1175823681e-22 1.0000000000e+00 1.0587911841e-22 1.0000000000e+00 5.2939559203e-23 1.0000000000e+00 2.6469779602e-23 1.0000000000e+00 1.3234889801e-23 1.0000000000e+00 6.6174449004e-24 1.0000000000e+00 3.3087224502e-24 1.0000000000e+00 1.6543612251e-24 1.0000000000e+00 8.2718061255e-25 1.0000000000e+00 4.1359030628e-25 1.0000000000e+00 2.0679515314e-25 1.0000000000e+00 1.0339757657e-25 1.0000000000e+00 5.1698788285e-26 1.0000000000e+00 2.5849394142e-26 1.0000000000e+00 1.2924697071e-26 1.0000000000e+00 6.4623485356e-27 1.0000000000e+00 3.2311742678e-27 1.0000000000e+00 1.6155871339e-27 1.0000000000e+00 8.0779356695e-28 1.0000000000e+00 4.0389678347e-28 1.0000000000e+00 2.0194839174e-28 1.0000000000e+00 1.0097419587e-28 1.0000000000e+00 5.0487097934e-29 1.0000000000e+00 2.5243548967e-29 1.0000000000e+00 1.2621774484e-29 1.0000000000e+00 6.3108872418e-30 1.0000000000e+00 3.1554436209e-30 1.0000000000e+00 1.5777218104e-30 1.0000000000e+00 7.8886090522e-31 1.0000000000e+00 3.9443045261e-31 1.0000000000e+00 1.9721522631e-31 1.0000000000e+00 9.8607613153e-32 1.0000000000e+00 4.9303806576e-32 1.0000000000e+00 2.4651903288e-32 1.0000000000e+00 1.2325951644e-32 1.0000000000e+00 6.1629758220e-33 1.0000000000e+00 3.0814879110e-33 1.0000000000e+00 1.5407439555e-33 1.0000000000e+00 7.7037197775e-34 1.0000000000e+00 3.8518598888e-34 1.0000000000e+00 1.9259299444e-34 1.0000000000e+00 9.6296497219e-35 1.0000000000e+00 4.8148248610e-35 1.0000000000e+00 2.4074124305e-35 1.0000000000e+00 1.2037062152e-35 1.0000000000e+00 6.0185310762e-36 1.0000000000e+00 3.0092655381e-36 1.0000000000e+00 1.5046327691e-36 1.0000000000e+00 7.5231638453e-37 1.0000000000e+00 3.7615819226e-37 1.0000000000e+00 1.8807909613e-37 1.0000000000e+00 9.4039548066e-38 1.0000000000e+00 4.7019774033e-38 1.0000000000e+00 2.3509887016e-38 1.0000000000e+00 1.1754943508e-38 1.0000000000e+00 5.8774717541e-39 1.0000000000e+00 2.9387358771e-39 1.0000000000e+00 1.4693679385e-39 1.0000000000e+00 7.3468396926e-40 1.0000000000e+00 3.6734198463e-40 1.0000000000e+00 1.8367099232e-40 1.0000000000e+00 9.1835496158e-41 1.0000000000e+00 4.5917748079e-41 1.0000000000e+00 2.2958874039e-41 1.0000000000e+00 1.1479437020e-41 1.0000000000e+00 5.7397185099e-42 1.0000000000e+00 2.8698592549e-42 1.0000000000e+00 1.4349296275e-42 1.0000000000e+00 7.1746481373e-43 1.0000000000e+00 3.5873240687e-43 1.0000000000e+00 1.7936620343e-43 1.0000000000e+00 8.9683101717e-44 1.0000000000e+00 4.4841550858e-44 1.0000000000e+00 2.2420775429e-44 1.0000000000e+00 1.1210387715e-44 1.0000000000e+00 5.6051938573e-45 1.0000000000e+00 2.8025969286e-45 1.0000000000e+00 1.4012984643e-45 1.0000000000e+00 7.0064923216e-46 1.0000000000e+00 3.5032461608e-46 1.0000000000e+00 1.7516230804e-46 1.0000000000e+00 8.7581154020e-47 1.0000000000e+00 4.3790577010e-47 1.0000000000e+00 2.1895288505e-47 1.0000000000e+00 1.0947644253e-47 1.0000000000e+00 5.4738221263e-48 1.0000000000e+00 2.7369110631e-48 1.0000000000e+00 1.3684555316e-48 1.0000000000e+00 6.8422776578e-49 1.0000000000e+00 3.4211388289e-49 1.0000000000e+00 1.7105694145e-49 1.0000000000e+00 8.5528470723e-50 1.0000000000e+00 4.2764235361e-50 1.0000000000e+00 2.1382117681e-50 1.0000000000e+00 1.0691058840e-50 1.0000000000e+00 5.3455294202e-51 1.0000000000e+00 2.6727647101e-51 1.0000000000e+00 1.3363823550e-51 1.0000000000e+00 6.6819117752e-52 1.0000000000e+00 3.3409558876e-52 1.0000000000e+00 1.6704779438e-52 1.0000000000e+00 8.3523897190e-53 1.0000000000e+00 4.1761948595e-53 1.0000000000e+00 2.0880974298e-53 1.0000000000e+00 1.0440487149e-53 1.0000000000e+00 5.2202435744e-54 1.0000000000e+00 2.6101217872e-54 1.0000000000e+00 1.3050608936e-54 1.0000000000e+00 6.5253044680e-55 1.0000000000e+00 3.2626522340e-55 1.0000000000e+00 1.6313261170e-55 1.0000000000e+00 8.1566305850e-56 1.0000000000e+00 4.0783152925e-56 1.0000000000e+00 2.0391576462e-56 1.0000000000e+00 1.0195788231e-56 1.0000000000e+00 5.0978941156e-57 1.0000000000e+00 2.5489470578e-57 1.0000000000e+00 1.2744735289e-57 1.0000000000e+00 6.3723676445e-58 1.0000000000e+00 3.1861838223e-58 1.0000000000e+00 1.5930919111e-58 1.0000000000e+00 7.9654595557e-59 1.0000000000e+00 3.9827297778e-59 1.0000000000e+00 1.9913648889e-59 1.0000000000e+00 9.9568244446e-60 1.0000000000e+00 4.9784122223e-60 1.0000000000e+00 2.4892061111e-60 1.0000000000e+00 1.2446030556e-60 1.0000000000e+00 6.2230152779e-61 1.0000000000e+00 3.1115076389e-61 1.0000000000e+00 1.5557538195e-61 1.0000000000e+00 7.7787690973e-62 1.0000000000e+00 3.8893845487e-62 1.0000000000e+00 1.9446922743e-62 1.0000000000e+00 9.7234613717e-63 1.0000000000e+00 4.8617306858e-63 1.0000000000e+00 2.4308653429e-63 1.0000000000e+00 1.2154326715e-63 1.0000000000e+00 6.0771633573e-64 1.0000000000e+00 3.0385816786e-64 1.0000000000e+00 1.5192908393e-64 1.0000000000e+00 7.5964541966e-65 1.0000000000e+00 3.7982270983e-65 1.0000000000e+00 1.8991135492e-65 1.0000000000e+00 9.4955677458e-66 1.0000000000e+00 4.7477838729e-66 1.0000000000e+00 2.3738919364e-66 1.0000000000e+00 1.1869459682e-66 1.0000000000e+00 5.9347298411e-67 1.0000000000e+00 2.9673649205e-67 1.0000000000e+00 1.4836824603e-67 1.0000000000e+00 7.4184123014e-68 1.0000000000e+00 3.7092061507e-68 1.0000000000e+00 1.8546030753e-68 1.0000000000e+00 9.2730153767e-69 1.0000000000e+00 4.6365076884e-69 1.0000000000e+00 2.3182538442e-69 1.0000000000e+00 1.1591269221e-69 1.0000000000e+00 5.7956346104e-70 1.0000000000e+00 2.8978173052e-70 1.0000000000e+00 1.4489086526e-70 1.0000000000e+00 7.2445432631e-71 1.0000000000e+00 3.6222716315e-71 1.0000000000e+00 1.8111358158e-71 1.0000000000e+00 9.0556790788e-72 1.0000000000e+00 4.5278395394e-72 1.0000000000e+00 2.2639197697e-72 1.0000000000e+00 1.1319598849e-72 1.0000000000e+00 5.6597994243e-73 1.0000000000e+00 2.8298997121e-73 1.0000000000e+00 1.4149498561e-73 1.0000000000e+00 7.0747492803e-74 1.0000000000e+00 3.5373746402e-74 1.0000000000e+00 1.7686873201e-74 1.0000000000e+00 8.8434366004e-75 1.0000000000e+00 4.4217183002e-75 1.0000000000e+00 2.2108591501e-75 1.0000000000e+00 1.1054295751e-75 1.0000000000e+00 5.5271478753e-76 1.0000000000e+00 2.7635739376e-76 1.0000000000e+00 1.3817869688e-76 1.0000000000e+00 6.9089348441e-77 1.0000000000e+00 3.4544674220e-77 1.0000000000e+00 1.7272337110e-77 1.0000000000e+00 8.6361685551e-78 1.0000000000e+00 4.3180842775e-78 1.0000000000e+00 2.1590421388e-78 1.0000000000e+00 1.0795210694e-78 1.0000000000e+00 5.3976053469e-79 1.0000000000e+00 2.6988026735e-79 1.0000000000e+00 1.3494013367e-79 1.0000000000e+00 6.7470066837e-80 1.0000000000e+00
有効数字がそれぞれ5桁で計算した結果を示せ．
maple
0.723657
- 0.723649
------------
maple
72365.7
- 1.23659
------------
maple
72365.7
- 0.001
------------
pythonにはroundという標準的な丸める関数がありますが，使用には注意が必要です．
print('0.4 =>', round(0.4))
print('0.5 =>', round(0.5))
print('0.6 =>', round(0.6))
0.4 => 0 0.5 => 0 0.6 => 1
print('0.5 =>', round(0.5))
print('1.5 =>', round(1.5))
print('2.5 =>', round(2.5))
print('3.5 =>', round(3.5))
0.5 => 0 1.5 => 2 2.5 => 2 3.5 => 4
print('0.05 =>', round(0.05, 1))
print('0.15 =>', round(0.15, 1))
print('0.25 =>', round(0.25, 1))
print('0.35 =>', round(0.35, 1))
print('0.45 =>', round(0.45, 1))
0.05 => 0.1 0.15 => 0.1 0.25 => 0.2 0.35 => 0.3 0.45 => 0.5
注釈 浮動小数点数に対する round() の振る舞いは意外なものかもしれません: 例えば、 round(2.675, 2) は予想通りの 2.68 ではなく 2.67 を与えます。これはバグではありません: これはほとんどの小数が浮動小数点数で正確に表せないことの結果です。
2. 組み込み関数 round() — Python 3.6.3 ドキュメント
さらに，pythonにはdecimalという，10進固定及び浮動小数点数演算用のライブラリがあります．この中のquantize関数で直感的な四捨五入が実現できます．こちらの方が良さそう．
from decimal import *
def pretty_p(result,a,b,operator):
print('context.prec:{}'.format(getcontext().prec))
print(' %20.14f' % (a))
print( '%1s%20.14f' % (operator, b))
print('-----------')
print( ' %20.14f' % (result))
n = 10
getcontext().prec = n
a=Decimal('0.723657').quantize(Decimal(10) ** -n)
b=Decimal('0.723649').quantize(Decimal(10) ** -n)
pretty_p(a-b,a,b,'-')
a=Decimal('0.723657').quantize(Decimal(10) ** -n)*100000
b=Decimal('0.123659').quantize(Decimal(10) ** -n)*10
pretty_p(a-b,a,b,'-')
a=Decimal('0.723657').quantize(Decimal(10) ** -n)*100000
b=Decimal('0.1').quantize(Decimal(10) ** -n)/100
pretty_p(a+b,a,b,'+')
context.prec:10
0.72365700000000
- 0.72364900000000
-----------
0.00000800000000
context.prec:10
72365.69999999999709
- 1.23659000000000
-----------
72364.46340999999666
context.prec:10
72365.69999999999709
+ 0.00100000000000
-----------
72365.70100000000093
n = 5
getcontext().prec = n
a=Decimal('0.723657').quantize(Decimal(10) ** -n)
b=Decimal('0.723649').quantize(Decimal(10) ** -n)
pretty_p(a-b,a,b,'-')
a=Decimal('0.723657').quantize(Decimal(10) ** -n)*100000
b=Decimal('0.123659').quantize(Decimal(10) ** -n)*10
pretty_p(a-b,a,b,'-')
a=Decimal('0.723657').quantize(Decimal(10) ** -n)*100000
b=Decimal('0.1').quantize(Decimal(10) ** -n)/100
pretty_p(a+b,a,b,'+')
context.prec:5
0.72366000000000
- 0.72365000000000
-----------
0.00001000000000
context.prec:5
72366.00000000000000
- 1.23660000000000
-----------
72365.00000000000000
context.prec:5
72366.00000000000000
+ 0.00100000000000
-----------
72366.00000000000000
同様に，0.81321/(23.173-23.094) を有効数字がそれぞれ5桁，4桁，3桁，2桁で計算した結果を示せ.(2018年度期末試験)
通常の解の公式を使った解と, 解と係数の関係(下記の記述を参照）を使った解とを出力するプログラムをpythonで作成すると以下の通りとなる． 解の有効数字が2種類の計算方法の違いでどう違うか， いくつかの精度(precision)で実行させた結果を使って解説せよ． (2017年度期末試験)
2次方程式 $ax^2+bx+c=0$の 係数$a, b, c$が特殊な値をもつ場合，通常の解の公式
$$ x = \frac {-b \pm \sqrt{{b}^{2}-4ac}}{2a} $$にしたがって計算するとケタ落ちによる間違った答えを出す．その特殊な値とは
$$ \sqrt{{b}^{2}-4ac} \approx |b| $$となる場合である．
ケタ落ちを防ぐには, $b > 0$の場合は,
$$ x_1 = \frac {-b - \sqrt{{b}^{2}-4ac}}{2a} $$として，ケタ落ちを起こさずに求め, この解を使って, 解と係数の関係より
$$ x_2 = \frac {c}{a\, x_1} \,\ldots(1) $$で求める．$b < 0$ の場合は，解の公式のたし算の方 $$ x_1 = \frac {-b + \sqrt{{b}^{2}-4ac}}{2a} $$ を使って同様に求める．
from numpy import sqrt
def solve_normal_formula(a,b,c):
x0=(-b-sqrt(b**2-4*a*c))/(2*a)
x1=(-b+sqrt(b**2-4*a*c))/(2*a)
return (x0,x1)
def solve_precise_formula(a,b,c):
x0=(-b-sqrt(b**2-4*a*c))/(2*a)
x1=c/(a*x0)
return (x0,x1)
b = 40
print(solve_normal_formula(1,b,2))
print(solve_precise_formula(1,b,2))
(-39.949937343260004, -0.050062656739996214) (-39.949937343260004, -0.050062656739996665)
from decimal import *
prec = 7 #13
b= '40'
print("\n b =", b)
print("\nprec=", prec)
getcontext().prec = prec
print(solve_normal_formula(Decimal('1'),
Decimal(b),
Decimal('2')))
print(solve_precise_formula(Decimal('1'),
Decimal(b),
Decimal('2')))
b = 40
prec= 7
(Decimal('-39.94994'), Decimal('-0.050065'))
(Decimal('-39.94994'), Decimal('-0.05006265'))
(1)式の導出を示しておく．解と係数の関係を導びくと $$ (x-x_1)(x-x_2) = x^2 -(x_1+x_2)x + x_1 x_2 \\ = \frac{a}{a}x^2 + \frac{b}{a}x + \frac{c}{a} $$ となる．(1)式は最後の定数項より， $$ \frac{c}{a} = x_1 x_2 $$ から導かれる．