Wheel Factoring を導入して少し高速化

昨日,就寝直前に発生していた障害が再現できない.InvertFuncで配列のインデックスが範囲を超えているというエラーだ.リリースモードでは起きないバグなのだろうか?⇒Nを極端に大きくした場合に発生しているようだ.N=2024202420242024で発生した.10桁でも起きる.9桁では起きない.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までは問題なく処理できるようになった.約19桁まで処理できるということになる.これより大きい数字を入力できるようにするためには,数値ボックスを止めてテキストボックスに変える必要がある.数値ボックスについているUpDownボタンは結構便利なので捨て難いものがある.しばらくは,このままでよいのではないか?BigIntegerは理論的には上限が存在しないので,テキストボックスならもっとずっと大きい数でも入力できるようになるのだが…

φの約数を求める計算にも上限が入っているようだ.TextBox2を使っている.ValueChangedで10^10という制限を設けている.⇒外しても上限の9223372036854775807までそこそこに動作する.この制限は撤廃でよいと思う.⇒いや,動作はしているが,表示が空だ.どこかで処理を打ち切っているのではないか?⇒確かに,対象数が大きくなると相当な時間が掛かるようになる.{1537228672809129301}ではほとんどハング状態だ.⇒Nが18桁を超えるとφの約数列が表示されなくなる.テキストが長過ぎるのだろうか?MaxLegthは32767だ.⇒確かに上限を超えている.45338ある.MAXTEXTLENGTHという定数に32767が入っている.ここまでで打ち切るのが無難かもしれない.

確かにTotientFuncの時間消費は大だ.いや,TotientFuncというより,FactoringSubではないか?FactoringSubはTotientFuncを呼び出している.やはり,TotientFuncで時間を喰っている模様だ.ここではTotientFuncを素数判定のために使っているのだが,もう少し効率的なアルゴリズムがあるのではないだろうか?以前はエラトステネスの篩を使っていたのだが,廃止してしまった.大きい数を連続して扱うときにはこれがあった方が有利なのではないか?

対象数Nの平方根まで持てばよいのだから,MAXARRAYSIZE = 268435455なら,72,057,593,501,057,025までが対象数になる.つまり,32桁までカバーできる.これは悪くないような気がする.エラトステネスの篩を一度に作るのではなく,必要に応じて拡張してゆけばよいのではないか?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なので,28がその周期なのだが,つまり,素数はこの周期の上にしか現れないということになっているのだが… どういうことになるのか?⇒※

※上の記述には誤りがある.{4, 2, 4, 4, 6, 2, 6}ではなくて,{4, 2, 4, 2, 4, 6, 2, 6}だ.従って,周期も28ではなく,30になる.この件に関しては上記URLの筆者に「編集の提案」という形で問い合わせを送った.この30という数はかなり興味深い.FBの数学・物理etcグループに投稿したAMSのDaily Epsilon 課題でまだ正解(の証明)を得ていない問題(下記URL↓)があり,その中に出てくる数字が30であるからだ.「どなたか30が最大という証明を考えてください」というコメントを付けてあるが,今のところ応答はない.

https://www.facebook.com/groups/2354748741306929/permalink/6215621771886254/

18桁の数123456789012345680に1をプラスした数の描画が完了するまでの時間を計測してみよう.5分経過したが,まだ停止しない… もう少し小さい数で測ってみよう.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}としなくてはならない.

この点を修正しても,まだ従来論理と完全一致しない.新版では大きい因数が先になるという現象が起きている.これはなぜだろう?⇒周期が7のままになっていた.⇒これで問題は解消した.多少は速くなることを期待してもよいのではないだろうか?⇒123456789012345677をテストしてみよう.⇒5分経過したが,まだ抜けてこない.旧版は止めて新版のロジックだけを走らせているのだが…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のとき,3x7を指定して,ゴミが出てしまう.

image

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

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

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

どうもまだ,間違っているところがある.1/7=0.142857142857…なので,U=142857でなくてはならないが,428571が出ている.固定部1,循環部6となっているが,固定部0でなくてはならない.
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となるのだが…

どうもこの不良は相当初期の頃から入り込んでしまっているようだ.上記のコメントは7週間前となっているので,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の数列が入っている.これから取り出しておこう.

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

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

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

image

晴れたので,布団を干しシーツを洗濯した

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

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

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

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

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

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

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

②も難易度はかなり高そうな気はする.N^p=KQ+Rであるとして,これからNを逆算するとすれば,N=p√(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)はnによって定まるが,ψ(n, K)はnと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が循環するようになる.かなり複雑な(実は単純な)動作だ.かなりいろいろと面白いことが分かってきた.レピュニットはnUが不規則な形になっている場合でもほとんど完全なレピュニットになる.例外は末尾に0が付いている場合だけだ.⇒オプションでテーブルのサイズを指定できるようにしたい.ただし,そこまでやるとすれば部分マトリックスを出せるようにする必要があるだろう.φ(K)も表示しておく必要がある.⇒対処した.

夏向きに部屋を模様替えする

風呂に入って洗濯した.気温が高かったので三度続けて洗濯ものを乾かすことができた.一度目はレースのカーテン,二度目は厚地のカーテン,三度目は予備のシーツその他.間仕切りに使っていた大きなポンチョは以前浴室にあった白いビニールのカーテンに交換した.東側の出窓を塞いでいた棚は壁際に移動して窓を開放したのでだいぶ風通しがよくなった.今年もまた猛暑になりそうだが,これでなんとかしのげるだろう.石油ストーブを片付けて靴箱も整理した.最近はあまり歩いてないので普段はしまっておくことにする.

PowerResidueFuncが生成する剰余数列と従来論理のFindOutCycleで出力する剰余数列とがまったく一致しない.FindOutCycleの方が正しいと思われるので,PowerResidueFuncの論理がどこかで間違っていることになる.ともかく,これらが一致するようにしなくてはならない.FindOutCycleでは最初にRTをmaxまで初期生成している.つまり,R(i) = BigInteger.ModPow(num, i, K)だ.これは剰余数列を生成する手順そのものなので,疑問の余地はない.PowerResidueFuncでは論理を逆数周期列と共有するためにあれこれ小細工が入っている.この工作が仇となっているものと思われる.オリジナルではR()は{0, 4, 2, 1, 4, 2, 1}のようになっている.一方,改訂版で{0, 1, 2, 0 4, 0, 0, 4}だ.

これはオリジナルでは発生順にR配列に格納しているのに対し,改訂版ではハッシュ表になっているためだ.別途テーブルを作って受け渡しするようにした.これで齟齬はなくなった.ただし,現行では循環が検出された時点で打ち切っているため,数列の長さはKよりも短い.⇒正しく動作するようになった.あとは,基本的に動作確認と使い勝手の向上だけだと思う.⇒まだ表示されていない項目があった.(N, B)だ.これはTextBox4に入っている.⇒DispInvertで表示しているのだが,呼ばれていない.⇒DUMPNUMSTRINGという制限が掛かっていた.

一文字変数名をいくつか変更した.N→n, @→#,#→@など.#という記号は剰余類群の位数の意味で使われることがある.#は剰余数列の周期だが,ほぼ位数に一致するので,こちらで使うのが適切ではないかと思う.これまで使っていた@は循環小数の循環節の桁数の方に回すことにした.この方がなんとなく感じがつかめるのではないかと思う.

出力保存に関してはもう少し手を入れる必要がある.⇒DispInvertを止めたのは表示に時間が掛かり過ぎるからではないか?(N, B)だけを別途出すようにした方がよさそうだ.マトリックスは260を上限としていたが,Int64の範囲で対応することにした.ただし,配列には上限があり,それ以上にテーブルを構築するのに無闇に時間が掛かるので,テーブルは128 x 128 までとすることにした.これで一応実用時間内で大概のことができるようになり,後はファイル保存だけになった.

従来はデスクトップに固定ファイル名のログファイルを出力するという方式だったが,ファイル名とパス名を指定できるようにすることにしたので,それを格納しておく方法が必要になってきた.レジストリに書き込むか,ないし構成ファイルに保存するということになるが,レジストリのIOはすべてC++で書いているので,ちょっとここではアクセスできない.構成ファイルは常備されているが,従来のINIファイルはすでに廃れていて,XMLのapp.configというのを使うことになっている.

構成ファイルをプロジェクトに追加するところまではできたが,(と言っても結構つまづいた)form1_loadで読み込もうとすると,

Configuration system failed to initialize

のようなエラーになってしまう.ネットで探してみたが,なかなか出てこない.最後に以下のURLに

In the App.config file, I deleted all the contents of the ‘system.diagnostics’ tag. I no longer have an error detected in the Settings.Designer.vb file when reading the setting.

というアドバイスがあったので試したところ,一発で解決した.

https://github.com/dotnet/project-system/issues/6784

PowerResidueFuncの論理は間違っている

昨日の修正でいくつかの項目が表示できなくなっている.まず,それらを補充しておこう.TextBox14に表示していた剰余周期数列が出なくなったのは,ResidueFuncからPowerResidueFuncに切り替えたためだ.とりあえず,ResidueFuncで行っていた出力を切り出して関数化しておこう.DispResidueCycleという関数を起こした.一応動いているが,どうも間違っているような気がする.⇒かなり厄介な話になってきた.ResidueFuncが生成する剰余列とPowerResidueFuncが出す剰余列はまったく別物だ.これでは使い物にならない.ということは,PowerResidueFuncの論理が間違っているということになる.

もう一つおまけに,FBの数学物理談話室の問題に応用しようとしたところ,Kの範囲が狭すぎて使い物にならないことが分かった.少なくともInt64の範囲までは扱えなくてはならない.また,K≧2となっているが,1でもよいことにしたい.というのは,単純にN^mを求めたい場合があるからだ.⇒Kの最大値をNと同じ9223372036854775807に設定して動作するようになったが,剰余周期数列の計算が収束しない.これは打ち切りでよいのではないか?

どうも完全に仕上がったのではないだろうか

一応循環単位数Uを表示できるようにはなっているが,まだだいぶ問題がありそうだ.というか,この過程で肝心の桁数の計算もおかしくなっている気配がある.総合的な点検が必要だ.

N=3,B=10のとき,Uが空になっている.固定部1循環部1という数字で,k循環節も/A:0.3&0/のようになっている.1/3=0.3333なのだから,/A:0.&3/でなくてはならないはずだ.str3の表示も0.3000…のようになっているので,間違いだ.つまり,固定部0,循環部1でなくてはならない.⇒とりあえず,強引に収めた.

N=7, B=10のとき,固定部5, 循環部6でUが700000となっている.これは明らかに間違いだ.1/7=0.1428571428571428…で循環節は142857でなくてはならない.つまり,固定部は0のはずだ.⇒IT(j)=0を見る,つまり,初項に一致した場合はfixed=0で抜けるようにした.

N=14, B=10でnUが9999990となるが,これはどういうことか?⇒循環節Uが714285となるのは正しい.714285✕14=9999990となるので,計算は合っている.N=15の場合も,U=6で6✕15=90となり,末尾に0が付いてくる.N=24の場合,U=6でRepunitは144となる.1/N=0.0416666666…なので固定部は.041の3桁,循環部は6で1桁となるから,U=6で正しい.従って,6✕24=144という計算は合っている.N=28のとき,U=571428となっている.これも数字は合っている.従って,571428✕28=15,999,984は正しい.

一見した限りでは正しく動作しているように思われる.B=10の場合には,①nUが99999…のようなレピュニット型になる場合,②99999…の末尾に0が付く場合,③それ以外の3パターンに区分できる.これらのパターンが発生する条件については,別途調べることにして先に進もう.

PowerResidueFuncで新旧不一致が発生している.N=7, B=10でValueChangedPro→ResidueFuncPro(n, K, False)を実行しているところだ.⇒明らかにこれは改訂版の方が間違っている.従来版では周期4だが,改訂版ではゼロになっている.⇒初項をRT(max)に格納するというのは余分だったようだ.この論理を外して動作するようになった.

どうも完全に仕上がったのではないだろうか?⇒MakeRecursionUnitというのはかなりコストの掛かる処理だ.Uが巨大数になってしまう.N=12345678のとき,桁数は335616もある.10進100桁でも十分大きいのでとんでもない数だ.⇒Nが10桁まではそこそこに出せるようになった.11桁ではさすがに時間が掛かる.現行ではMAXARRAYSIZE = &H7FFFFFFとしている.InvertFuncで処理する数の上限を決めた方がよい.⇒除数が10進で10桁を超える場合はエラーで戻るようにした.

これで12桁くらいまでは処理できるようになったが,FactoringSubで呼ばれるTotientFuncがネックになっている.この関数のアルゴリズムは十分効率的ではあるが,対象数が大きくなるとかなり厳しい.TotientFuncは素数判定のために使っているのだが,何かもっと効率的な方法はないものだろうか?⇒現状では12桁というのが限界のように思われる.過去(2023/04/28)には20桁のサンプルをテストしていたこともあるのだが…何が原因でこれほど重くなってしまったのか?⇒Bの大きさにも関係している.Bが素数である方が有利に作用する.

現在時間が掛かっている処理は,φの約数を取り出す処理なので上限を設けるか打ち切ってもよいのではないか?⇒あれこれやりくりしてかなり速くなった.こんなものではないだろうか?

循環単位数(recursion unit)Uは計算できるようになった

循環単位数(recursion unit)Uは計算できるようになった.これにnを掛けたnUがレピュニット(の整数倍)になることを確認する.この数はB進で表記したいのだが… notationにはnをB進表記した数が表示されている.⇒ConvertNum2Stringという関数があった.⇒計算はできるようになったが,話が少し違うような気がする.B=10の場合,N=3→9,6→36,7→999999,9→9,11→990,12→36,13→9999990,14→9999990,15→90,17→99999999999999990,18→90,19→9999999999999999990

N=7を除いて,すべて999…の末尾に0が付いている.中には6や12の場合のように,9の倍数ではあるが9の並びではない数字さえ出てくる.⇒明らかにこの末尾の0はゴミと考えられる.N=21, B=10のとき,MakeRecursionUnitはfixed=1, keta=6で呼び出される.このとき,QTには{0,4,7,6,1,9,0,0,}のような値が入っている.このfixedの1は冒頭の0のことで,そこから桁数分取ると476190になる.桁数は間違っていないと思うので,fixed=1というのが間違っているのだろう.つまり,476190ではなく047619なのだと思う.これは結局,InvertFunc → PowerResidueFunc がまた間違っていることを意味する.

ごく初期のバージョンを除いてすべて間違っている.正しく動作しているのは,「久留島喜内 2023-04-16」だけだ.このバージョンではまだRTはハッシュ化されていないので,R = RT(j)ならば,単純にfixed = jとしている.しかし,動作がかなり違うような気がする.RTには{1, 10, 16, 13, 4, 19}が入っていて,R=1を検索しているため,j=0でヒットしている.これに対し,現行版ではR=RT(0)を見ているため,R=10を検索する動作になっている.現行版ではRは{10, 16, 13, 4, 19, 1, 10}のように発生しているので,初項10から始まることになる.

旧版と現行版ではかなり大きな違いがある.旧版ではxがnを超えるまではパスする仕組みがある.いや,そもそもRT(0)が違うのではないか?確かにそれはあるが,RT(0)に1を書き込んでも変化しない.いずれにしてもRT(0)に対する書き込みには疑問がある.RTはハッシュテーブルだから,ハッシュ値が0になるということは当然あり得るのに,そこに先住者がいるというのはかなりおかしい.少し考え直す必要がある.

少なくとも現行版ではRT(0)を考慮する必要はない,というより,考慮するのは間違いだと言える.もし,初期値をどこかに格納する必要があるのなら,RT(0)ではなく,むしろRT(max)だろう.ここならどこからも届かない安全な置き場所だ.しかし,現行論理ではいきなり If R = RT(0) Then となる.かなりまずいと思う.

べき剰余数列と循環小数節の統合

べき乗剰余数列の周期を求める計算と循環小数の循環節長を求める計算ではどちらも剰余演算になるのでアルゴリズム的にはかなり類似したものになっている.これを統合することは可能か?というのが目下の課題だ.すでにかなりのところまで煮詰まってはいるのだが,まだ完全に解決したとは言い難い.強い類似性はあるが,相違する点も大きい.実際,真逆の計算なのだから違って当然なのだが,それを一つの計算で間に合わせようというところにかなりの無理がある.ただし,もしそれができれば,この2つの計算の「関係」がかかりクリアに見えてくるのではないかという期待がある.

PowerResidueFuncというnの逆数の循環周期を求める関数をべき剰余列にも適用するというのが目標だ.逆数1/nの循環周期では基数Bのべきを法nで除した剰余RにBを乗ずるという計算,つまりR_i=R_(i-1)*B mod H を反復することで周期を決定しているが,除算の対象を剰余Rではなく,R_i=B^i mod K のように切り替えることでほぼ同等の計算が実現できるという考え方だ.問題はRが反復していることの判定に微妙な差異があるという点だ.その前に一つよくわからない例外が発生した.

SeedTestClickで例外が発生する.K=7でi=24のとき,num0.valueにi値を書き込むところでオーバーフロー例外が発生する.しかし,書き込もうとしている値は24という小さな数であり,整数範囲を超えているなどのことは考えられない.何が起きているのか訳が分からないというのが実情だ.num0の上限値はInt64の最大値を指定しているので,9223372036854775807という十二分に大きな数字だ.⇒障害の真因は不明だが,最終的にこのパーツを作り直すことで問題解消した.

大体動くようになったが,nがKの倍数の場合に従来論理と一致しない現象が残っている.剰余数列は{0, 0}のようになるが,現行ではこれを固定部0,循環部1と解釈している.しかし,これは少しおかしいような気もする.循環小数としてみると,1/7を7進表記した場合には0.1となり,固定部1,循環部0となる.つまり,小数ではゼロが発生した時点で打ち切っている.剰余数列もそのような解釈でよいとすれば,固定部1:循環部0というのが正解になるような気がする.

たとえば,N=6,K=8の場合,剰余周期は{6*,4*,0,0}で固定部2,循環部1となっているが,これは固定部2,循環部0でよいのではないだろうか?だとすれば,剰余周期の表示も{6*,4*}のようになるはずだ.これは言い換えれば,統合版の論理の方が正しいということになる.もし,これでよいとすれば差し替えて終わりということになるのだが,従来論理を見直して一応修正しておこう.従来論理はResidueFuncから切り離して,FindOutCycleという独立の関数になっている.

FindOutCycleでは最初に剰余数列をストレートに生成し,末尾項から遡って末項の出現位置を探索,そこまでのステップを「周期」として確定した後,今度は初項から末項までドロップ項を探してカウントしているが,いくつかの点において問題がある.今の場合は,剰余数列には0項しか入っていないので,それを周期カウントに含めるのは問題がある.逆にドロップ項の検査では0項は対象に含まれていない(これは正当かもしれないが…)末項が0の場合は循環しないと断定してよいのだろうか?それは,当然だろう.従って,末項が0の場合は非ゼロ項が末項になり,そこから初項までが固定部ということになる.今の事例のようにすべてゼロ項の場合は固定部1ということになる.⇒対処した.

K=8 N=2のとき,新旧で固定部の不一致が起きた.現行では2,改訂版では3になっている.剰余数列は{2*,4*,0,}で末尾の0は循環には入らないから,固定部2桁というのが正しい.従って,改訂版の論理が間違っていることになる.しかし,これは末尾の0を固定部に入れてしまうという考え方もあるのではないか?というのは,完全に割り切れる場合には{0*}という固定部が生じるのは避けられないと考えられるからだ.言い換えると,末項の{0}は明示的に固定部に含めるという考え方でよいのではないか?もし,それが通用するのであれば,改訂版が正しいので,現行版を修正する必要があるということになる.FindOutCycleを修正してこの件は解決したが,まだ不一致が残っている.

K=12,N=2のとき固定部の不一致が生じる.改訂版では2,現行版では1になっている.循環部はどちらも2だ.剰余数列は{2*,4,8,4,}のようになるので,固定部は1というのが正しい.従って,現行版の方が正しく,改訂版は修正が必要ということになる.⇒一応動くようになったが,確信は持てない.もう少し整理する必要があるような気もする.ともかく,少し動かしてみよう.

どうも大体仕上がってしまったようだ.パラメータを一定範囲で変化させる包括検定を実施しても停止しないようになった.もう少し整理できそうな気もするが,先に進むことにしよう.試してみたいことが2つある.①小数循環節を整数化してその性質を調べる,②小数循環節に関連するレピュニット数を生成してその性質を調べる,の2点だ.まず,①から取り掛かることにしよう.⇒Recursion Unit U は10進数として表記する.(B進数として表記すると1/Nの循環節とダブってしまう).

レピュニット数の議論ではoriginatorという用語が出てくる.

B=1の場合の逆数の進数表記をどうするか?

N=7, B=1で起動するとInvertFuncで例外が発生する.QTなどの配列のRedimをPowerResidueFuncに移してしまったため,ここではQTにアクセスすることはできない.暫定的にQTのRedimはInvertFuncで実施することにしておく.これで例外の発生は抑制されるようになったが,問題は残っている.B=1の場合の逆数の進数表記をどうするか?という問題だ.2023/04/26には以下のような記述がある.

b=1というのはかなり特殊な場合で,ペアノの公理の後者関数のような形式になっている.このとき,b^ψ=1^1=1となり,いかなるnに対してもψ=1となると考えられるが,1進数の1/nは0.000…1 のような固定桁のみの数字になるため,k=0となり,ψと一致しない.gcd(n, b)=gcd(n, 1)=1でψが有効となるケースは基本的にψ=k=gcd(φ, k)になると考えられるが,b=1はそこでψ≠kとなる唯一のケースと考えられる.また,このとき,gcd(φ, k)=gcd(φ, 0)=φとなるため,三者がすべて異なる値を持つという特殊ケースになる.

とりあえず,これでよいということにしておこう.実装もそのようになっている.さて,このPowerResidueFuncを使い回しできるかどうか?が問題だ.つまり,この関数を使ってResidueFuncを実現できるかどうか?言い換えれば,InvertFuncとResidueFuncはどこが共通で,どこが違うのかを明らかにする必要がある.これらが一種の逆関数であることは明らかだが,剰余演算として実装されているところに共通点がある.しかし,もちろん,大きな違いがある.一方の対象がNであるとして,他方の対象は1/Nだ.これを共通処理化可能だろうか?

どうもかなり難しそうだ.その前に一度ResidueFuncをPowerResidueFuncのように抽象化してみてはどうか?そうすれば,もう少し見えてくるものもあるのではないか?

N=3, B=10で固定部1になっている.これはおかしい.1/3=0.333なのだから,固定部0でなくてはならない.どうも直近の修正で壊してしまったようだ.⇒PowerResidueFuncのループの入口にあったはずの「RT(0) = B Mod n」が消えてしまっている.RTを外部で初期化するように修正した際に削ってしまった模様だ.

ひとつ不良が出ている

ひとつ不良が出ている.n=7のとき1/nで固定部が1となっている.これは誤りだ.固定部:0でなくてはならない.10進数表記すると1/7=0.142857142857…で,循環節U=142857でnU = 142857*7 = 999999 = 9*111111となる.以下のURLによると,

An Efficient Factoring Algorithm by Repunit Number Method
http://www.aya.or.jp/~babalabo/repunit/index.html?fbclid=IwAR0HOaEACHvZ4MyWwXNk6_Dj4Tnr_BAJTjjLMKuBY3PyyhdLo61SbvqT3tY

[8′] 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.

[8′] は証明されているので定理だが,これを拡張した [9] Let m be an integer without factor 2 nor 5, and the recursion unit of m be U, then mU/9 is a repunit number Rn, i.e., mU=10n-1.も定理である.

今のケースではn=7は素数であるから,これを満たさなくてはならない.今回はB進⇔10進の相互変換ができるようになっているので,それだけで十分だろうと思ってレピュニット関係のロジックはすべて削ってしまったのだが,ちょっと早まったかもしれない.ともかく,どこかが間違っていることは確かなので調べてみよう.循環節はTextBox5に書き込まれている.この文字列はConvertNum2Stringで生成しているが,文字列自体には誤りはない.maxfixedに入っている固定桁数が間違っている.というか,maxfixedをここで使っているのが間違っている.

少なくとも初期バージョンでは正しい値が出力されている.⇒正しく出力されるのはバックアップの中で一番古い「久留島喜内 2023-04-16」しかない.それから一つ上がると固定桁:3など間違った値が出始める.この版ではDispParametorの引数で渡されたfixed値をそのまま使っている.現行版ではMakeRecurringDecimalで渡されるfixedの値がすでに誤っている.InvertFuncが返している値が1になっている.fixed = IT(j)でIT()には,{6, 2, 1, 4, 5, 3}が入っている.IT(0)には値が設定されていない.どうも,InvertFuncの論理全体を見直す必要が出てきた.

InvertFuncで剰余が一致したときの操作が間違っている.初項に戻る場合はつねに剰余が1になっているという推定が誤りだ.RT(0)には初項の剰余が格納されているので,それと比較する必要がある.これで多分問題は解決したのではないかと思う.repunitに関する検査を追加しておこう.どこに入れればよいか?レピュニットが関係してくるのは逆数検定に限られている.Prime TestにはDivideRepunitが残っている.

SeedTestを実行して,ValueChanged文字列 str3とstr4の不一致で停止した.10新表記の場合は両者が一致する必要がある.確かに小数第12位で不一致が発生している.

“0.0714285714285714″
“0.071428071428071428071428071428071428″

N=14, B=10で起きている.714285の方が正しい.つまり,str3は間違っている.MakeRecurringDecimalで生成するstr2にも欠陥がある.これしか入っていない.”/A:0.&071428/”.str3を生成しているのはMakeDecimalStringだ.この関数はQT配列から文字列を生成している.QTに入っている値は正しい.引数でfixedとketaを渡しているが,fixedに0が入っている.これは1でなくてはならない.上の修正で初項に一致した場合には強制的に1を代入しているが,これがまずいのだろう.

固定項には0の並びが含まれる場合がある.それをどうやって判定したらよいのだろう?QT(0)=0の場合はfixed=1となるようにして逃げた.とりあえずこれで動作するようになったが,これですべてカバーできるようになっているかどうかは疑問だ.どうもうまく行っていないようだ.n=6, B=10のとき,str3=0.11111111111111… のようになってしまう.R=4でIT(j)=1かつR=RT(0)が成立し,QT(0) が0でないのでfixed=0に設定している.1/6=0.16666…なので,固定部は1でなくてはならないのだが,誤動作している.

どうも,InvertFuncのアルゴリズムは仕上がっていないのではないだろうか?InvertFuncではQT, RT, ITの3つの配列を使っている.ここではB^i mod N を計算し,剰余が一致すれば周期が閉じたと判定して,固定部長と周期を割り出すことになっている.ループからの脱出条件は,①剰余がゼロになった,②一致する剰余を検出,③中止ボタンないしカウントオーバーの3つだ.問題が起きているのは,②の周期の終端に達して剰余が一致したという場合だ.ResidueFuncではドロップ項というのを見ているが,その操作が欠けているのではないだろうか?つまり,剰余が一致したというだけでは不十分なのではないか?

それを確認するためには,QTを見て,同一Qの発生時点まで進む必要があるように思われる.⇒できたようだ.ポイントは2つ.①QT(0)=0つまり,商が立たない場合には商が立つところまで進む必要があること,②異なる商が立っていた場合には,商が一致するところまで進む必要があること.これで,ドロップ項を見なくても固定部を正しく同定することができるようになった.この論理は剰余数列の方にも使えるのではないだろうか?もし,それが可能ならロスタイムを少しでも減少される効果はあるかもしれない.ともかく,一度バックアップを取っておこう.

N=14,B=16のとき,PrimeTestで停止した.fixed > 0 And i <> fixed + ketaが起きている.10進なら1/14=0.07142857142857…だが,16進なのでどういうことになるのか… i=4でfixed=6, keta=6になっている.fixedの6というのは明らかに誤りと思われる.多分,これは正しいのではないかと思う.str3=0.12492492492492492 だとすれば,固定部1,循環部3ということになる.1+3=4=iで計算も合う.

QT(0)=1で最初から商は立っている.Qは9なので,上からQTを探して,QT(3)=9を検出し,fixed=3としている.どこが間違っているのか?⇒QT(0)が非ゼロのとき,Qを探しに行ってオーバーランしている.i の範囲を超えたときには,fixed=1としなくてはならない.なぜなら,IT(j) = 1 And R = RT(0)であるから,初項から循環に入ることは間違いではないが,循環が完成するためにはRだけでなく,Qも一致する必要があるが,周期が短いときには,初項で発生するQがまだどこにも出現していない場合がある.そのような場合でも,初項のRが一致していることは初項の次の項から循環が始まることを意味しているからだ.

これでまず,穴は全部塞げたのではないかと思う.検算を兼ねてレピュニット関係の検定を導入しておくことは意味があると思う.まず最初に冒頭に出てきたnUというのを作ってみよう.nと除数pが互いに素の場合には循環周期Uをn倍したもの/p-1がレピュニット数になるというものだ.つまり,nのbベースの循環節をUとするとき,nとbが互いに素であれば,nU=n^(b-1)/(b-1).実装されているレピュニット関係の関数を集めておこう.まず,現行バージョンには以下がある.

  1. DivideRepunit ∑B^i mod n ある桁数のレピュニット数のnを法とする剰余
  2. RepuitFunc ∑B^i ある桁数のレピュニット数を求める

オリジナルソースを見てみたが,これですべてだ.つまり,大したことはやっていない.DivideRepunit はレピュニット数が桁数で割り切れるかどうかを見ているだけだ.これはちょっと後回しにして,InvertFuncとResidueFuncで同じ論理が使えるかどうかを見ておくことにしよう.まず,InvertFuncから主要ロジックを切り出してみよう.この処理の主眼は,nとbを指定して固定桁と循環部の長さを決定することだ.

N=7, B=2のとき,PowerResidueFuncで i <> fixed + keta のエラーになった.i=4 <> fixed=2+keta=3=5.N=7は2進表記では111だが,1/7は0.00100100100…のようになるので,fixed=2, keta=3が正しいように思われる.⇒ fixed ≠ IT(j) となる場合はあり得る.今の場合がそうだ.QT(0)=0の後に0が続く場合はそのようなことが起こる.⇒このような場合には(フラグを立てて)パスするようにした.

▲n=7, B=1のとき,ValueChangedで停止した.If QT(i) > Bというエラーだ.keta=0でfixed=6だ.B=1の場合はQT()は1の並びにならなくてはならないはずなのだが… ⇒従来論理ではEXITINVERTFUNCでQTの補充というのをやっているが,暫定的に止めてある.また,B=1, N=1の場合にはPowerResidueFuncを通さないでストレートにEXITINVERTFUNCに飛ぶようになっているので,QTが前の状態のまま残留しているのだろう.⇒どうも,これはかなり難しい問題だ.