ゼルコバの木テント村

まだ，午前９時前（AM08:55）だというのに，こんなキチガイじみた数字が出ている．このペースでカウントが伸びれば，日付が変わるころには確実に1000人を超えることになるだろう．

何か異変が起こったという感じなのだが…

そろそろ外部に持ち出せる雰囲気になってきたので，試しにEXEをVAIOにコピーして走らせてみたのだが，ウィンドウ枠が縦長過ぎてスクリーンに収まらない．対策を色々考えてみたが，画面の一番下のブロックを切り出して独立のコントロールパネルのようなもの仕立てるという結論になった．このパネルには主に５種の検定を実行するためのボタンがあるが，一般ユーザにはあまり関心がない部分と思われるので，却って使い易くなるのではないかと思う．このパネル上のボタンで有用なのは右端に置いてある，①Get φ，②Update，③Build Matrixの３つのボタンでできればこのボタンはメイン画面に置きたいところなのだが… まぁ，手元にリモコンを置いて操作するような感じになるのではないかと思う．ともかくやってみよう．⇒大体動作するようになった．

▲InvertTestで画面出力が出ない．⇒解消した．

▲Updateボタンでマトリックスが開き，BuildMatrixは無動作になっている．⇒修正した．

▲Ξの判定が間違っている．N=12347のとき，Ξは　0111111111111111111111111111111111111111111111111
1111111111111111111111111111111111111111111111111
1111111111111111111111111111111111111111111111111…

のような数になるが，ｎUは09999999999999999999999999999999999999999999999999
999999999999999 のように先頭にゼロが付いた数になっている．Uは 809, 913, 339, 272, 697, 821, 333, 117, 356, 442, 860, 613, 914, 311, 168, 704, 948, 570, 502, 956, 183, 688, 345, 347, 047, 865, 878, 351, 016, 441, 240, 787, 235, 765, 773, 062, … のような 10進で6173桁の数で，それを12347倍したものがｎUだから，５桁上がるとして6178桁になるはずだ．ｎUの先頭に入っている０というのは何かの間違いではないだろうか？

nUの数字列はConvertNum2Stringで生成されている．ToBCDはこんな巨大数を処理しきれないので，MaxKetaSuで打ち切っているが，バイト配列の先頭に0が入って，{0,9,9,9,…}のような内容になっている．数値をテキストに変換する場合には頭から割り込んできているはずだから，桁数が上限を超えていても処理できなくてはならない．⇒どうも怪しくなってきた．ToBCDでは対象数Nを頭からBで割り込んだ数を配列に格納し，それを逆順にした配列を返している．配列に格納されているのはN%Bの剰余Rだ．確かにこれはBで割った商を求めているのではなく，Bを進数に変換したときのバイト列を取っているのだから，正当な操作と言えるだろう．多分これしか方法はないものと思われる．

つまり，一度頭から尻尾までさばかなければ皿には盛り付けできないということだろう．ということは桁数がオーバーフローした時点でご破算にしなくてはならない．ToBCDでは上限がMAXTEXTLENGTH=32767だ．これでは到底歯が立たない．というか，そもそも整数値と桁数を混同しているように見える．ToBCDではともかく末尾桁まで計算しなくてはならないのだから，その限界を正確に見極める必要がある．ToBCDでは最初から作業配列のサイズを

Dim arraysize As BigInteger = MaxKetaSu

でMaxKetaSu = 1000で固定している．これはかなりまずい．作業用配列サイズをMAXARRAYSIZEまで拡大してみた．桁数は6178なので完全に処理できる．Ξの計算ではテキスト化したｎUではなく，ｎUの生値を使っているはずだから，この判定は正しいと言える．また，出力されたレピュニット数11111…も末尾までさばいた値から生成したものだから正しいと言える．表示桁数は1000桁になっているが，これはこれで十分だろう．ToBCDでオーバーフローした場合には，デバッグコンソールにメッセージを出してブレークするようになっているが，エラー復帰するようにした方がよい．間違った値が表示されると却って混乱の元になる．エラーコードを返すようにした．

UコードはNが６桁辺りで限界になってしまうが，桁数で１万桁台でも表示不能となっている．もう少し拡げることは可能なのではないか？MakeRecursionUnitはMaxOutput＝10000が上限になっている．上限をMAXARRAYSIZEに代えてみた．問題なく動作している．ただし．桁数が3万を超えると流石に重くなってくる．MAXTEXTLENGTH4 = 32767を上限としてみよう．⇒まぁ．そこそこの動きになった．

テキストボックスの名前の付け替えを行った．一応収束した．また，ClearParameterという関数を作ってValueChangedの冒頭で呼び出すようにした．これはすべてのテキストボックスを一旦クリアするためのものだ．これで更新されていないボックスが浮き出しになる．また，InvertFuncの冒頭でClearInverseを実行するようにした．

▲SeedTestではcycleとstripeを更新してない．また，PrimeTestではInvertFuncを呼び出しているが，パラメータはまったく更新していない．InvertTestでも同じだ．以前は動いていたような気がするのだが…

▲N=12347のとき，Updateボタンを押して例外が発生する．ClearInverseで例外が発生している．なぜだろう？notation.Text = “”というコードで例外が置きている．ただし，初回は問題なくパスしている．⇒パーツを作り直したら収まった．⇒原因は不明

PrimeTestではBの値が変化する以外なにも表示されない．動作はしているのだが，テキストが長過ぎて格納できないのだろう．いや，切り詰める処理は入っている．Uをテキストに変換したときの長さは6169で，テキストボックスの上限は32767だから十分入る大きさだ．⇒消去してから描画までの時間が長過ぎるためだ．描画した直後に消去されるような流れになっている．DispInvertFuncの冒頭でクリアするようにして見えるようになった．ただし，Bの因数リストとB進表記は空欄のままだ．⇒PrimeTestClickの中でDispNotationを実行するようにした．これでボックスBの中身はすべて出力されるようになった．

InvertTestにも同じ問題がある．⇒DispInvertFuncの中で実行した方が分かり易い．stripeはKを変えない限り一定だがNを変えればcycleは変化しなくてはならない．まずこれを見てみよう．⇒DispResidueCycleで出力されている．DispResidueCycleはResidueFuncProから呼び出される．呼び出しているのは，①ValueChangedPro，②modulo3.Click，③GoButton，④TestMatrixの4箇所だ．SeedTestではFactoringとInvertFuncは実行しているが，剰余計算はやっていない．modulo3というのは，現在ReisueTestと呼んでいるものだ．⇒名前を付け替えた．③のGoButtonはすでに廃止されている．

SeedTestのどこかでやけくそ時間が掛かっている．どこかでコントロールを失っている感じだ．カーソルは砂時計のままなので処理は完了していないのだが，ブレークで止められない．N=1234567．いや，処理は明らかにすでに終わっている．カーソルが戻っていないだけと思われる．⇒UseWaitCursorがオンになっていた．⇒ResidueTestとSeedTestではValueChangedProを実行するようにした．これでもれなく表示が更新されるようになった．いや，まだだいぶ落ちているような気がする．ResidueTestでは，1/nやUが更新されていない．ResidueTestでは逆数計算をやっていないからだ．これはそれでよいのではないか？

▲InvertTestを実行後の状態がおかしい．N=12347, B=10で因子が1,2,3,5 表記がDLHとなっている．前の状態が残ったままになっている．PrimeTestも同じだ．⇒Passflagが立った状態で値を戻しているためだ．passflagをリセット後に値を戻すようにした．

これで大体収まったのではないかと思う．

▲Getψ関数はψ値を確定するためのKを探索するルーチンだが，動作に不審がある．たとえば，N=12358でK=27のときψ=3で確定しているにも関わらず，Getψでk=３にジャンプするというのはおかしい．

▲各種の検定の結果をテーブルに出力する．マトリックスの異種文字数を出す．最初の項目は将来的課題とする．二番目は対処する．

▲出力文字数は制限しない（テキストボックスの容量まで）

現行ではnU/(B-1)の値にrepというラベルを付けているが，この数値はつねにレピュニットになるとは限らないので，正確とは言えない．多少うるさいが，直接nU/(B-1)と表記しておいた方が間違いがないのではないか？合わせて，rep値がレピュニッであるか否かを示すインジケータも表示しておきたい．昨日からEFAの読み直しを始めているが，レピュニット値を得られるか否かはｎの素因数分解に深く関係している．

ある数ｒがレピュニットであることをどのように判定すればよいだろう？レピュニットはBに依存するので計算にはBが現れなくてはならない．r=(B^e-1)/(B-1)だから，r*(B-1)+1してB^eに戻せるかどうかを試してみればよいのではないか？いまの場合，eは未知だから，X=r*(B-1)+1として，B^log(X)=Xとなることを確認すればよい．いや，これは間違いだ．B^(logX/logB)=Xだ．Bを底とする対数を直接求めることは可能だろうか？x=logX / lobBを計算してやればよい．⇒できた！かなりわかり易くなった．⇒ダメだ．誤動作している．

n=29, B=5のとき，r=11111111111111がレピュニットでないと判定されている．doubleの計算なので誤差が出てしまっているのだろうか？x=13.999999999999998という値だ．CTypeで変換しているが切り捨てになっているのではないだろうか？⇒そのようだ．Math.Round関数を使うようにしてみた．Roundの引数で四捨五入を行う小数位を指定できるが，15とした場合には落ちてしまう．15桁目は8で，14桁目までは9が入っている．引数を14として動作したので，これでよいことにする．整数値で検算しているので問題ないはずだ．

ダメだ．n=73, B=29でr=11111111111111111111111111111111111111111111111111

1111111111111111111111 の数字をレピュニットでないと判定している．71.999999999999986を71にまるめている．⇒第二引数を12まで落としてみた．これで(B^e-1)/(B-1)の値を完全に検査できる．この値には何か名前が付いていただろうか？この値は結構重要度が高いので，何かギリシャ文字を与えてもよいのではないだろうか？⇒φ，ψときて，次はχではないか？しかし，文字面ではｘと同じになってしまう．upsilonというのもある．nUの約数なのでΥもおもしろいと思うのだが，見た目がYになってしまう．Ξは横棒だが，一が３つ並んだ形状に見えるので，これを使ってみよう．つまり，上で「nU/(B-1)」と表記した値だ．Ξ=nU/(B-1)ということになる．⇒直感的で悪くないと思う．

このインジケータをチェックすると，ｎとBが互いに素であることと，Ξ値がレピュニット数であることは完全に同値であることが確認できる．これは大きな収穫だ．ダニエル・デジオジーとの対話の中で数字和（digit sum）の話が出てくる．どこか場所があれば出してみたいのだが… 数字根（digital root）というのもある．これは数字和を一つの数とみなしてその数字和を求め，のような操作を繰り返して一桁になったときの数字だ．ｎのB進表記を出すところでダンプ出力として出しておくことにしよう．DispInvertFunc→DispNotationで表示している．

▲Bを変えたとき，Bの因数とB進表記が変化しない．⇒対処した．

▲Bを変えたとき，何も外部出力されない．⇒暫定修正としてBottomLineを呼び出すようにした．

▲DispInvertFunc で循環節を重複して生成しているようだ．

DispInvertFunc→MakeRecurringDecimal→※→DispParametor2
DispInvertFunc | Inverse.Click→DispInvert→※

文字列の生成はMakeRecurringDecimalに統一し，DispInvertの中から呼び出せばよいのでは？いや，DispInvertFunc → DispInvertの呼び出しはDUMPNUMSTRINGの限られているから止めてよい．従って，

DispInvertFunc→MakeRecurringDecimal→※→DispParametor2
Inverse.Click→DispInvert→※

のようになる．つまり，DispInvertの論理をMakeRecurringDecimalの呼び出しに変えればよいのではないか？というか，DispInvertを廃止して，Inverse.Click→DispInvert→DispInvertFuncでもよいような気がするが… いや，そもそも，Inverse.Click→DispInvertFunc は実行されている．Inverse.ClickではDispInvertFuncの後に，DispInvertとDispRecurringDecimalを実行している．DispRecurringDecimalは出力画面に「循環小数」などの情報を出力している．これはDispInvertFunc に取り込んでしまった方が早い．⇒DispInvertは廃止した．

▲nRの計算が間違っているような気がする．K=14のとき，nR=14という値が表示されている．しかし，nRはn%KなのだからnR<14でなくてはならない．ｎRはKの剰余だが，少しだけ違うところがある．ｎR=0のときには，ｎR=Kとなるように調整してあるところが味噌だ．これによって，マトリックスの行番号とｎRが一致するようになっている．eRの場合も同様の調整が入っている．2023-06-07のログによれば，

「たとえば，K=4のとき，e=3, 6, 9ではeRは０になるが，マトリックスの指標としてみる場合には０という数字は現れないので，eR=0 → eR=K-1のように換算する必要がある．ｎR=0にも同様の問題はあるが，行の場合は行ヘッダが表示されないので，剰余が１列目に表示されているが，それが０になっているので，特に手当しなくても問題ないようには見える．⇒縦軸は番号Nとみなすべきだから，やはりゼロは適当ではないと思う．eRの場合と同様0→Kのように読み替えるべきだ．」

のようになっている．驚くべきことは，ｎRと数字根は完全に一致する．ただし，すべてのKについてそれが成り立っている訳ではない．

▲上の修正の影響と思われるが，Kを変更してもｎRが変化しなくなってしまった．いや，処理は抜けていないようだが，描画されないような状態になっている．何かタイミング的な問題があるのだろうか？Nが小さい範囲では問題は起きない．⇒ChangeCursorが入っていないため，処理が続いているのが見えていないことが原因と思われる．⇒対処した．

「ｎRと数字根は完全に一致」しているように見えたのだが，まったく起きなくなってしまった．少なくともK=45では一致していたと思うのだが…対象の数字を一桁小さくしたら一致するようになった．確かに数字根の計算には逆数演算が必要というか，そのカテゴリに含まれるが，逆数方面は上限がタイトなので動作しない場合があるのではないか？

▲べき乗の計算値やB進表記などをBottomLineで表示しているが，あまりよくないと思う．BottomLineは処理中に繰り返し呼び出されるので無駄な処理が入り過ぎる．DispParametor2に移してみよう．ダンプに紛れてしまうので読みづらい．ValueChangedの末尾に置いてみよう．

数字根がｎRと一致する条件がわかった．K=B-1だ．つまり，B=10のときK=9の剰余はBの数字根と一致する．B=10，K=45では一部一致する場合もあるが，不一致となる場合もある．また，次のようなことも言える．「KがB-1の倍数のときは，B^x％Kは一定値を取る．」いや，一定にならない場合もある．たとえば，B=10でK=27では{1, 10, 19}の３つの値に分離する．しかし，K=36なら一定だ．おそらくBの因数をp1*p2*…として，(B-1)*p1*p2…などであれば一定になるのではないだろうか？Kの素因数分解も必要かも知れない．B(B-1)というのも気になる数字だ．循環節の書式を下記のように変更した．

A: 0.005|4347826086956521739130

とても読み易くなった．

▲Kを変えてもe%K-1が変化しない．stripeは変化している．⇒e<Kならば値は変化しないと思う．e>Kになれば変化するようになる．剰余演算ではつねにこのようなことが起こる．つまり，Kより小さい値は剰余演算ではほとんど意味がない．⇒というか，Kの範囲内でしか通用しない．

▲Bを変化させたとき，Bの因数と表記が変化しない．Bの素因数を表示しているのはDispNotationだけだ．DispNotationを呼んでいるのはValueChangedProだけだ．⇒base2.ValueChangedイベントでDispNotationを呼び出すようにした．

▲Nを変更して（n,K)が更新されない．この値を更新しているのは，DispModPowだけだ．⇒いや，呼び出されているし，更新もされている．１や７などの数字で色も地味なので見落とされていただけだ．

▲横数列が尽数列になっても，原始根マークが付かない．K=11のとき，ψ＝K-1のときのみインジケータが点灯する．

Let p be a prime number which is neither 2 nor 5, and the recursion unit of p be U, then pU/9 is a repunit number Rn, i.e., pU=10n-1.

という命題を提示し，その証明を与えているが，その中でｐが素数でｐとeが互いに素のとき，ｂ^e-1≡0 mod p となるような最小の整数eが存在し，これはｐがb^n-1を割り切ることを意味するとしている．これは，フェルマーの小定理を言い換えただけのものであるから，正しい．従って，pU=b^e-1となるような最小の整数Uが存在するとし，これはUがpの循環単位であることを意味するとしている．

n=11, e=2, K=7, B=5のとき，＠＝５, U=284, nU=4444，rep=1111で，Uは11の循環単位である．このとき，5^5%11=284…1となるので，B^@%n=U…1が成立している．ここで，＠は循環節の桁数である．eという指数はべき剰余数列のψに相当する数と考えられるが，Q=264に相当する数がψないし#に出てくるような計算式を立てられるだろうか？#ないしψはKより小さいと考えられるので，相当大きなKが必要になる．nもかなり大きな数でなくてはならないだろう．K = 1111に設定するとかなり近い数字に近づく．Kが素数でないと大きな数が出てこない．ちょっと的が外れている感じ．

N^eの計算値を再利用するために出力用のテキストボックスを追加しようと思ったのだが，すでに十二分に画面が混み合っているので出力ペーンに書き出すことにする．EFAと実装が合わない点として，09 Mar 2005のメールには For example 21 is a composite but has a repunit number 111 with recursion unit 03 とあるが，これは道具箱の動作とあっていない．1/21=0.04761904761904761904761904…で，U=047619，nU=999999, rep=111111だ．これは明らかに誤認と思われる．このメールに記載されたリストでは，U=047619となっているので，何か勘違いしていたのではないだろうか？

EFAではn→r{m}(ｍ桁のレピュニット数)のとき，nをr{m}のオリジネータと呼んでいる．このとき，ｎの素因数分解とr{m}には明らかに何らかの相関がある．ｎは必ずしもレピュニット数のオリジネータになる訳ではないので，ｎから生成されるレピュニット類似の数をrep(n)と表記することにしよう．rep(n)がレピュニット数r(m)になる条件を調べることは意味があるように思われる．

▲Bの素因数分解も出す必要がある．レピュニットに関わりがある．

09 Mar 2005のメールには誤りがあるが，その下の方のリストは正しいので，注意深い読者なら多分気付いたことだろう．このメールの終わりの方に，ｐが素数であるための５条件というのを列挙している．この項目はチェックする必要がある．もし，これらが正しいとすると素数判定をかなり/ある程度高速化できる可能性がある．

(1) p has a repunit R_k=pU(p)/9
(2) k divides p-1
(3) for p > 3, (p,9)=1
(4) p divides R_k
(5) p and k are mutually prime

朝起きると，PCが落ちているというのは日課になっているが，今のところ何とか走ってくれている．これが起動できなくなったらどういうことになるのか？考えても仕方ないので考えないようにしている… VSを起動してkinaiを走らせたら，以下のパネルが出た．

一旦ソリューションを閉じて再起動で問題なくロードできた．昨日は新しいプロジェクトを起こして「AKS素数判定法」を実装しようとしていたのだが，Pythonコードからの移植なので動かすまでにはかなり掛かりそうだ．むしろ，Pythonコードを直接試してみた方が早いのではということになったので，下記URLからコードを借りることにしたのだが，この中にψ数を計算しているところがある．現行ではψ数が求められないというケースがしばしばあるが，ψ数を独自に計算するツールを一つ作っておいた方がよい．ψ(n, K)で適切なKを選択すれば任意のｎでψ数を計算できるはずなので，まず，それを実装してみることにする．

これはMatrix Testでやっていることに一番近いので，MatrixTestの隣りにボタンを置くことにしよう．（ボタンはすでに昨日作っている）

▲PsiFunctionで反例が出た．N=1234567926，K=5，e=44のとき，ψ値を検出できない事例が出た．⇒1234567926ではなく，1234567927だ．動作がおかしい．その後どこかへ跳んでしまった．画面が出てこない．直下のループに落ちていた．⇒このNはK=2でψを与えるようになっている．つまり，GCMが１であっても任意のKでψを取得できるとは限っていない．K=5の場合，剰余周期は{2,4,3,1}なので，べきを取ればR=1になる場合もあり得ると考えられる．⇒確かにeR=4ならR=1となる．ただし，ψ数はeに依存しない値なので，ψ不定という状態は変わらない．ちなみに，1234567927は素数である．

▲MatrixTestではφの約数列とψは更新されていない．⇒しかし，φの約数をどこかで（部分的に）書き込んでいる．⇒Kを更新した時点でイベントハンドラがValueChangedProを呼び出している．⇒K_TextChangedはpassflagでパスするようにした．出してもよいのだが，遅くなるだけだし，ここでは不要と思われる．

▲ψが決定されることとψに原始根マークが出るのは完全に一致しているように見える．それで正しいか？⇒ψ=K-1の場合を原始根と判定している．また，このときφ(K)=ψとなる．ψはつねにKよりも小さいのか？ψ＝#と考えられる．#<Kであるのだから当然ψ<Kとなるだろう．

ψを因数分解のツールとして使おうという意図があるのだが，ψが<Kのような小さい数であるとしたら，ほとんど役に立たないのではないだろうか？むしろ最大のべきを探すべきなのではないか？AKS素数判定では多項式の剰余というのを使っているが，これがネックになっているのではないだろうか？素数判定だけならむしろトーシェント関数で直接判定した方が効率がよいような気がする．Wheel Factoringは確かに効果がある．種になる素数を増やせばそれだけ効率は改善されるものと思われる．現在は{2,3,5,7}を種としているが，それを{2,3,5,7,11}にしたらどうなるかを試してみたい．下記URLでは100万までの素数を列挙するというプログラムを公開しているので，参考にした．

https://keisan.casio.jp/exec/user/1648462730

▲ResidueTestでは素因数などが変化しない．SeedTestでは変化するが，Qなどは不変．cycleなども変化しない．⇒検定種別により，注目している項目が異なるということのようだが，変化しないところは空欄にしておいた方がよい．PrimeTestではInvertが上限オーバーしているため，何も変化しない．FactoringSubなどは実行されているが，値が同一のため動作していないように見える．InvertTestではB進表記だけは動く．PrimeTestではそれすら変化しない．

{2,3,5,7}のホイールの円周は2x3x5x7=210だ．この円周上には48個の素数点がある．この素数点を抽出する簡単なプログラムを書いて，

{1,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,
89,97,101,103,107,109,113,121,127,131,137,139,143,149,151,
157,163,167,169,173,179,181,187,191,193,197,199,209}

のような数列を得た．小さい素数については，

{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59,
61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 121,
127, 131, 137, 139, 143, 149, 151, 157, 163, 167, 169, 173,
179, 181, 187, 191, 193, 197, 199, 209}

を素数リストに追加し，残る部分にホイールを適用するという構成にしたので，ホイールのスタート地点は211になる．これを導入して従来の円周30のホイールと差し替えた．W30の素数点は{4, 2, 4, 2, 4, 6, 2, 6}の８個だが，W210には48個の素数点が載っている．安定動作しているので，多分これで行けると思う．どの程度高速化しているかはあまりはっきりしないが，改善されていることは確かだと思う．

これで大体収まったと思うのだが，An Efficient Factoring Algorithm by Repunit Number Methodに書いてあることとはだいぶ乖離が発生してしまった．そもそもψ数が一致しない．定義は同じだと思うのだが，数字が合わない．また，逆数検定に循環単位Uというのを導入しているが，その辺りも書いてあることと実装にはかなりの隔たりがある．EFAではψ数をψ(b,n)のように表記し，以下のように定義している．

Given mutually prime numbers n and b, there exist nonnegative integer x<n and positive integers y and z such that z=x*b^y \n…x, that reads n divides x*b^y, the quotient z, and the remainder x, in other words, n*z+x=x*b^y. A set of {x,y,z} exists infinitely. u-length function y(b,n) is defined as the minimum such number y

乗法群の位相と呼ばれているord_n(a)の定義は現在使っているψ数とほぼ一致すると思う．⇒EFAに戻って読み直しを始めた．

▲BuildMatrixで例外が発生した．N=2684380，K=6，e=4．⇒マトリックスのパラメータ設定に関わるエラーだ．「DataGridView コントロールの SelectionMode が ColumnHeaderSelect に設定されているとき、列の SortMode を Automatic に設定することはできません。」⇒BuildPowerGridでグリッドを初期化している．冒頭で.SelectionMode = DataGridViewSelectionMode.RowHeaderSelectを実行して解決．

▲Kを減少方向に変化させて，ゼロ除算例外が発生した．DispModPowの中で起きているようだ．K=1のとき，K-1=0になる．⇒対処した．

e=2でKを変化させると，eRの値は，K=1で不定，K=2と3ではeR=０になる．この動作は正しいか？⇒どうも，話が合わないような気がする．⇒2023-06-04のログで，MakeStripePatternで剰余演算によりゼロ除算が発生したという報告がある．このときは，①除数がゼロになる場合は除算を実行しない．②べきが０になった場合にはK-1を代入する，という２つの対策を実施している．MakeStripePatternでは，i=1→Kで{x=i^pow mod K}の数列を生成している．べきpowはべきがKより大きい場合にはpow=beki mod K-1 で剰余を取っている．

たとえば，K=4のとき，e=3, 6, 9ではeRは０になるが，マトリックスの指標としてみる場合には０という数字は現れないので，eR=0→eR=K-1のように換算する必要がある．ｎR=0にも同様の問題はあるが，行の場合は行ヘッダが表示されないので，剰余が１列目に表示されているが，それが０になっているので，特に手当しなくても問題ないようには見える．⇒縦軸は番号Nとみなすべきだから，やはりゼロは適当ではないと思う．eRの場合と同様0→Kのように読み替えるべきだ．⇒これでやっと辻褄があう．行番号Nと列番号eには０が現れないようになった．行番号のKというのが0の並びになる行だ．eRもnRも完全に断続なしでシーケンシャルに1→K/K-1の数字が並ぶようになった．

▲なぜだろう？以前は行ヘッダ部でホバリングすると行番号が表示されていたのだが，出なくなった．⇒間違えて .Rows(i).HeaderCell.Value = i.ToStringを書き換えていたようだ．それだけでなく，ウィンドウをオープンにしたままでBuildMatrixを実行すると，やはり出なくなってしまう．⇒ToolTipTextに設定したら表示できるようになった．これしかないだろう．⇒いや，もしかしたら不要かもしれない．テーブルが選択されていない場合は，単にホバリングしだけでは出ないのかもしれない．⇒やはりそのようだ．単なる操作上の問題だった．

▲BuildMatrixのサイズボックスを数値ボックスに変える．⇒対処した．

▲マトリックスの列ヘッダをクリックして列選択できるようにしたい．⇒ColumnHeaderMouseClickとRowHeaderMouseClickを使えばいけるのではないか？⇒できた．楽勝だ！

BuildMatrixではTestMatrixを実行しないようにしたので，マトリックスに関する手当もほとんど完了だ．あとは，Kが大きいときに表示するマトリックスの構築だけだ．これはそれほど難しくない．マトリックスの左上点を決めてやれば，あとはこれまで通りの手順で生成できる．マトリックスの左上点はeRとnRでよいと思う．もしかすると，多少違和感があるかもしれないが，ともかくやってみよう．

▲シアン色に塗りつぶししている，φ，ψ，φ(K)の枠線が赤くなっているのはなぜだろう？⇒枠線がFixed3Dになっているためと思われる．FixedSingleに変えておいた．一部のパーツはFixed3Dになっている．

「マトリックスの左上点をeRとnRとする」というのは一案だが，(eR, nR)がマトリックスの中に入っている場合は，(1, 1)からの描画でもよいのではないか？その方が分かり易いし，使い勝手もよいような気がする．(eR, nR)が外れする場合には，左上に持ってくるのではなく，マトリックスの中央にくるような配置の方がよいのではないか？まず，この方向で修正を入れてみよう．

▲マトリックスのサイズ上限を128としているが，描画までに相当な時間が掛かっている．もっと小さくてもよいのではないか？大体，eRが60くらいでフルスクリーンになってしまう．128の半分の６４としてみよう．⇒これで，描画までに７秒くらいだ．こんなものではないだろうか？マトリックスには数字がぎっちり詰まっているだけだから，あまり大きくしてもほとんど意味がない．行選択ないし，列選択するためには，ヘッダ部が見えている必要もある．いや，ヘッダ部はスクロールしても常時見えるようにはなっている…

▲アプリ終了時にパラメータを保存し起動時にロードする．⇒対処済み

▲InvertTestでは「■MakeRecursionUnit Abort for Digit Overflow」というエラーメッセージしか出力されない．⇒確かにInvertTestでは何も出力していないようだ．⇒InvertTestは単純にInvertFuncを実行し，DispInvertFuncで結果を出力しているだけなので，InvertFuncがエラーを返した場合にはやることが何もない状態になる．⇒古いコードを参考にとりあえず，DispInvertとDispRecurringDecimalを復活させてみた．⇒どこかでハングしてしまっている模様だ．ketaが負の場合は実行しないようにしてとりあえず動いた．⇒DispInvertで出力する循環節の最大長を1000まで縮減するようにして，ぼぼ満足できる動作になった．

▲N=123467890ではPrimeTestもエラーメッセージだけになってしまう．InvertTestも同様だ．これはやむを得ないと見るしかないかもしれない．Int64 MAXの上限の手前でSeedTestするとほぼハング状態になる．まだ何か抜けきれていないところがあるようだ．Factoringだ．N={9223372036854775701}．そろそろターゲットをFactoringに向ける必要がでてきたようだ．⇒NのFactoringをやっているDispParametorを後回しにしてみたが，効果なしだ．φ(n)の素因数分解ですでに手間取っている．⇒しかし，これは何とか抜けてきた．ψ数の計算量を見てみよう．⇒Psi()を呼んでいるところがない．なぜだろう？使っているのはDumpMatrixだけだ．それでは，ψはどこで出しているのだろう？

Psi()ではなく，PsiFunctionが使われている．Factoringにも，FactoringSubとFactoringSub2という２つのバージョンがある．DispModePowをDispParameterに入れて，FactoringSubを外に追い出すことにしよう．DispParameterは２箇所で使われている．9223372036854775701 ＝ 3^4 x 11^2 x 41064384945901だ．おかしい．41064384945901はまだ分解できる．

41064384945901 = 293 x 40151484457 だ．40151484457もまだ分解できる．40151484457 = 4934119 x 28589579．どうも，これでもまだ終わりでない可能性はある．FactoringSub2はもしかすると，すでに捨ててしまっているかもしれない… いや，まだ残してあった．⇒数字のコピーで間違えていた．41064384945901ではなくて，941064384945901だった．その次の数字も間違っている．40151484457ではなくて140151484457だ．

いずれにせよ，現行のFactoringSubの動作は正しい．また，FactoringSub2はそれに比較すると話にならないくらい遅い．実際のところ使い物にならない．しかし，肝心のψ数は空欄のままだ．一つψ数を求めるボタンを作ってみよう．

N=268435472，K=6のとき，#2でcycle=2,4,2になるという事例を考えてみよう．この値はn^の設定に関わらず一定だ．べきn^を変えるとstripeは変化するが，cycleは変化しない．これはcycleがマトリックスの行でstripeがマトリックスの列なのだから当然だろう．べきを２に固定してnを変化させると，Rと#が変化し連動してcycleが変化する．これも当然の動作ではあるが… Rと行が一致しない場合をカラー表示すると，

• N=268435472　R=4 → 2,4,2,4
• N=268435473　R=3 → 3,3,3,3
• N=268435474　R=4 → 4,4,4,4
• N=268435475　R=1 → 5,1,5,1
• N=268435476　R=0 → 0,0,0,0
• N=268435477　R=1 → 1,1,1,1
• N=268435478　R=4 → 2,4,2,4
• N=268435479　R=3 → 3,3,3,3
• N=268435480　R=4 → 4,4,4,4
• N=268435481　R=1 → 5,1,5,1
• N=268435482　R=0 → 0,0,0,0

のようになり，R=1と4が特殊な事例になっている．R=2と5はこの系列には現れない．剰余数列には周期性があり，長さ６の周期{4,3,4,1,0,1}を繰り返している．この周期列は縦周期列に現れるものとまったく同じだ．つまり，べき＝2, 4, 6の縦数列にこのパターンが現れている．これもある意味で当然だ．縦周期というのはそのようなものとして定義されているのだから．マトリックスはKによって一意に決定される．また，n^e mod K＝ R の値はeの値に依らず一定である．

それでは，ある値ｎとeを指定したときの cycle がマトリックスのどの行に該当するのかを決定するパラメータを何と見ればよいのだろう．上のリストで見ると 1,1,1,1 のcycleは268435477に対応している．これはKに依存した値と思われるので，Kの剰余を取ってみよう．確かにR=1となっている．しかし，R=1だからと言って1,1,1,1になる訳ではない．268435476と268435482の間にはR=1が２つあり，その一つは5,1,5,1だ．基準はN=1であり，６の周期で変化しているのだから，N=1とN=268435472の剰余的な距離が分かればよい．N=1 → R=1のはずだから，N=268435472→R＝4との差分を見ればよいはずだ．

いや，上ではR の値はeの値に依らず一定としているが，間違っている．たとえば，N=268435472^2 mod 6≡4だが，268435472^3 mod 6≡2だ．実際，マトリックスの第２行を見ると，2,4,2,4の数列になっている．しかし，e=2でR=4となるのは2だけではない．4も同様にR=4となっているが，この行は4,4,4,4という数列を生成する．言い換えると，Rとeだけではcycleは決定できない．決め手はなんだろう？ψは不定になっているが，ψを決定しても決め手にはならないだろう．ψは存在すれば#に一致するはずだが，e=2でも3でも#の値は2で変わらない．26843547２^1のRを見ると，R=2だ．N=1ならR=1だから，その差は１．つまり，N=2≡268435472 mod Kと言える．

べき乗する前のNの剰余を見ておけばよいのではないだろうか？この剰余をNrとすると，確かにNrはマトリックスの行を指定しているように思われる．N^e mod K ≡ R の剰余Rはeの値によって変化するが，Nrの値はeに依らず（当然だが）不変である．いまNをｎと表記しているので，この数字はnRとラベル付けしておこう．これでマトリックスとの対応が疑問の余地なく一致するようになった．つまり，マトリックス上の行はｎRで一意に決定され，マトリックス上の列はe=n^で決定できる．⇒完璧だ．もはや疑問の余地はない．

このツールには５個の検定が含まれているが，検定範囲を指定するのではなく，検定の回数を指定するという方式に変えておくことにする．この方が設定し易いのではないかと思う．

▲SeedTestの１回分でInvertFuncが三度呼び出されているのはなぜか？①SeedTestClickからの呼び出し，②ValueChanged x 2 ．SeedTestClickの中では明示的にValueChangedを実行しているのでInvertFuncの実行が重複してしまう．①passflagが立っているときは，UpDown.ValueChangedでValueChangeProをパスするようにした．

▲MatrixTest中にｎRがまったく更新されないのはなぜか？この値はDispModPowで更新している．DispModPowはValueChangedから呼ばれているのだが…BuildMatrixでは基本的に画面の更新を行っていない．その後にTestMatrixが実行されるが，これはテーブル全体のテストでNの書き換えなどは行われない．⇒TestMatrixではResidueFuncProの中でNの書き換えを実施している．TestMatrixでResidueFuncProを実行したあと，DispModPowを呼び出すようにした．これでｎRは更新されるようになったが，ｎの素因数列は止まったままだ．これはDispParametorで出している．これが入ると少し重くなるとは思うが…

▲PrimeTestでPowerResidueFuncが逆数モードで２回呼び出されている．ValueChangedから実行されている．passflagは立っているのだが，実行されている．passflagは整数なので０と比較するようにして動作するようになった．

▲PrimeTestでNを268435400に設定してInverteFuncの処理範囲に入るようにしたが，UやｎUがまったく表示されていない．エラーも出ていないように見えるのだが… 同じ設定でInvertTestを実行したら，ハングしてしまった．⇒DispInvertで相当手間取っている．ValueChangedから呼び出しているDispInvertFuncはInvertTestとDispInvertを含んでいるので，そっくりこの処理を呼び出すだけでよいと思われるが，大きいNはほとんど処理しきれないというのが実情だ．しかし，DispInvertFuncではそこまで時間が掛かっていない．Uが出るためにはかなりNが小さくないと処理できないようだ．処理可能なｎの最大値は268435455だ．

しかし，InvertTestでは，n=268435400 B=10でハングしてしまう．なぜだろう？ToBCDで手間取っているようだ．しかし，この関数はDispInvert→ConvertNum2String→ToBCDで，DispInvertはDispInvertFuncに含まれているはずなのだが… おかしい．Uが出なくなってしまった．N=2684354は十分小さい数のはずなのだが… もっと小さい数にすれば出る．MakeRecursionUnitで桁数オーバーでアボートが起きている．⇒DispInvertFuncならUやｎUは表示されないが，処理は進む．おそらく，DispInvertFuncで使っているルーチンは途中で打ち切るようにできているのだろう．⇒MakeRecursionUnitがUを返してこない場合にはConvertNum2Stringは実行されない．つまり，土壷にはまらないようになっている．これは差し替えでよいだろう．

n=268435455のとき，digitは1342176になっている．MaxOutputは30000なのでそれよりもだいぶ大きい．QTのサイズは2684392あるのでUの計算はできるはずだ．⇒計算自体は単純な計算なので実行は可能だが，Uは生値なので途方もない大きさになる．やはり，ある限度で打ち切った方が現実的なのではないか？桁数で30000というのはそれだけでも無茶無茶大きな数字だ．やはり，これは打ち切るのが正解としておこう．1/nの循環節までは出せるのでこれで満足するしかないと思う．

InvertTestではBをB+50まで変化させているが，特定のBでやけに重くなるのはなぜだろう？重くなるBには 42, 46, 57がある．おそらくこれらは何らかの理由で循環節が長くなっているのだろう．n=2684354のとき，B=42で fixed=1 keta=27962だ．この桁数は上限を超えていないので，ｎU，RPも出力されているが，Uが空欄になっている．なぜだろう？B=43 では，fixed=0 keta=1342176でアボートしている．B=42が遅いのではなく，打ち切りになっていないために見かけ上遅くなっているというだけだ．46, 57の場合も同じだ．ダンプしてみたが表示しようもないほど大きな数だ．

U=30979015254076365607484929108097110855736447
23737560540653941886057714834360774845095407791
71918951268321512524684829108915145925588584343
70873530921047046611821998461775454883886861647
25702974448815854451091871481243706996114903023
35610469271871095106330135472657551521990096397
75521223776067953620107407306115688619234917853
24243166297078829776098227164689…これが延々と続く…

ともかく，Uという数は手に負えないほど大きくなるということが分かった．⇒テキストボックスに入り切る範囲で打ち切ることで描画できるようになった．これ以上のことはできない．これで処理できる場合には限界まで表示できるようになった．1/nもたとえば，B=47の場合，1342176桁もあるのだから，どこかで打ち切っているのだろう．

▲MatrixTestのとき，(n, B)と(φ(n), @)が更新されない．現行では，(n, B)はDispInvertFuncでしか更新されていない．というか，MatrixTestはBとはまったく関わりのないテストになっているはずだから，更新されないというのが正しいのではないか？この２つはBの枠内に移してしまおう．⇒これで剰余演算と逆数演算が完全に分離できた．

▲ResidueTestで画面の更新が止まったあと，かなりの時間砂時計のままになる．なぜか？おそらく，最後の辺りで桁数の大きな処理が入っているためだろう．27962桁というのが，n={2684354}, B=42で起きている．例のB=42, 46, 57だ．この桁数上限はもう少し小さい方がよい．10000桁まで減らしてみたら，だいぶ楽になった．これでよい．

▲PrimeTestはBを変化させるテストになっているが，なにも更新されない．なぜか？仮にUなどが上限で描画できなかったとしても，1/nは変化してもよいのではないだろうか？DispInvertFuncもDispInvertも実行されていないためだ．⇒DispInvertFuncを実行するようにした．

▲n=2684351, n^2=2, K=2, B=42のとき，DispInvertFuncの下記の位置でエラーが発生した．

extBox13.Text=U.ToString.Substring(0, TextBox13.MaxLength)

現行ではTextBox13.MaxLength=32767だが，U.ToStringの長さがそれよりも短いときにはエラーが発生する．⇒対処した．

▲BuildMatrixでパラメータを書き換えるのはまずい．MatrixTestは別としても，BuildMatrixでは画面を変化させない方がよい．TestMatrixをパスしてもよいのではないか？

stripeの横にn^e mod K を表示する．正しくは，n^e mod (K-1)ではないか？⇒そのようだ．eRはstripeの列番号に相当し，nRはcycleの行番号に対応する．これで完全にストーリーが完結した．

素数判定にb^p-1 ≡ 1 mod p を導入しようとしているところだが，その前にIsPrimeの代わりにTotientFuncを呼んでいるところをIsPrimeに切り替えておこう．⇒いや，すでに整理されている．FactoringSubなどで使われていたが，すでに廃止されている．

素因数分解などで i の範囲をmaxにしているところがあるのではないか？素因数を取るためには√maxまで検査すれば十分である．⇒いや，その辺りは抜かりなく対処している．FactoringSubでは max = BSqrt(n) としている．TotientFuncではそのような手当は行っていないが，(x / d < d)でレンジチェックしているので問題なさそうだ．多少とも効率を改善することになったかどうかはあまり，はっきりしない．

とりあえず，これでN=9223372036854775807の範囲（Int64）では万遍なく動作することになったと言えると思う．ただし，1/Nでは配列サイズの上限があるので，MAXARRAYSIZE = 268435455より大きいNでは循環節などの表示はできない．同様の理由でcycle, stripeが表示不能になる場合もある．ともかく，時間は掛かるが停止しないで動作できるようになったのではないか？

Nが10桁のときのSeedTestを実行しているところだが，さすがに遅い．ネックになっているのは素因数分解のところと思われるが，取り敢えず今のところ打つ手はない．というか，ψ数を使って因数分解する話があったはずなのだが… しばらくはこれを使って An Efficient Factoring Algorithm by Repunit Number Methodを読み直しながら遊んでみることにしよう．銀林浩「初等整数論入門」の読み直しもやらなくては…

▲N=1234567890，B=1でPrimeTestを実行して，下記のエラーが出た．

B=1のときは，固定部がN桁になるため，配列サイズオーバーが起きていた．上限を設定して回避するようにした．また，このエラーメッセージを出力ペーンにも出力するようにした．nを調整して，MAXARRAYSIZE=268435455以下に設定すれば，正常に描画できる．

▲N=268435455で上限オーバーは発生しなくなったが，1/nやU，ｎUなどが空欄のままになっている．特に，1/nは更新されずに古いままの値が残っているのはまずい．⇒今度は問題なく描画されている．⇒PrimeTestではなく，InvertTestだったかもしれない．いや，更新された．ただし，かなり重い．⇒描画は出ている．⇒base2.ValueChangedイベントが掛かって，ValueChangedをダブって実行していた．

▲ValueChangedの中でパラメータを変更するのは問題がある．⇒ValueChangedの中で変更しているのではない．Inverse_ClickではBをある範囲で変化させながらテストしているのでこのような動作になることは不可避だ．

▲B=10，K=1でMatrixTestを実行して，ゼロ除算例外が発生した．障害はResidueFuncProで起きている．DispResidueCycle → MakeStripePatternだ．べきを小さくするためにK-1の剰余を取っていた．剰余を取る際にはゼロ割りに注意する必要がある．

▲cycleの出力がおかしい．#が０になっているためだろう．N=268435468，B=10, K=7で，n^7．この値はPowerResidueFuncの戻り値だ．⇒除数上限オーバーが起きていた．PowerResidueFuncの戻り値が負のとき，TextBox５～19までは****を表示しているが，cyclelen，dropnum，TextBox14は手当てされていない．

PowerResidueFuncは一つのセッションで二度呼ばれる．①剰余と②逆数は同じ関数で計算している．引数が異なるので，結果は同じにならない．剰余計算の場合はKが除数となるので，Nが巨大数でも問題なく計算できる．実際，正しい答えが出ている．⇒誤動作していると思ったのは，Kを見落としていたためだ．

一つかなり奇妙なことに気付いた．N=268435472，K=6のとき，#2でcycle=2,4,2になる．しかし，マトリックスを見ると，これは２の行だ．剰余R=4で行番号と一致しない．

マトリックスの最左列はNを示していると考えられるので，N mod K がRと一致しなくてはならないと考えられるのだが… 実際，別の数字で確認すると，

• N=268435472　R=4 → 2,4,2,4
• N=268435473　R=3 → 3,3,3,3
• N=268435474　R=4 → 4,4,4,4
• N=268435475　R=1 → 5,1,5,1
• N=268435476　R=0 → 0,0,0,0
• N=268435477　R=1 → 1,1,1,1
• N=268435478　R=4 → 2,4,2,4
• N=268435479　R=3 → 3,3,3,3
• N=268435480　R=4 → 4,4,4,4
• N=268435481　R=1 → 5,1,5,1
• N=268435482　R=0 → 0,0,0,0

のようにすべての行が出揃っているのだが，微妙に一致しないところがある．剰余数列は６パターンあるが，剰余は0, 1, 4, 3, 4, 1, 0で１と４が２回重複出現している．これはおそらくべきが掛かっているためと思われるが，それがどのように作用しているのかはこれだけではまだ分からない．いまのテストの設定ではべき＝２となっている．

べきが変化するとstripeが影響を受けるが，おそらくcycleの方にもなんらかの影響があるものと思われる．いずれにしても，この小さなマトリックスの中で起きていることなので，解析は可能だろう．これはもちろん，Kが素数であるか否かなどの要因もからんでいるものと推定されるが，整理できると思う．動作も安定してきたので，そろそろリリースの準備に入ってもよいのではないだろうか？

▲いつの間にか保存ファイル名がtotient.txtからnamechanged.txtに変わっている．VBのデザイナー画面では「べき剰余数列.txt」と入っているが，app.configにはtotient.txtと記録されているので，totient.txtというのがデフォルトファイル名だ．namechangedというのは，どこかの時点で設定したファイル名で，この値はTextBox16に入力→リターンキーで変更すると，My.Settings.filename = filenameに格納される．アプリ起動時にはMy.Settings.filename から取り出され，TextBox16に格納している．「namechange.txt」という文字列は現行のファイル上には存在していない… 確かに，前にそのような名前を設定したことはあるが，それがどこから出てくるのかはナゾだ．

リリース版とデバッグ版は別アプリということになっているようだ．それぞれが異なるファイル名を保持できる．

昨日，就寝直前に発生していた障害が再現できない．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

1. n=a^b(a,b∈N,b>1)という形で表せるなら、nは合成数
2. o_r(n)>log^2(n)となるような最小の数 r を求める
3. a≤r なる a のいずれかが 1<(a,n)<nを満たす場合、nは合成数 （ここでの (a,n) は、aと n の最大公約数）
4. n≤r ならば、n は素数
5. 1≤a≤[sqrt(ϕ(r))logn] の範囲のaに対して(X+a)^n ≢ X^n+a(mod X^r−1, n)ならば、nは合成数
6. 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では因数分解しなくてはならないのだから，意味がないかもしれない．

マトリックスのサイズを指定できるようにした．上限は128×128．ただし，テーブル自体はつねにフルサイズで生成する（と言っても物理的な上限はある）テーブルが生成してあれば，部分テーブルとしてのマトリックスはいつでも再構築可能になるので将来的にも拡張し易いと思う．

▲K=7のとき，３ｘ７を指定して，ゴミが出てしまう．

第４列に余分な数字が出ている．値は間違ってはいないが，指定範囲を超えている．これは，前のデータが残っているためではないだろうか？⇒Columns.ClearとRows.Clearで全データを消去するようにした．

▲K=7のとき，べきが６でstripeが{1,1,1,1,1,1,1}になる．これは間違っている．末尾は０にならなくてはならないはずだ．実際，マトリックスはそうなっているのだが…stripeはTextBox18に書き込まれる．この文字列はMakeStripePatternで生成している．この値はテーブルから取り出したものではなく，その場で生成している．マトリックスはつねに表示されているとは限らないためだ．テーブルは行単位に生成されている．つまり，T(i, 1) = i Mod kで左端のセルを設定し，その値を乗じて横数列を生成している．stripeは数列を上端の１から始めて縦方向に，ModPow(i, pow, K)で値を決定している．

i=k-1のとき，つまりpow＝k-1のときはつねにpow Mod (K – 1)はゼロになるので，すべての項でi^０＝１となり，１が連続することになる．このテーブルには何か決定的な不備があるのだろうか？かなり由々しき事態になってきた… ⇒解決策が見つかった．べきがゼロになった場合は，元の値に戻って計算すればよい．実際は除数がK-1でゼロになっているので，べき＝K-1でよいはずだ．これで縦横計算が完全に接続した．これは結局，フェルマーの小定理で n^p-1≡1 mod p の現象に対応している．ただし，n=pの場合は例外で，p^p-1≡0 mod p になるということだろう．ｎ＝ｐでなくても，nがｐの倍数であれば，当然0になる．つまり，これが，ｎとｐが互いに素という条件が必要という理由だ．

どうもまだ，間違っているところがある．1/7=0.142857142857…なので，U＝142857でなくてはならないが，428571が出ている．固定部１，循環部6となっているが，固定部０でなくてはならない．
https://www.facebook.com/groups/2354748741306929/posts/6031847776930322/?comment_id=6050062105108889&reply_comment_id=6057553787693054

循環小数の固定部はTextBox12に入っている．これを更新しているのはDispParametor2だが，DispParametor2はこの値をInvertFuncからもらっていると見られるので，結局PowerResidueFuncの動作がおかしいということになる．トレースしてみよう．IT(j)=0ならば，fixed=0となる．また，R = RT(max)ならばfixed=0となる可能性もあるかもしれない．R = RT(max)のとき，invertでなければ，fixed = IT(j) – 1となるので，IT(j)=1ならばfiexe=0となるのだが…

どうもこの不良は相当初期の頃から入り込んでしまっているようだ．上記のコメントは７週間前となっているので，50日前とすると，4月10日前後だ．ということは，やはり，久留島喜内 2023-04-16の時期と思われる．それ以前のバージョンんは数能極形術に入っているので，この版をきっかけに模様替えしたのだろう．

久留島喜内 2023-04-16からInvertFuncを切り出して現行版に組み込み，並行検査できるようにした．データを比較してみよう．QTは{1,4,2,8,5,7,0}となっている．IT, RTはここからは見えない．旧版でもQTはまったく同じ数列が格納されている．InvertFuncが使っているRS()を大域に置いて比較できるようにした．RSは{0,3,2,6,4,5,1,3,0,0,…}という内容になっている．InvertOldでも同様にRT()を見える場所に置いた．内容は{1,3,2,6,4,5,0,0,0}だ．大きな相違点はRT(0)に1が入っているという点だ．

RSは最近追加された出力用配列だ．オリジナルはRTでこちらはハッシュ表になっている．問題はR=1が検出されていないという点だろう．これは，初項が書き込まれていないためと見られる．旧版ではR=1がつねにR(0)に格納jされているため，①が検出されると周期が確定するという作りになっている．新版にはこの機構が欠けている．新版では全検索しなくても既出が存在しないことは判定できるため，書き込みだけで終わってしまっている．対応するとすれば，R=1のときは，R=0と同様，直ちにブレーク動作に移るべきだろう．

動作するようになった．最後に一つ，剰余Rが一致しないという点が残っている．⇒InvertFuncは参照でRを返しているが，PowerResidueFuncは明示的にRを返していないが，RSにはRの数列が入っている．これから取り出しておこう．

▲ｎの因数の動作がおかしい．ｎを変化させてもしばらく変化しないことがある．ｎ=2357の前後で停滞する感じ．ｎの因数はTextBox8に格納されている．イベントが途中で消失しているのではないだろうか？valuechanged中は入力ボックスをディセーブルにするという機構を復活させて動作するようになった．⇒特にBを更新するのにはかなり時間が掛かる．⇒最後の結果は合っているようなので，ディセーブルは掛けないようにする．途中で同期が取れなくなる場合もあるが，実害はないと思われる．不安ならば，ValueChangedをクリックすればよい．

ValueChangedボタンでInvertFuncが実行されていない．リターンキーなら実行される．⇒modulo2.ClickとUpDown.ValueChangedで同じ処理を実行するようにした．stripeはKが変化したときしか変化していないが，これでよいのだろうか？⇒Kが固定の場合，ｎが変化するとcycleが変化し，べきが変化するとstripeが変化する．MatrixTestではべきは固定なので，Kが変化したときしか変化しない．べきを変化させるテストもおもしろいとは思うが．トータルの時間が相当大きくなってしまう．

どうもまだ，完全にコントロールできていない．

晴れてくれたので，布団を干してシーツを洗濯した．大物洗いはこれで大体完了だ．とりあえず，久留島喜内の道具箱は大体仕上がってのではないかと思う．出力ペーンに書き出している情報は未整理だが，ぼちぼち片付けることにする．ここまでできたら，一度リリースしてもよいのではないかと思うのだが… ともかく，ダンプをファイルに取ってざっとチェックしてみることにしよう．Matrix Testに少し時間が掛かり過ぎているような気もするが，ダンプの量が過大なのではないか？

▲saveボタンを押してエラーパネルが表示された．設定はデスクトップのtotient.txtに保存するようになっているが，そのようなフォルダは存在しないという．これはデフォルトパスをc:\desktopとしているためだろう．デスクトップは個別のユーザフォルダの中にあるので，ログインしないと場所を特定できない．EXEの置いてあるフォルダはアプリを起動すれば特定可能だが，場所はユーザによってまちまちであり，デフォルトにはなり得ない．⇒起動時は空のままとし，ロード時に登録がない場合にはEXEのフォルダを指定するというのでよいのではないか？

▲ResidueTestでは何も保存されていない．SeedTestは動作している．⇒書き込みはWriteStreamという関数でやっている．

▲Saveがオンの状態でファイルを削除すると，ファイルは削除されて，そのあと実行しても結果は保存されない．⇒いや，勘違いだと思う．設定はEXEのフォルダだった．デスクトップのファイルは無関係だから削除できる．開いているときには，ファイル削除でパネルが出る．

出力パネルへの出力とファイルへの書き込みを完全に同期させたので漏れの可能性はゼロになった．あとは検定結果の画面出力だけなのだが，いま一つもの足りない感じがあるので，もう少し遊んでみることにしよう．次のような項目が考えられる．①整数ｎがノントーシェントであるか否かの判定，②剰余を入力可能とし，逆算してｎを決定する問題．どちらも結構難しそうな気はするが…

①に関しては以下のようなことが知られている．(1)１はノントーシェントではない，(2)１以外のすべての奇数はノントーシェントである，(3)ノントーシェントであるような偶数は無数にある，(4)ほとんどの数はノントーシェントである，(5)ｐを素数とするとｐ－１はノントーシェントではない，(6)２ｐがノントーシェントであることと，２ｐ＋１が合成数であることは同値である，(7)4p がノントーシェントであることと、2p + 1, 4p + 1 がともに合成数であることも同値である，(8) p − 1 で表される数はノントーシェントではない，(9) (p − 1)p の形で表される矩形数もノントーシェントではない，(10) p − 1 で表される異なる数同士の積もノントーシェントにはならない．

かなりのことは分かっているが，すべての偶数についてノントーシェントであるか否かを判定する手順は確立していない．分からないものについては，不明とするしかないのではないか？ノントーシェントであることが何か重大な意味を持ち得るのか否か？それとも特に意味もなく散乱しているだけなのかは調べてみなければ分からないのではないだろうか？その意味ではやってみるだけの価値はあるような気はするのだが…

②も難易度はかなり高そうな気はする．N^p=KQ+Rであるとして，これからNを逆算するとすれば，N=ｐ√(KQ+R) R＜K を解かなくてはならない．Kを固定し，pとQは任意でよいとしてもNを決定するのはかなり難しい．Kが素数であれば，0～K-1までの数がすべて剰余として現れるので問題は比較的簡単になる．同様にN^1なら同じことが言えるので必ず解を持つと言える．逆に言うと，N^{k＋1}はN^1と同じ剰余を持つから同様の解を得られるが，これでは単純過ぎるような気がする．

▲ValueChangedを実行しても，剰余数列/ψ/#などが変化しない．どこかで止めているのではないか？⇒剰余数列はTextBox14に書き込まれる．これを書き込んでいるのは，DispResidueCycleとResidueFuncだ．どちらもsilentを見ているので，このフラグが落ちているのではないか？いや，動いている．リリース版を試しているのでバージョンが古かったのかもしれない．⇒いや，多分勘違いだと思う．φ(ｎ)はｎによって定まるが，ψ(n, K)はｎとKのいずれかが変化したときの変化するので，べき数を変えただけでは変化しない．

ψがKを変えないと変化しないのは分かるが，べき乗を変えても周期が変わらないのはなぜだろう？これはおかしいような気がする．剰余周期と言っているのは，1/nの循環節ではなく，べきの周期なのだが… いや，これはべき乗を変えても変化しないというのが正しい．そもそも剰余周期というのは，べきが変化したときの剰余数列だ．従って，Nを変えれば当然変化する．これはマトリックスで行を切り替える操作に相当する．逆に言うと，Cycleというボックスに表示されているのは，マトリックスの行に他ならない．であるとすれば，べき数に対応した垂直の剰余列を表示するということも考えられるのではないか？⇒いや，それが欠けているとすれば，仕様的には片手落ちと言うべきだろう．

▲べき乗の九九表のタイトルにKの値を明示する．⇒対処した．

上の①，②に追加して，③nU/(B-1)を表示するということも考えられる．これは普通はレピュニット数になると考えられるが，そうならない場合もある．そうならない場合を調べる手がかりになるのではないか？もう一つの拡張として，Kの値に関わりなくテーブルを生成することが考えられる．つまり，テーブルを生成するKの大きさに制限を設けないという意味だが，テーブル全体を表示するのではなく，ある範囲，たとえば100×100くらいのテーブルを生成するのなら実用速度で表示できる．指定したNとべきがテーブルの中央に配置されるようにテーブルの範囲を決定すればよい．これはGoogle Earthがやっていることだ．これができるとほとんど制限なしで使うことができるようになるのではないか？理想型かもしれないが，追求してみたい．

ともかく，できるところからやってみることにしよう．③は簡単に実装できた．「べき乗の九九表のタイトルにKの値を明示」も片付けた．①のノントーシェントはとりあえず後回しということにして，べきに対応した縦数列を出すことを考えよう．テーブルが構成されていれば，それから拾い出すだけも作れるが，テーブルがつねに存在するという訳ではないし，また，その内容と設定が一致するとも限らないので，独自に計算を実施しなくてはならない．これはほぼ，PowerResidueFuncをもう一つ書くのに等しいが，やる必要があるだろう．後戻りできないように，先に出力ボックスを作っておこう．

さて，この数列を何と呼べばよいか？行方向の数列，つまりべき剰余数列は仮にcycleというラベルを付けているのだが，マトリックスの縦数列はまだあまりはっきりした位置付けが見えてこない．この数列にも周期性はある程度認められるのだが，大きな特徴はK番目のセル，つまりK行目のセルはすべてゼロになるという点だ．この横一列のゼロというのは，素数の場合はK行にしか現れないが，合成数の場合には（おそらく素因数のべきが入っている場合）中間にも入る場合があり，このゼロ線で周期になる場合もあり，ならない場合もある．

この縦数列で特徴的なのは，回文になる場合があるという点だ．これらのことを考慮すると，ここでは固定部：循環部という切り分けは難しいと考えられるので，とりあえず，K個の値をすべて表示するということにしておく．このような事情から，暫定的にこのボックスのラベルはgradationとすることにした．グラデーションの幅はKであり，そこに虹のように分解された色が見えるというイメージだ．いや，gradationではやや文学的過ぎるような気もする．stripeというのはどうだろう？この方が適切な気がする．この方が収まりもよいし，直感的だ．

一応動作するようになったが，どうもよく分からないところがある．stripeにもcycleとはまた違う種類の周期があるように見える．⇒原因は分かった．べきも剰余を取る必要がある．それもKではなくてK-1だ．このような操作を入れると始めてK-1の周期でstripeが循環するようになる．かなり複雑な（実は単純な）動作だ．かなりいろいろと面白いことが分かってきた．レピュニットはｎUが不規則な形になっている場合でもほとんど完全なレピュニットになる．例外は末尾に０が付いている場合だけだ．⇒オプションでテーブルのサイズを指定できるようにしたい．ただし，そこまでやるとすれば部分マトリックスを出せるようにする必要があるだろう．φ(K)も表示しておく必要がある．⇒対処した．