2016年11月30日水曜日

予定表

徐々に年末が近付きつつある. 年末独特の心配事は他にもあるが, まずは来年の予定の記帳だ. 通常の手帳は12月(かせいぜい1月かそれも第1週だけとか)で終りであり, 次の会合が翌年の始めだったりすると, 書き込む場所がない.

また例年に送られてきて使っている企業の手帳も年末にならないと届かない.

そういうわけで, この数年, もしかしたら十数年は10月ころになると, 次のような簡単な予定表をA4の紙に印刷し, 手帳に挟んで利用していた.

下に示すのは2017年版である. 左上の01S, 02M, ...は1月1日(日曜), 1月2日(月曜), ... を表す. 右端の列の31Fは3月31日(金曜)に対応する.

1日分が非常に小さいが, この予定表を作るころは, 新春の予定も少いのであまり困ることはなかった. 絵から分るように日曜だけは太字にしてある.

もちろん新しい手帳が届くまでの使命で, 手帳が来るとそちらへ転記し, お役御免になる.



ところが, 世の中では皆さんがスマートフォンで予定を管理するようになり, 学生のころから毎年使っていた手帳は, 企業でも今年(2016)の分から遂に作るのを止めた.

今年はとりあえず市販の手帳に予定を記入してきたが, やはり使い易いとは思わなかった. そこでこの暮は, 来年の予定表を自作することにした. その方が楽しくもある.



上の写真がこうして出来た個人用予定表の表紙である. 右側が表紙, 左側が裏表紙である. 大きさはB7版になる. (こうして開いた状態はB6版である.)

裏表紙はすぐ分るように, 2017年の七曜表である. しかし, 最上段が左から1月2月3月, 次の段が4月5月6月, ... ときて, よく見ると5段あり, 最下段は実は2018年の1月2月3月なのである.

この予定表の特徴は, 2017年1月から2018年3月までになっていることだ. つまり3ヶ月のバッファーを準備している.

従って, 表紙の次は



のように2017年1月になり, 裏表紙の前は



のように2018年3月である.

1日分が狭いのではないかといわれそうだが, 最近は予定表に記入する事項も少くなり, たとえ沢山書き込む日があっても, 欄外に書いて線で結ぶリンク構造にすれば楽々収まる.

そういう次第で来年はこの予定表を活用しようと考えている.

いい遅れたが, 表紙の西暦年号の下は曜日(上段)と月齢(下段)を計算する定数(枠内は左から1月2月3月 | 4月5月6月 | 7月8月9月 | 10月11月12月)である.

どちらも私のブログの話題になったもので, m月d日の曜日は左欄外の4にm月の定数とdを足し, 7で除した法をとる. 月齢は左欄外の1に足し, 30で除した法をとる. 右欄外は2018年の定数である.

曜日の定数が1月と2月で赤になっているのは, 閏年には1を引くことの注意である. 閏年には欄外の定数も赤で表示するようにプログラムしてある.

こういう予定表を作るには紙の中央で綴じる必要がある. 中綴じホチキスを探したが, 「マックス多機能スケール ナカトジール」というものを見付けそれを購入して使った.

ただ, 上の写真で分るように, 丁度折り目のところで留めるのは結構むずかしい.

予定表を吐き出すPostScriptのプログラムを書くのは, これまでカレンダーのプログラムを何度も書いたことがあるので, 簡単である. 来年, 再来年の祝日はインターネットで探すとすぐに見付かる.

2016年10月5日水曜日

ループ検出

前回のこのタイトルのブログで, Gosper流のループ検出法の説明をしたが, これでループが発見できる理由があれだけではいまいちである.

私がどう納得しているかも書いておく必要があろう. xiが次々と置かれていく様子を添字だけで再掲する.

iT0 T1 T2 T3
1 1
2 1  2
3 3  2
4 3  2 4
5 5  2 4
6 5  6 4
7 7  6 4
8 7  6 4  8
9 9  6 4  8
10 9 10 4  8
これらの表中のxiに, 次のxがあるかどうかを調べるのだが, そこで添字の差だけの表にすると
iT0 T1 T2 T3
1 -1
2 -2  -1
3 -1  -2
4 -2  -3 -1
5 -1  -4 -2
6 -2  -1 -3
7 -1  -2 -4
8 -2  -3 -5  -1
9 -1  -4 -6  -2
10 -2  -1 -7  -3
つまり1つ前2つ前に同じものがあると, それらはT0で見付かり, 3つ前4つ前にあると, それらはT1で見付かり,...というわけだ. Tの添字が大きいほど, 繰り返し開始時のμの値の曖昧さが大きいことがこの表から判明する.

TAOCPでは演習問題3.1-7にこの話があり, その解答には

xnが表Tρnに格納されると, それはその後xn+1, xn+2, ..., xn+2ρn+1と比較される.

と記載してあるが, 上の説明の方が分り易いと私は思っている.

Gosperの元の記事はHAKMEM 132 (http://home.pipeline.com/~hbaker1/hakmem/flows.html#item132)にある.

2016年9月22日木曜日

ループ検出

数列xiが例えば32ビットの語の範囲の値しかとれなくて, xi+1=f(xi)を繰り返していくと, いつかは昔のある値が現れ, 以後同じ数列を辿ることになって, ループを回り始める.

このループの開始位置とループ長を求めるのがループ検出である.

そのひとつの方法にBill Gosperの考案したものがあり, あちこちに説明があるが, あまり真面目に読んだことはなかった. しかし気にはなっていたので, この数日プログラムを書いてみたりしながら, その仕掛けを調べてみた.

いまx1, x2, x3, ... が1,2,3,4,5,6,7,8,9.10,6,7,8,9,10,6,7,8,...であったとしよう.(x0から始める流儀もあり, また私はそれが好きだが, 今回の添字は1から始めることにした.) The Art of Computer Programming(TAOCP)の演習問題3.1,6の記述を借用すれば, 「x1,x2,...,xμ,..., xμ+λ-1は互いに相異なるが, n ≥μの時はxn=xnであるμとλがある.」

このμが開始位置, λがループ長で, 上の例ではμ=6, λ=5だ.

記憶場所も時間もふんだんにある時は, xの値をすべて順に記憶しておき, 新しいxについて, それを過去のすべてのxの値のそれぞれと比較する. そこに同じものがあれば, 前の値のあったところがループの開始場所であり, そこから今の値までの長さがループ長である.

しかし, こんなことはとてもやっていられないから, あの手この手を考えることになる. その一つにGosperの提案したものがある. 以下は添字などを私流に多少修正したものである.

これは簡単にいえば, 繰り返しが「現れるや否や」発見するのは諦める, 最新のxの1つ前, 2つ前, 4つ前, 8つ前, ... などを覚えておき, ループに入り同じ値が現れ始めたらなるべく早目にそれを発見しようという方針である.

「1つ, 2つ, 4つ, 8つ,...」から推察されるように, n個のxを調べた時は, log2 n個程度の配列Tを使う.

ところで, n>0を二進法で表した時, もっとも右の1の位置, もっとも左の1の位置をTAOCPに倣い, それぞれρ(n), &lambda(n)とする. ビットの位置は右端を0とする. (上のループ長の定義にもλがあって申し訳ない.)

ρ(1)=0, λ(1)=0; 
ρ(2)=1, λ(2)=1; 
ρ(3)=0, λ(3)=1; 
ρ(4)=2, λ(4)=2; 
ρ(5)=0, λ(5)=2; 
ρ(6)=1, λ(6)=2; 
ρ(7)=0, λ(7)=2; 
ρ(8)=3, λ(8)=3; 
ρの方は別名ルーラー関数(物差関数)という.



十進法の物差しとは様子が違い, ひと目盛りが半分ずつに分割されていて, インチの物差しが大体こういう感じである.

ρとλを計算する関数 rhoとlamは次のようだ.

ρは
(define (rho x)
 (define (r x rh)
   (if (= (modulo x 2) 0) (r (quotient x 2) (+ rh 1)) rh))
 (if (= x 0) 'error (r x 0)))

(map rho (a2b 1 16)) => (0 1 0 2 0 1 0 3 0 1 0 2 0 1 0)
またλは関数logがあれば
(define (lam n)
 (inexact->exact (floor (/ (log n) (log 2))))

(map lam (a2b 1 16))=>
(0 1 1 2 2 2 2 3 3 3 3 3 3 3 3)
ρのように自力で計算するには
(define (lam x)
 (cond ((= x 0) 'error)
       ((= x 1) 0)
       (else (+ (lam (quotient x 2)) 1))))
さて配列をT0, T1, T2, ... とする. (こちらは0から始める) そしてxiTρ(i)に順に置くのである.

xiが次々と置かれていく様子は次のようだ.

iT0 T1 T2 T3
1x1
2x1 x2
3x3 x2
4x3 x2x4
5x5 x2x4
6x5 x6x4
7x7 x6x4
8x7 x6x4 x8
9x9 x6x4 x8
10x7 x10x4 x8
注意すべきは, この表のうち見えているのは最後の行だけということである. つまりi=8の時には, x1, x2, x3, x5はすでに消えている.

さて, 上の例のようなループがあった時の同様な配列のx10=10までの変化の様子は,
iT0 T1 T2 T3
1 1
2 1  2
3 3  2
4 3  2 4
5 5  2 4
6 5  6 4
7 7  6 4
8 7  6 4  8
9 9  6 4  8
10 9 10 4  8
のようだ. ここでx11=6になった時, この配列は9, 10, 4, 8であり, この中に6はないから, 6を(ρ(11)=0だから)T0のところに置く. したがって配列は6, 10, 4, 8になる. 次, x12=7も配列にないから, 7をT2に置き, 配列は6, 10, 7, 8になる.

次はx13=8だが, i=8の時に作ったT3に8があるので, 循環が始まっていたのがわかる.

そこで問題はiの行のTkを置いたのは誰かということだ. 最初の表からxを削除した表は次のようだが
i\k 0  1  2  3
11
21 2
33 2
43 24
55 24
65 64
77 64
87 64 8
99 64 8
107 104 8
これで例えばi=7, k=2ならこの4を置いたのはi=4であったということが知りたい.

いろいろやってみると, iからiの下からk+1ビットをとったものの, 先頭の0,1を逆にしたものを引くといいらしいことが分かる.

iの下からk+1ビットをとるには, 2k+2-1で1のビット列を作ってマスクする.

i &(2k+2-1)

これを 2kでxorすると先頭のビットが反転する.

それをiから引くのである. Schemeで書くと
(define (m i k)
   (- i (fix:xor (fix:and i (- (fix:lsh 1 (+ k 1)) 1)) 
    (fix:lsh 1 k))))
であって, 実際に計算すると
(do ((i 1 (+ i 1))) ((= i 11))
(display (cons i 
 (map (lambda (k) (m i k)) (a2b 0 (+ (lam i) 1)))))
(newline))

(1 1)
(2 1 2)
(3 3 2)
(4 3 2 4)
(5 5 2 4)
(6 5 6 4)
(7 7 6 4)
(8 7 6 4 8)
(9 9 6 4 8)
(10 9 10 4 8)
となって上の表と一致する.

従ってi=13の時のx13>=8がi=12, k=3にあり, (m 12 3)=8が分ったから, λ=13-8=5であった.

プログラムを書いてやってみよう. 最初の数列xiは, μとλから
(define (f i mu lm)
 (if (< i mu) i (+ (modulo (- i mu) lm) mu)))

(map (lambda (i) (f i 6 5)) (a2b 1 19)) => 
(1 2 3 4 5 6 7 8 9 10 6 7 8 9 10 6 7 8)
と作る. プログラム全体と実行結果を下に示す.

(define (iloop i x)
 (define (kloop k)
  (if (> k (lam (- i 1))) (iloop i x)
  (begin 
   (if (= x (vector-ref tab k)) (list i k (m (- i 1) k))
     (kloop (+ k 1))))))
(vector-set! tab (rho i) x)
(display i) (display (take (+ (lam i) 1) (vector->list tab))) 
(newline)
(set! i (+ i 1))
(set! x (f i mu lm))
(kloop 0))

(define mu 6) (define lm 5)
(iloop 1 (f 1 mu lm))

1(1)
2(1 2)
3(3 2)
4(3 2 4)
5(5 2 4)
6(5 6 4)
7(7 6 4)
8(7 6 4 8)
9(9 6 4 8)
10(9 10 4 8)
11(6 10 4 8)
12(6 10 7 8)

=> (13 3 8)
つまりi=13のとき, k=3の所に同じものがあり, その時のiは8であった. 従ってλは13-8=5であった.

2016年8月13日土曜日

複素数用計算尺

複素数用の計算尺があるという話をブログに書いたのは, 2009年8月のことだ.

最近すこし図を描き直し, Webでも楽めるようにしたので, その話をしたい.



上の図は対数の複素数を示す. 点Aから右へ延びる実線は1より大きい正の実数軸. Aは(虚数部が0の)複素数1. Bは複素数2, ...で, 2,3,4,...と対数的に間隔が狭くなっている.

Aで実数軸と交わる実線の曲線は, 実数部が1で虚数部が何かの複素数である.

一方, 点Cから右へ延びる破線はiより大きい正の(?)虚数軸. Cは(実数部が0の)複素数i, D は複素数2i, ... で, こちらも2i,3i,4i,...と対数的に間隔が狭くなっている.

Cで虚数軸と交わる破線は, 実数部が何かで虚数部がiの複素数である.

従って点Eは, 複素数1+iを示す.

実はこの上の図は下の真数のGauss平面を対数化したもので, 真数の図では実数軸と虚数軸の交点に0があるが, 対数では0は無限に遠くにあるので, 実数軸と虚数軸が平行になってしまう. 対数の図で, 実数軸と虚数軸の間にある部分は, 真数のGauss平面の第一象限に対応する.



この図には描いてないが, 正の実数軸の反対側には負の実数軸があるわけで, それは最初の対数の図では, Cから右に伸びる虚数軸のさらに上に現る. Cから上へACの距離だけ行ったところに, 負の実数の出発点 -1があって, 正の実数軸と同様に=1, -2, -3,...が右に出てくるのである.

さらにその上というか, Aから出発する正の実数軸の下というかに, 負の虚数軸があるのだが, 図がごちゃごちゃするから省略してある.

省略はそれだけではない. 10, 10i -10などの右には, これと同じパターンがあって, 20,30,... 20i,30i,... -20,-30,...と続くのである.

また左には0.9, 0.8,...があるのは理解できるであろう. つまり上の図はほんの一部を示しているに過ぎない.

もっと広範に描いて見ると, 下のような図になる.



前置きが長くなった. この図を使って複素数の乗算をするにはどうするかを述べよう. 通常の計算尺では, 2つの対数尺A, Bが次の図の上のように上下に並んでいて,

2掛ける何かを計算したければ, 下の図のようにB尺の1をA尺の2に合わせ, 例えばB尺の3に相対するA尺の目盛り6を見て, 2×3は6なることを 知るのである.

対数でも同様で, 1+iに何かを掛けようとするなら, B尺に相当する赤い図の1(点A)をA尺に相当する黒い図のEに合わせ, B尺の何か(例えば1+i(E))に 相対するA尺の値(Dつまり2i)を知る.



(1+i)×(1+i)=12+i2+2iで, 1の2乗とiの2乗の和は消えるから2iで正解だ.

このwebプログラムは http://www.iijlab.net/~ew/complexsr.htmlに置いてある. B尺のAを置きたいところをクリックすると, A尺が薄い色になり, 濃い色のB尺の1がその場所へ移動するから, B尺の乗数に相対するA尺を読めばよい.

2016年7月7日木曜日

プラニメーター

何日か前に書いたプラニメーターのブログに機能をもうひとつ追加した.

シミュレーターの画面にsquare, circleなどクリックすべき箱があったが, その隣に無名の箱があった. 実はここには曲線で囲まれた図形の例を置くつもりであった.

それが出来たので, シミュレーターを更新した. そのcurveの箱をクリックすると次のような図になる. この水銀整流管のようなものは「公立はこだて未来大学」のロゴである. その面積を求めることができる.



startをクリックするとメーターに



134.179が現れる. 従って面積は(* 134.179 400)=>53671.6である.

ところでこの自動追尾の方法だが, もとのロゴのビットマップを次のように作り,

0000000001FFFFF800000000
000000003FFFFFFFC0000000
00000001FFFFFFFFF8000000
0000000FFFFFFFFFFF000000
0000003FFFFFFFFFFFC00000
000000FFFFFFFFFFFFF00000
000003FFFFFFFFFFFFFC0000
00000FFFFFFFFFFFFFFF0000
00003FFFFFFFFFFFFFFFC000
0000FFFFFFFFFFFFFFFFE000
その輪郭線をとりだすプログラムにかけると
0000000001FFFFF800000000
000000003E000007C0000000
00000001C000000038000000
0000000E0000000007000000
000000300000000000C00000
000000C00000000000300000
0000030000000000000C0000
00000C000000000000030000
00003000000000000000C000
0000C0000000000000002000
のようなものが得られる. この線を辿るデータを作るプログラムをSchemeで書いた. シミュレーターではそのデータを追いながらプラニメーターの先端を動かしているのである.

ところで始めのビットマップに1がいくつあるか数えてみた. 54003個あった. プラニメーターの結果は53671.6だったから, まぁ大体うまくいっていると思ってよかろう.

2016年7月4日月曜日

プラニメーター

小学生の時, 算数の時間にいろいろな図形の面積の求め方を習った. その折, 先生から三角形とか台形とかではない複雑な図形の面積はどう求めるか聞かれた.

私の当然の答は方眼紙に写しとり, 図形内部の1ミリ方眼を数えることであった. もうひとつの答は, 厚紙に写しとり, 図形を切り抜いて重さを測り, それと同じになる矩形の厚紙を作り, その面積を求めるのであった.

後者の答は象の重さを測るのに, 象を船に乗せ, その時の吃水線の位置を船端に印し, 次に同じ吃水線に船が沈むまで小石を船の乗せてその小石の総重量を測るという中国の子供の話にヒントを得ていたかもしれない.

その後大学生くらいになって, プラニメーター(面積計)なるものがあることを知った. しかし, 自分で使う機会はなかった.

東大工学部計数工学科に勤めるようになった時, その学科では以前の学生実験でプラニメーターを使うものがあったと聞いたが, すでにプラニメー ターやタイガー計算機を使う実験はなくなっていた.

プラニメーターは機械式アナログ計算機である. プラニメーターの図は当今ではインターネットで探すといくらでも見付かる.

その一つが次である.

とあるページから借用させていただいた.



これは詳しくいうと Amsler polar planimeter という形式のものだ.

使い方は右上のOは下にピン, 上に重りのあるもので, 机のある場所に重りで固定する.

そして右下のAのピンで測りたい図形の周囲を一巡する. それに従い, 2本の腕は回転したり腕の間の角度を変たりしなが動く. 腕の蝶番BのあたりにCという車輪がみえるが, この車輪の軸はBAの腕に平行で, 特定の方向の回転だけを蓄積する.

すると不思議にも車輪の移動した長さにBAの腕の長さを乗じたものが面積になるのである. 長さ×長さだからディメンジョン的には合っている.

今回はプラニメーターの原理はかつ愛し, 最近プログラムしたプラニメーターのシミュレーターのことを書きたい.

そのプログラムは
http://www.iijlab.net/~ew/planimetersim.html
に置いてある. そこへ行くと



のような図が現れる. 左下のcircleをクリックすると, 今度は



のような図が現れる.

左端の赤の点が固定点(最初の図のO)で, 四角や円の上の青の点が図形をトレースする点(同A)である. その間の2本の線が, プラニメーターの2本の腕で, 折れ曲っているところが腕の蝶番(B)に対応する. 折れ曲りの近くの直交している線分が, 計測用の車輪の位置のつもりだが, その回転が左上のメーターみたいなもにに示される.

たとえば四角の図を表示している場合, 左下のstartをクリックすると, 青の点が移動し, 計測が始まる. 一周すると停止する. その時のメーターの読みは399.99になっている.



この図の四角は一辺が400の正方形で, プラニメーターの二の腕(BA)の長さも400になっている. メーターの読みが約400なので, それに腕の長さ400を乗ずると, 正方形の面積160000が得られる.

計測点Aの移動はプログラムで制御しているとはいえ, 凄い精度である.

次に円をやってみよう. 円は直径が400, 半径200なので, 面積はπ×200×200である.

それを腕の400で割ると, メーターの読みはπの100倍になるはずである. 実際にやってみると314.156が得られて安心する.



しかし, 手動でプラニメーターを操作すると, もちろんこんな精度は得られない. だからといって, 正方形の場合, 計測点を移動する際に物差しを当ててずらしてはならないと学生実験の手引書に注意がある. つまり線から一方に外れているとその誤差がそのまま集積されるからで, 手動でふらふら動かす方が誤差が左右に揺れてキャンセルされるので望ましいと書いてあった.

2016年5月1日日曜日

Prefix adder

60年くらい前, 東大物理学科高橋研究室でパラメトロン計算機PC-1を作っていたころ, 36ビットの並列加算器で, 高速に繰上げ(carry)したいという問題があった.

後藤さんはパラメトロンの多数決素子を巧妙に使い, 桁数に対して対数時間で繰上げを検出する回路を考案して, それを取り入れた.

先日出版されたDigital Design and Computer Architecture Arm Editionを眺めていたら, 同じような考えによる, Prefix Adderという回路の話があった. 今回は それが話題である.





この上の図が16ビットの二進数Ai, Bi (i=0,1,...,15)と下から来る繰上げCinを足すときのCarry-Lookaheadの回路だ.

上の図の上段の白い四角, 中頃の黒い四角, 下段の角のとれた白い四角の説明が下の図である.

ある桁 i が, 下の桁 i-1 からの繰上げの有無に関係なく, 上の桁 i+1 へ繰上げを生じるとき, 繰上げをgenerateするということで, Gi=1, そうでないとき, Gi=0とする. AiとBiが共に1なら繰上げがあるから, Gi=AiBi.

ある桁 i が, 下の桁 i-1 からの繰上げを上の桁 i+1 へ伝えるとき, 繰上げをpropagateするということで, Pi=1, そうでないとき, Pi=0とする. AiとBiのいずれかが1なら繰上げがあるから, Pi=Ai+Bi.

そうすると, i からi+1へ出る繰上げCi

Ci=Gi+PiCi-1

になる. このGやPをPrefixという.

これは単一の桁についての話であったが, iからjの桁のブロックについても同じことがいえる. ただGi, Piの代わりにGi:j, Pi:jのように書く.

そうだと思って上の図を見てみると, 左下のあたりに黒四角で左下隅に14:-1と書いたものがある. つまり15の桁への下全体のブロックからの繰上げはG14:-1 (と, P14:-1 と, Cin) で決るということを示す.

GiとPiがAiとBiとで書けたように, Gi:jとPi:jがそのブロックを途中で分けたGi:kとGk+1:j, Pi:kとPk+1:jで同じように書くと, Gi:jはi,k間で繰上げが出るか, k+1,j間で繰上げが出て, i,k間で繰上げが伝わるかだから, Gi:j=Gi:k+Pi,kGk-1,j.

Pi:jはi,k間でもk-1,j間でも繰上げが伝わるかだから, Pi:j=Pi:kPk-1:j.

ここまではDigital Designの本に書いてあることだ. 上の図のように半分, 半分と区間を分るにはi,jからkをどう決めればよいかは, この本には書いてない. そこがこの本とThe Art of Computer Programmingの違うところである. TAOCPならk(i,j)が詳しく記述してある筈だ.

しかしそう難しいことではなくて, λx(x>0)をxを二進法で表したとき, 最も左にある1のビット番号とすると, (λ1=0, λ2=1,λ3=1, λ4=2,...)

k=2λ(i-j)+j

と簡単に得れれる. Schemeでシミュレートしてみたが, これで合っているようであった.

2016年3月1日火曜日

菱形六十面体

菱形六十面体についてはこれまで何度もこのブログに書いた. 夫々の菱形が黄金比であることは, MathWorldにそう書いてあったのを鵜呑みにしたきらいがある.

今回はそれを確かめてみたい.

下の図は正十二面体である. 中心を通る赤, 緑, 青の線はx, y, z軸を示す. それぞれの軸は正十二面体の30本ある稜線のうち6本の中央を通過する. 次の図との関係で偶数だけがついているが, その頂点を持つ正五角形で座標を考える.


頂点間の距離, 稜線0, 2の長さを2とすると, 正五角形の対角線8, 4の長さは2φ. 従って4の頂点の座標は(φ, φ, φ). それを基準に他の座標も分かって, 1は(φ2, -1, 0);2は(φ2, 1, 0); 6は(1, 0, φ2); 8は(φ, -φ, φ)である. (1のx座標は中心から0までの距離と, 中心から4までの距離が等しいことから計算できる.) この正五角形に5個の菱形を作ったのが下の図の菱形六十面体だ.



これをx軸の方向から撮影した写真がこれである.



写真では視差があるから, 注目している正五角形のついて正確な図にしたのが次の図だ.



ちょうど1のところをx軸がこちらへ向かっている. その両側へ1離れたところに02が来る. また1から真上に行ったところに中心10があるわけだが, その高さは現時点では未定.

さて菱形は真上から見れば菱形であるが, 少し横から見れば菱形に見える方向も皆無ではないものの, 一般には菱形に見えない. しかし平行四辺形であるのは確かである.

従って0,1,10,9の菱形, 1,2,3,10の菱形を考えると, 0,91,3の辺は1,10の辺と平行である. また9,100,1と平行なども分かる.

すると8,75,40,1などに平行である. 一方 10,5, 6,7の菱形は対角線の方向から見ているから菱形のままであり, 10,55,6の長さは等しい. 5, 6の高さの差(φ2-φ=1)は既知だから, 10の高さも確定し, 計算するとφ-1=1/φである. 5,4の幅が1だから, 5のy座標も1/φと決まる.

ところで残りは奥行きである. 2, 3 の辺は, 横から見ると, 正面図の6, 5の辺に相当していることが分かる. 従って3の奥行きが分かり, 3, 4 がz軸に平行なことも分かる. 従って10, 5, さらに5, 6もz軸に平行だったことも分かった.

つまり正面から見た菱形10,5, 6,7は真の形であった.

という次第で菱形の縦横比は1:φであったわけである.

2016年2月24日水曜日

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する.

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



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

2016年2月15日月曜日

月齢カレンダー

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

それではその年の月齢定数をCとすると, その年のmd日の近似的な月齢を (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の年が回って来る.

2016年2月14日日曜日

菱形六十面体

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) a3とあり, また菱形六十面体の項にはその体積が4√(2(5+√5))a3とあるから, 体積的にも確かにそうである.

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

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



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

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

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

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



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



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

全体
北極
赤道
南極

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

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

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

2016年2月3日水曜日

菱形六十面体

菱形六十面体についてこのブログで何度も取り上げた. その最初, 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)同じ値になった. めでたしめでたし.

2016年1月25日月曜日

曜日の計算

久し振りに曜日の計算の話だ. 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 だから水曜であった.