ゼルコバの木テント村

どうも，やはり精度 2^30 以上ではメモリ不足が起きるのは避けられないようだ．

最初にFpr,atter というのをインストールすることを求められたのでインストールした．新しい版では，コードを実行する前に，先行してプレコンパイルというのが実行されるようだ．

Precompiling VSCodeDebugger finished.
   1 dependency successfully precompiled in 35 seconds
   Activating new project at `C:\Users\babalabo\.julia\environments\v1.12`

2^29までは順調に進んだが，2^30でクラッシュしてしまった．状態は以前と変わらない．計算速度は２割くらい速くなったような気はする．

Precision set to 2^29=536870912 bits, which is approximately 161614248 decimal digits.
Θ=4.881152430408162405204287101960529897794735314100e-161614249,　N = 536870912
sinΘ=1.533459261642224750278082996492416767270795121478e-161614248
pai=3.141592653589793238462643383279502884197169399375e+00, time = 509 seconds
n = 1, tanp4=1.000000000000e+00, diff=4.881152430408e-161614249, time = 2678 seconds
◎　Convergence achieved at n = 1 precis = 2^29 = 536870912 keta = 161614248
mul.c:832: MPFR assertion failed: (mpfr_uprec_t) bq + cq <= ((mpfr_prec_t) ((((mpfr_uprec_t) -1) >> 1) – 256))

Julia のπ定数は以前として間違っている．ただし，REPL では正しい値が出力されている．

pi＝3.141592653589793115997963468544185161590576171875e+00

julia> BigFloat(pi)
3.14159265358979323846264338327950288419716939937510582097494459230781640628619

フェイスブックのページを閲覧中，広告に載っていた「続きを読む」ボタンに引っかかって危うくウィルスに感染するところだった．

カスペルスキーがブロックしてくれたので，何とか難を逃れることはできたが．．．危ないところだった．．．

どうも，Julia の動作がおかしい．標準整備されているπ定数の値がかなり狂っている．

NET  ：3.141592653589 793238462643 383279502884 197169399375 105820974944
     592307816406 286208998628 034825342117

電卓：3.141592653589 793238462643 3832795

Julia ：3.141592653589 793115997963 468544185161 590576171875

自家：3.141592653589 793238462643 383279502884 197169399375

Julia のπは１６桁目からすでに間違っている．この値は精度を設定しても変わらない．だとすると，アルゴリズムの中で自家製πと比較したときに差異が出なくてはならないのだが．．．どうなっているのだろう？

これは推測だが，πという数値はシンボルとしても機能しているのではないか？つまり，πに近い数値は丸めてπとみなすという動作が入っているような気がする．その証拠に，それぞれを２で割って差分を取ったところ，明らかに差異が現れるようになった．⇒これで検証できるようになった．というより，Julia のπはおかしな値を持っているので，検証には使えないことが分かった！！！REPLなら正しい値が取れる．

julia> BigFloat(pi)
3.141592653589793238462643383279502884197169399375105820974944592307816406286198

前は正常動作していたような気がするので，再インストールしてみよう．⇒バージョンが上がって，v1.10 から v1.12 になった．事態の好転を期待できるだろうか？いくらかでも改善されていればうれしい，．．

Precision set to 536870912 bits, which is approximately 161614248 decimal digits.
base=-1.000000000000000000000000e+00 + 0.000000000000000000000000e+00 i
Θ=4.881152430408162405204287101960529897794735314100e-161614249,　N = 536870912  
z=1.000000000000000000000000e+00 + 1.533459261642224750278083e-161614248 i       
sinΘ=1.533459261642224750278082996492416767270795121478e-161614248
pai=3.14159265358979323846264338327950e+00, time = 699 seconds

base=-1.000000000000000000000000e+00 + 0.000000000000000000000000e+00 i
Θ=4.881152430408162405204287101960529897794735314100e-161614249,　N = 536870912  
z=1.000000000000000000000000e+00 + 1.533459261642224750278083e-161614248 i       
sinΘ=1.533459261642224750278082996492416767270795121478e-161614248
pai=3.14159265358979323846264338327950e+00, time = 699 seconds

Θ=4.881152430408162405204287101960529897794735314100e-161614249,　N = 536870912  
z=1.000000000000000000000000e+00 + 1.533459261642224750278083e-161614248 i       
sinΘ=1.533459261642224750278082996492416767270795121478e-161614248
pai=3.14159265358979323846264338327950e+00, time = 699 seconds
z=1.000000000000000000000000e+00 + 1.533459261642224750278083e-161614248 i       
sinΘ=1.533459261642224750278082996492416767270795121478e-161614248
pai=3.14159265358979323846264338327950e+00, time = 699 seconds

sinΘ=1.533459261642224750278082996492416767270795121478e-161614248
pai=3.14159265358979323846264338327950e+00, time = 699 seconds

pai=3.14159265358979323846264338327950e+00, time = 699 seconds

n = 1, tanp4=1.000000000000e+00, diff=4.881152430408e-161614249, time = 3427 seconds
Θ=6.986524479602259580995891220200500532802060184779e-80807125,　N = 268435456 

いや，どうも動作がおかしい．もう一度やり直してみよう．ターミナル出力に前回のログが残っているので，分かりづらくなっている．最後の方を見ると，pai の計算に699秒～１２分，tan(π/4)の計算に3427秒～５７分掛かっているが，計算は完了しているように見える．diff=4.88e^161614249 に対し，θ=6.99e80807125 だから，diff < θ が成立している．その後のプリント文でメモリ不足が起きている可能性はある．precision(BigFloat)を出力しようとしている．どうもこれがまずかったのではないか？

修正して再実行してみたが，やはりメモリ不足が起きている．最後にN = 268435456までは出力されているので，どうもブレークできていないようだ．ブレークの条件は，diff < θ で，いずれも BigFloat だ．

  Activating project at `C:\Users\babalabo\.julia\environments\v1.10`
Precision set to 536870912 bits, which is approximately 161614248 decimal digits.
base=-1.000000000000000000000000e+00 + 0.000000000000000000000000e+00 i
Θ=4.881152430408162405204287101960529897794735314100e-161614249,　N = 536870912
z=1.000000000000000000000000e+00 + 1.533459261642224750278083e-161614248 i
sinΘ=1.533459261642224750278082996492416767270795121478e-161614248
pai=3.14159265358979323846264338327950e+00, time = 664 seconds
n = 1, tanp4=1.000000000000e+00, diff=4.881152430408e-161614249, time = 3266 seconds
Θ=6.986524479602259580995891220200500532802060184779e-80807125,　N = 268435456
ERROR: LoadError: OutOfMemoryError()

diff = 4.881152430408e-161614249, θ=4.88115243040816240e-161614249でほぼ同じだ．等号が入っていないためだろうか？やり直してみよう．結果が出るまでに一時間は掛かる．．．

！解けた！3614＋677＝4291秒＝７２分．pai を計算するまでなら，１１分だ．桁数は161614248桁，１兆６千億桁だ．多分もう一桁上がってもそれほど大きくはならないだろう．

  Activating project at `C:\Users\babalabo\.julia\environments\v1.10`
Precision set to 536870912 bits, which is approximately 161614248 decimal digits.
base=-1.000000000000000000000000e+00 + 0.000000000000000000000000e+00 i
Θ=4.881152430408162405204287101960529897794735314100e-161614249,　N = 536870912
z=1.000000000000000000000000e+00 + 1.533459261642224750278083e-161614248 i
sinΘ=1.533459261642224750278082996492416767270795121478e-161614248
pai=3.14159265358979323846264338327950e+00, time = 677 seconds
n = 1, tanp4=1.000000000000e+00, diff=4.881152430408e-161614249, time = 3614 seconds
◎　Convergence achieved at n = 1 precis = 536870912 keta = 161614248
  *  Terminal will be reused by tasks, press any key to close it.

精度 2^30を仕掛けてみたら，エラーになってしまった．

keta::BigInt = BigInt(floor(log10(2)*precis)) でエラーになっている．2^30=1073741824だ．プログラムを整理して，任意精度を１～３０まで連続テストするように書き直してテストしたところ，やはり３０で落ちている．2^29 までは完了しているのだが，2^30 ３０冒頭で落ちているようだ．keta::BigInt = BigInt(floor(log10(2)*precis))  が通っていない．assertion failed だが，メモリ不足は起きていないようだ．限界と言ってしまえばそれまでだが．．．

円周率πの近似値１億桁の壁が突破できない．８千万桁まではクリアしているが，一億桁を超えられない．比較検証用にJulia 備え付けπの基準値を用いれば多分パスするはずだが，ある種の循環論法になってしまうのでやりたくない．基準π値を用いない方法として，Julia の標準三角関数にπを代入して確認する方法を試している．sinπ=0，cosπ=1，tan(π/4)=1となるから，これらと一致すれば，得られたπ値は（指定精度の近似値として）正しいと言える．tanで試してみた結果は，

Precision set to 536870912 bits, which is approximately 161614248 decimal digits.
Θ = 2.382564904887951073216169781732674520415196125559e-323228497,N=1073741824
z = 1.000000000000000000000000e+00 + 7.485048401896851550940689e-323228497 i
sinΘ = 7.485048401896851550940688646230092416056892404614e-323228497
pai = 3.1415926535897932384626433832795028841971693993751058209749445923e+00
n=1, Θ=2.382564904887951073216170e-323228497 tanp4=1.000000000000000000000000e+00

tanp4=tan(pai/4) とBigFloat(1)を比較しているのだが，一致しない．多分，tan4pに微小な端数が残っているのだろう．１とtan(pai/4)の差分Δがゼロに等しいことの代わりに，Δ<θ で代用することは考えられる．このロジックではθ～０とみなされ．これまでにも，この方法は試みてはいるのだが，そこまでゆく前に極微小量の計算でクラッシュしてしまっていた．tanに切り替えたのはそれを回避するためだが，もしかすればうまくゆく可能性はあるような気がする．

ダメだ．Δ=1-tan(pai/4)として，Δ > θ になってしまう．Δ=4.88e-161614249 で．θ=2.38e^323228497．Δは設定された精度で１だが，θは擬似的な０なのでそれよりもずっと小さい．これは避けられないと思う．base～１なので，これと比較すればよいのでは？いや，BigFloat(1)と比較しているのだから，それもおかしい．どうしたらよいのだろう？打つ手がなくなってきた．．．

現行とは逆にθをより小さくする方向で動かしたらどうなるか？いや，現行ではNはｎに従って増加するようになっているのだから，θは小さくなっているのではないか？

ジアン角単位系を用いて，円周率πの近似値１億６千万桁の計算することに成功した．計算式履いたってシンプルなもので，

π ～ Im( (-1)^ε )) / ε for extreamly small ε．（ε はジアン角単位系の極めて小さい角度）

を単純に算術演算しているだけだ．Julia には BigFloat という任意長精度の浮動小数点型あるので，実装は容易だ．計算の所要時間は2180秒，36分で．これが速いのか遅いのかはよく分からないが，EXEならもっと速くなることは確かなので，やってみることにしよう．

Julia には定数πが「無理数」として備わっている．これは多分，必要に応じて，任意サイズのπが使えるという意味だろう．これを使えば，現在のコードをもう少し改善できる．比較参照用データは不要になるし，もっと高速でΘを小さくできる．現在はΘに 2^-n という小数を与えているが，これを 10^-n とすればずっと小さくなるし，計算精度も測りやすくなる．実装は簡単なのでやってみることにしう．

久しぶりに開発環境に戻ってきた．「円周率πの近似値計算」をやるつもりなのだが，どうもさっぱり勝手が見えなくなってしまった．まず，言語の選定から始めなくてはならない．条件は複素関数とBigFloatが使えることだ．C++, Julia, Python はいすれもこれらの機能を持っている．これは最終的にはツールとして使えるようにしたいので，その意味ではC++しか選択肢はないのだが，C++の複素数は64ビットまでしか使えない．他の２つは任意長DOUBLEが使えるので，とりあえず，Juliaを使うことにした．Juliaはこれまでにもかなり使い込んでいる．

バックアップにはフォルダが１０個残っている．しかし，有り物を動かそうとしてもうまく動かない．Eドライブ（開発支援）には５個フォルダがあり，その中の「蜜蜂天国」から入ったのだが，ソースファイルが見当たらない．⇒フォルダを開いたら見えるようになった．ワークスペースを指定したのだが．．．別途フォルダを選択する必要があるようだ．

実行したら，Revise というモジュールをインストールするように促されたので実行したら，次のようなワーニングが出た．

julia> using Revise
[ Info: Precompiling Revise [295af30f-e4ad-537b-8983-00126c2a3abe]
WARNING: using IR.SSAValue in module LoweredCodeUtils conflicts with an existing identifier.
WARNING: using IR.SlotNumber in module LoweredCodeUtils conflicts with an existing identifier.

どうも余分なことをしてしまったようだ．確かJuliaはソースを修正しても直ちに反映しないような動作になっていたような気がするが，それを改善してすぐ動くようにするパッケージと思われる．⇒再起動して動くようになった．．．少し感覚が戻ってきた．見覚えのあるこんな画面が出てきた．いや，すごいな．こんなことまでやってる．

このGIFアニメはいつまでも終わらない．．．121KBというそれほど大きなファイルではないのだが．．．再帰関数というフォルダには何本か数値計算的なものが入っているので，ここから開始することにしよう．何をやろうとしているのかはよく分からないが，BigFloatが使われているので，目的に近いような気がする．

F = 6.0 / (2.0 + F) という再帰関数を実行し，A=1.64575131106459059から，(A+1)^2=6.9999999999999999 という値が出力されている．pre::Int = 256 となっているので，256ビット整数演算をやっているようだ．「階乗の末尾連続ゼロ.jl は動かなかった．「y＝x^2-3/2」でもエラーが起きている．まぁ，とりあえず動いているので，なんとかなるだろう．

黒木のべき乗和公式からファウルハーバーのべき乗和公式への変換を実現するという課題をそろそろ片付けておきたい．大筋は大体見えてきたのだが，まだ細部で煮詰まっていないところがある．

その前に一つ気になる点として，黒木のべき乗和公式の塁和が j=2→k+1 になっているという点がある．文献ではここは，j=0→k-1になっているというところだ．つまり，∑ {k=0→k-1} binom(k+1, j) P_j(n) のところが，黒木の公式では ∑ {k=2→k+1} binom(k+1, j) F_{k+1-j}(n) となっている．インデックスのシフトは加算の個数を変えないので，それ自体は安全だが．．．

どうも，この式自体が最初から間違っていたような気がする．この公式を導出したところ（記事）を再点検する必要がありそうだ．⇒いや，もしかすると間違っていない可能性もある．j=2のとき，k+1-j=k-1，j=k+1のとき，k+1-j=0 となるので，P_j に関しては一致している．しかし，二項係数の部分は合わないのではないか？(k+1, 2)と(k+1, k-1)，および，(k+1, k+1)と(k-1, 0) は一致するだろうか？後者が一致することは間違いないとして，前者はどうか？

m=k+1-j とすれば，m=0のとき，j=k+1，m=k-1のとき，j=2 となるから，数字的には合っているのではないか？逆に言えば，m=k+1-j とすれば, j=2のとき，m=k-1，j=k-1のとき，m=0 となるから，通常の公式と一致する．この，「k+1-j」がそのまま出てくる式も結構見かけるので，そのままでもよいという気もしないでもない．この式は2024年の7月から公開されているので，間違っていたとすればとんだ恥さらしだが，とりあえず，数式的には合っていると考えてよいのではないか？