Skip to content

Instantly share code, notes, and snippets.

@ohno

ohno/Big.ipynb Secret

Created August 24, 2022 17:32
Show Gist options
  • Save ohno/e4d656aec716a71a741fff7ae23b9058 to your computer and use it in GitHub Desktop.
Save ohno/e4d656aec716a71a741fff7ae23b9058 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
Big.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Juliaで任意精度演算に入門してみた\n",
"\n",
"Juliaでは任意精度演算が標準でサポートされています. 特別なライブラリなども不要なあたりに次世代の計算科学を担っていく覚悟と矜持を感じますね. バックエンドではBigFloat型には[the GNU MPFR library](https://www.mpfr.org/)が, BigInt型には[the GNU Multiple Precision Arithmetic Library (GMP)](https://gmplib.org/)が使われているようです."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 環境"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Julia Version 1.7.1\n",
"Commit ac5cc99908 (2021-12-22 19:35 UTC)\n",
"Platform Info:\n",
" OS: Windows (x86_64-w64-mingw32)\n",
" CPU: Intel(R) Core(TM) i7-4650U CPU @ 1.70GHz\n",
" WORD_SIZE: 64\n",
" LIBM: libopenlibm\n",
" LLVM: libLLVM-12.0.1 (ORCJIT, haswell)\n"
]
}
],
"source": [
"versioninfo()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 基本的な使い方\n",
"\n",
"デフォルトでは, 整数はInt64型, 小数はFloat64型になります. Float64型はいわゆる倍精度浮動小数点数というものです. `b = big(2)`とすると, BigInt型になります."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(a) = Int64\n",
"typeof(b) = Float64\n",
"typeof(c) = BigInt\n",
"typeof(d) = BigFloat\n",
"a = 2\n",
"b = 2.0\n",
"c = 2\n",
"d = 2.0\n"
]
},
{
"data": {
"text/plain": [
"2.0"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a = 2\n",
"b = 2.0\n",
"c = big(2)\n",
"d = big(2.0)\n",
"\n",
"@show typeof(a)\n",
"@show typeof(b)\n",
"@show typeof(c)\n",
"@show typeof(d)\n",
"\n",
"@show a\n",
"@show b\n",
"@show c\n",
"@show d"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"標準で用意されているほとんどの関数は, 引数にBigIntまたはBigFloatを渡すと, 戻り値にはBigIntまたはBigFloatが返ってきます. 例えば$\\sqrt{2}$を計算してみましょう."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(a) = Float64\n",
"typeof(b) = BigFloat\n",
"a = 1.4142135623730951\n",
"b = 1.414213562373095048801688724209698078569671875376948073176679737990732478462102\n"
]
},
{
"data": {
"text/plain": [
"1.414213562373095048801688724209698078569671875376948073176679737990732478462102"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a = sqrt(2)\n",
"b = sqrt(big(2))\n",
"\n",
"@show typeof(a)\n",
"@show typeof(b)\n",
"\n",
"@show a\n",
"@show b"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`sqrt(big(2))`の戻り値を[こちらのサイト](https://apod.nasa.gov/htmltest/gifcity/sqrt2.1mil)と比較すると\n",
"\n",
"- 1.414213562373095048801688724209698078569671875376948073176679737990732478462102\n",
"- 1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850...\n",
"\n",
"となっており, 最後の1桁だけズレているものの, 78桁は一致しています. すごい. これはまだまだ序の口です."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 有効桁数\n",
"\n",
"[有効桁数は任意の長さに変えられる](https://docs.julialang.org/en/v1/base/numbers/#Base.MPFR.setprecision)とこのことです. まず, デフォルトでは"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(a) = Float64\n",
"typeof(b) = BigFloat\n",
"precision(a) = 53\n",
"precision(b) = 256\n",
"sqrt(a) = 1.4142135623730951\n",
"sqrt(b) = 1.414213562373095048801688724209698078569671875376948073176679737990732478462102\n"
]
},
{
"data": {
"text/plain": [
"1.414213562373095048801688724209698078569671875376948073176679737990732478462102"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a = 2.0\n",
"b = big(2.0)\n",
"\n",
"@show typeof(a)\n",
"@show typeof(b)\n",
"\n",
"@show precision(a)\n",
"@show precision(b)\n",
"\n",
"@show sqrt(a)\n",
"@show sqrt(b)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"です. 例えばこれを`setpresition(BigFloat, 53)`とすると, 仮数に53ビット使われますので, BigFloat型がFloat64型と同じ有効桁数になります. 浮動小数点数の仮数・基数・指数については[Wikipedia](https://ja.wikipedia.org/wiki/%E6%B5%AE%E5%8B%95%E5%B0%8F%E6%95%B0%E7%82%B9%E6%95%B0)を参照してください."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(a) = Float64\n",
"typeof(b) = BigFloat\n",
"precision(a) = 53\n",
"precision(b) = 53\n",
"sqrt(a) = 1.4142135623730951\n",
"sqrt(b) = 1.4142135623730951\n"
]
},
{
"data": {
"text/plain": [
"1.4142135623730951"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"setprecision(BigFloat, 53)\n",
"\n",
"a = 2.0\n",
"b = big(2.0)\n",
"\n",
"@show typeof(a)\n",
"@show typeof(b)\n",
"\n",
"@show precision(a)\n",
"@show precision(b)\n",
"\n",
"@show sqrt(a)\n",
"@show sqrt(b)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`setpresition(BigFloat, ?)`で仮数の長さを変更すると, 関数の戻り値の桁数も変わります. グローバルに変更されるので, 並列計算の途中で変えるといった使い方はやめた方がよいです. 例えば 2^13bit = 8192bit = 1024B = 1KB を仮数部に使って$\\sqrt{2}$を先ほどよりも正確に求めてみましょう."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"sqrt(big(2)) = 1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462309122970249248360558507372126441214970999358314132226659275055927557999505011527820605714701095599716059702745345968620147285174186408891986095523292304843087143214508397626036279952514079896872533965463318088296406206152583523950547457502877599617298355752203375318570113543746034084988471603868999706990048150305440277903164542478230684929369186215805784631115966687130130156185689872372352885092648612494977154218334204285686060146824720771435854874155657069677653720226485447015858801620758474922657226002085584466521458398893944370926591800311388246468157082630100594858704003186480342194897278290641045072636881313739855256117322040245091227700226941127573627280495738108967504018369868368450725799364729060762996941380475654823728997180326802474420629269124859052181004459842150591120249441341728531478105803603371077309182869314710171111683916581726889419758716582152128229518488472089694633862891562882765952635140542267653239694617511291602408715510135150455381287560052631468017127402653969470240300517495318862925631385188163478001569369176881852378684052287837629389214300655869568685964595155501644724509836896036887323114389415576651040883914292338113206052433629485317049915771756228549741438999188021762430965206564211827316726257539594717255934637238632261482742622208671155839599926521176252698917540988159348640083457085181472231814204070426509056532333398436457865796796519267292399875366617215982578860263363617827495994219403777753681426217738799194551397231274066898329989895386728822856378697749662519966583525776198939322845344735694794962952168891485492538904755828834526096524096542889394538646625744927556381964410316979833061852019379384940057156333720548068540575867999670121372239475821426306585132217408832382947287617393647467837431960001592188807347857617252211867490424977366929207311096369721608933708661156734585334833295254675851644710757848602463600834449114818587655554286455123314219926311332517970608436559704352856410087918500760361009159465670676883605571740076756905096136719401324935605240185999105062108163597726431380605467010293569971042425105781749531057255934984451126922780344913506637568747760283162829605532422426957534529028838768446429173282770888318087025339852338122749990812371892540726475367850304821591801886167108972869229201197599880703818543332536460211082299279293072871780799888099176741771\n"
]
},
{
"data": {
"text/plain": [
"1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462309122970249248360558507372126441214970999358314132226659275055927557999505011527820605714701095599716059702745345968620147285174186408891986095523292304843087143214508397626036279952514079896872533965463318088296406206152583523950547457502877599617298355752203375318570113543746034084988471603868999706990048150305440277903164542478230684929369186215805784631115966687130130156185689872372352885092648612494977154218334204285686060146824720771435854874155657069677653720226485447015858801620758474922657226002085584466521458398893944370926591800311388246468157082630100594858704003186480342194897278290641045072636881313739855256117322040245091227700226941127573627280495738108967504018369868368450725799364729060762996941380475654823728997180326802474420629269124859052181004459842150591120249441341728531478105803603371077309182869314710171111683916581726889419758716582152128229518488472089694633862891562882765952635140542267653239694617511291602408715510135150455381287560052631468017127402653969470240300517495318862925631385188163478001569369176881852378684052287837629389214300655869568685964595155501644724509836896036887323114389415576651040883914292338113206052433629485317049915771756228549741438999188021762430965206564211827316726257539594717255934637238632261482742622208671155839599926521176252698917540988159348640083457085181472231814204070426509056532333398436457865796796519267292399875366617215982578860263363617827495994219403777753681426217738799194551397231274066898329989895386728822856378697749662519966583525776198939322845344735694794962952168891485492538904755828834526096524096542889394538646625744927556381964410316979833061852019379384940057156333720548068540575867999670121372239475821426306585132217408832382947287617393647467837431960001592188807347857617252211867490424977366929207311096369721608933708661156734585334833295254675851644710757848602463600834449114818587655554286455123314219926311332517970608436559704352856410087918500760361009159465670676883605571740076756905096136719401324935605240185999105062108163597726431380605467010293569971042425105781749531057255934984451126922780344913506637568747760283162829605532422426957534529028838768446429173282770888318087025339852338122749990812371892540726475367850304821591801886167108972869229201197599880703818543332536460211082299279293072871780799888099176741771"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"setprecision(BigFloat, 2^13)\n",
"@show sqrt(big(2))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[こちらのサイト](https://apod.nasa.gov/htmltest/gifcity/sqrt2.1mil)と比較すると, 最後の1桁だけズレていますが, 小数点以下2466桁まで正しく計算できていることがわかります. 最後の数桁をページ内検索すると確かめられます. ちなみに, `setprecision(BigFloat, 2^20)`として同様の計算を行うと最後の1桁だけズレますが小数点以下315652桁まで正しく計算できます."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"デフォルトの`setprecision(BigFloat, 256)`に戻しておきます."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"precision(big(2.0)) = 256\n",
"sqrt(big(2.0)) = 1.414213562373095048801688724209698078569671875376948073176679737990732478462102\n"
]
},
{
"data": {
"text/plain": [
"1.414213562373095048801688724209698078569671875376948073176679737990732478462102"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"setprecision(BigFloat, 256)\n",
"@show precision(big(2.0))\n",
"@show sqrt(big(2.0))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 型の規則\n",
"\n",
"それぞれの型の和を考えるとき, Int＋Float＝Floatとなり, 64＋Big＝Bigとなります.\n",
"\n",
"| | Int64 | Float64 | BigInt | BigFloat | \n",
"| ---------- | -------- | -------- | -------- | -------- | \n",
"| ＋Int64 | Int64 | Float64 | BigInt | BigFloat | \n",
"| ＋Float64 | Float64 | Float64 | BigFloat | BigFloat | \n",
"| ＋BigInt | BigInt | BigFloat | BigInt | BigFloat | \n",
"| ＋BigFloat | BigFloat | BigFloat | BigFloat | BigFloat | \n",
"\n",
"それぞれの型の積を考えるときも同じく, Int×Float＝Floatとなり, 64×Big＝Bigとなります.\n",
"\n",
"| | Int64 | Float64 | BigInt | BigFloat | \n",
"| ---------- | -------- | -------- | -------- | -------- | \n",
"| ×Int64 | Int64 | Float64 | BigInt | BigFloat | \n",
"| ×Float64 | Float64 | Float64 | BigFloat | BigFloat | \n",
"| ×BigInt | BigInt | BigFloat | BigInt | BigFloat | \n",
"| ×BigFloat | BigFloat | BigFloat | BigFloat | BigFloat | \n",
"\n",
"検証："
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(a) = Int64\n",
"typeof(b) = Float64\n",
"typeof(c) = BigInt\n",
"typeof(d) = BigFloat\n",
"a = 2\n",
"b = 2.0\n",
"c = 2\n",
"d = 2.0\n",
"typeof(a + a) = Int64\n",
"typeof(a + b) = Float64\n",
"typeof(a + c) = BigInt\n",
"typeof(a + d) = BigFloat\n",
"typeof(b + a) = Float64\n",
"typeof(b + b) = Float64\n",
"typeof(b + c) = BigFloat\n",
"typeof(b + d) = BigFloat\n",
"typeof(c + a) = BigInt\n",
"typeof(c + b) = BigFloat\n",
"typeof(c + c) = BigInt\n",
"typeof(c + d) = BigFloat\n",
"typeof(d + a) = BigFloat\n",
"typeof(d + b) = BigFloat\n",
"typeof(d + c) = BigFloat\n",
"typeof(d + d) = BigFloat\n",
"typeof(a * a) = Int64\n",
"typeof(a * b) = Float64\n",
"typeof(a * c) = BigInt\n",
"typeof(a * d) = BigFloat\n",
"typeof(b * a) = Float64\n",
"typeof(b * b) = Float64\n",
"typeof(b * c) = BigFloat\n",
"typeof(b * d) = BigFloat\n",
"typeof(c * a) = BigInt\n",
"typeof(c * b) = BigFloat\n",
"typeof(c * c) = BigInt\n",
"typeof(c * d) = BigFloat\n",
"typeof(d * a) = BigFloat\n",
"typeof(d * b) = BigFloat\n",
"typeof(d * c) = BigFloat\n",
"typeof(d * d) = BigFloat\n"
]
},
{
"data": {
"text/plain": [
"BigFloat"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a = 2\n",
"b = 2.0\n",
"c = big(2)\n",
"d = big(2.0)\n",
"\n",
"@show typeof(a)\n",
"@show typeof(b)\n",
"@show typeof(c)\n",
"@show typeof(d)\n",
"\n",
"@show a\n",
"@show b\n",
"@show c\n",
"@show d\n",
"\n",
"@show typeof(a+a)\n",
"@show typeof(a+b)\n",
"@show typeof(a+c)\n",
"@show typeof(a+d)\n",
"\n",
"@show typeof(b+a)\n",
"@show typeof(b+b)\n",
"@show typeof(b+c)\n",
"@show typeof(b+d)\n",
"\n",
"@show typeof(c+a)\n",
"@show typeof(c+b)\n",
"@show typeof(c+c)\n",
"@show typeof(c+d)\n",
"\n",
"@show typeof(d+a)\n",
"@show typeof(d+b)\n",
"@show typeof(d+c)\n",
"@show typeof(d+d)\n",
"\n",
"@show typeof(a*a)\n",
"@show typeof(a*b)\n",
"@show typeof(a*c)\n",
"@show typeof(a*d)\n",
"\n",
"@show typeof(b*a)\n",
"@show typeof(b*b)\n",
"@show typeof(b*c)\n",
"@show typeof(b*d)\n",
"\n",
"@show typeof(c*a)\n",
"@show typeof(c*b)\n",
"@show typeof(c*c)\n",
"@show typeof(c*d)\n",
"\n",
"@show typeof(d*a)\n",
"@show typeof(d*b)\n",
"@show typeof(d*c)\n",
"@show typeof(d*d)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"下記の3つの積だけを見ると分かりやすいですが, 同じBigFloatであっても, Float64×BigFloatでは有効桁数はFloat64に引っ張られるので, 精度を保ちたい場面ではBigFloatだけで計算を行うことをお勧めします."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(a) = Float64\n",
"typeof(b) = BigFloat\n",
"a = 1.4142135623730951\n",
"b = 1.414213562373095048801688724209698078569671875376948073176679737990732478462102\n",
"a * a = 2.0000000000000004\n",
"a * b = 2.000000000000000136716173153238459361596582123216045722514805262001562413306807\n",
"b * b = 1.999999999999999999999999999999999999999999999999999999999999999999999999999983\n"
]
},
{
"data": {
"text/plain": [
"1.999999999999999999999999999999999999999999999999999999999999999999999999999983"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a = sqrt(2)\n",
"b = sqrt(big(2))\n",
"\n",
"@show typeof(a)\n",
"@show typeof(b)\n",
"\n",
"@show a\n",
"@show b\n",
"\n",
"@show a*a\n",
"@show a*b\n",
"@show b*b"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 応用例：階乗\n",
"\n",
"Int64型で扱える最大の整数は9223372036854775807です."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typemax(Int64) = 9223372036854775807\n"
]
},
{
"data": {
"text/plain": [
"9223372036854775807"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@show typemax(Int64)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"階乗$n!$は$n$に対して爆発的に大きな値を取るため, Int64型の限界がすぐに見れます."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1!\t=\t1\n",
"2!\t=\t2\n",
"3!\t=\t6\n",
"4!\t=\t24\n",
"5!\t=\t120\n",
"6!\t=\t720\n",
"7!\t=\t5040\n",
"8!\t=\t40320\n",
"9!\t=\t362880\n",
"10!\t=\t3628800\n",
"11!\t=\t39916800\n",
"12!\t=\t479001600\n",
"13!\t=\t6227020800\n",
"14!\t=\t87178291200\n",
"15!\t=\t1307674368000\n",
"16!\t=\t20922789888000\n",
"17!\t=\t355687428096000\n",
"18!\t=\t6402373705728000\n",
"19!\t=\t121645100408832000\n",
"20!\t=\t2432902008176640000\n"
]
},
{
"ename": "LoadError",
"evalue": "OverflowError: 21 is too large to look up in the table; consider using `factorial(big(21))` instead",
"output_type": "error",
"traceback": [
"OverflowError: 21 is too large to look up in the table; consider using `factorial(big(21))` instead",
"",
"Stacktrace:",
" [1] factorial_lookup",
" @ .\\combinatorics.jl:19 [inlined]",
" [2] factorial",
" @ .\\combinatorics.jl:27 [inlined]",
" [3] top-level scope",
" @ .\\In[11]:2",
" [4] eval",
" @ .\\boot.jl:373 [inlined]",
" [5] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)",
" @ Base .\\loading.jl:1196"
]
}
],
"source": [
"for n in 1:60\n",
" println(\"$(n)!\\t=\\t\", factorial(n))\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"このように通常の整数だと$21!$が計算できずにエラーが出てしまいますが, 以下のように` factorial(n)`を` factorial(big(n))`に置き換えることでより大きい数値を計算できるようになります."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1!\t=\t1\n",
"2!\t=\t2\n",
"3!\t=\t6\n",
"4!\t=\t24\n",
"5!\t=\t120\n",
"6!\t=\t720\n",
"7!\t=\t5040\n",
"8!\t=\t40320\n",
"9!\t=\t362880\n",
"10!\t=\t3628800\n",
"11!\t=\t39916800\n",
"12!\t=\t479001600\n",
"13!\t=\t6227020800\n",
"14!\t=\t87178291200\n",
"15!\t=\t1307674368000\n",
"16!\t=\t20922789888000\n",
"17!\t=\t355687428096000\n",
"18!\t=\t6402373705728000\n",
"19!\t=\t121645100408832000\n",
"20!\t=\t2432902008176640000\n",
"21!\t=\t51090942171709440000\n",
"22!\t=\t1124000727777607680000\n",
"23!\t=\t25852016738884976640000\n",
"24!\t=\t620448401733239439360000\n",
"25!\t=\t15511210043330985984000000\n",
"26!\t=\t403291461126605635584000000\n",
"27!\t=\t10888869450418352160768000000\n",
"28!\t=\t304888344611713860501504000000\n",
"29!\t=\t8841761993739701954543616000000\n",
"30!\t=\t265252859812191058636308480000000\n",
"31!\t=\t8222838654177922817725562880000000\n",
"32!\t=\t263130836933693530167218012160000000\n",
"33!\t=\t8683317618811886495518194401280000000\n",
"34!\t=\t295232799039604140847618609643520000000\n",
"35!\t=\t10333147966386144929666651337523200000000\n",
"36!\t=\t371993326789901217467999448150835200000000\n",
"37!\t=\t13763753091226345046315979581580902400000000\n",
"38!\t=\t523022617466601111760007224100074291200000000\n",
"39!\t=\t20397882081197443358640281739902897356800000000\n",
"40!\t=\t815915283247897734345611269596115894272000000000\n",
"41!\t=\t33452526613163807108170062053440751665152000000000\n",
"42!\t=\t1405006117752879898543142606244511569936384000000000\n",
"43!\t=\t60415263063373835637355132068513997507264512000000000\n",
"44!\t=\t2658271574788448768043625811014615890319638528000000000\n",
"45!\t=\t119622220865480194561963161495657715064383733760000000000\n",
"46!\t=\t5502622159812088949850305428800254892961651752960000000000\n",
"47!\t=\t258623241511168180642964355153611979969197632389120000000000\n",
"48!\t=\t12413915592536072670862289047373375038521486354677760000000000\n",
"49!\t=\t608281864034267560872252163321295376887552831379210240000000000\n",
"50!\t=\t30414093201713378043612608166064768844377641568960512000000000000\n",
"51!\t=\t1551118753287382280224243016469303211063259720016986112000000000000\n",
"52!\t=\t80658175170943878571660636856403766975289505440883277824000000000000\n",
"53!\t=\t4274883284060025564298013753389399649690343788366813724672000000000000\n",
"54!\t=\t230843697339241380472092742683027581083278564571807941132288000000000000\n",
"55!\t=\t12696403353658275925965100847566516959580321051449436762275840000000000000\n",
"56!\t=\t710998587804863451854045647463724949736497978881168458687447040000000000000\n",
"57!\t=\t40526919504877216755680601905432322134980384796226602145184481280000000000000\n",
"58!\t=\t2350561331282878571829474910515074683828862318181142924420699914240000000000000\n",
"59!\t=\t138683118545689835737939019720389406345902876772687432540821294940160000000000000\n",
"60!\t=\t8320987112741390144276341183223364380754172606361245952449277696409600000000000000\n"
]
}
],
"source": [
"for n in 1:60\n",
" println(\"$(n)!\\t=\\t\", factorial(big(n)))\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"このように, BigInt型を使うことでInt64型の最大値9223372036854775807を超えて, $21! = 51090942171709440000$を計算することができました. では一体, どこまで大きな数を扱えるのか各自で試してみましょう. $1000!$の例を示します."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1000!\t=\t402387260077093773543702433923003985719374864210714632543799910429938512398629020592044208486969404800479988610197196058631666872994808558901323829669944590997424504087073759918823627727188732519779505950995276120874975462497043601418278094646496291056393887437886487337119181045825783647849977012476632889835955735432513185323958463075557409114262417474349347553428646576611667797396668820291207379143853719588249808126867838374559731746136085379534524221586593201928090878297308431392844403281231558611036976801357304216168747609675871348312025478589320767169132448426236131412508780208000261683151027341827977704784635868170164365024153691398281264810213092761244896359928705114964975419909342221566832572080821333186116811553615836546984046708975602900950537616475847728421889679646244945160765353408198901385442487984959953319101723355556602139450399736280750137837615307127761926849034352625200015888535147331611702103968175921510907788019393178114194545257223865541461062892187960223838971476088506276862967146674697562911234082439208160153780889893964518263243671616762179168909779911903754031274622289988005195444414282012187361745992642956581746628302955570299024324153181617210465832036786906117260158783520751516284225540265170483304226143974286933061690897968482590125458327168226458066526769958652682272807075781391858178889652208164348344825993266043367660176999612831860788386150279465955131156552036093988180612138558600301435694527224206344631797460594682573103790084024432438465657245014402821885252470935190620929023136493273497565513958720559654228749774011413346962715422845862377387538230483865688976461927383814900140767310446640259899490222221765904339901886018566526485061799702356193897017860040811889729918311021171229845901641921068884387121855646124960798722908519296819372388642614839657382291123125024186649353143970137428531926649875337218940694281434118520158014123344828015051399694290153483077644569099073152433278288269864602789864321139083506217095002597389863554277196742822248757586765752344220207573630569498825087968928162753848863396909959826280956121450994871701244516461260379029309120889086942028510640182154399457156805941872748998094254742173582401063677404595741785160829230135358081840096996372524230560855903700624271243416909004153690105933983835777939410970027753472000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000\n"
]
}
],
"source": [
"n = 1000\n",
"println(\"$(n)!\\t=\\t\", factorial(big(n)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" また, 少なくとも$10000!$までは計算できました. もはや正しいのかわからないですね. 検証法を調べた方や思いついた方がいたら教えてください. なお, BigIntではBigFloatとは違って, 勝手に有効桁数がが変わっているのではないかと思われますが, こちらも検証方法がわかりませんでした."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 応用例：円周率\n",
"\n",
"[こちらのサイト](https://www.tstcl.jp/ja/randd/pi.php)に載っている値と比較すると, デフォルトの桁数のBigFloat型で展開した円周率πは小数点以下75桁までは正しく, 最後の2,3桁がズレます."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"big(π) = 3.141592653589793238462643383279502884197169399375105820974944592307816406286198\n",
"exact = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348...\n"
]
}
],
"source": [
"@show big(π)\n",
"println(\"exact = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348...\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"もちろん, 精度を上げることはできますが, どんな桁数にしても最後の2,3桁はズレると思った方が良いです. 下記の例では最後の1桁だけズレていますが, 小数点以下2465桁までは合っています."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"BigFloat(π, precision = 2 ^ 13) = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420198938095257201065485863278865936153381827968230301952035301852968995773622599413891249721775283479131515574857242454150695950829533116861727855889075098381754637464939319255060400927701671139009848824012858361603563707660104710181942955596198946767837449448255379774726847104047534646208046684259069491293313677028989152104752162056966024058038150193511253382430035587640247496473263914199272604269922796782354781636009341721641219924586315030286182974555706749838505494588586926995690927210797509302955321165344987202755960236480665499119881834797753566369807426542527862551818417574672890977772793800081647060016145249192173217214772350141441973568548161361157352552133475741849468438523323907394143334547762416862518983569485562099219222184272550254256887671790494601653466804988627232791786085784383827967976681454100953883786360950680064225125205117392984896084128488626945604241965285022210661186306744278622039194945047123713786960956364371917287467764657573962413890865832645995813390478027590099465764078951269468398352595709825822620522489407726719478268482601476990902640136394437455305068203496252451749399651431429809190659250937221696461515709858387410597885959772975498930161753928468138268683868942774155991855925245953959431049972524680845987273644695848653836736222626099124608051243884390451244136549762780797715691435997700129616089441694868555848406353422072225828488648158456028506016842739452267467678895252138522549954666727823986456596116354888\n"
]
},
{
"data": {
"text/plain": [
"3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420198938095257201065485863278865936153381827968230301952035301852968995773622599413891249721775283479131515574857242454150695950829533116861727855889075098381754637464939319255060400927701671139009848824012858361603563707660104710181942955596198946767837449448255379774726847104047534646208046684259069491293313677028989152104752162056966024058038150193511253382430035587640247496473263914199272604269922796782354781636009341721641219924586315030286182974555706749838505494588586926995690927210797509302955321165344987202755960236480665499119881834797753566369807426542527862551818417574672890977772793800081647060016145249192173217214772350141441973568548161361157352552133475741849468438523323907394143334547762416862518983569485562099219222184272550254256887671790494601653466804988627232791786085784383827967976681454100953883786360950680064225125205117392984896084128488626945604241965285022210661186306744278622039194945047123713786960956364371917287467764657573962413890865832645995813390478027590099465764078951269468398352595709825822620522489407726719478268482601476990902640136394437455305068203496252451749399651431429809190659250937221696461515709858387410597885959772975498930161753928468138268683868942774155991855925245953959431049972524680845987273644695848653836736222626099124608051243884390451244136549762780797715691435997700129616089441694868555848406353422072225828488648158456028506016842739452267467678895252138522549954666727823986456596116354888"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@show BigFloat(π, precision=2^13)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[ライプニッツの公式](https://ja.wikipedia.org/wiki/%E3%83%A9%E3%82%A4%E3%83%97%E3%83%8B%E3%83%83%E3%83%84%E3%81%AE%E5%85%AC%E5%BC%8F)を使って, 円周率$\\pi$の計算を実装してみましょう.\n",
"\n",
"$$\n",
"\\pi = \\sum_{n=0}^\\infty \\frac{4(-1)^n}{2n+1}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0\t4.0\n",
"1\t2.666666666666667\n",
"2\t3.466666666666667\n",
"3\t2.8952380952380956\n",
"4\t3.3396825396825403\n",
"5\t2.9760461760461765\n",
"6\t3.2837384837384844\n",
"7\t3.017071817071818\n",
"8\t3.2523659347188767\n",
"9\t3.0418396189294032\n",
"10\t3.232315809405594\n",
"1000\t3.1425916543395442\n",
"2000\t3.1420924036835256\n",
"3000\t3.1419258758397897\n",
"4000\t3.1418425911015158\n",
"5000\t3.1417926135957908\n",
"6000\t3.1417592924821482\n",
"7000\t3.141735490326666\n",
"8000\t3.1417176379662446\n",
"9000\t3.1417037523562383\n",
"10000\t3.1416926435905346\n",
"∞\t3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628...\n"
]
}
],
"source": [
"x = 0.0\n",
"\n",
"for n in 0:100000\n",
" x += 4 * (-1)^n / (2*n+1)\n",
" if n in union(0:10, 1000:1000:10000); println(\"$n\\t\", x); end\n",
"end\n",
"\n",
"println(\"∞\\t3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628...\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"収束する気配が見られないので実装が間違っているのかと思いましたが, [Wikipedia](https://ja.wikipedia.org/wiki/%E3%83%A9%E3%82%A4%E3%83%97%E3%83%8B%E3%83%83%E3%83%84%E3%81%AE%E5%85%AC%E5%BC%8F)にも収束が遅いことが書いてありました. [少し調べた](https://ja.wikipedia.org/wiki/%E5%86%86%E5%91%A8%E7%8E%87%E3%81%AE%E6%AD%B4%E5%8F%B2)ところ, \n",
"2022年6月9日にGoogleの技術者・岩尾エマはるかさんが100兆桁の円周率を計算したときには, 計算はチュドノフスキー級数, 検証はBPP公式で行ったらしいことがわかりました. チュドノフスキー級数は実装が面倒くさそうなので[BPP公式](https://xn--w6q13e505b.jp/formula/bbp.html)\n",
"\n",
"$$\n",
"\\pi = \\sum_{n=0}^{\\infty} \\frac{1}{16^n}\\left(\\frac{4}{8n+1} - \\frac{2}{8n+4} - \\frac{1}{8n+5} - \\frac{1}{8n+6}\\right)\n",
"$$\n",
"\n",
"で計算します."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0\t3.1333333333333333\n",
"1\t3.1414224664224664\n",
"2\t3.1415873903465816\n",
"3\t3.1415924575674357\n",
"4\t3.1415926454603365\n",
"5\t3.141592653228088\n",
"6\t3.141592653572881\n",
"7\t3.141592653588973\n",
"8\t3.1415926535897523\n",
"9\t3.1415926535897913\n",
"10\t3.141592653589793\n",
"11\t3.141592653589793\n",
"12\t3.141592653589793\n",
"13\t3.141592653589793\n",
"14\t3.141592653589793\n",
"15\t3.141592653589793\n",
"16\tInf\n",
"17\tInf\n",
"18\tInf\n",
"19\tInf\n",
"20\tInf\n",
"∞\t3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628...\n"
]
}
],
"source": [
"x = 0\n",
"for n in 0:20\n",
" x += ( 4/(8*n+1) - 2/(8*n+4) - 1/(8*n+5) - 1/(8*n+6) ) / 16^n\n",
" println(\"$n\\t\", x)\n",
"end\n",
"println(\"∞\\t3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628...\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"ライプニッツの公式と比べると収束が早すぎて驚きますね. n=16以降で発散してしまっていますが, 恐らく`/ 16^n`の部分でFloat64として計算してしまいInfになってしまったものと思われます. 一応全ての項のnをbig(n)としておきましょう. 定数の方も不安ですが, 恐らく自動的にBigFloatになると思われます."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"n\tπ\n",
"0\t3.133333333333333333333333333333333333333333333333333333333333333333333333333315\n",
"1\t3.14142246642246642246642246642246642246642246642246642246642246642246642246642\n",
"2\t3.141587390346581523052111287405405052463875993287757993640346581523052111287413\n",
"3\t3.141592457567435381837004555057293394007389950594818748976963987105974935589556\n",
"4\t3.141592645460336319557021222442381831727406617979907186696980654491063373309567\n",
"5\t3.141592653228087534734378035536204469558528012197801934814422303215851016442912\n",
"6\t3.141592653572880827785240761895898484239065603786606461624339948713011338401089\n",
"7\t3.14159265358897270494077776717018944697112048981182286063349689713907385057788\n",
"8\t3.141592653589752275236177868398102225795024633409061087027733521554726399222753\n",
"9\t3.141592653589791146388776965910347414779015888488996772587067242388086742316342\n",
"10\t3.141592653589793129614170564041344858816452676296281615895622361248723158812323\n",
"11\t3.141592653589793232711292261930077163422606275435901151635298427143347507662168\n",
"12\t3.141592653589793238154766322501863827762609260414389714560723777440657408452123\n",
"13\t3.141592653589793238445977501940281666096938425156252904675060680768386424254539\n",
"14\t3.141592653589793238461732482037982486800056278143046732780578758091752793513451\n",
"15\t3.141592653589793238462593174670682882792683045699610435502181039811336565883137\n",
"16\t3.141592653589793238462640595138128445061235672871301070791890279329645172273369\n",
"17\t3.141592653589793238462643227424822458237094279625505676929921156072080520241485\n",
"18\t3.141592653589793238462643374515761485970237552267559842751936453791382429710676\n",
"19\t3.141592653589793238462643382784091514246623611329334708720045257099553183225298\n",
"20\t3.141592653589793238462643383251362615881909316518417908555365030283214294098105\n",
"30\t3.141592653589793238462643383279502884197157502154596455091303765610235512627314\n",
"40\t3.141592653589793238462643383279502884197169399375105814747586586693114454191957\n",
"50\t3.141592653589793238462643383279502884197169399375105820974944592304140986208374\n",
"60\t3.141592653589793238462643383279502884197169399375105820974944592307816406286233\n",
"70\t3.141592653589793238462643383279502884197169399375105820974944592307816406286233\n",
"80\t3.141592653589793238462643383279502884197169399375105820974944592307816406286233\n",
"90\t3.141592653589793238462643383279502884197169399375105820974944592307816406286233\n",
"100\t3.141592653589793238462643383279502884197169399375105820974944592307816406286233\n",
"∞\t3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628...\n"
]
}
],
"source": [
"x = big(0.0)\n",
"println(\"n\\tπ\")\n",
"for m in 0:100\n",
" n = big(m)\n",
" x += ( 4/(8*n+1) - 2/(8*n+4) - 1/(8*n+5) - 1/(8*n+6) ) / 16^n\n",
" if n in union(0:20, 20:10:100); println(\"$n\\t\", x); end\n",
"end\n",
"println(\"∞\\t3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628...\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"最後の2桁はズレてしまうようですが, ここまで綺麗に計算できると気持ちが良いですね. 余力のある方は有効桁数を増やしたり, より効率的な[Bellard公式](https://xn--w6q13e505b.jp/formula/bbp.html)で計算してみましょう.\n",
"\n",
"$$\n",
"\\pi = \\frac{1}{64} \\sum_{n=0}^{\\infty} \\frac{(-1)^n}{1024^n}\n",
" \\left(\n",
" -\\frac{32}{4n+1} - \\frac{1}{4n+3} + \\frac{256}{10n+1}\n",
" -\\frac{64}{10n+3} - \\frac{4}{10n+5} - \\frac{4}{10n+7}\n",
" + \\frac{1}{10n+9}\n",
" \\right)\n",
"$$\n",
"\n",
"\n",
"**演習**　[BPP公式](https://xn--w6q13e505b.jp/formula/bbp.html)と[Bellard公式](https://xn--w6q13e505b.jp/formula/bbp.html)それぞれを用いて円周率を100桁決定しなさい. その際に必要な反復回数を比較しなさい."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 消費メモリ\n",
"\n",
"デフォルトでのBigFloat型のメモリ少量を見積もります. 今のPCのメモリは8～16GBくらいでしょうから, 1行列あたりせいぜい1GBが限度だと思われます. とすると, メモリの制約からすれば, Float64では12000×12000, BigFloatでは6000×6000程度の行列の取り扱いが限界だと思われます. メモリよりも行列を用意する際の計算速度の方が問題になるかもしれませんね."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\tFloat32\t\tFloat64\t\tBigFloat\t\t\n",
"i\ti×i×4\t\ti×i×8\t\ti×i×32\n",
"1\t 4.000 B\t 8.000 B\t 32.000 B\n",
"10\t400.000 B\t800.000 B\t 3.200 KB\n",
"100\t 40.000 KB\t 80.000 KB\t320.000 KB\n",
"1000\t 4.000 MB\t 8.000 MB\t 32.000 MB\n",
"2000\t 16.000 MB\t 32.000 MB\t128.000 MB\n",
"3000\t 36.000 MB\t 72.000 MB\t288.000 MB\n",
"4000\t 64.000 MB\t128.000 MB\t512.000 MB\n",
"5000\t100.000 MB\t200.000 MB\t800.000 MB\n",
"6000\t144.000 MB\t288.000 MB\t 1.152 GB\n",
"7000\t196.000 MB\t392.000 MB\t 1.568 GB\n",
"8000\t256.000 MB\t512.000 MB\t 2.048 GB\n",
"9000\t324.000 MB\t648.000 MB\t 2.592 GB\n",
"10000\t400.000 MB\t800.000 MB\t 3.200 GB\n",
"11000\t484.000 MB\t968.000 MB\t 3.872 GB\n",
"12000\t576.000 MB\t 1.152 GB\t 4.608 GB\n",
"13000\t676.000 MB\t 1.352 GB\t 5.408 GB\n",
"14000\t784.000 MB\t 1.568 GB\t 6.272 GB\n",
"15000\t900.000 MB\t 1.800 GB\t 7.200 GB\n",
"16000\t 1.024 GB\t 2.048 GB\t 8.192 GB\n",
"17000\t 1.156 GB\t 2.312 GB\t 9.248 GB\n",
"18000\t 1.296 GB\t 2.592 GB\t 10.368 GB\n",
"19000\t 1.444 GB\t 2.888 GB\t 11.552 GB\n",
"20000\t 1.600 GB\t 3.200 GB\t 12.800 GB\n",
"30000\t 3.600 GB\t 7.200 GB\t 28.800 GB\n",
"40000\t 6.400 GB\t 12.800 GB\t 51.200 GB\n",
"50000\t 10.000 GB\t 20.000 GB\t 80.000 GB\n",
"60000\t 14.400 GB\t 28.800 GB\t115.200 GB\n",
"70000\t 19.600 GB\t 39.200 GB\t156.800 GB\n",
"80000\t 25.600 GB\t 51.200 GB\t204.800 GB\n",
"90000\t 32.400 GB\t 64.800 GB\t259.200 GB\n",
"100000\t 40.000 GB\t 80.000 GB\t320.000 GB\n"
]
}
],
"source": [
"using Printf\n",
"\n",
"function KMG(byte)\n",
" if byte > 10^9\n",
" return @sprintf(\"%7.3f GB\", byte/10^9)\n",
" elseif byte > 10^6\n",
" return @sprintf(\"%7.3f MB\", byte/10^6)\n",
" elseif byte > 10^3\n",
" return @sprintf(\"%7.3f KB\", byte/10^3)\n",
" else\n",
" return @sprintf(\"%7.3f B\", byte)\n",
" end\n",
"end\n",
"\n",
"x = convert(Float32, 1.0)\n",
"y = convert(Float64, 1.0)\n",
"z = convert(BigFloat, 1.0)\n",
"println(\"\\t\", typeof(x), \"\\t\\t\", typeof(y), \"\\t\\t\", typeof(z), \"\\t\\t\")\n",
"println(\"i\\ti×i×\", sizeof(x), \"\\t\\ti×i×\", sizeof(y), \"\\t\\ti×i×\", sizeof(z))\n",
"for i in union(1,10,100,1000:1000:20000, 30000:10000:100000)\n",
" println(i, \"\\t\", KMG(i^2*sizeof(x)), \"\\t\", KMG(i^2*sizeof(y)), \"\\t\", KMG(i^2*sizeof(z)))\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"仮数部分だけで256bit = 32Bなのに, 全体で32Bというのは間違いでは？という気もします. また, `sizeof()`はBigFloat型の行列に対して機能しないようなので, コミュニティで質問した方が良さそうな気がしています. 試しに`sizeof(big.(rand(1,2)))`と実行してみてください. Float64での見積もり結果が返ってきます."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(A[1, 1]) = Float64\n",
"typeof(B[1, 1]) = BigFloat\n",
"sizeof(A[1, 1]) = 8\n",
"sizeof(B[1, 1]) = 32\n",
"typeof(A) = Matrix{Float64}\n",
"typeof(B) = Matrix{BigFloat}\n",
"sizeof(A) = 32\n",
"sizeof(B) = 32\n"
]
},
{
"data": {
"text/plain": [
"32"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = rand(2,2)\n",
"B = big.(A)\n",
"\n",
"@show typeof(A[1,1])\n",
"@show typeof(B[1,1])\n",
"@show sizeof(A[1,1])\n",
"@show sizeof(B[1,1])\n",
"\n",
"@show typeof(A)\n",
"@show typeof(B)\n",
"@show sizeof(A)\n",
"@show sizeof(B)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 固有値問題\n",
"\n",
"Schrödinger方程式の近似解法として固有値問題を解くとき, 固有値だけではエネルギーしかわからないので, 固有ベクトルから期待値も計算したいわけです. そこで, BigFloat型で固有値と固有ベクトルの両方を求める方法を考えていきましょう.\n",
"\n",
"多倍精度で数千×数千の行列の固有値と固有ベクトルを全て求めるのは（定量的にどの程度の問題になるかわかりませんが）ロスが大きいと思われるので, せいぜい下から20,30個程度の固有値と固有ベクトルだけ求められれば良いです. Juliaの標準ライブラリがLAPACKのラッパーであることは有名ですが, [LAPACK.jl](https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/lapack.jl#L5122)を見る限りでは, 全ての固有値と固有ベクトルを求める[DSYGVD](https://www.nag-j.co.jp/lapack/dsyev.htm)（[DSYGV](https://www.nag-j.co.jp/lapack/dsyev.htm)の分割統治法版）のラッパーは実装されているものの, 一部の固有値・固有ベクトルを求める[DSYGVX](https://www.nag-j.co.jp/lapack/dsyevx.htm)のラッパーはそもそも用意されていないようです.（お前が書けやという話ですが）\n",
"\n",
"同様にBigFloat型でも, [GenericLinearAlgebra.jl](https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl)で固有値が計算できるようですが, 一部の固有値の計算はできません. この問題については[いろいろと検討された](https://discourse.julialang.org/t/eigenvalues-and-extended-precision/68119/9)ようですが, 解決していないようです. 固有ベクトルに関しては求めるルーチンがあるだけマシということで, よしとしましょう. 問題は固有ベクトルです. そもそも固有ベクトルを求める関数が見つからないので[柳澤 優香『任意精度演算を用いた反復改良による数値計算手法』(2019)](https://doi.org/10.11540/bjsiam.29.1_4)などを参考に実装する予定です.\n",
"\n",
"本当は[中田真秀](http://nakatamaho.riken.jp/)さんの[MPLAPACK](https://github.com/nakatamaho/mplapack)をJuliaから呼び出せると最高ですが, インターフェースがありません. Juliaのコミュニティの方でも[議論されたことはある](https://discourse.julialang.org/t/mplapack-multiprecision-blas-lapack/70514)ようですが, 現状では開発されている様子はありません.\n",
"\n",
"また, スレッドの最後に[Nemo.jl](https://fredrikj.net/blog/2018/07/high-precision-linear-algebra-in-julia-bigfloat-vs-arb/)の方が行列積が高速だということが紹介されていますが, 固有値問題への応用例が見つかりません. [Arb本体の方では実装された](https://fredrikj.net/blog/2018/12/eigenvalues-in-arb/)とのことですが, 使い方がわからなかったので今回はGenericLinearAlgebra.jlでの計算例を紹介します."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 厳密固有値\n",
"\n",
"ここでは対称行列`A`の固有値を求めていきます. ベンチマークとして, 固有値が全て整数になるような対称行列を考えました.固有値は$-5,-1,8$であり, [Wolfram Alphaによる検証](https://www.wolframalpha.com/input?i=eigenvalues+%7B%7B1%2C2%2C6%7D%2C+%7B2%2C0%2C2%7D%2C+%7B6%2C2%2C1%7D%7D)を利用したり, $|A-\\lambda\\boldsymbol{I}|=0$を計算することで検証できます."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"det(big.(A) - 8I) = 0\n",
"det(big.(A) + 5I) = 0\n",
"det(big.(A) + 1I) = 0\n"
]
},
{
"data": {
"text/plain": [
"0"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using LinearAlgebra\n",
"\n",
"A = [ 1 2 6;\n",
" 2 0 2;\n",
" 6 2 1]\n",
"\n",
"@show det(big.(A) - 8*I)\n",
"@show det(big.(A) + 5*I)\n",
"@show det(big.(A) + 1*I)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 近似固有値\n",
"\n",
"LinearAlgebra.jlによる[固有値問題](https://zenn.dev/ohno/articles/cb10dc5b3f5bbc)の計算例についてはこちらを参照してください. "
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"eigvals(Symmetric(A)) = [-5.000000000000001, -0.9999999999999996, 8.0]\n"
]
},
{
"data": {
"text/plain": [
"3-element Vector{Float64}:\n",
" -5.000000000000001\n",
" -0.9999999999999996\n",
" 8.0"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using LinearAlgebra\n",
"\n",
"A = [ 1 2 6;\n",
" 2 0 2;\n",
" 6 2 1]\n",
"\n",
"@show eigvals(Symmetric(A))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"上記の例をBigFloatに対応させます. [GenericLinearAlgebra.jl](https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl/tree/d16065d74e18b2633e8479cb9a9928f1e786300f)を使うと`eigvals()`がBigFloatの行列にも対応します."
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"eigvals(big.(A)) = BigFloat[-5.0, -1.0, 7.999999999999999999999999999999999999999999999999999999999999999999999999999862]\n"
]
},
{
"data": {
"text/plain": [
"3-element Vector{BigFloat}:\n",
" -5.0\n",
" -1.0\n",
" 7.999999999999999999999999999999999999999999999999999999999999999999999999999862"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# using Pkg\n",
"# Pkg.add(\"GenericLinearAlgebra\")\n",
"using GenericLinearAlgebra\n",
"\n",
"A = [ 1 2 6;\n",
" 2 0 2;\n",
" 6 2 1]\n",
"\n",
"@show eigvals(big.(A))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"LinearAlgebra.jlでは`Symmetric()`など行列の種類を渡す仕様でしたが, GenericLinearAlgebra.jlでは機能しませんでした. 自動で判定されているのでしょうか. `Symmetric{BigInt, Matrix{BigInt}}`を渡すと機能しません."
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(Symmetric(A)) = Symmetric{Int64, Matrix{Int64}}\n",
"typeof(Symmetric(big.(A))) = Symmetric{BigInt, Matrix{BigInt}}\n",
"typeof(big.(Symmetric(A))) = Matrix{BigInt}\n"
]
},
{
"data": {
"text/plain": [
"Matrix{BigInt} (alias for Array{BigInt, 2})"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@show typeof(Symmetric(A))\n",
"@show typeof(Symmetric(big.(A)))\n",
"@show typeof(big.(Symmetric(A)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## まとめ\n",
"\n",
"- BigInt型, BigFloat型の基本的な使い方と型のルールについて確認しました.\n",
"- Int64型では計算できない階乗を実際に計算できることを示しました.\n",
"- JuliaでBigInt型を用いて円周率を計算する例を示しました.\n",
"- BigInt型の行列のメモリ消費量を見積もりました.\n",
"- BigFloat型の行列に対する固有値問題について議論しました."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 参考文献\n",
"\n",
"- https://docs.julialang.org/en/v1/manual/integers-and-floating-point-numbers/#Arbitrary-Precision-Arithmetic\n",
"- https://docs.julialang.org/en/v1/base/numbers/#BigFloats-and-BigInts\n",
"- https://jutho.github.io/KrylovKit.jl/dev/man/eig/#KrylovKit.geneigsolve"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.7.1",
"language": "julia",
"name": "julia-1.7"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.7.1"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment