Zenn
🍜

Juliaで任意精度演算に入門してみた

2022/08/25に公開
Julia
Jupyter Notebook
#
任意精度演算
tech

Juliaでは任意精度演算が標準でサポートされています. 特別なライブラリなども不要なあたりに次世代の計算科学を担っていく覚悟と矜持を感じますね. バックエンドではBigFloat型にはthe GNU MPFR libraryが, BigInt型にはthe GNU Multiple Precision Arithmetic Library (GMP)が使われているようです.

環境

入力
versioninfo()
出力
Julia Version 1.7.1
Commit ac5cc99908 (2021-12-22 19:35 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-4650U CPU @ 1.70GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, haswell)

基本的な使い方

デフォルトでは, 整数はInt64型, 小数はFloat64型になります. Float64型はいわゆる倍精度浮動小数点数というものです. b = big(2)とすると, BigInt型になります.

入力
a = 2
b = 2.0
c = big(2)
d = big(2.0)

@show typeof(a)
@show typeof(b)
@show typeof(c)
@show typeof(d)

@show a
@show b
@show c
@show d
出力
typeof(a) = Int64
typeof(b) = Float64
typeof(c) = BigInt
typeof(d) = BigFloat
a = 2
b = 2.0
c = 2
d = 2.0

標準で用意されているほとんどの関数は, 引数にBigIntまたはBigFloatを渡すと, 戻り値にはBigIntまたはBigFloatが返ってきます. 例えば\sqrt{2}を計算してみましょう.

入力
a = sqrt(2)
b = sqrt(big(2))

@show typeof(a)
@show typeof(b)

@show a
@show b
出力
typeof(a) = Float64
typeof(b) = BigFloat
a = 1.4142135623730951
b = 1.414213562373095048801688724209698078569671875376948073176679737990732478462102

sqrt(big(2))の戻り値をこちらのサイトと比較すると

  • 1.414213562373095048801688724209698078569671875376948073176679737990732478462102
  • 1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850...

となっており, 最後の1桁だけズレているものの, 78桁は一致しています. すごい. これはまだまだ序の口です.

有効桁数

有効桁数は任意の長さに変えられるとこのことです. まず, デフォルトでは

入力
a = 2.0
b = big(2.0)

@show typeof(a)
@show typeof(b)

@show precision(a)
@show precision(b)

@show sqrt(a)
@show sqrt(b)
出力
typeof(a) = Float64
typeof(b) = BigFloat
precision(a) = 53
precision(b) = 256
sqrt(a) = 1.4142135623730951
sqrt(b) = 1.414213562373095048801688724209698078569671875376948073176679737990732478462102

です. 例えばこれをsetpresition(BigFloat, 53)とすると, 仮数に53ビット使われますので, BigFloat型がFloat64型と同じ有効桁数になります.（確かにそうですが, 基数と指数が結果に影響しないか確かめてください. いずれにせよ型によって別のアルゴリズムが使われる場合は結果が異なることがあります. 例えばLinearAlgebra.jlはBogFloatに対応していないので, GenericLinearAlgebra.jlが使われ, 誤差の入り方も異なります.） 浮動小数点数の仮数・基数・指数についてはWikipediaを参照してください.

入力
setprecision(BigFloat, 53)

a = 2.0
b = big(2.0)

@show typeof(a)
@show typeof(b)

@show precision(a)
@show precision(b)

@show sqrt(a)
@show sqrt(b)
出力
typeof(a) = Float64
typeof(b) = BigFloat
precision(a) = 53
precision(b) = 53
sqrt(a) = 1.4142135623730951
sqrt(b) = 1.4142135623730951

setpresition(BigFloat, ?)で仮数の長さを変更すると, 関数の戻り値の桁数も変わります. グローバルに変更されるので, 並列計算の途中で変えるといった使い方はやめた方がよいです. 例えば 2^13bit = 8192bit = 1024B = 1KB を仮数部に使って\sqrt{2}を先ほどよりも正確に求めてみましょう.

入力
setprecision(BigFloat, 2^13)
@show sqrt(big(2))
出力
sqrt(big(2)) = 1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462309122970249248360558507372126441214970999358314132226659275055927557999505011527820605714701095599716059702745345968620147285174186408891986095523292304843087143214508397626036279952514079896872533965463318088296406206152583523950547457502877599617298355752203375318570113543746034084988471603868999706990048150305440277903164542478230684929369186215805784631115966687130130156185689872372352885092648612494977154218334204285686060146824720771435854874155657069677653720226485447015858801620758474922657226002085584466521458398893944370926591800311388246468157082630100594858704003186480342194897278290641045072636881313739855256117322040245091227700226941127573627280495738108967504018369868368450725799364729060762996941380475654823728997180326802474420629269124859052181004459842150591120249441341728531478105803603371077309182869314710171111683916581726889419758716582152128229518488472089694633862891562882765952635140542267653239694617511291602408715510135150455381287560052631468017127402653969470240300517495318862925631385188163478001569369176881852378684052287837629389214300655869568685964595155501644724509836896036887323114389415576651040883914292338113206052433629485317049915771756228549741438999188021762430965206564211827316726257539594717255934637238632261482742622208671155839599926521176252698917540988159348640083457085181472231814204070426509056532333398436457865796796519267292399875366617215982578860263363617827495994219403777753681426217738799194551397231274066898329989895386728822856378697749662519966583525776198939322845344735694794962952168891485492538904755828834526096524096542889394538646625744927556381964410316979833061852019379384940057156333720548068540575867999670121372239475821426306585132217408832382947287617393647467837431960001592188807347857617252211867490424977366929207311096369721608933708661156734585334833295254675851644710757848602463600834449114818587655554286455123314219926311332517970608436559704352856410087918500760361009159465670676883605571740076756905096136719401324935605240185999105062108163597726431380605467010293569971042425105781749531057255934984451126922780344913506637568747760283162829605532422426957534529028838768446429173282770888318087025339852338122749990812371892540726475367850304821591801886167108972869229201197599880703818543332536460211082299279293072871780799888099176741771

こちらのサイトと比較すると, 最後の1桁だけズレていますが, 小数点以下2466桁まで正しく計算できていることがわかります. 最後の数桁をページ内検索すると確かめられます. ちなみに, setprecision(BigFloat, 2^20)として同様の計算を行うと最後の1桁だけズレますが小数点以下315652桁まで正しく計算できます.

デフォルトのsetprecision(BigFloat, 256)に戻しておきます.

入力
setprecision(BigFloat, 256)
@show precision(big(2.0))
@show sqrt(big(2.0))
出力
precision(big(2.0)) = 256
sqrt(big(2.0)) = 1.414213562373095048801688724209698078569671875376948073176679737990732478462102

型の規則

それぞれの型の和を考えるとき, Int＋Float＝Floatとなり, 64＋Big＝Bigとなります.

Int64 Float64 BigInt BigFloat
＋Int64 Int64 Float64 BigInt BigFloat
＋Float64 Float64 Float64 BigFloat BigFloat
＋BigInt BigInt BigFloat BigInt BigFloat
＋BigFloat BigFloat BigFloat BigFloat BigFloat

それぞれの型の積を考えるときも同じく, Int×Float＝Floatとなり, 64×Big＝Bigとなります.

Int64 Float64 BigInt BigFloat
×Int64 Int64 Float64 BigInt BigFloat
×Float64 Float64 Float64 BigFloat BigFloat
×BigInt BigInt BigFloat BigInt BigFloat
×BigFloat BigFloat BigFloat BigFloat BigFloat

検証：

入力
a = 2
b = 2.0
c = big(2)
d = big(2.0)

@show typeof(a)
@show typeof(b)
@show typeof(c)
@show typeof(d)

@show a
@show b
@show c
@show d

@show typeof(a+a)
@show typeof(a+b)
@show typeof(a+c)
@show typeof(a+d)

@show typeof(b+a)
@show typeof(b+b)
@show typeof(b+c)
@show typeof(b+d)

@show typeof(c+a)
@show typeof(c+b)
@show typeof(c+c)
@show typeof(c+d)

@show typeof(d+a)
@show typeof(d+b)
@show typeof(d+c)
@show typeof(d+d)

@show typeof(a*a)
@show typeof(a*b)
@show typeof(a*c)
@show typeof(a*d)

@show typeof(b*a)
@show typeof(b*b)
@show typeof(b*c)
@show typeof(b*d)

@show typeof(c*a)
@show typeof(c*b)
@show typeof(c*c)
@show typeof(c*d)

@show typeof(d*a)
@show typeof(d*b)
@show typeof(d*c)
@show typeof(d*d)
出力
typeof(a) = Int64
typeof(b) = Float64
typeof(c) = BigInt
typeof(d) = BigFloat
a = 2
b = 2.0
c = 2
d = 2.0
typeof(a + a) = Int64
typeof(a + b) = Float64
typeof(a + c) = BigInt
typeof(a + d) = BigFloat
typeof(b + a) = Float64
typeof(b + b) = Float64
typeof(b + c) = BigFloat
typeof(b + d) = BigFloat
typeof(c + a) = BigInt
typeof(c + b) = BigFloat
typeof(c + c) = BigInt
typeof(c + d) = BigFloat
typeof(d + a) = BigFloat
typeof(d + b) = BigFloat
typeof(d + c) = BigFloat
typeof(d + d) = BigFloat
typeof(a * a) = Int64
typeof(a * b) = Float64
typeof(a * c) = BigInt
typeof(a * d) = BigFloat
typeof(b * a) = Float64
typeof(b * b) = Float64
typeof(b * c) = BigFloat
typeof(b * d) = BigFloat
typeof(c * a) = BigInt
typeof(c * b) = BigFloat
typeof(c * c) = BigInt
typeof(c * d) = BigFloat
typeof(d * a) = BigFloat
typeof(d * b) = BigFloat
typeof(d * c) = BigFloat
typeof(d * d) = BigFloat

下記の3つの積だけを見ると分かりやすいですが, 同じBigFloatであっても, Float64×BigFloatでは有効桁数はFloat64に引っ張られるので, 精度を保ちたい場面ではBigFloatだけで計算を行うことをお勧めします.

入力
a = sqrt(2)
b = sqrt(big(2))

@show typeof(a)
@show typeof(b)

@show a
@show b

@show a*a
@show a*b
@show b*b
出力
typeof(a) = Float64
typeof(b) = BigFloat
a = 1.4142135623730951
b = 1.414213562373095048801688724209698078569671875376948073176679737990732478462102
a * a = 2.0000000000000004
a * b = 2.000000000000000136716173153238459361596582123216045722514805262001562413306807
b * b = 1.999999999999999999999999999999999999999999999999999999999999999999999999999983

応用例：階乗

Int64型で扱える最大の整数は9223372036854775807です.

入力
@show typemax(Int64)
出力
typemax(Int64) = 9223372036854775807

階乗n!はnに対して爆発的に大きな値を取るため, Int64型の限界がすぐに見れます.

入力
for n in 1:60
    println("$(n)!\t=\t", factorial(n))
end
出力
1!	=	1
2!	=	2
3!	=	6
4!	=	24
5!	=	120
6!	=	720
7!	=	5040
8!	=	40320
9!	=	362880
10!	=	3628800
11!	=	39916800
12!	=	479001600
13!	=	6227020800
14!	=	87178291200
15!	=	1307674368000
16!	=	20922789888000
17!	=	355687428096000
18!	=	6402373705728000
19!	=	121645100408832000
20!	=	2432902008176640000
OverflowError: 21 is too large to look up in the table; consider using `factorial(big(21))` instead

このように通常の整数だと21!が計算できずにエラーが出てしまいますが, 以下のように factorial(n)を factorial(big(n))に置き換えることでより大きい数値を計算できるようになります.

入力
for n in 1:60
    println("$(n)!\t=\t", factorial(big(n)))
end
出力
1!	=	1
2!	=	2
3!	=	6
4!	=	24
5!	=	120
6!	=	720
7!	=	5040
8!	=	40320
9!	=	362880
10!	=	3628800
11!	=	39916800
12!	=	479001600
13!	=	6227020800
14!	=	87178291200
15!	=	1307674368000
16!	=	20922789888000
17!	=	355687428096000
18!	=	6402373705728000
19!	=	121645100408832000
20!	=	2432902008176640000
21!	=	51090942171709440000
22!	=	1124000727777607680000
23!	=	25852016738884976640000
24!	=	620448401733239439360000
25!	=	15511210043330985984000000
26!	=	403291461126605635584000000
27!	=	10888869450418352160768000000
28!	=	304888344611713860501504000000
29!	=	8841761993739701954543616000000
30!	=	265252859812191058636308480000000
31!	=	8222838654177922817725562880000000
32!	=	263130836933693530167218012160000000
33!	=	8683317618811886495518194401280000000
34!	=	295232799039604140847618609643520000000
35!	=	10333147966386144929666651337523200000000
36!	=	371993326789901217467999448150835200000000
37!	=	13763753091226345046315979581580902400000000
38!	=	523022617466601111760007224100074291200000000
39!	=	20397882081197443358640281739902897356800000000
40!	=	815915283247897734345611269596115894272000000000
41!	=	33452526613163807108170062053440751665152000000000
42!	=	1405006117752879898543142606244511569936384000000000
43!	=	60415263063373835637355132068513997507264512000000000
44!	=	2658271574788448768043625811014615890319638528000000000
45!	=	119622220865480194561963161495657715064383733760000000000
46!	=	5502622159812088949850305428800254892961651752960000000000
47!	=	258623241511168180642964355153611979969197632389120000000000
48!	=	12413915592536072670862289047373375038521486354677760000000000
49!	=	608281864034267560872252163321295376887552831379210240000000000
50!	=	30414093201713378043612608166064768844377641568960512000000000000
51!	=	1551118753287382280224243016469303211063259720016986112000000000000
52!	=	80658175170943878571660636856403766975289505440883277824000000000000
53!	=	4274883284060025564298013753389399649690343788366813724672000000000000
54!	=	230843697339241380472092742683027581083278564571807941132288000000000000
55!	=	12696403353658275925965100847566516959580321051449436762275840000000000000
56!	=	710998587804863451854045647463724949736497978881168458687447040000000000000
57!	=	40526919504877216755680601905432322134980384796226602145184481280000000000000
58!	=	2350561331282878571829474910515074683828862318181142924420699914240000000000000
59!	=	138683118545689835737939019720389406345902876772687432540821294940160000000000000
60!	=	8320987112741390144276341183223364380754172606361245952449277696409600000000000000

このように, BigInt型を使うことでInt64型の最大値9223372036854775807を超えて, 21! = 51090942171709440000を計算することができました. では一体, どこまで大きな数を扱えるのか各自で試してみましょう. 1000!の例を示します.

入力
n = 1000
println("$(n)!\t=\t", factorial(big(n)))
出力
1000!	=	402387260077093773543702433923003985719374864210714632543799910429938512398629020592044208486969404800479988610197196058631666872994808558901323829669944590997424504087073759918823627727188732519779505950995276120874975462497043601418278094646496291056393887437886487337119181045825783647849977012476632889835955735432513185323958463075557409114262417474349347553428646576611667797396668820291207379143853719588249808126867838374559731746136085379534524221586593201928090878297308431392844403281231558611036976801357304216168747609675871348312025478589320767169132448426236131412508780208000261683151027341827977704784635868170164365024153691398281264810213092761244896359928705114964975419909342221566832572080821333186116811553615836546984046708975602900950537616475847728421889679646244945160765353408198901385442487984959953319101723355556602139450399736280750137837615307127761926849034352625200015888535147331611702103968175921510907788019393178114194545257223865541461062892187960223838971476088506276862967146674697562911234082439208160153780889893964518263243671616762179168909779911903754031274622289988005195444414282012187361745992642956581746628302955570299024324153181617210465832036786906117260158783520751516284225540265170483304226143974286933061690897968482590125458327168226458066526769958652682272807075781391858178889652208164348344825993266043367660176999612831860788386150279465955131156552036093988180612138558600301435694527224206344631797460594682573103790084024432438465657245014402821885252470935190620929023136493273497565513958720559654228749774011413346962715422845862377387538230483865688976461927383814900140767310446640259899490222221765904339901886018566526485061799702356193897017860040811889729918311021171229845901641921068884387121855646124960798722908519296819372388642614839657382291123125024186649353143970137428531926649875337218940694281434118520158014123344828015051399694290153483077644569099073152433278288269864602789864321139083506217095002597389863554277196742822248757586765752344220207573630569498825087968928162753848863396909959826280956121450994871701244516461260379029309120889086942028510640182154399457156805941872748998094254742173582401063677404595741785160829230135358081840096996372524230560855903700624271243416909004153690105933983835777939410970027753472000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

また, 少なくとも10000!までは計算できました. もはや正しいのかわからないですね. 検証法を調べた方や思いついた方がいたら教えてください. なお, BigIntではBigFloatとは違って, 勝手に有効桁数がが変わっているのではないかと思われますが, こちらも検証方法がわかりませんでした.

応用例：円周率

こちらのサイトに載っている値と比較すると, デフォルトの桁数のBigFloat型で展開した円周率πは小数点以下75桁までは正しく, 最後の2,3桁がズレます.

入力
@show big(π)
println("exact   = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348...")
出力
big(π) = 3.141592653589793238462643383279502884197169399375105820974944592307816406286198
exact   = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348...

もちろん, 精度を上げることはできますが, どんな桁数にしても最後の2,3桁はズレると思った方が良いです. 下記の例では最後の1桁だけズレていますが, 小数点以下2465桁までは合っています.

入力
@show BigFloat(π, precision=2^13)
出力
BigFloat(π, precision = 2 ^ 13) = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420198938095257201065485863278865936153381827968230301952035301852968995773622599413891249721775283479131515574857242454150695950829533116861727855889075098381754637464939319255060400927701671139009848824012858361603563707660104710181942955596198946767837449448255379774726847104047534646208046684259069491293313677028989152104752162056966024058038150193511253382430035587640247496473263914199272604269922796782354781636009341721641219924586315030286182974555706749838505494588586926995690927210797509302955321165344987202755960236480665499119881834797753566369807426542527862551818417574672890977772793800081647060016145249192173217214772350141441973568548161361157352552133475741849468438523323907394143334547762416862518983569485562099219222184272550254256887671790494601653466804988627232791786085784383827967976681454100953883786360950680064225125205117392984896084128488626945604241965285022210661186306744278622039194945047123713786960956364371917287467764657573962413890865832645995813390478027590099465764078951269468398352595709825822620522489407726719478268482601476990902640136394437455305068203496252451749399651431429809190659250937221696461515709858387410597885959772975498930161753928468138268683868942774155991855925245953959431049972524680845987273644695848653836736222626099124608051243884390451244136549762780797715691435997700129616089441694868555848406353422072225828488648158456028506016842739452267467678895252138522549954666727823986456596116354888

ライプニッツの公式を使って, 円周率\piの計算を実装してみましょう.

\pi = \sum_{n=0}^\infty \frac{4(-1)^n}{2n+1}
入力
x = 0.0

for n in 0:100000
    x += 4 * (-1)^n / (2*n+1)
    if n in union(0:10, 1000:1000:10000); println("$n\t", x); end
end

println("∞\t3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628...")
出力
0	4.0
1	2.666666666666667
2	3.466666666666667
3	2.8952380952380956
4	3.3396825396825403
5	2.9760461760461765
6	3.2837384837384844
7	3.017071817071818
8	3.2523659347188767
9	3.0418396189294032
10	3.232315809405594
1000	3.1425916543395442
2000	3.1420924036835256
3000	3.1419258758397897
4000	3.1418425911015158
5000	3.1417926135957908
6000	3.1417592924821482
7000	3.141735490326666
8000	3.1417176379662446
9000	3.1417037523562383
10000	3.1416926435905346
∞	3.14159265358979323846264...

収束する気配が見られないので実装が間違っているのかと思いましたが, Wikipediaにも収束が遅いことが書いてありました. 少し調べたところ,
2022年6月9日にGoogleの技術者・岩尾エマはるかさんが100兆桁の円周率を計算したときには, 計算はチュドノフスキー級数, 検証はBPP公式で行ったらしいことがわかりました. チュドノフスキー級数は実装が面倒くさそうなのでBPP公式

\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)

で計算します.

入力
x = 0
for n in 0:20
    x += ( 4/(8*n+1) - 2/(8*n+4) - 1/(8*n+5) - 1/(8*n+6) ) / 16^n
    println("$n\t", x)
end
println("∞\t3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628...")
出力
0	3.1333333333333333
1	3.1414224664224664
2	3.1415873903465816
3	3.1415924575674357
4	3.1415926454603365
5	3.141592653228088
6	3.141592653572881
7	3.141592653588973
8	3.1415926535897523
9	3.1415926535897913
10	3.141592653589793
11	3.141592653589793
12	3.141592653589793
13	3.141592653589793
14	3.141592653589793
15	3.141592653589793
16	Inf
17	Inf
18	Inf
19	Inf
20	Inf
∞	3.14159265358979323846264...

ライプニッツの公式と比べると収束が早すぎて驚きますね. n=16以降で発散してしまっていますが, 恐らく/ 16^nの部分でFloat64として計算してしまいInfになってしまったものと思われます. 一応全ての項のnをbig(n)としておきましょう. 定数の方も不安ですが, 恐らく自動的にBigFloatになると思われます.

入力
x = big(0.0)
println("n\tπ")
for m in 0:100
    n = big(m)
    x += ( 4/(8*n+1) - 2/(8*n+4) - 1/(8*n+5) - 1/(8*n+6) ) / 16^n
    if n in union(0:20, 20:10:100); println("$n\t", x); end
end
println("∞\t3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628...")
出力
n	π
0	3.133333333333333333333333333333333333333333333333333333333333333333333333333315
1	3.14142246642246642246642246642246642246642246642246642246642246642246642246642
2	3.141587390346581523052111287405405052463875993287757993640346581523052111287413
3	3.141592457567435381837004555057293394007389950594818748976963987105974935589556
4	3.141592645460336319557021222442381831727406617979907186696980654491063373309567
5	3.141592653228087534734378035536204469558528012197801934814422303215851016442912
6	3.141592653572880827785240761895898484239065603786606461624339948713011338401089
7	3.14159265358897270494077776717018944697112048981182286063349689713907385057788
8	3.141592653589752275236177868398102225795024633409061087027733521554726399222753
9	3.141592653589791146388776965910347414779015888488996772587067242388086742316342
10	3.141592653589793129614170564041344858816452676296281615895622361248723158812323
11	3.141592653589793232711292261930077163422606275435901151635298427143347507662168
12	3.141592653589793238154766322501863827762609260414389714560723777440657408452123
13	3.141592653589793238445977501940281666096938425156252904675060680768386424254539
14	3.141592653589793238461732482037982486800056278143046732780578758091752793513451
15	3.141592653589793238462593174670682882792683045699610435502181039811336565883137
16	3.141592653589793238462640595138128445061235672871301070791890279329645172273369
17	3.141592653589793238462643227424822458237094279625505676929921156072080520241485
18	3.141592653589793238462643374515761485970237552267559842751936453791382429710676
19	3.141592653589793238462643382784091514246623611329334708720045257099553183225298
20	3.141592653589793238462643383251362615881909316518417908555365030283214294098105
30	3.141592653589793238462643383279502884197157502154596455091303765610235512627314
40	3.141592653589793238462643383279502884197169399375105814747586586693114454191957
50	3.141592653589793238462643383279502884197169399375105820974944592304140986208374
60	3.141592653589793238462643383279502884197169399375105820974944592307816406286233
70	3.141592653589793238462643383279502884197169399375105820974944592307816406286233
80	3.141592653589793238462643383279502884197169399375105820974944592307816406286233
90	3.141592653589793238462643383279502884197169399375105820974944592307816406286233
100	3.141592653589793238462643383279502884197169399375105820974944592307816406286233
∞	3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628...

最後の2桁はズレてしまうようですが, ここまで綺麗に計算できると気持ちが良いですね. 余力のある方は有効桁数を増やしたり, より効率的なBellard公式で計算してみましょう.

\pi = \frac{1}{64} \sum_{n=0}^{\infty} \frac{(-1)^n}{1024^n} \left( -\frac{32}{4n+1} - \frac{1}{4n+3} + \frac{256}{10n+1} -\frac{64}{10n+3} - \frac{4}{10n+5} - \frac{4}{10n+7} + \frac{1}{10n+9} \right)

演習　BPP公式とBellard公式それぞれを用いて円周率を100桁決定しなさい. その際に必要な反復回数を比較しなさい.

消費メモリ

デフォルトでのBigFloat型のメモリ少量を見積もります. 今のPCのメモリは8～16GBくらいでしょうから, 1行列あたりせいぜい1GBが限度だと思われます. とすると, メモリの制約からすれば, Float64では12000×12000, BigFloatでは6000×6000程度の行列の取り扱いが限界だと思われます. メモリよりも行列を用意する際の計算速度の方が問題になるかもしれませんね.

入力
using Printf

function KMG(byte)
    if byte > 10^9
        return @sprintf("%7.3f GB", byte/10^9)
    elseif byte > 10^6
        return @sprintf("%7.3f MB", byte/10^6)
    elseif byte > 10^3
        return @sprintf("%7.3f KB", byte/10^3)
    else
        return @sprintf("%7.3f B", byte)
    end
end

x = convert(Float32, 1.0)
y = convert(Float64, 1.0)
z = convert(BigFloat, 1.0)
println("\t", typeof(x), "\t\t", typeof(y), "\t\t", typeof(z), "\t\t")
println("i\ti×i×", sizeof(x), "\t\ti×i×", sizeof(y), "\t\ti×i×", sizeof(z))
for i in union(1,10,100,1000:1000:20000, 30000:10000:100000)
    println(i, "\t", KMG(i^2*sizeof(x)), "\t", KMG(i^2*sizeof(y)), "\t", KMG(i^2*sizeof(z)))
end
出力
  Float32		Float64		BigFloat		
i	i×i×4		i×i×8		i×i×32
1	  4.000 B	  8.000 B	 32.000 B
10	400.000 B	800.000 B	  3.200 KB
100	 40.000 KB	 80.000 KB	320.000 KB
1000	  4.000 MB	  8.000 MB	 32.000 MB
2000	 16.000 MB	 32.000 MB	128.000 MB
3000	 36.000 MB	 72.000 MB	288.000 MB
4000	 64.000 MB	128.000 MB	512.000 MB
5000	100.000 MB	200.000 MB	800.000 MB
6000	144.000 MB	288.000 MB	  1.152 GB
7000	196.000 MB	392.000 MB	  1.568 GB
8000	256.000 MB	512.000 MB	  2.048 GB
9000	324.000 MB	648.000 MB	  2.592 GB
10000	400.000 MB	800.000 MB	  3.200 GB
11000	484.000 MB	968.000 MB	  3.872 GB
12000	576.000 MB	  1.152 GB	  4.608 GB
13000	676.000 MB	  1.352 GB	  5.408 GB
14000	784.000 MB	  1.568 GB	  6.272 GB
15000	900.000 MB	  1.800 GB	  7.200 GB
16000	  1.024 GB	  2.048 GB	  8.192 GB
17000	  1.156 GB	  2.312 GB	  9.248 GB
18000	  1.296 GB	  2.592 GB	 10.368 GB
19000	  1.444 GB	  2.888 GB	 11.552 GB
20000	  1.600 GB	  3.200 GB	 12.800 GB
30000	  3.600 GB	  7.200 GB	 28.800 GB
40000	  6.400 GB	 12.800 GB	 51.200 GB
50000	 10.000 GB	 20.000 GB	 80.000 GB
60000	 14.400 GB	 28.800 GB	115.200 GB
70000	 19.600 GB	 39.200 GB	156.800 GB
80000	 25.600 GB	 51.200 GB	204.800 GB
90000	 32.400 GB	 64.800 GB	259.200 GB
100000	 40.000 GB	 80.000 GB	320.000 GB

仮数部分だけで256bit = 32Bなのに, 全体で32Bというのは間違いでは？という気もします. また, sizeof()はBigFloat型の行列に対して機能しないようなので, コミュニティで質問した方が良さそうな気がしています. 試しにsizeof(big.(rand(1,2)))と実行してみてください. Float64での見積もり結果が返ってきます.

入力
A = rand(2,2)
B = big.(A)

@show typeof(A[1,1])
@show typeof(B[1,1])
@show sizeof(A[1,1])
@show sizeof(B[1,1])

@show typeof(A)
@show typeof(B)
@show sizeof(A)
@show sizeof(B)
出力
typeof(A[1, 1]) = Float64
typeof(B[1, 1]) = BigFloat
sizeof(A[1, 1]) = 8
sizeof(B[1, 1]) = 32
typeof(A) = Matrix{Float64}
typeof(B) = Matrix{BigFloat}
sizeof(A) = 32
sizeof(B) = 32

固有値問題

Schrödinger方程式の近似解法として固有値問題を解くとき, 固有値だけではエネルギーしかわからないので, 固有ベクトルから期待値も計算したいわけです. そこで, BigFloat型で固有値と固有ベクトルの両方を求める方法を考えていきましょう.

多倍精度で数千×数千の行列の固有値と固有ベクトルを全て求めるのは（定量的にどの程度の問題になるかわかりませんが）ロスが大きいと思われるので, せいぜい下から20,30個程度の固有値と固有ベクトルだけ求められれば良いです. Juliaの標準ライブラリがLAPACKのラッパーであることは有名ですが, LAPACK.jlを見る限りでは, 全ての固有値と固有ベクトルを求めるDSYGVD（DSYGVの分割統治法版）のラッパーは実装されているものの, 一部の固有値・固有ベクトルを求めるDSYGVXのラッパーはそもそも用意されていないようです.（お前が書けやという話ですが）

同様にBigFloat型でも, GenericLinearAlgebra.jlで固有値が計算できるようですが, 一部の固有値の計算はできません. この問題についてはいろいろと検討されたようですが, 解決していないようです. 固有ベクトルに関しては求めるルーチンがあるだけマシということで, よしとしましょう. 問題は固有ベクトルです. そもそも固有ベクトルを求める関数が見つからないので柳澤 優香『任意精度演算を用いた反復改良による数値計算手法』(2019)などを参考に実装する予定です.

本当は中田真秀さんのMPLAPACKをJuliaから呼び出せると最高ですが, インターフェースがありません. Juliaのコミュニティの方でも議論されたことはあるようですが, 現状では開発されている様子はありません.

また, スレッドの最後にNemo.jlの方が行列積が高速だということが紹介されていますが, 固有値問題への応用例が見つかりません. Arb本体の方では実装されたとのことですが, 使い方がわからなかったので今回はGenericLinearAlgebra.jlでの計算例を紹介します.

厳密固有値

ここでは対称行列Aの固有値を求めていきます. ベンチマークとして, 固有値が全て整数になるような対称行列を考えました.固有値は-5,-1,8であり, Wolfram Alphaによる検証を利用したり, |A-\lambda\boldsymbol{I}|=0を計算することで検証できます.

入力
using LinearAlgebra

A = [ 1  2  6;
      2  0  2;
      6  2  1]

@show det(big.(A) - 8*I)
@show det(big.(A) + 5*I)
@show det(big.(A) + 1*I)
出力
det(big.(A) - 8I) = 0
det(big.(A) + 5I) = 0
det(big.(A) + 1I) = 0

近似固有値

LinearAlgebra.jlによる固有値問題の計算例についてはこちらを参照してください.

入力
using LinearAlgebra

A = [ 1  2  6;
      2  0  2;
      6  2  1]

eigvals(Symmetric(A))
出力 3-element Vector{Float64}
-5.000000000000001
-0.9999999999999996
 8.0

上記の例をBigFloatに対応させます. GenericLinearAlgebra.jlを使うとeigvals()がBigFloatの行列にも対応します.

入力
# using Pkg
# Pkg.add("GenericLinearAlgebra")
using GenericLinearAlgebra

A = [ 1  2  6;
      2  0  2;
      6  2  1]

eigvals(big.(A))
出力 3-element Vector{BigFloat}
-5.0
-1.0
 7.999999999999999999999999999999999999999999999999999999999999999999999999999862

LinearAlgebra.jlではSymmetric()など行列の種類を渡す仕様でしたが, GenericLinearAlgebra.jlでは機能しませんでした. 自動で判定されているのでしょうか. Symmetric{BigInt, Matrix{BigInt}}を渡すと機能しません.

入力
@show typeof(Symmetric(A))
@show typeof(Symmetric(big.(A)))
@show typeof(big.(Symmetric(A)))
出力
typeof(Symmetric(A)) = Symmetric{Int64, Matrix{Int64}}
typeof(Symmetric(big.(A))) = Symmetric{BigInt, Matrix{BigInt}}
typeof(big.(Symmetric(A))) = Matrix{BigInt}

まとめ

  • BigInt型, BigFloat型の基本的な使い方と型のルールについて確認しました.
  • Int64型では計算できない階乗を実際に計算できることを示しました.
  • JuliaでBigInt型を用いて円周率を計算する例を示しました.
  • BigInt型の行列のメモリ消費量を見積もりました.
  • BigFloat型の行列に対する固有値問題について議論しました.

参考文献

この記事の元になったノートはGistにアップロードしてあります.
https://gist.github.com/ohno/e4d656aec716a71a741fff7ae23b9058

Discussion