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の年が回って来る.