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