MathematicaのCal

いろいろなプログラム言語でCalのプログラムを書いてみた. 今回はMathematicaである.

Mathematicaというと, 一発処理のような印象である. 例えば Prime[25] と入力すると, 25番目の素数 97が出力される. しかしなんとかプログラムを書いてみようと思い, 例によってCalを書いてみた. 結局いろいろなノウハウを覚えつつ試行錯誤の結果完成したのが次のプログラムである. その見どころを説明しよう.



0行目. 仮引数はyr_のように後に下線をつける. 本体で使うときはその下線はいらない. システムの関数手続きは大文字で始めるが, ユーザの関数は小文字で始める. 本体の定義は:=の後に書く.

1行目. Module[{ の次にローカル変数を書く. ローカルな関数名も書く.

2,3行目. printd[d_]は, その月の最後の日付をzとして, 1≤d≤zの時は出力し, それ以外は3文字の空白にする. さらにdが1桁なら前に2文字の空白, 2桁なら1文字の空白をつける.

4,5行目. prline[d_]は, dからd+6までをprintdで文字列にしたものを, StringJoin(<>を使う)し, Printで出力する.

6,7行目. 年と月を出力. また曜日名を出力.

8,9行目. うるう年なら, ローカル変数lを1に, そうでなければ0にする.

10行目. 各月の朔日の曜日を決めるため, 次の年の初めまでの日数を計算.

11行目. 1,2月の日数と, 月初めから年末までの日数(mod 7)を設定.

12,13行目. 3月以降の月に上のような値を計算する.

14行目. 朔日の曜日wを確定する.

15行目. dを第1週の日曜に相当する日(1-w)から月末日まで, 7ずつ増やしてprlineする.

こうして出来たのが次だ.



なれない言語でプログラムを書くのはしんどい.

月齢カレンダー

このブログの2008年12月に月齢計算のことを書いた. それは計算機を使って月齢を計算する自己流アルゴリズムであったが, その昔, 私は高校生の頃, 先輩に教えられて毎年月齢定数というのを使って月齢を概算していた.

それではその年の月齢定数をCとすると, その年のm月d日の近似的な月齢を (C + m + d) mod 30 で得るものであった.

それで時々年の初めころ, 今年の月齢定数はいくつかなとほどほどに計算して記憶していたこともあった.

ところでインターネットであちこち見ていたら堀さんの式というのに出会った. http://www.asj.or.jp/geppou/archive_open/1968/pdf/19680704.pdf

詳しくは元記事を見ていただくことにして, 大体は月齢定数と同様の趣旨であるが, mの代わりに f(m)を足すのである. この値が1月から順に0,2,0,2,2,4,5,6,7,8,9,10なのである. 6月から後はf(m) = m - 2なので, その辺は昔の月齢定数と同様だが, 年初がすこし違っていて, この式の方がより正確になるらしい. 元記事のタイトル「鬼鬼西」は1月から6月までの値の覚え方である.

Cの方は, (((y - 11) mod 19) × 11) mod 30 だそうで, 本年2016年なら 2005 mod 19 = 10なので, 110 mod 30 = 20が今年の月齢定数になる. 元記事の掲載された1968年はこの定数がちょうど0だったのでタイムリーな発表であった.

f(2)=2だから (20+2+8) mod 30 = 0 となり, 2月8日が新月つまり春節であった. 3月の満月をxとすると, 20 + 0 + x = 15 (mod 30) だからx = 25となって, 今年の復活祭は3月27日である. また3月9日には部分日蝕があるが, その月齢は 20 + 0 + 9 = 29 だ.

なおこの定数は 2005 mod 19 = 10 から後8年はこの値が1ずつ増え, 従って定数も毎年11ずつ増えるから初めから計算しなおすことはなく, 来年2017年は1, 2018年は12, 2019年は23,...と分かる. こうして..., 7, 18と来たら, その次は29ではなく0とする(なる?)のだ. 次回0なのは2025年である. 1968年, 1987年, 2006年, 2025年と19年ごとに定数が0の年が回って来る.

菱形六十面体

2010年5月17日のこのブログの続きである. そこにはWolfram MathWorldのRhombic Hexecontahedronの項に 20 golden rhombohedra can be combined to form a solid rhombic hexecontahedron. と書いてあるのが気になっていて, そのrhombohedronの頂点の立体角を計算し, たしかに20個あると全立体角, 4πステラジアンになると書いた.

golden rhombohedra は菱形六面体とでも訳そうか, 6つの面が黄金菱形になっているものである.

Wolfram MathWorldのGolden Rhombohedronの項は, これには鋭角のものと鈍角のものがあるというが, ここで使うのは鋭角の方である. この項には1辺の長さaの鋭角黄金六面体の 体積は1/5√(10-2√5) a³とあり, また菱形六十面体の項にはその体積が4√(2(5+√5))a³とあるから, 体積的にも確かにそうである.

ざっと考えてみると, 菱形六十面体の母体である正十二面体には頂点が20あるから, 多分その頂点の数と菱形六面体が対応しているであろうと推察できるが, あまり自信はない.

昨今, 3Dプリンタが使えるようになったので, ちょっと作ってみようという気になった.



上の写真の上は20個の菱形六面体が出来たところ, 下はその5個を使い, 菱形六十面体の内, 正十二面体の一面に相当する面が完成した所である.

この写真では分らないが, 3Dプリンタで作ったものは, 存外 外面がざらざらで, 接着するのが大変であった. 次回はもっと時間がかかっても, もう少し精密な面にしてみたい.

さて, 次の5個で同様に正十二面体の反対側の面に相当する面も作る. この両方の5個組を北極と南極とすれば, 残りの10個で赤道に相当する部分を構成することになる.

それには10個の内2個ずつを接着して2個組を5組作り, それを一周するように繋げるのである. これは仲々面倒で, どの面とどの面が合さるのか, いろいろこねくり回さなければならない. まぁそうこうする内に固定出来た. 下の写真は接着中の写真で, 小型のクランプやしゃこ万力で押えているところである. 左に転がっている2個は北極と南極の5個組である.



出来上がった赤道部分のこのようである.



次の図はPostScriptで描いた, 菱形六十面体と上でいう北極, 赤道, 南極の図である. 手前右上から光が当っているように影を付けてみたが, 効果はいまいちであった.

全体
北極
赤道
南極

相対的な位置が分るように, 外接する正十二面体と, 中心を通るx,y,z軸も描いてある.

これで見ると20個すべての菱形六面体が中心に接しているのが分る. つまり赤道の組は中央が厚さ0である.

3Dプリンタで作ったモデルにはこのような精度は残念ながらない. でもこんな実験が出来るのも3Dプリンタのお蔭である.

菱形六十面体

菱形六十面体についてこのブログで何度も取り上げた. その最初, 2010年5月16日のに「角D'B'E'は黄金菱形なので, 黄金比φ=1.61に対して2cot^-1(φ)=63.435°である. その一端を持ち上げると, 角DBEは少しずつ広がり, 72°になるわけだ. その時の傾き角を計算すると, 驚いたことに31.7175°であった. 」と書いている.


長い間 気になっていたが, 先頃 証明した. 図を見てほしい.

図の下は破線のAD'B'E'が縦1, 横φの黄金比の菱形である. 実線のADBEは角Aが72度で縦が1の菱形である. この横の対角線の長さを以下eとする. 黄金比の方の角Aの半分は31.7175°である.

図の上の直角三角形は底辺が下の実線の菱形の横の対角線, 斜線が黄金比の対角線で, これを描いてみると角Aがやはり31.7175°になるのである. この角を以下θとする.

証明したいのは, 上の31.7175°, つまりcos θ=e/φが下の31.7175°, つまりtan 2/φに等しいということである. 一連の計算は下の通り.


φの値はよく知られている. tan 36°はcos 36°がφ/2 から計算できる.

(0)上の図のθが欲しいので, atanの式を書く. 以下atanの引数を計算する.

(1) φとeの値を代入する. (2)根号の中, 第1項を展開し, 第2項の分母を有理化する. (3)根号の中を通分し, 1/eの値を書く.

(4)根号の中の分子を計算する.

(5)根号の中, 分母分子を2倍すると, 分母分子とも開平出来た.

一方 1/φも計算すると(6)(7)同じ値になった. めでたしめでたし.

曜日の計算

久し振りに曜日の計算の話だ. 2011年5月25日の私のブログ(曜日の計算)に

「21世紀中の 2000+y年m月d日の曜日は, (y+floor(y/4)+h(m)+d+4) mod 7 を計算するとそれが曜日になる.」

ただし「うるう年の1月と2月はh(m)の値から1を引く必要がある.」

とあり,

m 1月 2月 3月 4月 5月 6月 7月 8月 9月10月11月12月
h   2   5   5   8   3   6   1   4   7   2   5   0

の表がある.

今日の話題は上の式の「+4」と直ぐ上のh(m)である. 21世紀は+4だと書いたが, では他の世紀ではどうなるか. 20世紀に当て嵌めて計算してみると, その時代には+5であった.

我々の生きている20, 21世紀中は, この+5と+4を覚えておけばいいが, もう少し一般的に出来ないかと考えてみた. グレゴリオ暦は4世紀で曜日が元に戻るから, 年を100で割った整数部をさらに4で割った剰余に対してこの値を決めればよいわけだ.

1900年からは(剰余が3で)+5
2000年からは(剰余が0で)+4
もっとやってみると
2100年からは(剰余が1で)+2
2200年からは(剰余が2で)+0
だから, この剰余の0,1,2,3に対応して +4,+2, 0,-2と覚えておけばよいことが分った.

剰余  0  1  2  3
b    +4 +2  0 -2

上の式は (y+floor(y/4)+h(m)+d+b(c)) mod 7 になる.

h(m)の表は, 1月の2, 2月の5は覚えるとして, 後は奇数の月も偶数の月も5回ずつあり, 偶数月は12-m, 奇数月は3,5,7月は8-m, 9,11月は16-mと覚えるのがよさそうだとこの頃は思っている.

使ってみよう. グレゴリオ暦に改暦された1582年10月15日の曜日は (82+20+2+15-2) mod 7=117 mod 7=5 だから金曜. 日本でグレゴリオ暦に改暦された1873年1月1日の曜日は (73+18+2+1+0) mod 7=94 mod 7=3 だから水曜であった.

菱形六十面体を織る

前回のこの題名のブログの終の方に「ちょっとややこしいのが糊しろの部分で, 正六面体のように簡単ではない. その説明はまたの機会にしたい.」と書いたその部分の説明だ.

下の絵を見てほしい. 線がごちゃごちゃ描いてあるが, 前回の最後の図の各帯が菱形六十面体のどの面を形成しているかを示している.



例えばAの帯は, 真中少し下のAの回りのオレンジ色の線に対応している. Aの, 時計でいえば1時の方向に0の面があり, そこから右下へ進んでαの辺から外へ出, 少し下のαの辺から展開図の戻り, 一つおきに1,2,3の面を通り, θから再び外へ出, 上のθの4から戻ってくることを示している.

緑で示すBの帯は, 5,6,7,8,9だが, 5はIの五角形の下の方にあり, 右上に進み, γから出て, 右下のγから戻り, Bを取り巻くように6,7,8,9が見える.　Bの場合, 外へ出るのは1回である.

KとLは, 外へ出ないから遥かに簡単である.

ところでこれらの閉曲線には, その色で塗りつぶした三角形と, 白抜きの三角形が書き込んである. 実はそこが糊しろである. 白抜きの方が出発点. 線上にある三角の頂点の方へ進む. これが前回の帯の絵の上端に相当する. Aで言えば0の上の白い部分(下の49と重なる部分)であり, 塗りつぶしの方が下端になる. 線上の三角の頂点の方から下端に到る.

これを眺めると, K, Lは別として, 他の帯は上の五角と下の五角を貫いていて, A,B,C,D,Eは下では円弧, 上では直線であり, F,G,H,I,Jは逆になっている. K以外の糊しろはすべて上の五角の方にある.

つまり下の五角は出来上がった正十二面体の下半分の鉢を構成し, 上の五角は上半分を構成していたわけだ.

底の方から組み立てて, 糊しろが込み合う上面で組み立てが完了するような構造になっていた.

だから組み立て方は, A,B,C,D,Eをほぼ中心で重ね合わせ, その周囲をKの帯で固定し, その辺からF,G,H,I,Jも加わって下の鉢を作り, 中段を過ぎて上の周囲の面に至り, それをLで補強した後, F,G,H,I,Jを上の面に相互に差し込むのである. (糊しろといっても糊をつけて貼るのではない.)

元々の帯が巧みに設計されていて, この図を描くにも, 円弧と円と直線の基本の曲線群を用意してしまえば, たちまちにして完成した次第だ.

EDSACのプログラム技法

EDSACのプログラムを集積したウェブページがある. (詳しくはそこの "To browse the collection of software, click this link."をクリックする.) ここには私の書いたエラトステネスの篩も置いてあって嬉しい(21行目). そこにDavid J. Wheelerの素数のプログラムがあった(15行目). Martin Campbell-Kellyの書いたEDSACの文献(Programming the EDSAC: Early Programming Activity at the University of Cambridge, IEEE Annals of the History of Computing, Vol 1, No. 2.)によると, EDSACには「EDSACの本」で有名なイニシアルオーダーの前にもう一つのイニシアルオーダーがあって, 1949年5月6日の運転開始の日にはその方が使われていた. 同年6月22日にケンブリッジで開催された計算機会議でEDSACが公開され, その時, Wilkesの書いた平方の表とWheelerの書いた素数のプログラムのデモが行われた.

Wheelerは1927年生まれだから, 22歳のWheelerが当時どういうプログラムを書いていたかには非常に興味がある. これはもう読むしかない.

ところがその内容たるや

[Primes]
T107SO92SO93SO94SS5ST6SO95SO95ST7S
A96SR4SS97SE42SL4SA97SG45SS98SG100ST7S
A97SA4ST97SH97SN97SL64SL64SA96SE39ST7S
A96SU1SA4ST96SA99ST97SS88ST7SH91SA1S
E72SV91SS89SE71SA89STLOSH90SV1SL4S
TLA7SA98SG67SA6SA4SG36SE33SP2SP500S
JSP16S#S@S&S!SP2LP1LPLP1L
T7SA99ST97SA4SA96ST96SE39S

だけなのだ. これを読む手がかりはここにある.

読んでみると, やはり面白いプログラムであった. 特にそのプリントルーチンには感心した. 以下にその解説をしたい.





31番地 T107S. このプログラムは最初のイニシアルオーダーInitial Order1で読み込むが, 読み込まれるプログラムが106番地まであるなら, 最初をT107Sとするのであった. つまりイニシアルオーダーには格納命令があり, 最初はT31Sになっていて, 31番地から命令を格納する. つまり, T31S, T32S, ..と変わっていく. そして格納作業の終りは, 格納命令から31番地にある命令を引き, それが正であるかを見る. T107Sを引くと0(正)になるので31番地へ進んでくる. このT命令はアキュムレータのクリア用にも見えるが, すでに0になっている.

32番地は92番地にある#S(figure shift)を出力する. 以下のプログラムが数字を出力するので, その準備.

33,34番地では93, 94番地にある@S(cr), $S(lf)を出力. 次のS5Sはイニシアルオーダーの5番地にあるP5S(10*2^-16)を引き, -10を作り, 6番地(m[6])に入れる. 6番地は -10,-8,-6,..と変わり, 1行に5個の素数を印刷する.

37, 38番地. !S(space)を2個出力. (このプログラムのあちこちにある)T7Sは(67番地のものを除いて)アキュムレータをクリアする.

このプログラムでは素数の候補pから3,5,...と奇数(d)を引き(以後 dを減数という), 0になるかどうかを見ている. (前にも書いたようにEDSACには除算命令はない.) pは96番地, dは97番地にあって, それぞれ5と3に初期化されている.

40から46番地はそのpからdを引く. ただ初めは16*dを負になるまで引き, 次にdを正になるまで足す. これでpをdで割った剰余が得られる. 実際には16*dを引くのではなく, pを4ビット右シフトして引き, 4ビット左シフトして足す. (おおざっぱにいうと単にdを引くより, 16倍加速されている. こういう剰余の取り方はEratosthenesの篩でも使える.)

47番地. EDSACには零ジャンプがないので, 正になったところで1を引き, 負なら零であったということで, 合成数の判定をして100番地へ進む.

剰余が0でなければ, 50,51,52番地でdを2増やし, 53番地からH97Sで乗算レジスタにdを置き, N97S, L64S, L64Sで-d^2を作る. それにpを足し, 減数の二乗が素数候補より小さければ39番地へ戻り,　次のテストに掛かる.

減数が大きくなって素数と判定出来ると, 60番地から先で素数の出力にとり掛かる.

U1Sで素数を1番地に入れ, A4Sでpを2増やして63番地で96番地に入れておく.

素数は4桁で出力するつもりで, 66番地で-4を7番地のカウンタへ入れる.

68番地から76番地は1桁の出力で, 1番地に1234があるとして説明すると, 68番地H91Sで乗算レジスタに2^-11を入れる(32*2^-16=2^-11).

アキュムレータに1234を足す. 72番地へ飛んで89番地にある1000を引く. アキュムレータは234(*2^-16)になり, 正なので71番地へ行く. V91Sだからアキュムレータは234*2^-16+2^-22になる.

再び1000を引くと今度は負になり, 74番地で1000を足し, アキュムレータは234*2^-16+2^-22に戻る.

これを0番地の長語に入れる. つまり左半分の1番地は234, 右半分の0番地はファンクション部に1となる. そこでo命令は0番地を出力するから1が印字される(EDSACの文字コードは2012年9月28日のこのブログ参照).

77番地から80番地は234を10倍する. H90Sで乗算レジスタに10*2^-4を置き, V1Sでアキュムレータは234*10*2^-4になり, 次の左シフトで2340になり, 1番地に戻される.

81番地からは7番地のカウンタを1足し, まだ4桁出力していなければ, 67番地へ戻る.

次は2340だから, 71番地のループは2回実行され, 右半分の方は2*2^-11になって2が印字され, 1番地は3400になる.

4桁出力すると84番地へ来て, 6番地のカウンタにA4Sで2を足し, 5個未満なら36番地, 5個終わっていたら33番地へ戻る.

一方合成数の時は100番地へ進み, 101, 102番地でdを3に戻し, 103,104,105番地でpを2増やして39へ戻る.

この解釈が正しければ, 期待される出力は

  0005  0007  0011  0013  0019
  0019  0023  0029  0031  0037
...

なる筈であった. サイボウズ・ラボの西尾君がEDSACのシミュレータを持っているというので, 試してもらったらその通りであったとメイルが来た.

このプログラムの疑問は素数候補が5から始まることである.

減数を3から始めると素数候補の3は割り切れるから, 3は素数にならない. 3を素数にするには2から割り始めなければならないからなのだ.