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

昨日から走らせている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の場合しかカバーしていない.包絡線を見るにはある程度の本数を引かないと見えないので,前者からスタートするしかないのではないか?コンパスの軸を直線と交換し,重ね書きできるようにすればある程度見えるようになるのでなないだろうか?

やってみよう.

網目タイツ包絡線図

FBグループの「数学物理etc(談話室)」で和算に関わる議論が続いている.その一環として「宿題Ⅲ楕円コンパス直線の1次方程式をθで偏微分することでアステロイド曲線の式が得られることを示せ.もし,そのような式を得ることができないとすれば,その理由を述べよ.」と言う課題を提示した.自己回答しなくてはならないが,少し下調べすることがある.①楕円と楕円コンパス線が接している,つまり,楕円コンパス線が楕円の接線になっていることを確認する.②その接点が包絡線上の点であること確認するという2点だ.

FishnetTightsEnvelope.wxmという網目タイツ包絡線図を出しているmaximaのソースコードがあるので,ここから始めることにする.楕円を出しているコードも必要だ.どうもmaximaのファイルを後から探すのは結構大変な気がする.一々開いて走らせなくてはならないし,走らせても動かないものもあったり…結局,FBの投稿に戻ってそこで探すようになったりする.こちらの方が図が見えているので探し易い.いま探しているのは「宿題Ⅱ:アルキメデスの楕円コンパスの拡張」というタイトルで掲載しているものだ.

ArchimedesCompass

Envelope5.gifというファイルが残っていないということは,作れなかったということだろうか?おそらく,もう少し精度を落として別のファイルで保存したのだろう.pitch=64で成功したことがあるような気がするのだが…pitch=32では保存できた.この位で満足するしかないのではないだろうか?楕円コンパスの包絡線に関する研究の端緒となった画像(楕円コンパス.png)を掲載しておこう.このファイルは拡張子がPNGになっているが,どういう経緯でそういうことになったのだろう?これを出力したwxmは残っているだろうか?これは,多分multiplot_modeをオンにして描画していると思う.

楕円コンパス

これと似た図にLineEnvelopeというのがある.下図はいま,LineSegment.wxmの出力を保存したものだ.上の曲線群は楕円で,下の曲線群は直線となっているところが異なる.これらの2つのベースはいずれも楕円コンパスだが,楕円曲線の包絡線と直線群の包絡線が同じ形状(多分一致しているはずだ)をしているという根拠を突き止めたい.

LineEnvelope

この図の線分を延長したものが網目タイツのラジアルな直線なので,この2つで比較してもよいのだが,最終的には動かしてみたいので,元になるアニメが出力できるwmxを探しているところだ.ただし,楕円曲線に関しては多分フルスペックのアニメは作っていないと思う.このwxmでもアニメを作っている形跡はあるが,出力まで行ったかどうかはわからない.楕円を出しているアニメは多分これ↓だけだろう.

Archimedes

このファイル(Archimedes.gif)は多分,アルキメデスの楕円.wxmかその前後の版で作っているはずだ.

玉露の粉茶がまだ残っていた!

お酒が切れてしまったので,業スーに買い出しにゆかなくてはならないのだが,すでに3時.欠品を点検してみたところ,まだしばらくは持ちそうなので,Wellciaでお酒とタバコだけ補充するという方針に切り替えた.4L25度の焼酎がWellciaでは業スーより100~200円くらい高くなるが,全体の出費は押さえられるのでこの方がベターだ.Wellciaではお茶を補給できないのだが難点だが(粉茶はないし普通のお茶は高い),この前(敬老の日)に市から支給された地域通貨(ネギー)で市内で購入した玉露の粉茶というのがまだほとんど手付かずで残っていた.

maximaのコードを書くのも一種のプログラミングなので,やはり,ログを付けながらでないと捗らない.二重ループを構成する必要があるが,まず,内側のシャフトがスライドしながら回転するところを作ってみよう.wikiにあるアルキメデスの楕円コンパスの模倣だ.⇒maximaの具合がかなり悪い.保存された「アルキメデスの楕円.wxm」を開いたところ,ほとんど壊れていて白紙状態になっている.「楕円コンパス.wxmx」は問題なく開けた.ここから始めることにしよう.

このサンプルでは,

parametric(k*cos(θ), (1-k)*sin(θ), θ, 0, 2*%pi)]

で楕円を描いている.とりあえず,この座標をX軸とY軸に射影して2点を結べば楕円コンパルのロッド部分になるはずだ.このサンプルはdraw2dで描画しているが,引数リストの中でFORループを使えるだろうか?というか,多分drawに書き換えないと動かないと思う.アニメのコードから始めた方がよさそうだ.⇒drawではgr2dで描画している.makelistはそれ自体ループになっているが,その中でFORループを使うことはできるのではないか?⇒makelistのループはkというパラメータで楕円の半径を制御している.

このkという値は,シャフト上の点の位置を示すものだが,それを延長したものが本来のアルキメデスコンパスだ.シャフト両端の座標は半径1の円周上の点をXY軸に移したものだから,kの値に関わりなく描画できなくてはならない.つまり,makelistの内側にループを作る必要がある.とりあえず,makelistをコメントアウトして,画面出力できるように改造しておこう.⇒どうも,maximaは思ったより難物だ.カンマのある/なしなどで動作しないのはよいとしても,安定がすこぶる悪い.「refusing」というぶっきら棒なメッセージが出たあと,メインの編集画面がほとんど壊れてしまっている.

maxima で遊ぶ

  1. アニメGIFに保存
    file_name = “zzz”,
    terminal  = ‘animated_gif,
  2. 縦横比を一致させる
    *proportional_axes=xy
  3. 表示範囲
    xrange = [-0.1, 1.1], yrange = [-0.1, 0.5],
  4. グリッド。y 軸の目盛
    grid = true, ytics = 0.2,
    xaxis = true, yaxis = true,
  5. 線幅,色
    line_width = 1,
    color = purple,
  6. terminal = screen (デフォルト), png, pngcairo, jpg, gif, eps, eps_color, epslatex, epslatex_standalone, svg, canvas, dumb, dumb_file, pdf, pdfcairo, wxt, animated_gif, multipage_pdfcairo, multipage_pdf, multipage_eps, multipage_eps_color, aquaterm
  7. dimensions=[600, 500] 出力端末サイズ 1/100cm単位の長さ
  8. scenes: drawで連続出力に用いるシーン配列オブジェクト


zzz

oscillator_phasespace_potential

Archmedes

とりあえず,ここまでできた.

Arichimedes

しばらくテント村にご無沙汰しているが,今日の来訪者カウントは23時現在で,776を記録している.こうなると,近い将来にも来訪者/日が1000を超える可能性が出てきた.

image

カンブリア爆発?もしかするとわたしの誕生日に関係していたのかも?FBでは80個を超えるお祝いをもらっている.それにしても,更新が止まっているテント村のカウントが上がっているというのはナゾだ.