進捗をみるためにプログレスバーを導入した

プログレスバーを導入したのだが,どうも動作が思わしくない.初期起動で少しだけ進んだ状態のまま止まってしまう.⇒進捗率6%で抜けてしまっている.⇒ループから横道に抜けていた.基数の情報(約数)などを知りたいのだが,常時画面に出すのではなく,Invert出力でよいのではないか?表示されるべき項目には,基数:素因数分解,nとのGCM,bのレピュニット数など.

CASIOのサイトでφ関数が計算できる.その結果と出力が一致していない.n=1234567のとき,CASIOではφ=1224720となっているが,KINAIでは41088が出力されている.Wolframでも同じ数を出しているので,こちらが間違っているのだろう.φの計算はFukuzo氏に紹介されたものをそのまま使っているつもりなのだが…オリジナルのコードは残っているだろうか?⇒オリジナルのコードに戻してみたが,結果は同じだ.1234567は1, 127, 9721, 1234567しか約数を持っていない.従って,1と1234567を除けば,127, 9721, 127×9721の3つしかない.いや,これらの倍数で1234567より小さいものも含まれるから,実数はもっと大きい.確かに41088という数字は小さ過ぎるような気がする.

オリジナルのコードは新潟大の竹内研究室からのコピーだし,最近のものは,奥村晴彦著『c言語による「アルゴリズム辞典」』が出典だ.2つとも間違っているというのはかなり考えづらい.加工の仕方が悪いのだろうか?確かに古いコードでLongという整数型を使っているのは問題かもしれない.それもあるが,PFactorをsieveを使わないように改造したことも影響している.PFactorは小さく約数から返すようになっているが,逆に大きい約数から返すようにしないと動作しないはずだ.オリジナルコードは廃棄してしまったが,sieveを使うようにするか,ないし,単純で非効率なコードに書き換えるしかない.

奥村のコードは少なくとも整数範囲では正しく動作しているので,こちらを改造して対応することにする.⇒なぜだろう?今度は正しい値がでるようになった.修正箇所としては,パラメータを表示している論理をInvertFuncの外に出したといことしかしていないのだが…いや,原因は分かった.InvertFuncは無茶苦茶時間の掛かる処理でまだφを計算するところまで進んでいなかったというだけの話だ.プログレスバーを取り付けてみたが,数分実行しても値はゼロからまったく変化していない.⇒カウンタを10000で割った剰余の1/100を出すようにして動作するようになった.10000点が1サイクルに当たる.

この論理の高速化は可能だろうか?この計算は循環周期が短ければかなり簡単に終わる場合もあり,周期がnに近いような場合にはとてつもなく時間が掛かってしまう.これは避けられない.⇒10000を1000に短縮した.この方が計算が速そうに感じるので心理的効果はある.これだと変化が速いので,10000点くらいまでは我慢できそうだ.いや,せいぜい5000点くらいかもしれない.値の更新では途中で打ち切ることにする.3000くらいで十分なのではないだろうか?⇒いや,もっと短くてもよい.大きい数字の場合にはInvertを使うという仕様だ.

逆数周期表示を見易くするために,記号の前後に空白文字を入れてみたが,コピーするとき切れてしまうので,空白を入れないように変更した.また,この行に固定部と循環部の桁数を出すようにした.これらの値はInvertFuncが最後まで実行されないと与えられない.処理が中断された場合には0, 0 となる.

だいぶいい感じになってきた

だいぶいい感じになってきた.特に,B進数⇔10進数の相互変換が自在にできるようになったのがうれしい.これは相当役に立つと思う.循環小数は基数が決まらないと意味がないので,基数を変化させることができることはその辺りを調査するときに使うツールとしては必須の機能だ.ここまでやると,Invertを実行しなくても,逆数の循環単位などを常時表示できるようにしたくなってくる.整数Nの逆数である循環小数Dの固定部Xと循環部Yを2つの整数とみなすと,XのB進数表記の桁数をx,Yの桁数をyとすれば,0≦x, y<N,0≦X<B^x,0≦Y<B^yで,

D = X / b^x + Y * ρ(x, y) (1)

ここで ρ は ρ(x, y) = Σ{k=0→∞}(1 / {B^y}^k) / b^x のような無限級数.(1)式はどんな小数にも適用できる等式だが,桁数 x, y が不定というところが厄介なところだ.4093という素数をB=11で逆数を取ったところ,循環単位の桁数が2046になった.ρ(x, y) は無限等比級数で,r=1/B^y < 1 により収束する.Σ{n=1→∞}a_n = a_1/(1 – r) .これが正しいとすると,1 – r という数がどのようなものであるのかが気になる.この値を計算するには,B^y という数を決定し,その逆数を取って,それと1の差分を計算する必要がある.どうもこの当たりが鍵になりそうな気がするので,実装を試みることにしよう.

電卓としても使えるようにn^pの計算で剰余だけでなく計算数値も出すようにした.循環小数値をX&Yの形式で表記することを考えてみる.たとえば,252_10であれば,00&390625のようになる.つまり,1/252 = 00&390625(10)と表記できる.この表記では頭のゼロは落とせないので,/00/390625/10とするのがよいのではないだろうか?あるいは,10/00/390625の方がベターかもしれない.⇒このような値を与えられたとき,これから1/252という元の数字を導くことができるだろうか?

InverseClickでエラーが発生した.n=255, b=12のとき,固定桁 fixed=1, 循環桁 k=16で,1+16*2=33桁の数字列を取り出したあと,zerosup = keta – nnstr.Length が –4 という負値になったため, New String(“0″c, zerosup)でエラーが発生した.nnに入っている値は,{1608573608807851864064416092523532}という大きなもので,b^keta / n の値だが,ConvertNum2String でテキストに変換したとき,桁数が37しか取れなかった.

ConvertNum2Stringが誤動作しているようだ.間にNULLキャラが入っているように思われる.いや,”A8”とか,”BA3”のような複数文字がつながった部分がある.⇒エンコードの問題だ.UTF8を使っていたが,UNICODEに変更して動作するようになった.

System.Text.Encoding.Unicode.GetString(BitConverter.GetBytes(c))

SeedTestでWRONG TOTIENTが出ているが,これは何だろう?「pk <> k And tot > 1」という理由だ.tot はφ関数,kは循環単位長,pkはtotとkのGCMだ.SeedTestは1~Nの範囲で検定を行い,整数の属性をチェックしている.pk=2, k=0, tot=2だ.対象数は4で基数は10.⇒ツールで常時φやψが見えるようにしておこう.

InverseClickで「R <> b And R <> 1」で停止した.R=0, b=20, b=128.⇒今度は,R=4, b=10, n=12 で停止した.InvertFuncではまずx=b=10として,x\n=10\12=0…10.次にx=x*b=10*10=100として,100\12=8…4,次に,x=R=4で,x*b=4*10=40\12=3…4.このRは既出だから,循環が閉じて復帰している.この計算は正しい.つまり,循環が閉じるときのRはbや1に限った話ではない.従って,ここで停止する必要はないと言ってよい.

ファイルを開こうとして,競合エラーが発生する.これはおもしろくない.⇒エラートラップを仕掛けてresumeするようにした.⇒ファイルに保存をオプションとする.

InverseClickで「ERROR:Count cannot be less than zero. (Parameter ‘count’)」というエラーになった.zerosupが-1になっている.上で出ていた症状だ.nnstr.Length=2でketa=1のため,差分がマイナスになってしまう.nnstrには”10”という文字列が入っている.nnも10だ.fixed=1でkが0であるため,keta=1+2*0で1になっている.n=1という数が入っている.n=1というのは,1/nが1になるというかなり特殊な数字だ.”循環小数=0.” + nnstrとあるように,ここでは1/nが1より大きくなることを想定していない.⇒暫定的にn=1の場合はエラーを回避して抜けるようにした.

▲φφというパラメータをチェックする必要があるだろうか?この値はφ(φ())で,nのトーシェント関数のトーシェントだから,nを変えない限り変化しない.

▲ψがゼロになる場合がある.このようなときは,循環桁数=φになる.通常は,φ≧ψでψ=循環桁数になっているのだが…これは調べておく必要がある.どうもBの約数と関係があるようなので,Bの約数も表示できるようにするとよいのだが…

ToBCDで例外が発生している.以下の計算式で発生する.

Dim byteCount As BigInteger = Math.Log(val) / Math.Log(b) + 1

n=255と小さな数字だが,逆数を整数化した数は巨大数になる.3921568627450980392156862745098 ToBCDでは冒頭で配列サイズを決めるためにMath.Logを使っている.対数計算しないように方式変更した.⇒最初にInt16.MaxValueサイズの配列を確保し,これを使って直接n\bの剰余(BCD値)を格納して桁数をカウントするようにした.

だいぶ整ってきた

だいぶ整ってきた.

image

対象数値ないし除数,基数を変えると自動的にボックス内の数値を更新するようになっている.除算を実行して余りがある場合にはGCMは1より大きく,割り切れるときにはGCMは除数そのものになる.notationの欄ではB進表記した数を入出力できるようにしたいのだが,まだ実装されていない.画面下部の剰余数列のブロックではkないしnをある範囲で変化させて周期数列を生成しているが,単点のn^p mod k という値を取り出したい場合もあるので,個別に計算できるようにした.ボタンの数は大幅に減って上段では,①PrimeTest, ②SeedTest, ③Invert と④Clear だけになった.これらの処理の内容はまだかなり整理が必要だが,大枠は完全に固まったといってよい.

対象数値をLongの範囲とするという制限は撤廃することにした.どうしても制限が必要な場合は出てくるかもしれないが,原則として無制限としておきたい.重くなる原因としてダンプ文字列が長くなることがあるかもしれない.コラッツのときは文字数制限してある分量を超えると捨てるようにしていたが,かなり効果があった.ただし,ファイル出力を選べばダンプをファイルに保存できるようになっていたので,その機能は必要かもしれない.出力パネルはもう少し横幅があった方がよい.

InvertFuncはPrimeTest, SeedTestからも呼び出している.ダンプの量が多いので,Invertコマンドの場合以外はこれらの値をダンプしないようにした.ただし,ファイルには出力する.逆数の循環数列を設定パネル上に表示して Invertボタンを廃止することも考えたが,対象数が4桁を超えると相当な桁数になってしまうので,現実的ではない.詳細はInvertボタンで表示し,PrimeTestやSeedTestでは隠蔽するというのが妥当だ思う.ディテールはまだまだ整備しなくてはならないが,最後の関門 notationへの変換・逆変換をやっつけてしまうことにしよう.

Convert.ToStringは基数が2, 8, 10, 16の場合しか扱わない.それ以外では例外が発生する.これは自前で作るしかない.ToBCDという関数も思ったような動作にならないので,丸ごと書き換えた.⇒動くようになった.かなりおもしろい.基数の範囲は2~62とした.アルファベットのA~Zとa~zに別々の値を与えて,最大でzzzzzzz…のような数まで表示できる.あとは,これを表記から数値に変換する処理を組み込むだけになった.ただし,こちらはユーザ入力から数値に変換するので,エラーに対処する必要がある.

いや,このエラーは無視してもよい.特に動作上のエラーにはならない.たとえば,2進数列に2という数字が紛れ込んだ場合には,2×2の値がその桁にあるとして計算するだけだ.10進に戻って,もう一度2進で表示すると正しい2進数列になっている.

16という数が入力できない.(6をBSで消去して再入力しようとして)1はBのレンジに入っていないため,弾かれてしまうようだ.⇒1を範囲に加え無処理で抜けるようにした.

N進数への変換・逆変換はほぼ仕上がったようだ.基数は最終的に2~64まで対応ということにした.62でzまで使い切るので,64の場合には{と|が余分に使われる.もし,その進数の範囲を超える文字が入力された場合には,直ちに変換して正しい数字列に書き換わるので,その文字が使えないことを指摘しなくてもよいだろう.あと,残った課題は,①3種の検定の出力を整理して使い易いものにする.②巨大数でハングしないように適切な例外処理を入れる,の2点だけになった.そのあとは,このツールをどう活用するか?というフェーズに入るが,おそらく因数分解の効率化というところから入ることになるのではないかと思う.相当強力なツールになったので,役に立つと思う.

3種テストの方はすべて出力をファイルに保存できるようになっているが,剰余の方は残っている.こちらも保存できるようにした方がよい.ただし,こちらはあまり長いダンプは出ないので,画面からコピーするだけでもよいかもしれない.ファイルを開いて追加出力できればよいのだが…惜しむらくは,文字が少し小さ過ぎるという点ではないだろうか?いまさらフォントサイズを変えるという訳にもゆかないが,もう一回り大きければだいぶ楽だったのではないだろうか?⇒完成間近なのでアイコンを見つけなくてはならない.⇒いまいちパッとしないが,これ「 image 」 でいくことになった.

昨日から走らせている検定がまだ終わらない

昨日から走らせている1281のInvertがまだ終わらない.ψ関数を書き換えて,対数を使って近似値を求め,厳密値を再計算して一致することを確認するという論理にして見たのだが,やけくそ時間が掛かるようになってしまった.この方式では b^x ≡ 1 mod n を与える最小の x を求めるのに,b^x=ny+1 の両辺の対数を取って,x=log(ny+1)/log(b) {y=1→b^n}を計算しているが,まだ,x=10のところを延々とやっている.これで見ると,既存方式の効率の良さが分かる.まぁ,これは見込みないので停止して次に進むことにしよう.⇒やはり,打ち切りボタンは必要だ.⇒元の論理を復活させたら,瞬時に処理完了した.

現行のpsi関数は対象数値が大きくなると時間を消費するものになる.主要な原因は剰余テーブルが大きくなってスキャンに時間が掛かるようになるためだ.Rの大きさはまちまちで順序不定のため頭から全数走査している.もう少し効率的な管理法はないものだろうか?一番簡単なのはハッシュを使うことではないか?ヒットしない場合には多少余分のコストは掛かるが,何もしないよりはましであるような気がする.いや,Rは重複のないユニークな数であるはずだし,R<Nなのだから,単純にT[R]=Rとしておけばよい.これなら相当早くなるはずだ.

いや,それだけでは不十分だ.ψ=i – j で計算されているが,このやり方ではjを決定できない.どこかにこの値を保持しておく必要がある.もう一つ配列を作るか,構造体を作って一緒に格納するか?⇒もうひとつ配列を作ることにした.⇒すでに書き込まれたRというのが出てきた.これはかなりおかしい.対象nは1281で出力ではpsi=60となっている.⇒このテーブルは生成時に初期化されていないのではないだろうか?いや,初期化しても同じだ.

計算は合っている.周期列の桁数は60でψは60になっている.同じRが現れているのにゲームオーバーにならないのはなぜだろう?これらは,冒頭のi=1, 2, 3 で書き込まれた分だ.FORループでは,x=x*bを実行して,bのべきを生成しているが,x≧nになるまでは除算は実行されない.XとRを交換した場合も同じだ.つまり,除算はつねにx≧nに達した時点で実行されるようになっている.多分論理的にはこれで正しいのだろう.問題は重複が発見されたときのインデックスを書き込むかどうか?という点だ.ψは最小である必要があり,そのためにはiから減ずるjは大きい方がよいので,多分書き込むというのでよいのだろう.

「すでに書き込まれたR」が発生するのは,つねにループの冒頭部分だ.これはこれでよいのではないかと思う.⇒修正で相当な効果があるかと思ったのだが,実質ほとんど変わらなかった.12813をテストして,どちらの方法でも約30秒くらい掛かる.テーブルのサイズは2134くらいなので,CPUパワーからすると問題にならないのかもしれない.もう一桁増やして128131を試してみたが,有意差は出なかった.どちらも1分5秒位掛かっている.numが128131に対し,i=3468なので,相対的にはやはりそれほどのものではないということだろうか?まぁ,処理の主要部分は除算なので影響は出ないのかも知れない.処理速度に大差がないとしたら,どちらの論理が優位というのも難しい.まぁ,走査しなくてよいという点では改訂版の方に多少分があるような気もするが…

どうも状態がかなりおかしくなっている.BigInteger.ToStringでおかしな動作が起きている.ゼロが書き込めなくなっていたが,それだけでなく,何も書けない(つねに2が出力される)ような状態になっている.応急措置として,VS2019の修復インストールを試みているところだが,訳が分からない.designerを直接いじるのはまずいということは知っているが,テキストボックスを削除したとき,名前が残っていて誤動作が起きるので削除操作している.しかし,それだけではこれほど具合が悪くならないと思うのだが… 単純ミスだった.すぐ後に置いたGCMの計算でテキストボックスに書き込みしていた.

入力ボックスでnの値を変えると,自動的にφ,ψの値が更新されるようにしようと思ったのだが,少し無理がある.nの桁数が大きくなるとφないしψの計算に多大な時間が掛かるようになり,打ち切りボタンが必要になってくるので,やはり,スタートボタンは必要だ.せめて,φとψの並行計算というのは実施できないだろうか?FactoringSubでは素数判定にφ関数を使っているので,数が大きくなると動作が止まってしまう.

発掘した素数検定ツールも動き始めている

かなり仕事がはかどった.発掘した素数検定ツールも動き始めている.できることが大幅に広がった気がする.ひょっとすると「効率的な素因数分解」まで進めるかもしれない.いまやっている「べき乗の剰余数列の周期性」と「1/nの小数部の循環周期」に繋がりのあることも見えてきた.よく分からないのはPrimeTestとSeedTestだ.おそらくPrimeTestでは素数判定アルゴリズムの検証を試みているのではないかと思われるが,SeedTestのところはいまのところまったく不明だ.出力も雑多な記号の羅列のようになっていて,作成意図がまったく分からない.

先に進む前にいくつか片付けておきたい.まず,ファイルをもうひとつ作って,そこに不要になったコードを移しておこう.廃棄してしまえばそれまでだが,後日何か参照する必要が出てくる可能性もあるので温存しておいた方が無難だ.ファイル名はごみ箱としておけばよい.おそらく,現行版ではまだSieveを使っているはずだが,廃止することにしたので,PrimeTableも破棄ということになる.グローバル変数もごみ箱に移して,何かのときには復活できるようにしておいた方がよい.

PrimeTableはPFactorで使われている.PFactorは最小の約数を返す関数だ.確かに,素数はかなり疎らに存在するから全数検査というのは割が悪い.TotientFuncの論理を使って書き直してみよう.①最初に2をチェックしてその後は奇数だけを検査,②√xの範囲を検査としたのでそこそこ動作するだろう.PrimeTableはPrimeTestSubでも使われている.PrimeTableは素数テーブルなので,これがないと素数判定をランタイムで実行しなくてはならない.⇒代替としてTotient関数を使うしかない.φ値がx-1ならxは素数とみなされる.これしかないのではないだろうか?あまり割がよいとは言えないがやむを得ない.

PrimeTestSubではすべての偶数をチェック範囲に入れている.これはムダなことだ.修正が必要.⇒やってみよう.⇒修正完了した.20614132のInverseを実行したらループから抜けてこない.確かにかなり大きい数ではあるが…この数は,4.85104102370160-10^8だが.周期<20614132であるとしても相当大きなものになりそうだ.このPCでは計算できないかもしれない.⇒おそらく,これはダンプに途方もない時間が掛かっているのだろう.スクロールが止まっているので様子は見えないが,数万あるいはそれ以上の桁数のダンプを出すのは重過ぎる.⇒確かにそのようだ.ダンプを停めたら軽く動いた.

20614132程度の数なら一瞬で因数分解は終わる.これでSieveを使わずに素数判定できるようになった.現在,出力用のテキストボックスが3つあるが,出力はすべてOutボックスに出すようにしてもよいのではないか?その方が画面も簡素になる.あるいは,値を設定すると自動的にφとψを計算して表示するようにしてもよい.その方がおもしろいのではないか?1/nの値を出すのもおもしろい.NotationはOutに出力するのではなく,トグルでBベースと10進を切り替えできるようにするのがよいのではないか?そうすれば,repunitというボタンは廃止できるし,ketaも廃止できるだろう.valueの値を書き換えるという仕様はよくない.

このためにはB進数の文字列⇔数値変換ができなくてはならない.VBにConvert.ToStringという関数がある.ただし,ToStringはInt32だ.この逆変換としてConvert.ToInt32というのがある.入力ボックスに入力される値としてはInt32以上を禁止しても実用上は差し支えないのではないか?φとψを常時画面に表示するようになれば,PrimeTestは不要になる.PrimeTestの目的は φ mod ψ ≡ 0 の検証であると考えられるので,例外が発生した場合にだけ停止するようにしておけばよい.

Inverseの出力ももう少しシンプルなものにしたい.Inverseでは対象数の小数点表示を出すべきだ.周期列とφないしψの関係がくっきり出るような表示が必要だ.一番分かりづらいのがSeedTestだが,狙っているのは冪剰余列でやっているようなところではないのだろうか?べき剰余列でも,φなどを表示するべきだ.冪剰余の連続実行ボタンで,InverseやSeedTestを実行できるようにすることも考えられる.チェックがオンになっている機能を実行するようにすればよい.

Convert.ToStringはかなり動作に制限がある.引数が8ビット,つまり,1バイトの数を指定した進数に変換するというだけだ.つまり,数を一度BCDに変換して格納しておく必要がある.逆数を常時表示しておきたいのだが,表示桁数が足りない.doubleで15桁,decimalでも29桁しかない.InverseFuncで出力を停めているが,ループから抜けてこない.配列サイズオーバーが発生してしまう.とりあえず,対象数値をInt32に限定しないと無理なようだ.Int32に限定してもまだオーバーフローしてしまう.7FEFFFFF=2,146,435,071が限界ということのようだ.⇒1/nを10^m倍して固定小数点化したものを文字列に変換して出力するようにした.これなら桁数がいくら増えても表示できる.

一つおもしろいことに気付いた.n=7として,1/7を計算すると,{142857}という循環節が現れる.これを数とみなしてその逆数1/142857を取ると,この分数の循環単位は{000007}となる.つまり,{142857}と{000007}はある種の逆数のような関係になっている.ψ(n, b)=循環単位桁数であることはほぼ間違いないいが,これらを計算しているコードを見ると,ほぼ事実上同じコードになってしまっている.これでは同じ値になるのは当然だ.両者はそれぞれ独立の概念として定義され,それに従ってコーディングしているつもりだが,実質的には同一のものに別の名前を与えていただけなのだろうか?

ψは,ある基数bのもとで,b^x=1 mod n となるような最小の数として定義されている.循環単位桁数はその言葉通りの意味だ.ψの計算コードが間違っているのだろうか?どうもその気配が濃厚だ.定義によれば,B^x=1 mod n となる数 x を探しているはずなのに,循環が発生する時点での値を返している.どうも,まだ読み切れていないのだろうか?PrimeTestは単点テストと思っていたが,Bの範囲を指定してリスト出力する多点検定だった.⇒処理を中断する打ち切りボタンが必要だ.

356をInvertしてfixedの値が間違っている.前方の0がカウントされていない.⇒対処した.

べき乗の剰余数列の周期性について

2005年に書いた An Efficient Factoring Algorithm by Repunit Number Method に関係するソースコードを発掘したので,手を入れて動かしてみたい.これはVBで書かれていたもので,VB6を再インストーしてみたが,動作しない.VS2017である程度ファイルを読み込めたので,使えそうな部分だけをコピーして蘇生させてみた.とりあえずなんとか動くようになった.これを仕上げておけば,手元に置くツールとして何かと有用なのではないかと思う.剰余演算や因数分解,φ関数,GCD,単位反復数などが出てくる.使い方の分からないところもあり,よく動いていないところもある.

  1. Modulo 動作しない
  2. PrimeTest 動作しない
  3. Repunit 動作しない
  4. Division 除算を実行して,商と剰余を出力
  5. Totient φ(n)を計算しているはずだが,間違っている
  6. SeedTest エラーが発生する
  7. Factoring 因数分解のはずだが,何も表示されない
  8. Inverse エラーになる
  9. GCM 2数の最大公約数を表示

まず,動作しない機能を動かすところから始めよう.

Moduloは,num0, base, Factorの入力値を取り出して1~Factorまでの剰余数列を画面出力している.これは結局べき剰余でやっていることと同じなので統合してしまうことにしよう.⇒並列動作できるようにしてみたが,両者がやっていることはかなり違う.Moduloの場合は,y*b^p mod k でbは固定,pを1~Factorまで変化させている.つまり n^p mod k で見ると,p^0, p^1, … p^F mod k を出力している.

Modulo では,yが0になると抜けるようになっているので,k=10でb=10の場合は,0を出力しただけで終了してしまう.⇒ようやく,話が合い始めた.結局のところやっていることは同じだ.つまり,ModuloボタンとGoボタンは当機能と言ってよい.Moduloは捨ててよいだろう.

image

この時期はまだ,周期性というのに気付いていないので,ともかくFで指定した回数だけべきを繰り返すという作りになっているが,現行版ではkで指定した範囲のべきを生成して,その中の一周期分を取り出している.やはり,こちらの方が分かり易いと思う.新版では,nとkの範囲を指定して,n x k のテストを一括実行できるようになっている.一応ここまでできたので,バックアップを取っておこう.⇒失敗した.VSのフォルダで作業していた.開発用ドライブに移動する.

PrimeTest というのがあるが,少し大きそうなので,先にRepunitを見てみよう.これはどういう計算を行っているのだろう?たとえば,22 # 5 →  245411のような数が出力される.入力された b, psi という2つの値から,Σ b^i {i=0→psi-1}という計算をやっている.この値はnum0に入るが,剰余計算では除数(法)に当たる.Decimalというボタンがあって,このボタンを押すと,指定した基数Mで対象数のM進数表記を出力するようになっている.先にこちらを片付けよう.

この機能はアクションとするより,むしろ状態として扱った方がよいのではないか…つまり,たとえば,16進表示を指定すると,画面のすべての数字が16進になるというようなイメージだ.まぁ,そこまでやる必要があるのかないのかもわからないので,とりあえずは,動くようにするというところまでとしておこう.ここで,とりあえず,longを一括bignumberに変換しておこう.⇒BigNumberではなく,BigIntegerだった.⇒定数はBigIntegerにできない.整数除算の\もBigIntegerでは使えない.何を使えばよいのだろう?整数除算はBigInteger.Devide, 乗算はMultiply,BigInteger.Parseなどという関数もある.

quotient = BigInteger.DivRem(dividend, divisor, out remainder)

で,商と余りが同時に取り出せるようだ.^はPowを使わなくてはならない.PowerModという関数が入っているが,これはBigIntegerのModPowと同じものなのではないだろうか?いや,どうも様子が違う.奇数の場合には特殊な操作をしたり,平方してからmodを取ったり… この関数を呼び出しているのは,ModuloMultiとDivideRepunitだ.どりあえず,現行のままBigInteger対応修正だけ入れておこう.BigIntegerでは論理演算はできない.これらに関わる変数は暫定的にLongとしておいた.BigItegerではStrという関数も使えないようだ.⇒Decimalは通常10進と解釈されるので,Notationに改めた.Notationでは16進の場合などアルファベットを使うようにしたいのだが,先に進もう.

Divisionという機能は除算を実行して,商と余りを計算するだけなので,とりあえず現状のままとする.Totientというボタンではφ関数を計算しているが実行時エラーが起きる.Int(n/p)で割り切れているか否かを見ているためだ.この関数は丸ごと差し替えることにする.廃棄する前に比較して動作チェックしておこう.既存関数はLongで動かすようにした.⇒かなり悪い.どうも,PFactorという関数が間違っているようだ.⇒PrimeTableが初期化されていないのではないだろうか?これをやっているのは,Sieveという関数だ.ファイルをロードしたときに実行されている.このテーブルは再定義されるべきではない.Form_Loadイベントが掛かってこないため,sieveが実行されていない.⇒修正した.

MaxPrimeは10^9と定義されているが,実際はその平方根の範囲のテーブルしか作っていない.BigNumberはどこまで大きくなるか分からないので,あらかじめテーブルを作っておくというのは実際的ではない.基本的にPrimeTableを使っている関数は廃止すべきだろう.⇒とりあえず,いましばらく論理だけは置いておこう.⇒Factoringは取り敢えず問題なく動作しているようだ.PrimeTestSubを呼び出して,結果をパネルに出力する.素数の場合には何も出力せず,Factorにそのままの値を書き込む.この処理は取り敢えず不要な気がする.最終的な仕様としては,値を書き込んだ時点で自動実行するようにすべきだろう.

PrimeTestSubという関数は4箇所から呼び出されている.①FactoringClick,②PrimeTestClick,③SeedTestClickx2.Factoringで最後の値を出力しないのはやはり,まずいと思う.たとえば,246をFactoringすると,2,3が出力され,41がFactorに入るが混乱する.やはり,すべての素因数を列挙すべきだろう.となると,1も入れなくてはならない.⇒末尾に「,」が出ないようにする.⇒PrimeTestが動くようになった.verboseというグローバル変数で素因数分解の出力を抑制し,それ以外の付加情報を出力する.

GCMはGCMfuncという自前の関数を使っている.これはaとbのどちらか小さい方yが大きい方xを割り,その余りRを取って,y→x, R→yのように転換してRがゼロになるまでループするというロジックだ.おそらく,これは最速と思われるが,BigIntegerには,GreatestCommonDivisorという関数があるので,今後はこれを使うことにする.SeedTestが動作しないのは,リスト要素が明示的に追加されていないためと思われる.VB6ではこれで動いていたはずなのだが…⇒List1はSeedTestClickで参照されているだけで,内容はパネルに出力されるようになっているから,実質的には不要と思われるので,廃止することにする.

SeedTestClickは内部でPrimeTestSubを呼び出して,その結果を使っている.しかし,画面が乱れるので,PrimeTestSubでは出力させたくない.⇒大体動くようになったのではないかと思う.あとは,数値をチェックして動作を確認するだけだ.一度バックアップを取っておこう

もう一つ機能があった.Inverseというやつだ.これが分からない.⇒それでも一応動いているようだ.n=14のとき,

0, 10
8, 4
3, 4
n=12, K=0, ψ=1, R=0, φ=4, φφ=2

のようなダンプが出る.inverseというのは1/nのことを指している.そして,ここで問題にしているのは,有理数1/nの小数部における「周期列の長さ」である.周期列の長さとは,小数位の桁数のことであるから,その数字列の基数,つまり,x進数という表記形式が決まらなければ決定できない.ψ(b, n) のbはこの意味の基数である.ψ(b, n) はつねにn より小さい.the length y(b,n) of a recursion unit of 1/n is always smaller than n. fixed というパラメータが出てくるが,これは周期性のない固定部の長さのことだろう.この意味では Inverse では,まず,その進数表記に基づいたその数を表示すべきだろう.

Inverse で出てくる商と剰余のリストは何を示しているのだろう?確かにこれは1/nの計算ステップになっている.たとえば,128 base 10の場合は,0.0078125となるが,固定部の長さは7で,Qは{0, 0, 7, 8, 1, 2 5}だ.Rは{10, 100, 104, 16, 32, 64}となっているのだが… この関数の中では対数計算を行っている.pow = log(R) / log(b)という値を計算し,pow=b^powとしている.RTという配列を持っていて,そこにRという値を格納している.xの初期値は1dでそれをnで割る操作を反復しているのだが,x*bを計算して1桁上げてから整数除算を実行しているので,つねにその位の計算が実行されていることになる.

このアルゴリズムでは同じRが出現したときに周期列が完成したと判定し,そこで処理を打ち切ることになる.これはまったく正しいと思う.こういう式が出て来る.

n*U(b,n)=r*(by(b,n) -1)

U(b, n)は recursion unit と呼ばれる数で,rは最初の被除数かつ,最後に現れる剰余とされる.n=21, b=10 の場合,剰余数列は,{1, 10, 16, 13, 4, 19, 1}のようになり,ψ(10, 21)=6,U(10,21)=47619,r=1である.Inverseは大体動くようになったが,まだ動かしていなかったものがある.repunitだ.いや,もしかしたらやってたかな?設定値が欠けていると計算できなくなるようだ.⇒対処した.今度はFactoringが動かなくなっている.あれ,おかしい.動いた.どこを押してたのだろう?いや,やはり何か動作不良なところがある.SeedTestを実行するとおかしくなるような気がする.どこかでout.textを消去する余分な動作が入っているのではないだろうか?ともあれ,かなり整ってきた.

image

日によって来訪者カウントにムラがある

日によって来訪者カウントにかなりのムラがある.今日はまだ午前10時なのに,465もある.現在の延カウントは106662だ.ともかく,楕円コンパスの包絡線の宿題を片付けなくてはならない.「直線と楕円包絡線の2」というwxmがあるので,ここから始めよう.これはすでにマルチプロットになっているので,あとはGIFアニメに落とすだけだ.

image

と言っても,この作業がまだルーチン化していない.いまのところdrawwdで直接描画しているので,これを terminal::animated_gifに落とさなくてはならない.

数学物理etc談話室のディスカッション

数学物理etc談話室のディスカッションが活発になってきているので,これに時間を取られて他のことがほとんどできない状態になっている.テーマがどんどん広がってゆくので手が付けられない.いまも,①コンウェイの定理のところで,コンウェイの六角形でL=a+b+cになる条件を求めるという問題,②加納正司氏の提出した4桁数字根の問題,③界斜等円図から派生した楕円包絡線とコンパス線の関係など,課題が山積みになっている.④QNO.77(2023/03/18)の回答にも一部疑問が残っているので,どこかでクリアしておかなくてはならない.

④はWolfram alphaでも解けないので,かなり難しい問題だ.T氏の回答はおそらくどこかからコピペしてきたものと思われるので,氏に質問しても答えられないのではないかと思う(ので黙っている).置換積分ないし部分積分をやっているはずだが,どの部分を置換しているのかが突き止められない.maximaでも解けなかった.

maxima専用フォルダにファイルを整理

どうも作業の進捗がはかばかしくない.やることは決まっていて,それほど難しいことをしている訳ではないのだが,アニメを作るところで,というかそれ以前のところでもたついている.一つにはファイル管理,つまりライブラリアンの問題があるような気がする.Your Daily Epsilon of Math Calendar のフォルダで作業しているのだが,各種のファイルが溢れかえってカオス状態になっている.ファイルの命名規則lが定まっていないという問題もある.これまでのところは仕方ないとしても,何か整理された分類法を考える必要がある.

  1. アニメGIFには必ずタイトルを付けること
  2. GIFタイトルとファイル名は同名とすること
  3. GIFアニメはPhotoにアップロードとしてアルバムとして管理
  4. wxmファイルは常時バックアップが必要だが,そろそろ日付管理が必要
  5. 日付管理するとなれば,ログを残さなければ意味をなさない
  6. maximaで作業する専用のディレクトリを作った方がよいのでは?
  7. 安定版には追番を振らないファイル名を付ける

maxima専用フォルダを作るというのはよいと思うが,どうやって管理するかが問題だ.できればテーマごとに分離した方がよいのだが,必ずしも確立したテーマが存在する訳ではない.しかし,ファイル名を決めるためにはまず,テーマを確定することが必要だ.その場限りのまちまちの名前を付けておくと,後日検索するのに手間取ることになる.ともかく,いま,YDEMCフォルダに入っているwxmとgifファイルを別フォルダに移してみよう.いま開いているwxmは3つある.閉じる前に名前をメモしておかないと忘れてしまう.

  1. 包絡線図3 「Ellipse-Lines」 楕円→星芒線を出すアニメ 3x3
  2. 直線と楕円包絡線2 「Line and Ellipse Envelope」 楕円,直線→星芒線 2x2
  3. 楕円・直線コンパス アルキメデスの楕円コンパスのアニメ 3x3 最後にエラーが出る Improper argument: gr2d 描画は完了している

もし,wxmのコーディングが「開発」の一端であるとすれば,作業環境もFireBirdに移すべきではないか?どっちみち,ゼルコバの木と同時平行作業はできないのだから,FireBirdで作業した方が効率がよい.そうすれば,Vultureの画面も外部とのコミュニケーション用に空けておくことができる.GIFを投稿するときには,FireBird→Vultureにコピー→投稿でよい.Vultureは256GBのSDDを持っているので,そこをベースとしておけば,バックアップの手間も省けるし,Cドライブ容量の節約にもなる.⇒とりあえず,SDDにmaximaというフォルダを設置した.

EPSILON 23年3月 からGIFファイル20本,wxmファイル54本をVultureのmaximaに移動した.しかし,maximaのような不安定なコードを開発環境で実行することにも不安はある.maximaの暴走は場合によっては外部アプリないしシステムにまで影響する可能性を排除できない.BlackBirdはいまのところ手ぶらだが,ここで作業するというのもかなり無理がある.いいアイディアだと思ったのだが… maximaを開発機で実行するのはやはり,問題が大き過ぎると考えざるを得ない.ネットアクセスの部分をBlackHawkに移管できないだろうか?できないことはないような気がするのだが… 現在一番安定が悪いのが,このマシンだ.

BlackBirdは久々に立ち上げたので,カスペルスキーの定義データベースの更新にかなり時間が掛かりそうだ.一応まだ動いていないので,外部アクセス用としてはしばらく使えるだろう.もし,うまくゆけば,毎朝のルーチンになっているYour Daily Epsilon of Math Calendarの投稿もこのマシンで実行されることになり,メリハリがしっかりしたものになる.⇒BlackHawkをSNSアクセス用に使うというのはよいが,遅いのが難点だ.Cドライブの空きは2GBあるから,まだなんとかなる範囲なのだが… Vultureではログ出力を確認するだけになるので,モニターはwxmだけで占有できる.上記リストの項目2は動いているので,アニメ化するだけだ.またおかしくなっている.何もエラーは出していないのに,画面がバタバタしている.

楕円包絡線の原理が少し分かりかけてきた.楕円族の中で星芒線に接しているのは一部の楕円だけであり,それ以外のものは交差しているか,ないし完全にその外部にある.その一方直線族はすべての直線が星芒線に接している.コンパス線はペン先の移動に応じて異なる径の楕円を描くので,複数の楕円が1本のコンパス線を共有しているが,ある径以上の楕円は星芒線の描画に貢献していない.Line and Ellipse Envelopeを見ると,径が1より大きい楕円では扁平率が変化していない.これはなぜだろう?少なくとも複数の扁平率の楕円が混合しないと星芒線は現れない.最近作ったアニメの概要を拾っておこう.

  1. 楕円とコンパス線の包絡図 楕円とコンパス線の包絡図 短径0.5,長径1.5の楕円を描画するコンパス線のマルチアニメ 終りの方になると描画が極端に遅くなる これは公開している網目タイツ包絡線とほぼ同一
  2. アルキメデスの楕円2 最初のセクションは動かない→修正して動くようにした.Arichimedes2.gifを出力,短径≦1の範囲で扁平率を変えて楕円を描画.もう一つのセクションもほぼ同じ.こちらの方がアニメが滑らかに動く,この差はkのピッチの違いだけだ.
  3. アルキメデス Arichimedes.gifを出力 上と同じ その他,いくつかバリエーションがある
  4. 楕円コンパス(GIF) 半径0.5の円とコンパス軸部のアニメ 2番目のセクションはペン先まで含めたコンパスで楕円を描くアニメ 2x2 GIFは生成していない
  5. 楕円コンパス3 星芒線に含まれる楕円の静止画を出力
  6. 楕円コンパス4 楕円コンパス(GIF)と同じ,2番目のセクションも同じだが,スピードが速い
  7. 楕円コンパス6(GUI) アルキメデス楕円コンパスの安定版 makelistの2重ループでArchimedesCompass.gifを出力
  8. 楕円コンパス包絡図 楕円コンパス(GIF)と同じ 2番目も同じ
  9. 楕円コンパス包絡図2 「Line Segment Envelope」 半径0.5の円とコンパス軸部のマルチアニメ これで見ると,包絡線に寄与しているのはコンパス軸部だけであることが分かる SⅡは壊れている
  10. 楕円とコンパス線 SⅠは動作しない SⅡはEllipse-Compass Lines.gifを出力,このアニメはアルキメデスの楕円2と同じ 
  11. 楕円と楕円コンパス線の接点 「楕円とコンパス線の包絡図」短径0.5,長径1.5の楕円とそのコンパス線のマルチアニメ.楕円とコンパス線の包絡図と同じ pitch 64で最後まで線が飛ばない appendの中でmakelist
  12. 直線と楕円包絡線 円の中にゆっくり楕円を描くアニメ アルキメデスの楕円2と同じ
  13. 直線と楕円包絡線 「Line Segment Envelope」 直線と楕円包絡線2と同じ 
  14. 直線と楕円包絡線 「Line and Ellipse Envelope」 直線と楕円包絡線2と同じだが,中央に半径1の青い円を出している 終盤の描画は遅い
  15. 包絡線図 短径2,長径3未満の青い楕円と直径1と0.5の円,青のコンパス線のアニメ 動点が2つある
  16. 包絡線図2 「Ellipse-Lines」包絡線図と同じ 短径は2,長径は3未満 3x3 動きが遅いので時間が掛かる
  17. 包絡線図3 「Ellipse-Lines」マルチプロットアニメ 楕円のマルチアニメ 3x3 最後にエラーが出る
  18. 包絡線図4 「Ellipse-Lines」multiplotとmakegifの切り替えができるが仕上がっていない 何も表示しないでエラーで停止する

楕円と楕円コンパス線の接点

さて,始めよう.アルキメデスのコンパスが描く楕円の包絡線とコンパスの軸部の直線が描く包絡線が同じものであることを示すというのが課題だ.状況を具体的なアニメで示すことが一番分かり易いので,その方向に向かっている.⇒エラーが出ている.

WARNING:Couldn’t write to #<SB-SYS<FD-STREAM for “descriptor 744” {1003B33583}>:パイプを閉じています。
Trying new stream.WARNING:Couldn’t write to #<SB-SYS<FD-STREAM for “descriptor 724” {1004D235B3}>:パイプを閉じています。Trying new stream.

既存コードを流用しているだけなので,問題なく動くはずなのだが… else文を停めたら動き始めた.どうも,else文が実行されてしまっているようだ.構文的には間違っていないように見えるのだが…

楕円とコンパス線

この図で見ると楕円と楕円コンパス線はまったく接していないように燃える.楕円自体が星芒線の領域からはみ出している.楕円の描画が間違っているのだろうか?この楕円はparametric(d*cos(θ), (1-d)*sin(θ), θ, 0, 2*%pi)という式で描画している.長軸の半径を与えるdは定数で1.5だ.コンパスのペン先の同点を描画してみればすぐに分かる.「楕円コンパス5.wxm」がGIFなしの楕円コンパスだ.

ここでは,kを-1~2の間でコンパスの状態を描画している.kというのは,ペン先の位置を示すものと考えられる.ブルーの点がペン先で楕円はparametric(k*cos(θ), (1-k)*sin(θ), θ, 0, 2*%pi)で描画している.このkは上記のdに相当するから,一致している.どうもわからなくなってきた.下のアニメで見ると楕円とコンパス線は完全に交差していて,接するようなところは見られない.

Envelope5

これはどういうことか?おそらく,楕円の包絡線と直線の包絡線は形状的には類似しているが,スケールがまったく異なるのだろう.楕円の包絡線のアニメはまだ作っていなかったのではないだろうか?それをやって見る必要がある.楕円を描いているアニメは現在のところ2つある.一つはArchimedesCompass.gifで,もう一つはArchimes.gifだ.2つとも昨日のログに掲載している.前者は長径が0.5から2までをカバーしているが,後者は長径が0.5~1の場合しかカバーしていない.包絡線を見るにはある程度の本数を引かないと見えないので,前者からスタートするしかないのではないか?コンパスの軸を直線と交換し,重ね書きできるようにすればある程度見えるようになるのでなないだろうか?

やってみよう.