昨日，就寝直前に発生していた障害が再現できない．InvertFuncで配列のインデックスが範囲を超えているというエラーだ．リリースモードでは起きないバグなのだろうか？⇒Nを極端に大きくした場合に発生しているようだ．N=2024202420242024で発生した．10桁でも起きる．９桁では起きない．arraysizeに{268435455}という数が入っている．この値はMAXARRAYSIZE=&HFFFFFFFという定数だ．⇒ketaに負値が入っている．「上限オーバー」というエラーが起きている．
PowerResidueFuncの冒頭で除数が10^9を超えるとエラー復帰するようになっている．⇒PowerResidueFuncの戻り値（-1）をRS()配列に格納しようとして例外が発生する．⇒エラー復帰した場合にはRS配列に代入しないようにした．また，処理をパスする上限をMAXARRAYSIZEに変更した．I上限オーバーして計算値が得られないときは，表示欄に*****を出すようにしているが，RP（レピュニット数）の欄には出ていない．TextBox19を使っている．⇒対処した．(φ(n), @)にも*****が出ているが… 上限オーバーした場合にはPowerResidueFuncはまったく動作していないので，＠値つまり，循環節の周期も計算されていない．
TotientFuncにも上限が設定されている．ここでは上限は10^18になっている．この数字はどこから出てきたのか？Nが19桁になるとこの上限に抵触するようになる．⇒試みにこの制限を撤廃してみたが，特に問題は起きなかった．おそらく，時間が掛かる処理のときにブレークを掛けたらたまたまこの処理に当たったのではないだろうか？
N=9223372036854775807までは問題なく処理できるようになった．約１９桁まで処理できるということになる．これより大きい数字を入力できるようにするためには，数値ボックスを止めてテキストボックスに変える必要がある．数値ボックスについているUpDownボタンは結構便利なので捨て難いものがある．しばらくは，このままでよいのではないか？BigIntegerは理論的には上限が存在しないので，テキストボックスならもっとずっと大きい数でも入力できるようになるのだが…
φの約数を求める計算にも上限が入っているようだ．TextBox2を使っている．ValueChangedで10^10という制限を設けている．⇒外しても上限の9223372036854775807までそこそこに動作する．この制限は撤廃でよいと思う．⇒いや，動作はしているが，表示が空だ．どこかで処理を打ち切っているのではないか？⇒確かに，対象数が大きくなると相当な時間が掛かるようになる．{1537228672809129301}ではほとんどハング状態だ．⇒Nが１８桁を超えるとφの約数列が表示されなくなる．テキストが長過ぎるのだろうか？MaxLegthは32767だ．⇒確かに上限を超えている．45338ある．MAXTEXTLENGTHという定数に32767が入っている．ここまでで打ち切るのが無難かもしれない．
確かにTotientFuncの時間消費は大だ．いや，TotientFuncというより，FactoringSubではないか？FactoringSubはTotientFuncを呼び出している．やはり，TotientFuncで時間を喰っている模様だ．ここではTotientFuncを素数判定のために使っているのだが，もう少し効率的なアルゴリズムがあるのではないだろうか？以前はエラトステネスの篩を使っていたのだが，廃止してしまった．大きい数を連続して扱うときにはこれがあった方が有利なのではないか？
対象数Nの平方根まで持てばよいのだから，MAXARRAYSIZE = 268435455なら，72,057,593,501,057,025までが対象数になる．つまり，３２桁までカバーできる．これは悪くないような気がする．エラトステネスの篩を一度に作るのではなく，必要に応じて拡張してゆけばよいのではないか？AKS素数判定法というのがある．
https://qiita.com/srtk86/items/aaaf7225d9861512a63d
- n=a^b(a,b∈N,b>1)という形で表せるなら、nは合成数
- o_r(n)>log^2(n)となるような最小の数 r を求める
- a≤r なる a のいずれかが 1<(a,n)<nを満たす場合、nは合成数 （ここでの (a,n) は、aと n の最大公約数）
- n≤r ならば、n は素数
- 1≤a≤[sqrt(ϕ(r))logn] の範囲のaに対して(X+a)^n ≢ X^n+a(mod X^r−1, n)ならば、nは合成数
- nは素数
ここでo_r(n)は位数と呼ばれるもので，n^e≡1 mod r を満たす最小の自然数eを表す．（つまり，我々の用語でψ数）このアルゴリズムに関する論文があった．Manindra Agrawal，Neeraj Kayal，Nitin Saxenaという三人のインド人が書いたもので，2002年頃には着手していたようだが，2005年に大幅に改訂している．これはわたしが，An Efficient Factoring Algorithm by Repunit Number Methodを書いていた時期に重なる．彼らがわたしの書いていたものから着想を得た可能性はかなり高いと思う．Theory Edgeのメンバーであったかどうかまでは分からないが，入っていたとしても不思議ではない．
https://www.cse.iitk.ac.in/users/manindra/algebra/primality_v6.pdf
素数判定の高速アルゴリズムにはこの他にも「ミラービン素数判定法」などというのがあるが，こちらは確率的な判定なので100%と言えないところが難点だ．AKS法は論文のタイトルが「PRIMES is in P」となっているように，多項式時間アルゴリズムであることが際立った特徴で，O(log^10.5(n)）であるとされている．まだ，この論文は読んでいないが，わたしが構想していたことが実現されているような気がする．というか，ほとんど完全に被っていると思う．
上記のURLにはPythonで書いたコードが掲載されているが，「10桁いかないぐらいの素数判定にも数分単位の時間を取る」とあるので，実用的にはほとんど意味がないような気がする．我々のアルゴリズムなら，20桁ちかくの数を2,3分で因数分解し，約数列を出力するところまでできる．おそらく，このアルゴリズムのネックはトーシェント関数の計算ではないかと思う．それを計算するのにトーシェント関数の公式を生で使っている… いずれにしても，このアルゴリズム（AKS）はじっくり調べなくてはならないので，ともかく実装してみることにしよう．
素数判定を高速化するテクニックの一つとして，Wheel factorizationというのがある．{2, 3, 5}を既知の素数であるとして，i=7から始めて，{4, 2, 4, 4, 6, 2, 6}という周期で数値を加算しながら（つまり飛ばしながら）判定するという方法だ．FactoringSubではtotient関数で素数判定しているが，これを使えば，その手間をかなり削減できるだろう．ところで，この周期は何を意味しているのだろう？4+2+4+4+6+2+6=28なので，２８がその周期なのだが，つまり，素数はこの周期の上にしか現れないということになっているのだが… どういうことになるのか？⇒※
※上の記述には誤りがある．{4, 2, 4, 4, 6, 2, 6}ではなくて，{4, 2, 4, 2, 4, 6, 2, 6}だ．従って，周期も28ではなく，30になる．この件に関しては上記URLの筆者に「編集の提案」という形で問い合わせを送った．この３０という数はかなり興味深い．FBの数学・物理etcグループに投稿したAMSのDaily Epsilon 課題でまだ正解（の証明）を得ていない問題（下記URL↓）があり，その中に出てくる数字が３０であるからだ．「どなたか３０が最大という証明を考えてください」というコメントを付けてあるが，今のところ応答はない．
https://www.facebook.com/groups/2354748741306929/permalink/6215621771886254/
１８桁の数123456789012345680に１をプラスした数の描画が完了するまでの時間を計測してみよう．５分経過したが，まだ停止しない… もう少し小さい数で測ってみよう．12345678901234567+1→あっという間に終わってしまう．さらに+1すると今度はまた，やけくそ掛かるようになった．Wheel法というので得られる利得は7/28=1/4なのでそれほど大きな差は出ないかもしれないが，ともかくやってみよう．
顕著に速くなった．間違っていなければよいが… どうやって検算すればよいだろう．⇒とりあえず，従来論理を別関数として並行実行して比較してみよう．⇒問題なさそうだ．⇒いや，まるきりテストになっていない．前の論理がそっくり残っていた．速くなったと思ったのは錯覚かもしれない．N=12345678901234569ではこれまでと同じくらい掛かる．つまり，いつまで経っても停止しない．完全に元の論理に戻してみると，n=12345678901234576まではそこそこに動作している．比較のために設けたFactoringSub2は呼び出されていなかった．
どうも，このWheel factorizationというのはチョンボ臭い．確かに小さいNに対しては問題なく動作しているが，Nが大きくなると間違いが出てくる．⇒Wheel factorizationというメソッドが間違っているのではなく，上記のURLの記事に誤りがある．紹介されているアルゴリズム自体は間違っていないのかもしれないが，「{2,3,5}が既知の素数であるとし、i=7から探索を開始する。このとき、{4,2,4,4,6,2,6}の周期で数値を加算しながら素数かどうかを確認できる。」というのは誤りだ．正しくは{4,2,4,2,4,6,2,6}としなくてはならない．
この点を修正しても，まだ従来論理と完全一致しない．新版では大きい因数が先になるという現象が起きている．これはなぜだろう？⇒周期が７のままになっていた．⇒これで問題は解消した．多少は速くなることを期待してもよいのではないだろうか？⇒123456789012345677をテストしてみよう．⇒５分経過したが，まだ抜けてこない．旧版は止めて新版のロジックだけを走らせているのだが…14分経過，まだ解けない．この数は7×17636684144620811のように因数分解できる．7で割ることに失敗しているとは考え辛いので，17636684144620811が素数であるという判定に手間取っているのだろう．p-1=17636684144620810は，2×5×1763668414462081のように因数分解される．ということは，1763668414462081の素数判定に手間取っているのだろうか？
TotientFuncの中でもFactoringSubと同様に素数で割り込んでゆく処理が入っているので，上記のWheel法が適用できそうだ．現行ではそのような作りにはなっていないが，TotientFuncからTotientFuncを再帰的に実行したらどういうことになるだろう？却って遅くなってしまうのだろうか？いや，もちろん素数であるか否かを判定するよりも (x Mod d ≡ 0)を見る方が早いのではないだろうか？
もし，そうであるとすれば，FactoringSubの中でTotientFuncを実行するのは意味がないということにならないだろうか？素数であることを確認する以前に割り込んでしまった方が簡単なのではないか？⇒いや，TotientFuncを実行した方が速いようだ．なぜだろう？⇒どちらとも言えないというより，TotientFuncを使う方が遅いかもしれない．123456789012345671で計測してみよう．⇒この問題は結構難しいのかもしれない．123456789012345670では起動から表示までに11秒掛かる．TotientFuncを使っても使わなくてもほとんど変わらない．
素数判定の代わりにN^(p-1)≡1 mod pを見ればよいのではないだろうか？⇒FactoringSubではTotientFuncを使うのを止めた．小さい方から割って行っているので，素数であるか否かの判定は不要である．割れればそこで割り込むだけだ．TotientFuncがp-1であっても素数でない疑似素数というのが存在するので，TotientFuncを素数判定に使うのは適切ではない．いや，それも間違っている．フェルマーテストでは疑似素数を素数と認定することがあり得るが，素数以外でφ＝p-1となることはあり得ない．従って，totient関数は素数判定に使って問題ない．
totient関数をWheelを使って書き換えてみよう．⇒実装した．まだ，動作確認していないが，Nの素因数列が不定になっている．元の論理に戻しても同じだ．この値はTextBox8に入っている．⇒ここにも上限がかかっていた．10^10以上では計算をパスしている．とりあえず，制限を撤廃した．totient関数の改訂版と旧版を比較して動作確認しておこう．
▲MatrixTest中stripeが全く変化していない．少なくともKが変動すれば変化があってもよいのではないか？⇒この値はTextBox18だ．ResidueFuncとDispResidueCycleで更新している．DispResidueCycleはテスト中に実行されていて更新が掛かっている．⇒Kが変化するタイミングで更新されている．問題ないと思う．周りが活発に変化しているのに，ほとんど変化がないので止まっているように見えるのだろう．
totient関数の改訂版は一応テストをパスしたので，差し替えてよいだろう．123456789012345671も時間はかなり掛かるが出せるようになった．この数は素数である．一度バックアップを取っておこう．⇒FactoringSubでまだTotientFuncが残っていた．
▲123456789012345671はかなり時間が掛かるが，素数のインジケータが点灯してからなおしばらく掛かるのはなぜだろう？素数であることが判明した時点でほとんどの計算は完了しているような気がするのだが… 逆数の計算なら時間は掛かってもおかしくないが，上限を超えているためパスしているはず… PrimeNというチェックボックスだ．DispParametorという関数で表示している．GetDivisorsがその後で実行されている．ここがネックになっている．
pが素数であれば，pと互いに素であるような任意の整数aにおいて，a^(p-1)≡1 mod pが成立する．従って，pと互素であるようなaを選択してa^(p-1)≡1であれば合成数であると結論することができる．これはかなり効率に寄与するのではないかと思われるので，IsPrimeで採用することにしよう．2x3x5=30＝a，ないし，2x3x5x7=210＝aで検査することが考えられる．これはFactoringSubでも使えるのではないか？現行ではWheel Factoringを使って，3～7までの検査とそれ以降を分離しているので，前半が終わったところにこの検査を入れてやればかなり高速化する可能性はある．いや，Factroingでは因数分解しなくてはならないのだから，意味がないかもしれない．