前回の記事は有理数計算のサンブルだったので、今回は多倍長整数の計算。
有理数の計算で、オーバーフローに対するガードがないことを不満点としてあげたが、既約分数のべき乗を計算したりすると分母や分子がLongの範囲を超えてしまうのは仕方がない。4バイトとか8バイトといったサイズの制約を超えるためには、配列A(i)を用意して Σ A(i)・N^i という形で整数を表現すればいいだけなので簡単だ。Nを 2^15 以下にしておけば掛け算する時の一時的なオーバーフローの心配もない。
今回もデータ構造をクラス化したりせず、単純な配列で値を保持した。関数一覧にある機能はだいたい明らかだと思う。
b1 = str2bigInt("0xABCDEF0123456789") ' 文字列から構築
?bigInt2str(b1)
12379813738877118345
b2 = bigInt_pow(b1, 22) ' べき乗
?bigInt2str(b2)
1095747356866942874670236654359969216276597683447788174955606986733509423670022225232131983743614090882854465815739052044839779267369899822789127883550951615180272664979938646672434148953351414478944827697602226184465600916188807253315306160536630641445850300617883726624620685918980495021734095162050648727094223906841315178756746911585634217866565736702066664882279583761024961817642917680220796189606573581695556640625
' 階乗 1000!
fac1000 = foldl1(p_bigInt_mult, mapF(p_long2bigInt, iota(1,1000)))
fac1000s = bigInt2str(fac1000)
?fac1000s
402387260077093773543702433923003985719374864210714632543799910429938512398629020592044208486969404800479988610197196058631666872994808558901323829669944590997424504087073759918823627727188732519779505950995276120874975462497043601418278094646496291056393887437886487337119181045825783647849977012476632889835955735432513185323958463075557409114262417474349347553428646576611667797396668820291207379143853719588249808126867838374559731746136085379534524221586593201928090878297308431392844403281231558611036976801357304216168747609675871348312025478589320767169132448426236131412508780208000261683151027341827977704784635868170164365024153691398281264810213092761244896359928705114964975419909342221566832572080821333186116811553615836546984046708975602900950537616475847728421889679646244945160765353408198901385442487984959953319101723355556602139450399736280750137837615307127761926849034352625200015888535147331611702103968175921510907788019393178114194545257223865541461062892187960223838971476088506276862967146674697562911234082439208160153780889893964518263243671616762179168909779911903754031274622289988005195444414282012187361745992642956581746628302955570299024324153181617210465832036786906117260158783520751516284225540265170483304226143974286933061690897968482590125458327168226458066526769958652682272807075781391858178889652208164348344825993266043367660176999612831860788386150279465955131156552036093988180612138558600301435694527224206344631797460594682573103790084024432438465657245014402821885252470935190620929023136493273497565513958720559654228749774011413346962715422845862377387538230483865688976461927383814900140767310446640259899490222221765904339901886018566526485061799702356193897017860040811889729918311021171229845901641921068884387121855646124960798722908519296819372388642614839657382291123125024186649353143970137428531926649875337218940694281434118520158014123344828015051399694290153483077644569099073152433278288269864602789864321139083506217095002597389863554277196742822248757586765752344220207573630569498825087968928162753848863396909959826280956121450994871701244516461260379029309120889086942028510640182154399457156805941872748998094254742173582401063677404595741785160829230135358081840096996372524230560855903700624271243416909004153690105933983835777939410970027753472000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
' 個々の数字をひとつずつ切り離して配列化
fac1000m = mapF(p_mid(fac1000s), zip(iota(1, Len(fac1000s)), repeat(1, Len(fac1000s))))
digit = Array("0", "1", "2", "3", "4", "5", "6", "7", "8", "9")
' 0,1,2,3,4,5,6,7,8,9それぞれの数字の出現する回数を調べる
printM mapF(p_count_if(lambdaExpr(p_equal,1,ph_1), fac1000m), digit)
472 239 248 216 229 213 231 217 257 246
' "0"が断然多くてそれ以外は出現回数にそれほど差はなかった
（bigInt2str 関数による表示はとてつもなく遅かったが改善させた。）
関数一覧
| 関数 | 内容 | |
|---|---|---|
| long2bigInt | LongからbigIntを生成 | |
| bigInt2double | bigIntからDoubleを生成 | |
| double2bigInt | DoubleからbigIntを生成 | |
| bigInt_log | bigIntの対数 | |
| bigInt_sgn | bigIntの符号 | |
| bigInt_base | bigIntの基数(2 ^15) | |
| bigInt_abs | bigIntの絶対値 | |
| bigInt_plus | bigIntの加算 | |
| bigInt_minus | bigIntの減算 | |
| bigInt_mult | bigIntの乗算 | |
| bigInt_divide_mod | bigIntの除算 (商とMod) | |
| bigInt_divide | bigIntの除算 | |
| bigInt_mod | bigIntのMod | |
| bigInt_pow | bigIntのベキ乗 | |
| bigInt2str | bigIntからStringへの変換（10進表示） | |
| str2bigInt | StringからbigIntへの変換 | |
| bigInt_equal | bigIntの比較 (a = b) | |
| bigInt_not_equal | bigIntの比較 (a <> b) | |
| bigInt_less | bigIntの比較 (a < b) | |
| bigInt_less_equal | bigIntの比較 (a <= b) | |
| bigInt_greater | bigIntの比較 (a > b) | |
| bigInt_greater_equal | bigIntの比較 (a >= b) | |
| bigInt_gcd | 最大公約数 |
Double値からbigIntを生成する関数 double2bigInt は正確な値を求めるものではない。この逆の操作である bigInt2double はオーバーフローに対するガードはしていないので、Double の最大値である ±e+308 を超えている場合の挙動はわからない。有理数の分母と分子をこれで表現する対応はまだやっていない。
多倍長整数 x を配列 m で表現する内部構造は以下の通り。
| 配列要素 | m(0) | m(1) | m(2) | ・・・ | m(N+1) |
|---|---|---|---|---|---|
| 内容 | ±基数B | A(0) | A(1) | ・・・ | A(N) |
このとき、
x = ±A(0) + A(1)・B^1 + A(2)・B^2 + ・・・ + A(N)・B^N
である。基数Bは 2 ^ 15 = 32768 固定だが、10進表示の関数 bigInt2str を実装する際に一時的に B = 10000 のものを作るときがある。
「±基数B」 とあるのは、基数の値と数値そのものの符号を兼ねていて、たとえば 32770 と -32770 の表現はこうなる。
32770 → Array(32768 , 2, 1) ・・・ 2 + 1 * 32768 ^ 1
-32770 → Array(-32768 , 2, 1) ・・・ -(2 + 1 * 32768 ^ 1)
今回もHaskellとは全然関係ない話題だし、bigIntの実装にもほとんど使っていない。
VBAソースコードは misc_ratio.bas にある。
VBAHaskellの紹介 その20（有理数の計算）
VBAHaskellの紹介 その1（最初はmapF）
ソース https://github.com/mYmd/VBA