プラニメーター

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

シミュレーターの画面に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だったから, まぁ大体うまくいっていると思ってよかろう.

プラニメーター

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

私の当然の答は方眼紙に写しとり, 図形内部の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が得られて安心する.



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

Prefix adder

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

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

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





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

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

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

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

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

C_i=G_i+P_iC_i-1

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

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

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

G_iとP_iがA_iとB_iとで書けたように, G_i:jとP_i:jがそのブロックを途中で分けたG_i:kとG_k+1:j, P_i:kとP_k+1:jで同じように書くと, G_i:jはi,k間で繰上げが出るか, k+1,j間で繰上げが出て, i,k間で繰上げが伝わるかだから, G_i:j=G_i:k+P_i,kG_k-1,j.

P_i:jはi,k間でもk-1,j間でも繰上げが伝わるかだから, P_i:j=P_i:kP_k-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でシミュレートしてみたが, これで合っているようであった.

菱形六十面体

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

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

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


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



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



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



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

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

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

すると8,7や5,4も0,1などに平行である. 一方 10,5, 6,7の菱形は対角線の方向から見ているから菱形のままであり, 10,5と5,6の長さは等しい. 5, 6の高さの差(φ²-φ=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:φであったわけである.

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プリンタのお蔭である.