2013年12月11日水曜日

Christopher StracheyのGPM

Christopher Stracheyという名前を初めて見たのは1959年にパリで開催された第1回のIFIPの会議の報告書で, そのpp.336-341に彼の
Time sharing in large, fast computers
という論文があった.

要するに割込みの機能を使ってのマルチプログラミングの提案であるが, 当時我々東大物理の高橋研のパラメトロン計算機でも, すでに割込みと並列計算の実験を始めていたので, 同じようなことは, 洋の東西で同じ頃に気附くものだと思った.

そのStracheyに出逢ったのは, 1971年4月 英国Warwick大学で開かれたIFIP WG2.2の会議でであった. この会議には主査のMike Woodger(NPL)を始め, Edsger Dijkstra, Willem van der Poel, Simulaを開発したOle-Johan Dahlなど有名な計算機科学者が大勢いた. その会議中 パーティーが主催者のJohn Buxtonの超古い家で開催され, みなが車に分乗して行くことになった. その時昔の小さいMiniに乗せてくれたのが大柄なStracheyであった. 「Eiitiは後に入れ」と私を後部座席に押し込み, 助手席に乘ったBrian Randellと機関銃のような英語で話し合っていた.

さて最近のIEEE Annals of the History of Computing(July-September 2013)を眺めていたら, 何年か前に私が情報処理学会誌に「Wilkes先生を悼む」を執筆した時にWilkes先生の写真を 送ってくれたDavid HartleyさんがCPLのことを書いていて, その参考文献を辿っているうちに, CPLのアセンブラをGPM(General Purpose Macrogenerator)で作ったという話から, GPMを開発したStracheyがComputer Journalに寄せた論文
C. Strachey, A general purpose macrogenerator, Computer Journal Vol.8, No.3, 225-241
を探しあてた.

CPLはCambridge Plus Londonとも思えるが, Combined Programming Languageのアクロニムである. この言語はその後Martin RichardsがMITにいるころ, BCPLに単純化して開発し, さらにベル研究所に行ってC言語になった, ということを覚えている人はもう少ないに違いない.

ところでそのmacrogeneratorである.

macroというからにはマクロ定義とマクロ呼出しがあるわけで, aというマクロをb~1dと定義すると, aがマクロ名, b~1dがマクロ本体で, 本体のうちbとdは定数, ~1は第1引数である. 一方, $a,c;がマクロ呼出しで, (もとの論文では先頭は§であるがASCII文字にないので, 以下では$を使う.) マクロ名の本体に対して, cが第1実引数となって, 呼出しの結果はbcdとなる.

マクロ呼出しは $マクロ名,引数1,引数2,...; の形. このマクロ名は引数0としても使うことが出来る.

マクロ呼出しは $ と ; で, 定数は < と > で囲まれているから, マクロ名や引数の中に再帰的に書くことが出来る. その場合, 局所的なマクロ呼出しは呼出しの結果, 局所的な定数は, 一番外の < と > を外した定数本体でその部分を置き換える. 置き換えて出来た結果の文字列で外側のマクロを呼び出す.

GPMではマクロ定義もマクロ呼出しの形で, 上のマクロaの定義は
$def,a,<b~1d>;
のように書く. このマクロ呼出しは空の文字列を結果として返す.

マクロ本体は, この例のように, 本体を評価されたくない時は < と > で囲む.

マクロ呼出しの入れ子の例は次の通り:

引数にマクロ呼出しがある例
$def,a,<b~1d>;
と定義し,
$a,$a,c;;
と呼び出すと,
bbcdd
が返る.

マクロ名にマクロ呼出しがある例
$def,a,<b~1d>;$def,bcd,<b~1c~2d>;
と定義し
$$a,c;,e,f;
と呼び出すと
bdcfd
が返る.

Stracheyの論文にはsucとsuccessorというマクロの例がある.
$def,suc,<$1,2,3,4,5,6,7,8,9,10,$def,1,<~>~1;;>;
$def,successor,
 <$~2,$def,~2,~1<,$suc,>~2<;>;$def,9,<$suc,>~1<;,0>;;>;
sucは$suc,4;のように呼び出す. すると本体の1というマクロを呼び出すが, 引数の評価中にマクロ1が定義される. 引数~1が4なのでこのマクロの本体は~4になり, 1のマクロから5がとられて4の次の5が返る.

successorは 2桁の数の次を返す. $successor,3,4;だと3,5が, $successor,2,9;だと3,0が返るものだ.

3,4の場合はsuccessorの本体が
$4,$def,4,3,$suc,4;;$def,9.$suc,3;,0;;
となりマクロ4を呼び出そうとする. その$に対応する;は一番右にあるので, その前にマクロ4と9を定義する.
$def,4,3,$suc,4;;
だからマクロ4は3,5になり
$def,9,$suc,3;,0:
だからマクロ9は4,0になる. この状態でマクロ4を呼び出すから3,5か返る.

1の桁が9だと, マクロ9が2回定義され, その状態でマクロ9を呼び出すから後で定義した繰り上げのあるマクロの結果が返るのである.

この例にあるように, 本体評価中に定義されたものは, 評価が終わると定義のリストから削除されることに注意しよう.

Stracheyの元の論文にはCPLで書いたGPMの実装が載っている. CPLには再帰呼出しなどないから, スタックの上にリンクをつなげたりしていて, あまり読みたいとも思わない.

以下は私がSchemeで実装したものである. StracheyのGPMには四則演算の機能もあるが, そういうのはGPMの本質的なものではないから, 割愛した.

(define ch '())
(define (gpm str env args)
 (let ((outs "") (index 0))
  (define (getch)
   (set! ch (string-ref str index))
   (set! index (+ index 1))
   (if (or (char=? ch #\space) (char=? ch #\newline)) (getch)))
  (define (readquote)
   (let ((outs ""))
    (define (rq)
      (getch)
      (cond ((char=? ch #\>) outs)
            ((char=? ch #\<)
             (set! outs
              (string-append outs "<" (readquote) ">"))
             (rq))
            (else
             (set! outs (string-append outs (string ch)))
             (rq))))
    (rq)))
  (define (readstring)
   (let ((outs ""))
    (define (rs)
     (if (= index (string-length str)) outs
      (begin (getch)
       (cond ((char=? ch #\,) (cons ch outs))
             ((char=? ch #\;) (cons ch outs))
             ((char=? ch #\<)
              (set! outs (string-append outs (readquote)))
              (rs))
             ((char=? ch #\~)
              (getch)
              (set! outs (string-append outs
               (list-ref args (- (char->integer ch) 48))))
              (rs))
             ((char=? ch #\$)
              (set! outs (string-append outs (readmacrocall)))
              (rs))
             (else
              (set! outs (string-append outs (string ch)))
              (rs))))))
    (rs)))
  (define (readmacrocall)
    (let ((actuals '()))
     (define (rm)
      (let* ((result (readstring))
             (ch (car result)) (s (cdr result)))
       (cond ((char=? ch #\,)
              (set! actuals (cons s actuals)) (rm))
             ((char=? ch #\;)
              (set! actuals (cons s actuals))
              (macrocall (reverse actuals))))))
     (rm)))
  (define (macrocall actuals)
   (cond ((string=? (car actuals) "def")
          (set! env (cons (cdr actuals) env)) "")
         (else
          (let ((def (assoc (car actuals) env)))
           (apply (cons (cadr def) actuals))))))
  (define (apply expr)
   (gpm (car expr) env (cdr expr)))
  (define (gloop)
   (let ((result (readstring)))
    (cond ((string? result) (string-append outs result))
          ((char=? (car result) #\,)
           (set! outs (string-append outs (cdr result) ","))
           (gloop)))))

  (gloop)))


上の例を実行してみよう.
(define str "$def,a,;$a,c;")
(gpm str '() '()) => bcd

(define str "$def,a,;$a,$a,c;;")
(gpm str '() '()) => bbcdd

(define str "$def,a,;$def,bcd,;$$a,c;,e,f;")
(gpm str '() '()) => becfd

(define str "$def,suc,<$1,2,3,4,5,6,7,8,9,10,$def,1,<~>~1;;>;
$suc,4;")
(gpm str '() '()) => 5

(define str "$def,suc,<$1,2,3,4,5,6,7,8,9,10,$def,1,<~>~1;;>;
$def,successor,
 <$~2,$def,~2,~1<,$suc,>~2<;>;$def,9,<$suc,>~1<;,0>;;>;
$successor,3,4;,$successor,2,9;")
(gpm str '() '()) => 3,5,3,0
ついでに, 論文にあったsumは
$sum,α,β,γ;
で2桁の十進数α,βに1桁の数γを足すもので, $3,4,2;は以下のように3,6になる.
(define str "$def,suc,<$1,2,3,4,5,6,7,8,9,10,$def,1,<~>~1;;>;
$def,successor,
 <$~2,$def,~2,~1<,$suc,>~2<;>;$def,9,<$suc,>~1<;,0>;;>;
$def,sum,<$s,~1,~2,0,
 $def,s,<$~3,$def,~3,<$s,>$successor,~1,~2;
  <,>$suc,~3;<;>;$def,>~3<,~1<,>~2;;>;;>;
$sum,3,4,2;")
(gpm str '() '()) => 3,6
これは難しいからマクロ呼出しのトレースをしてみる.
00 (def suc $1,2,3,4,5,6,7,8,9,10,$def,1,<~>~1;;)
01 (def successor $~2,$def,~2,~1<,$suc,>~2<;>;
  $def,9,<$suc,>~1<;,0>;;)
02 (def sum $s,~1,~2,0,$def,s,<$~3,$def,~3,<$s,>
  $successor,~1,~2;<,>$suc,~3;<;>;$def,>~3<,~1<,>~2;;>;;)
03 (sum 3 4 2)
04 (def s $~3,$def,~3,<$s,>$successor,~1,~2;<,>
  $suc,~3;<;>;$def,2,~1<,>~2;;)
05 (s 3 4 0 )
06 (successor 3 4)
07 (def 4 3,$suc,4;)
08 (def 9 $suc,3;,0)
09 (4 )
10 (suc 4)
11 (def 1 ~4)
12 (1 2 3 4 5 6 7 8 9 10 )
13 (suc 0)
14 (def 1 ~0)
15 (1 2 3 4 5 6 7 8 9 10 )
16 (def 0 $s,3,5,1;)
17 (def 2 3,4)
18 (0 )
19 (s 3 5 1)
20 (successor 3 5)
21 (def 5 3,$suc,5;)
22 (def 9 $suc,3;,0)
23 (5 )
24 (suc 5)
25 (def 1 ~5)
26 (1 2 3 4 5 6 7 8 9 10 )
27 (suc 1)
28 (def 1 ~1)
29 (1 2 3 4 5 6 7 8 9 10 )
30 (def 1 $s,3,6,2;)
31 (def 2 3,5)
32 (1 )
33 (s 3 6 2)
34 (successor 3 6)
35 (def 6 3,$suc,6;)
36 (def 9 $suc,3;,0)
37 (6 )
38 (suc 6)
39 (def 1 ~6)
40 (1 2 3 4 5 6 7 8 9 10 )
41 (suc 2)
42 (def 1 ~2)
43 (1 2 3 4 5 6 7 8 9 10 )
44 (def 2 $s,3,7,3;)
45 (def 2 3,6)
46 (2 )
行0,1,2はマクロ定義, 3が呼出しだ. sumを呼び出すと$s,3,4,0とマクロを呼び出す準備をするが ;はまだないので$def,s,...をやってsを定義する.

この際殆どの引数は< >に囲まれているが, 右から15文字目くらいの~3はsumの最後の引数2になる. 4行目で見るとおり. ここでsを呼び出す($s,3,4,0;).

sの本体は$0の呼出しをしようとするが$def,0,があるので0を定義する(16行目). 3,5は$successor,~1,~2;で, 1は$suc,~3;で作られる. それを作るのがトレースの15行目までだ. 次にマクロ2を定義する(17行目).

やっと0が呼び出せ(18行目), Ss,3,5,1;を呼び出す(19行目).

今度はマクロ1と2を同じように定義し, $s,3,6,2;を呼び出す(33行目).

するとマクロ2と2が定義され(44, 45行目), 後から定義された2の本体3,6が返るのである.

万歳. うまくいっているようだ.

2013年11月8日金曜日

Illiacのブートストラップ

アメリカイリノイ州のIllinois大学で1950年頃に作られた計算機がIlliacである.

この計算機は最初のプログラムの入れ方が面白かった. 最近 同僚が神田の古書店で, 電気通信研究所から放出された図書の山から「M1プログラム作製法」を見つけてきたというので, 見せてもらった.

早速プログラムを最初に入れるブートストラップのところを読んだ. 今回はそれをちょっと書いてみたい.

1959年頃の計算機は, 内部は二進法であっても, 機械語のプログラムは十進数で書く. 相対番地も使える. それを読み込み, 二進に直し, 相対番地を絶対番地に直す短いプログラムが用意してあった. イギリスCambridge大学のEDSAC計算機ではそれをイニシアルオーダーといった. アメリカIllinois大学のIlliacではD.O.I.(Decimal Order Input)といった. 東大のPC-1はEDSACがお手本だったのでやはりイニシアルオーダーといった. これらのプログラムは始めから二進法で作ってある.

これらの二進十進変換プログラムの読込み法がそれぞれの計算機でいろいろ工夫された. EDSACではイニシアルオーダーが電話交換機用のロータリースイッチに配線され, スタートボタンを押すとそれが回転して読み込まれる.

PC-1では二進法でパンチした紙テープの読込み方が決められ, それに従ってパンチされた.

Illiacの方法はこれらにくらべてユニークであった. ブートストラップという. つまり靴の紐のことで, 編み上げの靴紐をしばるように順々に完成していくのである.

Illiacは二進法40ビットの計算機である. 命令は1語の左と右に20ビットずつ2つ収める.

20ビットの命令は左8ビットが機能部, 右12ビットが番地部で, 機械語ではそれを十六進法で表す. ただしIlliacの十六進法は10,11,...がA,B,..ではなく, K,S,N,J,F,Lという不思議な文字であった. (king sized number just for loveと覚える.)

ブートストラップの命令は次のようだ.

0 80 028 40 001
1 80 028 40 002
2 19 026 26 000
1 80 028 40 000
0 L4 001 40 001
1 80 028 40 0F6


左端の0,1,2はその右の語が格納される番地である.

ブートストラップが始まるとシーケンスコントロールが0になり, 最初の80028 40 001が0番地と命令対レジスタに入る.

その次の行からがテープにパンチされている. 途中の空白や改行は見易いように入れてあるので, 実際にはない.

命令80はアドレスで示す数の1/4桁をアキュムレータに読めであり, 右命令の40はアキュムレータを番地へ格納せよである. 従って十六進の28は40だから10桁, つまり次の命令80 028 40 002が読み込まれ, 1番地へ格納される.

すると直ぐこの命令が実行されて, 19 026 26 000 が2番地に入る.

命令19はアキュムレータに1/2を置き, それを番地部の数だけ右シフトするので, 1/2は左端の符号ビットの次の位置が1で, 38ビット右シフトするから, アキュムレータの右端に1を置くことになる. 26は番地部の左命令にジャンプ.

そこで0番地の命令対を再び実行し, 80 028 40 000 が1番地に入り, 続いてこれを実行するからL4 001 40 001が0番地に入る.

L4は番地部の数をアキュムレータに加算する命令である.

従って1番地80 028 40 000は80 028 40 001になり, この命令により最後の80 028 40 0F6が1番地に入る. 現時点では

0 L4 001 40 001
1 80 028 40 0F6
2 19 026 26 000
となって2番地を実行しようとしている.

アキュムレータに1を置き0に戻り1番地の右の番地を1殖やして次の10桁を読み格納するというループが廻り始める. だから続く命令はF7, 通常の十六進ではe716=23110から順に格納される.

M1の記憶装置は256語だったから, 231から255まで25語に50命令が入ることになる. Illiacは1024語だったかも知れないから, そうなら1番地右の定数が違っている.

IlliacのD.O.I.はCambridge大学から來ていたDavid Wheelerが作ったといわれる. このブートストラップもWheelerの考案かも知れない.

こういうブートストラップが作れたのは, Illiacに10桁読むという命令があったからで, EDSACやPC-1では1文字読む命令しかないのでこういうことはできない. EDSACそれも記憶場所に読み込むから, それをまたアキュムレータに読み出しで処理することになった. PC-1はどうせそうなるからというので, 読み込む先はアキュムレータであった.

最初のプログラムの読込み方は, その後ミニコン時代にも話題になった. 大方は正面のパネルにスイッチを使って始めの何語かを手でいれる方式であった.

ミニコンのメーカーからその入れ方が指示されるが, それにあきたりないユーザーは別の入力シーケンスを考え, どれだけ短くなったかを競い合った. 最近はそういう楽しみはなくなった.

2013年10月28日月曜日

数学体験館

この10月にオープンした東京理科大学の数学体験館を見た. パラメトロン計算機や微分解析機や沢山の機械式計算器のある近代科学資料館の地下に ある. 直接地下に行ける入口も作られた. 入場無料. 館長は秋山仁さん. 開館時間は12時(土曜は10時)〜16時. 日曜・月曜・祝日は休み.

実はそれを見るのではなく, 近代科学資料館の所用があっていったら, 体験館はまだ30分は開いていると聞き, いそいで見て回ったので, そのうち改めてゆっくり見学したい.

写真を撮影してはいけないというので, ここに示す図はちらしをスキャンしたものである.





ちらしの絵から面白そうなものを拾うと


これは(1/4)+(1/4)2+(1/4)3+...=1/3の証明である.

正三角形が4つの相似な正三角形に分割され, 上の1つが更に4つに分割され,...た図である. 例えば右下の青の正三角形は全体の1/4であり, 上の同じ大きさの正三角形は全体の1/4であって, その右下の青の正三角形はその1/4だから全体の(1/4)2. このようにして青の正三角形は上の式の左辺になり, 青は赤や黄とともに3色で全体になるから, 青の面積は1/3というのである.

しかし展示は図だけなので, 証明は自分で考えなければならない.



これは最小公倍数・最大公約数算出器. 名前だけみると素晴しい計算機だが, 素因数分解は先に出来ていて, 例になっている90は2の球が1個, 3の球が2個, 5の球が1個で並んでおり, 24の方は2の球3個と3の球1個が並んでいる.

左右に分類器があって, 球をいれるとそれが分類されて出てくるから, それぞれの球の最大個数, 最小個数がわかれば公倍数や公約数が得られる.

これはどうもなぁという装置である.



球の体積の公式を導く話だ. 右にあるスイカをサッカーボールのように筋を入れて, 球の中心との母線で錐体に切り出す. 左にあるようにスイカの錐体が沢山できる. 底面は球の一部だがどんどん細い錐体にすると底面はほとんど平面になり, その錐体の体積は底面積に高さ(つまり球の半径r)を掛けた1/3である. 錐体全部の底面積はスイカの表面積だから4πr2. 従ってスイカの体積は4πr3/3である.

表面積を知っている人なら体積も知っていそうなものではあるが.



三平方スライド. 円のなかに直角三角形と各辺を1辺とする正方形が出来ている. 左の図は小と中の正方形に青と黄の板が入っているところ. 円板を回転すると青と黄の板が分割されスライドして大の正方形にうまく収まることをしめす.

昔ボストンのサイエンスミュージアムに同じ趣向だが水が入っているのを見たことがあった.



二項分布パチンコ. あまりよくは見えぬが右上が坂の上になっていて, パチンコの玉が沢山入っている. 出口を開けると一斉に左下に流れ落ちるが, 途中に通路をランダムにする釘が打ってあり, それによってパチンコ玉は下の区画のどれかに落ちる. その形が二項分布に見えるという実験である.

なにかの本でも見たような気がする.

まぁこういう実験道具が沢山おいてあって, 数学がこういうふうに体験できると分かって面白い.

最初のちらしのフロアーマップにあるように, 立派な数学工房がある. 体験館の実験具の多くはどうもそこで試作されたらしい. お願いすると使わせても貰えるらしい. なにか作ってみるものはないだろうか.

計算機屋からみると, ユークリッドの互除法, エラトステネスのふるい, ハノイの塔のようなアルゴリズム関係の展示もあればいいのにと思うが, すでに展示場は一杯である.

2013年10月17日木曜日

微分解析機

前回の微分解析機のシミュレータの続きである.

実物の微分解析機には出力装置があり, 2つの変数軸をxとyに入れると解の図形が得られる.

MIT Schemeにあるscheme graphicsの最も簡単な使い方は
(define mydevice (make-graphics-device 'x))
で描画の装置を宣言する. 装置の描く場所の正方形は左右がx軸, 上下がy軸であり, 左端はx=-1, 右端はx=1, 下端はy=-1, 上端はy=1である. この座標系で(x0,y0)から(x1,y1) まで直線を引くには
(graphics-draw-line device x0 y0 x1 y1)
とする.

そこで前回のブログのsinz, coszを描くには
(define (draw device step x0 y0 x1 y1)
 (if (> step 0) (begin
  (graphics-draw-line device x0 y0
   (stream-car x1) (stream-car y1))
  (draw device (- step 1) (stream-car x1) (stream-car y1) 
    (stream-cdr x1) (stream-cdr y1)))))
を使う. (draw ..)が描画の関数で, deviceに対してstep回描画する. x0, y0は解の曲線の始点; x1, y1はstreamである.

一方 駆動する方は, 座標を描いた後
(for-each (lambda (y) (graphics-draw-line mydevice -1 y 1 y))
  '(-0.75 0 0.75))
(for-each (lambda (x) (graphics-draw-line mydevice x -1 x 1))
  '(-0.75 0 0.75))

(draw mydevice 6282 0 0.75
 (scale-stream sinz 0.75) (scale-stream cosz 0.75))
のようになっている. stepは積分がzの1/1000で進行し, 2πで一周するから6282ステップ. 次の0と0.75はxとyの初期値(つまり鉛筆を置く位置)で, 1だと画面すれすれになるので0.75倍にしてみた.

次のsinzとcoszが出力軸の変数で, こちらも0.75にギアダウンしている. そういうわけで, この呼び方も出力装置をシミュレートしているといえるであろう.

こうして描いたサークルテストの図が次である.



円の始点と終点の真上でちょっと食い違っているようにも見えるが....

次はezを描く. 描きたい範囲を-1と1の間に変換しなければならない. 下の図を見てほしい. 左と上が実際の座標で右と下がGraphicsの座標である. これによるとx軸は0.4倍すればよい. y軸の方は1が-0.2になっているから0.4倍してから0.6を引く.



従ってプログラムは次のようだ. 始めのfor-eachは座標を描く.

(for-each (lambda (y) (graphics-draw-line mydevice -1 y 1 y))
  '(-0.6 -0.2 0.2 0.6))
(for-each (lambda (x) (graphics-draw-line mydevice x -1 x 1))
  '(-0.8 -0.4 0 0.4 0.8))
 (scale-stream z 0.4) 
  (add-streams (scale-stream ez 0.4)
    (scale-stream ones -0.6)))
回数の1386は丁度上の端に相当するxの値がlog 4=1.386であることによる.



こういう図を描くと, e-zで左半分も描きたくなる.
(define e-z (integrator (delay (scale-stream e-z -1))
  (delay z)
 1.0))
を用意し,
(draw mydevice 2500 0 -0.2
 (scale-stream z -0.4) 
  (add-streams (scale-stream e-z 0.4) 
   (scale-stream ones -0.6)))
を追加すれば, 左も描ける.



2013年10月16日水曜日

微分解析機

8月30日のブログに「これらのシミュレーションにはSchemeのstream処理が適していると思うが, まだ手が付かずにいる.」と書いた.

やっとどうやらシミュレーションが出来るようになったので, 今回はその話題である.

SICPの206ページのintegral関数のように披積分関数を遅延にするのが味噌である. SICPの例題では独立変数がdtだけだが, 微分解析機では「回す」値になにが入るかわからないので, こちらもstreamに対応させなければならない.

独立変数に方は回転角の差が必要になるので, stream-cadrからstream-carを引いたものに, 披積分関数のstream-carを掛けることになる.

まずintegratorは次のようだ.
(define (integrator delayed-integrand delayed-variable 
 initial-value)
 (define (mult-streams a b)
  (cons-stream (* (stream-car a)
            (- (stream-car (stream-cdr b)) (stream-car b)))
    (mult-streams (stream-cdr a) (stream-cdr b))))
 (define int
  (cons-stream initial-value
   (let ((integrand (force delayed-integrand))
         (variable (force delayed-variable)))
   (add-stream (mult-streams integrand variable) int))))
int)
integratorが貰う引数は遅延になっており, intの計算のなかでforceを使って強制する.

次はシミュレーションに必要な補助関数である. これらはSICPにあるものだ.
(define (scale-stream stream factor)
  (stream-map (lambda (x) (* x factor)) stream))
(define (add-streams s1 s2)
 (stream-map + s1 s2))
(define (integers-starting-from n)
 (cons-stream n (integers-starting-from (+ n 1))))
(define integers (integers-starting-from 1))
さてこれらが用意できたら, 8月14日のブログにあった関数からやってみる. ブログの図を見ながらプログラムを眺めてほしい.

最初はsin zとcos z. 独立変数はzである. zは整数のstreamで, 1/1000にして使う. sin zのπ/4, つまり45度での値をstream-refで見ることにする.
(define z (scale-stream integers 0.001))
(define sinz (integrator (delay cosz) (delay z) 0))
(define cosz (integrator (delay (scale-stream sinz -1))
  (delay z) 1.0))
(display (stream-ref sinz 785)) ;=>.7071024791302126
値としてはよさそうである.

次はez.
(define z (scale-stream integers 0.001))
(define ez (integrator (delay ez) (delay z) 1.0))
(stream-ref ez 1000) ;=>2.716923932235898
z2とz3は次のようにする.
(define z (scale-stream integers 0.001))
(define izdz (integrator (delay z) (delay z) 0))
(define z2 (scale-stream izdz 2))
(define iz2dz (integrator (delay z2) (delay z) 0))
(define z3 (scale-stream iz2dz 3))
(stream-ref z3 500) ;=>.12499950000000003
(stream-ref z2 500) ;=>.25049999999999994
0.5の3乗が0.125, 2乗が0.25ならまぁ問題はない.

tan zは多少手強いが,
(define ones (cons-stream 1 ones))
(define z (scale-stream integers 0.001))
(define tanz (integrator (delay 1+tanzsq) (delay z) 0))
(define itanzdtanz (integrator (delay tanz) (delay tanz) 0))
(define tanzsq (scale-stream itanzdtanz 2))
(define 1+tanzsq (add-streams tanzsq ones))
(stream-ref tanz 785) ;=>.9979528633948679
(stream-ref tanz (quotient 3142 6)) ;=>.5761873674775257
テストしたのはtan π/4とtan π/6で, 1/√ 3=.5773502691896258だからこれも合格である.

1/zとlogezはzを1から積分する. z=1での初期値はどちらも0である.
(define z (scale-stream integers 0.001))
(define logez (integrator (delay 1/z) (delay z) 0))
(define 1/z (integrator (delay (scale-stream 1/z -1)) 
   (delay logez) 1))
(stream-ref 1/z 500) ;=>.6664863143078377
(stream-ref logez 1718) ;=>.999948126092027
1/zの500での値は1/1.5なので0.6666.... logezの1718はloge2.718だから1でいいわけである. zの値は1からのはずだが, 回転角の差を使うから, integersで構わない.

微分解析機には結果を描画する出力装置もついているから, 次はscheme-graphicsを駆動するようなプログラムにしたい.

2013年9月25日水曜日

微分解析機

微分解析機でバックラッシュを減らす仕掛けには, フロントラッシュ以外にも, ラッシュロック(lashlock)なるものがあった.

微分解析機では独立変数軸で回転されるディスクは, 被積分関数軸の値によって前後に動かされる. つまりディスクの載っている台座が, 被積分関数軸に繋がる送りネジ(lead screw)に嵌めたナットと一緒になっていて, 移動するわけだ.

被積分関数の値は結果の積分値に大きく影響するので, この送りネジの工作精度は重要で, ケンブリッジの微分解析機では, 送りネジの作製に微分解析機全体の1/10のコストがかかったという.

ここでもバックラッシュを減らす方法が検討されている. 東京理科大学に保存されている微分解析機にもラッシュロックがちゃんとあった.

下の図でハッチのあるのが送りネジである. それに左の大きいナットと右の小さいナットが嵌めてあり, 大きいナットの上(破線の上)に台座が固定されている. だから送りネジはこのナットを経て台座を動かす.



2つのナットの間には, 太い線で示すバネが挟んであり, 兩ナットを離そうとしている.

小ナットには縁に何か所が切れ目があり, 大ナットから突出した棒の先の爪が差し込めて, 小ナットの回転を防いでいる.

つまり送りネジのネジ山を両端に押してバックラッシュをなくそうとしているのである. なるほどすごい仕掛けだ.

Crank本によると, この仕掛けもトルクアンプと同じく, ベツレヘムスチール社のNiemanが発明したそうだ.

2013年9月11日水曜日

微分解析機

微分解析機ではバスシャフト, クロスシャフトに積分機, トルクアンプ, 入力卓, 出力卓, 加算器などが接続してある. しかし佐々木本によると, 「ガタ補正装置」というものもあったらしい.

これは一体なにかと思う一方, 理科大の微分解析機になにやら不思議な機器がついていた. Crank本で判明したのだが, その不思議なものこそガタ補正装置であった. 英語では「フロントラッシュ」という.

微分解析機の最大の泣き所はトルクアンプでのバックラッシュであったらしい.

バックラッシュとは, 通常の平歯車ではこう起きる. 次の図のように2つの平歯車MとNがかみ合っていて, Mが駆動する側, Nが駆動される側とする.



Mが右回転していると, MとNのかみ合う部分では, Mの歯の下がNの歯を上から押し下げてNを左回転させる. Mの歯の下とNの歯の上は常に接しているから連続的に駆動できる.

Mが停止し逆回転を始めると, 今度はMの歯の上がNの歯の下を押し上げることになるが, こちらには隙間があって, MがNに当たるまでNは逆回転を始めずに停止している. 従ってMとNの回転角の関係を図に描くと, 回転が逆になると違う関係になる. これがバックラッシュである.

トルクアンプでも逆回転が始まると, 逆側のドラムの回りの糸が締まるまでは出力が始まらず, しばらく不感時間があってバックラッシュ現象が生じる.

このバックラッシュを出来るだけ早くキャンセルしようとするのがフロントラッシュである. 名前からして結構ふざけているがなかなかの機構である.

次がフロントラッシュのメカの分解図である. Cambridge大学のWilkes先生のところで撮った写真だと断りがあった.



フロントラッシュはバスシャフトの間に入れ, バスシャフトの回転数を加速する機能を持つ.

右が入力シャフトで, その左に両隣のシャフトへの引掛りをもったベルトに挟まれたドラムがある. ベルトはドラムが勝手に回らないよう軽くブレーキをかけている. ドラムに入力シャフトが通っているが, 入力シャフトとドラムはシャフトの回りではフリーになっている.

ドラムの左は遊星歯車. そして左端に内歯車とそれと一体になった出力シャフトがある.

ドラムには右方にペグが2本. この挿入位置は円周上に点在し, 場所は変更し調整出来る. 一方入力シャフトにはペグにぶつかるピンがある.

通常は入力シャフトがピンでドラムのペグを押し, 同時回転し, 従って遊星歯車も同時回転し, 入力シャフトと出力シャフトは同時回転する. しかし入力シャフトが逆転すると, ピンがペグから離れ, もう一方のペグに当るまで, ドラムはバンドの摩擦で停止する.

ドラムと入力シャフトの回転速度に差があると, 遊星歯車の出力は加速され, 出力シャフトも少し早く回り, やがてピンがもう一方のペグに当たって通常の回転数に戻るという仕掛けである.

この辺の事情を図にしたのが下だ.



左の図で横軸が入力, 縦軸が出力である. DからAは入力と出力が一体となってある方向に回転しているところとする. Aで入力が逆転し始めるても, 出力はすぐには追従しないからBまでは出力は動かぬ. BからCで出力は入力にすこし遅れて逆転する.

右の図はフロントラッシュの場合で, 横軸は入力シャフトのまだ先にある真の回転角である. Bまでは逆転直後で入力軸が遅れているところ. やがて遊星歯車が回転して出力は真の回転角に追いつく様子を示す. 逆側のC, Dでも事情は同じである.

なんとしてもバックラッシュを防ごうとう努力の後が見えるではないか.

2013年8月30日金曜日

微分解析機

初等関数を積分機を使って作りだす方法で, 前回のブログの続きである.

tan z

tan zはπ/2で無限大になるから要注意

d tan z/dz =1+tan2だから
tanz=∫ 1+tan2z

tan2のつもりのtanzsqなる変数を用意し

tanz = ∫ (1 + tanzsq) dz
tanzsq = 2∫ tanz d tanz

tanzは1+tanzsqで「動かし」, zで「回す」.
tanzsqはtanzで「動かし」, tanzで「回し」2倍する.

図は次のようだ.



3本のバスシャフトにまたがる[・Σ・]のような記号は加算器で, ・のバスシャフトの値の和がΣのバスシャフトに出力される.

1/z, logez

logez = ∫ 1/z dz
1/z = -∫ 1/z d logez だから

logezは1/zで「動かし」, zで「回す」.
1/zは-1/zで「動かし」, logezで「回す」.

図は次に示す.



これらのシミュレーションにはSchemeのstream処理が適していると思うが, まだ手が付かずにいる.

2013年8月14日水曜日

微分解析機

微分解析機のことを調べようとすると, 参考になる文献は
城憲三「計算機械」共立全書
佐々木達治郎「計算機械」河出書房
くらいであった. 微分解析機を使われたり, 設計されたりした渡辺勝先生に先頃お目にかかったら, John CrankのThe Differntial Analyser(1947年刊)を読むべきだといわれた.

私が在職していた東大工学部計数工学科の図書室にあったので, さっそく読み始めた. まことに実用的な本であった. 英国のCambridge大学にあった微分解析機に基いて豊富な説明があった.



そこの微分解析機も, 積分機の円板の方を被積分関数に応じて動かす方式なので, 独立変数の入力は円板を`rotate'する(回す); 被積分関数の入力は積分機を`displace'する(動かす)と使い分けていて実に理解しやすい.

従って積分機の図も次のようで, 積分機を押し上るような矢印が被積分関数, 円板の横の一部網かけのドライブの軸が独立変数, 円板中央のT字状の軸が積分出力である. それらに繋がる縦の軸がクロスシャフト. 横の軸がバスシャフトである.



積分機を2台使ってsin, cosを発生させるのは, どこにでも例題があるが, もっと他の基本関数の発生法もなかなか面白かったので, 今回はそれが話題である.

本番で使うときは, 係数処理をしなければならないが, とりあえずどの関数は積分機何台で出来るかという一覧表があった. その中から:

sin z, cos z

d sin z/dz = cos z
d cos z/dz = -sin zだから
sin z = ∫ cos z dz,
cos z = ∫ -sin z dz

cos zで「動かし」, zで「回す」と sin zが得られ,
-sin zで「動かし」, zで「回す」と cos zが得られる.

図は上に示した通り. バスシャフトとクロスシャフトの交点が丸で囲んであるのは, クロスシャフトの回転を通常と反対にすることを示す.

z=0からπまでをプログラムでシミュレートしてみると(積分区間を1024分割した):

(define k 1024)
(define z 0) (define sinz 0) (define cosz 1) 
(define pi (* 4 (atan 1)))
(define dz (/ pi k))
(do ((i 0 (+ i 1))) ((> i k))
(set! z (+ z dz))
(set! sinz (+ sinz (* cosz dz)))
(set! cosz (- cosz (* sinz dz))))

(display (list z sinz cosz)) =>
(3.1423596439839487 -7.670673988518897e-4 -.9999994116371386)

ez

d ez/dz = ezだから
ez = ∫ezdz

ezで「動かし」, zで「回す」と, ezが得られる.

次の図のようだ.



z=0から1までをプログラムでシミュレートしてみると(積分区間を1024分割した):

(define k 1024)
(define z 0)
(define ez 1)
(define dz (exact->inexact (/ 1 k)))

(do ((i 0 (+ i 1))) ((> i k))
(set! z (+ z dz))
(set! ez (+ ez (* ez dz))))

z => 1.0009765625
ez => 2.719609006545993

z2, z3

z2 = 2 ∫z dz
z3 = 3 ∫z2 dz

zで「動かし」, zで「回し」, 2倍するとz2が得られ,
z2で「動かし」, zで「回し」, 3倍するとz3が得られる.



バスシャフト間を繋ぐ[1 2]や[1 3]は, 回転をそれぞれ2倍, 3倍する歯車を示す. z3はこのように積分機2台が必要である.

シミュレーションは不要であろう.

他の関数は次の機会にしよう.

この本を読んだ収穫のひとつが独立変数は「回す」, 被積分関数は「動かす」と言い分けることであった. 若き日のWilkes先生を含め, Cambridge大学の微分解析機のグループが, rotate, displaceを使うと意思の疎通が捗ると知り作業していたらしいことが想像できる.

大学院生のころよく読んだ海外誌がBSTJ(Bell System Technical Journal)であった. C.E.Shannonの通信理論の論文もBSTJに掲載された. BSTJの電話交換機の論文では, Sheといえば交換手, Heといえば加入者であって, 便利な使い分けだと感心したが, こういうこともコミュニケーションで大事なノウハウであろう.

2013年7月30日火曜日

面積計を使う調和解析器

Henriciの考案した調和解析器は, 前回のブログの関数と距離計の読みの図の少し上にあるように, 下の距離計は ∫ cos(θ)dy ∝ B1, 左の距離計は ∫ sin(θ)dy ∝ -A1を計測していた.

一般のnについては,

nAn=-1/π∫sin(nθ)dy,
nBn=1/π∫cos(nθ)dy,

を計測する. これは普通に書いてある

An=1/π∫0y cos(nθ)d&theta,
Bn=1/π∫0y sin(nθ)d&theta,

から次のようにして得られる. また式が多いのでTexで書くことにする.


積分の範囲は前回の図に示すように, 一周期分である.

2013年7月27日土曜日

面積計を使う調和解析器

厳密にいうとこれは面積計は使わないが, Olaus Henriciが考えた似たような仕掛けの調和解析器の話だ. 下の図は, Henriciの考案したものを, チューリッヒの計測器メーカーCoradiが1890年ころに作ったものである.



このどっしりした機械を, 例えば前回のブログにあったような, 解析したい図形の上に置く. 枠の横に長い方をx軸に平行にする. この枠は3個の車輪で支えられている. 1個は手前中央に見える幅広のD. 残りの2個は枠の内側の両端でE. これらの車軸がx軸と平行なので, この枠はy方向にしか動かない.

枠を前後に動かすと, Eが回転し, 2個のEを繋ぐシャフトも回転し, 中央に上半分が見える車輪Cもy方向の移動距離に比例した分回転する.

車輪Cの真上に上下方向の回転軸があり, その下に枠に収まった球のようなものが見える. この球は下の車輪に接していて, やはりx軸と平行な軸で回転する. その球を支えている枠は, 上下方向の回転軸と同時に球の周囲を回転出来る.

回転させるのは, 回転軸の上に見えるプーリーHとそれに接するワイヤーである.

大きい枠の左には, ワイヤーがついていて左右に動きそうな枠Wが見える. この枠の手前の下にFなるポインターがあり, このポインターで関数の曲線を左から右に追跡するのである.

関数曲線の上下で枠全体が前後に動きながら, ポインターで右方向へ追跡すると, ワイヤーが動き, 従ってプーリーも回転して, 球を支える枠も回転する仕掛けである.

さて下の図が球を支える枠の構造だ.



中央のSが球で, yの矢印は枠がy方向に動くときに球が回転する方向を示す. xはポインターが右へ, つまりx方向へ移動すると枠が回転する方向を示す.

球にはこの図で見るように, 下と左に球の回転量を計る距離計のようなものが接している(この辺がいちおう面積計っぽいところだ).

この図のような位置で球が前後方向へ回転すると, 下の距離計は積算するが, 左のは動かない.

しかしプーリーの回転に従って枠も時計廻りにθだけ回転すると, 下の距離計はyの回転量のcos(θ)倍を計測し, 左のはsin(θ)倍を計測するようになる.

この枠はxが0から2πまで移動するときに丁度1回転するようになっている.

したがって下の距離計の積算は∫cos(θ)dy, 左の距離計は∫sin(θ)dy になり, それぞれB1, -A1に関係する値が得られる.

プーリーの直径を変えると, 枠の回転は2倍にも3倍にもなり, ほかの係数も求めることができる.

なんか不思議な機構だが, プログラムでシミュレートしてみた. お馴染のf=A1 cos(θ)+A2 cos(2θ)+B1 sin(θ)+B2 sin(2θ) である.

下の図は黒線がfを表す. 枠を2πで1回転したとき, つまりn=1の時の下の距離計の値をB1(青), 左のをA1(赤)で示す. またn=2の時, 下のをB2(緑), 左のをA2(橙)で示す.



単一の三角関数の場合もやってみると下のようになる. f=A1cos(θ), A1=1 がこれだ.



f=B2 sin(2θ), B2=1も示す.



振幅0に対応する距離計は最後に0になり, それ以外の距離計は係数に従った値が読み取れることが分る.

ソナグラフの存在しないころ, この解析器で音響の周波数解析をした人がいたらしい.

2013年7月26日金曜日

面積計を使う調和解析器

前回のブログ(2013年7月18日)の続きだ. あの奇妙な形の面積で係数が得られる理由はこうだ.

関数とA1の絵をもう一度見てみよう.



上の図のPの座標は(θ, f(θ))である. これを円柱に巻きつけて, zero-edgeが中央にくるように見たのが下の図である.

下の図で, 枠の中心を原点としたPの座標を(x,y)とすると,
x=sin(θ)
y=f(θ)
である.

x=g(θ), y=f(θ)でθ=0から2πで閉曲線になるとき, 閉曲線で囲まれた図形の面積は
S=∫f(θ)g'(θ)dθ
だから, A1の面積は(fをθ, 2θだけの関数として)
S=∫(A1 cos(θ)+A2 cos(2θ)+B1 sin(θ)+B2 sin(2θ)) cos (θ)dθ
で A1∫cos(θ)cos(θ)dθ=A1π 以外は消えるからこの面積の1/πが A1 になる.

こんな説明でいいだろうか. ところで, 積分が怪しいとか忘れたとかのときは, WolframAlphaにお伺いするとよい.

integrate cos(x) cos(x) dx from 0 to 2 pi => π
integrate cos (x) cos (2x) dx from 0 to 2 pi => 0
integrate cos (x) sin (x) dx from 0 to 2 pi => 0

など. 便利である.

2013年7月18日木曜日

面積計を使う調和解析器

世の中には面白いことに気付く人がいるものだ.

下の図は適当に描いた周期関数である. 関数は赤, 緑, 青, 橙の各点を通り, 赤に戻る. 左端の赤い縦線はθ=0のところで, zero-edgeという.



これが透明な紙に描いてあるとし, 幅が丁度2πだから, 半径1の, これも透明な円柱に巻き付けたとする. AがA'に, BがB'に出会う. (下の図の上の部分は左がまだ巻く前, 右が巻いたところ.)



巻いた円柱は, zero-edgeが上に來るように寝かせる. それが右だ. 4つの点は先ほどと同じで, 赤は真上にあり, 曲線は右へ進んで緑の点で縁に来る. そこから青を経て橙の点までは円柱の下の面を進み, 橙からは上の面を通ってzero-edgeに戻る.

下の部分は, 上の関数の図を横に2倍に延ばし, 右は前と同じ径の円柱にぐるぐると2度巻き, zero-edgeが上になるように置いたものである. 赤から右向きに出発, 裏の中央に緑が来る. 再び表側に来て青の点を通り, また裏で橙を通って赤に戻ることを示す.

このようにして出来た透明な円柱を, 上や左から見たのが次の図である.



A1は1回巻きを上から見たもの(zero-edgeが中央にある), B1は1回巻きを左から見たもの(zero-edgeが右端にある), A2, B2は2回巻きをそれぞれ上と左から見たものである.

そこで面積計を取り出し, 各々の閉曲線の面積を計測する. 面積計算のプログラムは時計回りに追跡した時に正の値を返す. A1, A2は赤から緑へ行くのが時計回りで, ほぼ時計回りだから, 面積は正になり, B1, B2は赤から緑へ行くのが反時計回りだから, 面積は負になる.

これらの図を描くプログラムで, 途中の点の座標を記録しておいて, 面積を計算すると,
A1: 6301.53662
A2: 12601.5176
B1: -6365.42822
B2: -12728.2783

これらの名前から想像できるように, この面積は最初の周期関数のフーリエ係数の何倍かになっているのである. (最初の曲線はA1=A2=B1=B2=0.5だった.)

もっと先の方の係数An, Bnを求めるには, 元の関数の図を横にn倍に延ばし, 円柱にn回巻きにし, 上や左からみた面積を計測すればよい.

こんなややこしい図はやめて, 単一の三角関数で実験してみよう.

Cos θの場合. A1=1(他は0)でテスト


A1: 12668.2227
A2: -127.295372 ≈ 0
B1: 0.499938339 ≈ 0
B2: 1.99811935 ≈ 0

A1は円柱の半径を半径とする円である. これらの図は2πを400(point)にとって描いたので, 円柱の半径wは200/π.
正確な値は(200/π)2×π=12732.3954

Sin 2θの場合. B2=1(他は0)でテスト


A1: -0.99832052 ≈ 0
A2: -1.99973631 ≈ 0
B1: 0.00795254577 ≈ 0
B2: -25460.5547

このB2と前のA1で, 図の形は同じ円なのに, 面積が違うのは, B2の方は円周を2度回っているので, 面積が2倍になったのである

似たようなテストをCos 2θ, Sin θでもやって, 円の面積を計算して置き(Sin 2θ, Cos θの値と符号が違うだけだが), 最初の関数のテストで得られた面積を割って見ると

(define a1 6301.53662) ;最初の関数の面積
(define a2 12601.5176)
(define b1 -6365.42822)
(define b2 -12728.2783)

(define cosa1 12668.2227) ;単一関数の面積
(define cosa2 25333.3145)
(define sinb1 -12731.8721)
(define sinb2 -25460.5547)

(/ a1 cosa1) => .49742862666915383 ;最初の関数の係数
(/ b1 sinb1) => .4999601134855886
(/ a2 cosa2) => .4974286961147543
(/ b2 sinb2) => .49992148442861695
となって, 元の係数が全部0.5であったことが判明する.

ところで, このブログでは, 平面図や側面図を描いたので, 面積計で計測出来たわけだが, 円柱に関数の図を巻き付けただけでは, 面積計は使えないのではという疑問は当然だ. 実はそういうことで, アイディアはいいのだが, 装置としては実現されなかった.

Cliffordという人が言い出したそうで, Proceedings of the London Mathematical Societyのvol.v(1890年頃)に載っているらしい.

2013年7月13日土曜日

高速フーリエ変換

先日, 調和解析の器械のことを調べていて, フーリエ変換をやってみる必要があった.

データの個数の少い例なので, プログラムは簡単に書ける. 書きなが, 学生のころ(1953年ころ)に計算させられたことを思い出した. あの時に比べるといまは極楽だな.

ところで高速フーリエ変換(FFT)というのがあって, 私は30年も前に, 高橋秀俊先生の追悼文を「コンピュータソフトウェア」誌に寄稿したときに, 高橋先生がFFTのプログラミングに熱中されていたことにも触れた. (この記事はCiNiiで探すと見つかる.)

当時の雑誌を取り出してみると, そのプログラムが掲載されている. もとは東大大型計算機センターのライブラリプログラムだからFortranで書かれれていたのを, 私が当時使っていたFranz Lispに書き直したものであった. 添字の名前などにはいかにもFortranという気分が残っているが.

いまさらFranz Lispでもあるまいと, Schemeに書きあらためた. 面白いプログラムなので, ちょっとその説明を書いてみたい.

最初からそのプログラムである.

(define (fft n)
 (let* ((n1 (/ n 2))
        (n11 (/ n1 2))
        (s (make-vector n11))
        (x (make-vector n))
        (a (make-vector n1))
        (b (make-vector n1)))

  (define (concat x i)  ;記号sと1をもらい, 記号s1を作る.
   (string->symbol
    (string-append (symbol->string x) (number->string i))))
  (define (inits n11)   ;配列sの初期化
   (do ((i 0 (+ i 1))) ((= i n11))
    (vector-set! s i (concat 's i))))
  (define (initx n)     ;配列xの初期化
   (do ((i 0 (+ i 1))) ((= i n))
    (vector-set! x i (concat 'x i))))
  (define (restorex)    ;a,bをxに戻す
   (do ((i 0 (+ i 1))) ((= i n1))
    (vector-set! x i (vector-ref a i))
    (vector-set! x (+ i n1) (vector-ref b i))))
  (define (printx)      ;配列xを出力
   (newline) (do ((i 0 (+ i 1))) ((= i n))
   (display i) (display " ")
   (display (vector-ref x i)) (newline)))

  (inits n11)    ;sの初期化
  (initx n)      ;xの初期化
  (printx)
  (do ((n2 n (/ n2 2))) ((= n2 1))
   (let ((n21 (/ n2 2)))
    (do ((i 0 (+ i 1))) ((= i n21))
     (let ((i1 (+ i n21)))
      (vector-set! a i
       (list (vector-ref x i) '+ (vector-ref x i1)))
      (vector-set! b i
       (list (vector-ref x i) '- (vector-ref x i1)))
      (cond ((< (+ i n21) n1)
       (let ((i2 (+ i n11)))
        (vector-set! a i2 (vector-ref x (+ i n1)))
        (vector-set! b i2 (vector-ref x (+ i1 n1)))
        (let ((i3 i) (i6 (+ i n1)))
         (do ((j n21 (+ j n21))) ((= j n11))
          (set! i3 (+ i3 n2))
          (set! i6 (- i6 n21))
          (let* ((i31 (+ i3 n1)) (i4 (+ i3 n21))
           (i41 (+ i4 n1)) (i5 (+ i j)) (j1 (- n11 j))
           (ap (list
            (list (vector-ref x i4) '* (vector-ref s j1))
           '- (list (vector-ref x i41) '* (vector-ref s j))))
           (bp (list
            (list (vector-ref x i41) '* (vector-ref s j1))
           '+ (list (vector-ref x i4) '* (vector-ref s j)))))
           (vector-set! a i5 (list (vector-ref x i3) '+ ap))
           (vector-set! b i5 (list (vector-ref x i31) '+ bp))
           (vector-set! a i6 (list (vector-ref x i3) '- ap))
           (vector-set! b i6
            (list '- (vector-ref x i31) '+ bp))))))))))
     (restorex) (printx)))))

(fft 16)    ;n=16で起動

もとのライブラリの説明には
「大きさnのデータX0,X1,...,Xn-1(nは2の巾乗)のFourier変換. Am=∑i=0n-1Xicos(2πmi/n) (m=0,1,...,n/2), Bm=∑i=0n-1Xi sin(2πmi/n) (m=1,2,...,n/2-1)をすべてのmについて計算する. cosineの最後の変換はプログラム上 sineの最初の変換の場所に入っている. すなわちB0は常にゼロのために, B0の場所には求めたAn/2が入ってくる.」
とある. このような記述も懐しい. なおプログラムの作成は1966年8月22日だ.

sとxの初期化の後のdoループで, n2をn,n/2,n/4,...とlog2n回まわし, その次のdoループでiを0,1,2,...と回すからn log nのオーダーのアルゴリズムなのが分かる

このプログラムはフーリエ係数を計算するのではなく, なにを計算しているかを示すのが目的であった.

だから出力はこんな具合だ. 1回目と2回目のxの内容を示す.
0 x0    (x0 + x8)
1 x1    (x1 + x9)
2 x2    (x2 + x10)
3 x3    (x3 + x11)
4 x4    (x4 + x12)
5 x5    (x5 + x13)
6 x6    (x6 + x14)
7 x7    (x7 + x15)
8 x8    (x0 - x8)
9 x9    (x1 - x9)
10 x10   (x2 - x10)
11 x11   (x3 - x11)
12 x12   (x4 - x12)
13 x13   (x5 - x13)
14 x14   (x6 - x14)
15 x15   (x7 - x15)
最初はxにx0, x1,...,x15のような文字が入っている. 1回目の処理でxは(x0 + x8)のように変る. 計算用のプログラムでなら,
(vector-set! a1 i (list (vector-ref a i) '+ (vector-ref a i1)))
(vector-set! a1 i (+ (vector-ref a i) (vector-ref a i1)))
と直して, x0+x8の値にするところである.

xの配列はさらに以下のようになる. s1, s2, s3はsin(π/8), sin(π/4), sin(3π/8)である. sinの表は第1象限でしか保存しない.

0 (((x0 + x8) + (x4 + x12)) + ((x2 + x10) + (x6 + x14)))
1 (((x1 + x9) + (x5 + x13)) + ((x3 + x11) + (x7 + x15)))
2 ((x0 - x8) + (((x2 - x10) * s2) - ((x6 - x14) * s2)))
3 ((x1 - x9) + (((x3 - x11) * s2) - ((x7 - x15) * s2)))
4 ((x0 + x8) - (x4 + x12))
5 ((x1 + x9) - (x5 + x13))
6 ((x0 - x8) - (((x2 - x10) * s2) - ((x6 - x14) * s2)))
7 ((x1 - x9) - (((x3 - x11) * s2) - ((x7 - x15) * s2)))
8 (((x0 + x8) + (x4 + x12)) - ((x2 + x10) + (x6 + x14)))
9 (((x1 + x9) + (x5 + x13)) - ((x3 + x11) + (x7 + x15)))
10 ((x4 - x12) + (((x6 - x14) * s2) + ((x2 - x10) * s2)))
11 ((x5 - x13) + (((x7 - x15) * s2) + ((x3 - x11) * s2)))
12 ((x2 + x10) - (x6 + x14))
13 ((x3 + x11) - (x7 + x15))
14 (- (x4 - x12) + (((x6 - x14) * s2) + ((x2 - x10) * s2)))
15 (- (x5 - x13) + (((x7 - x15) * s2) + ((x3 - x11) * s2)))
n=16の場合は次が最後の表示である. 0行目は係数のA0, 1行目はA1, 2行目はA2などである.

A0は(x0+x1+...+x15)/8だが, このプログラムは1/nや2/nの計算はさぼっているので, その辺りは注意が必要だ. 要するにxにcosやsinを掛けて足すところまでにしか関心を持たない.
0 ((((x0 + x8) + (x4 + x12)) + ((x2 + x10) + (x6 + x14))) +
   (((x1 + x9) + (x5 + x13)) + ((x3 + x11) + (x7 + x15))))
1 (((x0 - x8) + (((x2 - x10) * s2) - ((x6 - x14) * s2))) +
   ((((x1 - x9) + (((x3 - x11) * s2) - ((x7 - x15) * s2)))
    * s3) -
    (((x5 - x13) + (((x7 - x15) * s2) + ((x3 - x11) * s2)))
    * s1)))
2 (((x0 + x8) - (x4 + x12)) +
   ((((x1 + x9) - (x5 + x13)) * s2) -
    (((x3 + x11) - (x7 + x15)) * s2)))
3 (((x0 - x8) - (((x2 - x10) * s2) - ((x6 - x14) * s2))) +
   ((((x1 - x9) - (((x3 - x11) * s2) - ((x7 - x15) * s2)))
    * s1) -
   ((- (x5 - x13) + (((x7 - x15) * s2) + ((x3 - x11) * s2)))
    * s3)))
そう思って見ると0は確かにx0からx15までの和である.

2も割り合いと易しい. xの列にcosの曲線を思いだし, cos(π/4)=sin(π/4)=s2に注意しながら, 下のようにs2を掛けるわけだが,
x0 x1 x2  x3 x4  x5 x6 x7 x8 x9 xa  xb xc  xd xe xf
 1 s2  0 -s2 -1 -s2  0 s2  1 s2  0 -s2 -1 -s2  0 s2
たしかにそう出来ている.

1についてもやってみる.

x0 x1 x2 x3 x4  x5  x6  x7 x8  x9  xa  xb xc xd xe xf
1  s3 s2 s1  0 -s1 -s2 -s3 -1 -s3 -s2 -s1  0 s1 s2 s3
一方式の方は, 1行目の(x0 - x8)はOK.

(((x2 - x10) * s2) - ((x6 - x14) * s2))はs2の相方でこれもOK.
2行目の
(((x1 - x9) + (((x3 - x11) * s2) - ((x7 - x15) * s2))) * s3)
の(x1 - x9)はs3が掛る.

3行目の
- (((x5 - x13) + (((x7 - x15) * s2) + ((x3 - x11) * s2))) * s1)))
の(x5 - x13)はs1が掛る.

まだ残っているのは,

(x3 - x11)*(s2 * s3 - s2 * s1)と
- (x7 - x15)*(s2 * s3 + s2 * s1)だ.

s2=sin(π/4)=cos(π/4),
s3=sin(3π/4)=cos(π/8)だから

s2*s3-s2*s1=sin(π/4)*cos(π/8)-cos(π/4)*sin(π/8)
=sin(π/4-π/8)=sin(π/8)=s1.
同様にしてs2*s3+s2*s3=s1なので,

(x3-x11)*s1, -(x7-x15)*s3
というわけだ.

2回目の結果も注意に値する.

0行目, 2行目に注意すると
0 (((x0 + x8) + (x4 + x12)) + ((x2 + x10) + (x6 + x14)))
2 ((x0 - x8) + (((x2 - x10) * s2) - ((x6 - x14) * s2)))
4 ((x0 + x8) - (x4 + x12))
6 ((x0 - x8) - (((x2 - x10) * s2) - ((x6 - x14) * s2)))

0はx0+x2+...+xeだからA0,

2は
x0 x2 x4  x6 x8  xa xc xe
 1 s2  0 -s2 -1 -s2  0 s2
だからA1,

4は
x0 x2 x4 x6 x8 xa xc xe
 1  0 -1  0  1  0 -1  0
だからA2,

のように偶数番目のxに対するフーリエ変換になっている. また奇数行目は奇数番目のxに対するフーリエ変換になっている.

上述のように, 記号の代入をやめて実際に計算すると, 数値計算用のフーリエ変換プログラムに改造できる.

xにsin(4πi/16)つまり2瘤のsin関数を置いてn=16で計算すると, B2がご覧のように8で, 他が0になる. このB2を2/n倍すると, 通常の係数1が得られる.
0 -6.432490598706545e-16
1 -8.378483910846401e-16
2 -2.0670557090385995e-15
3 1.1182945470482387e-15
4 2.4492935982947064e-16
5 -6.284358273892975e-16
6 5.97479550061776e-16
7 1.3277071107435814e-15
8 1.133107779529596e-15
9 -2.127518923182123e-16
10 8.000000000000002
11 1.1929584103349277e-15
12 0.
13 7.030996906759864e-16
14 1.7763568394002505e-15
15 2.771068273407289e-16
高速フーリエ変換のテストはSinカーブや鋸歯状波でやってみるのが一番である.

2013年7月6日土曜日

面積計を使う調和解析器

コンピュータやFFTの無かった時代, 調和解析には機械式の道具が活躍した. 調和解析器(harmonic analyzer)といえばKelvin卿のを思いだすが, もっと違う方式のものを最近知った.

G.U.Yule, "On a Simple Form of Harmonic Analyser," Proceedings of the Physical Society of London, XIII, 1894 に発表されたもので, 私はO. Mader, "Ein einfacher harmonischer Analysator mit beliebiger Basis," Elecktrotechnishen Zeitschrift, 1909を読み, どういうものかが分かった.

R.K.Otnes, "Notes on Mechanical Fourier analyzers"には下のような図が出ている.



左のパンタグラフのようなものは上下2台の面積計である. 右がその仕掛けで, 右上に円板がいくつか見えるのは, 求める係数により, 取り替える歯車である.

Maderの論文の図を書き直したのがこれだ. 計測の出発位置を示す.



下の方の曲線が計測すべき周期関数で, x=0からx=aまでが1周期である.

上のL字形の部分は台車(Wagen)で上下に動く. その台車の右の腕の先にKという軸があり, 梃子(Winkelhebel)FKSがその軸で回転できる. Kのx座標は下の曲線のa/2で, この時のKのy座標をh0とする.

梃子の下の腕の長さはmで, その先のポインタ(Fahrstift)Fで曲線を追跡する.

長さlの反対の腕の先にはローラーSがあり, その上下で左の歯軌条(Zahnstange)ZZ'が一緒に上下する.

台車にはDを軸とする半径Rの歯車(gezahnte Scheibe)があり, 台車に対するZZ'の相対的上下移動に従って回転する. 出発位置でのDの座標を(b,g)とする. この歯車には, 中心からrだけ離れた左と上に孔PsとPcがあって, それに面積計のポインタを接続する. つまりFが曲線を追跡し, 一周して来ると台車や梃子が動き, 歯車も上下に移動しながら回転し, その歯車の1点の描く面積を計算するのである. RはFがx=0からx=aまで移動した時に丁度n(係数の次数)回転するように出来ている. (製品によってはPsとPcの歯車が別になっていたらしい. 最初の図で同じ大きさの円板が2つずつあるのはその例だ.)

次数nに従ってRは変わるから, それぞれの歯車に対応して異なる軸用の孔が用意されている.

さて, Fを曲線の(x,y)まで移動した時の図が下だ.



Kの座標はy+hだ. hはFがx=0にあるときはh0だが, xがa/2になるまで増えるとhもh0より増え, a/2を過ぎると減ってx=aでh0になる.

ここからは式が多いのでTexで書く.





という次第で, Fls, Flcが求まるとBn, Anが得られる.

Maderの論文の最後に例があった. 周期が360ミリの鋸歯状波である.

下の図でA(0,0)からB(360,200)へ増加し, BからC(360,0)へ降下する波形だ. この図ではFはAにある. 歯車の赤丸はPs, 青丸はPcを示す. 最初の係数A0はcos 0=1を掛けて足すから, 面積計だけを使う. 面積は36000平方ミリだから平均の高さA0=100ミリである.



m=360ミリ, l=180ミリ, n=1, R=al/2πmn=28.6479ミリ, rもRと同じにした. 従ってK=90ミリである. (このKは梃子の軸のKとはもちろん違う.)

ABを20分割, BCを10分割, CAを20分割した各点での挺子の位置と歯車とPs, Pcの位置を示したのが次の図である.



これではよく見えないから, Ps, Pcの位置だけ取り出すとこうなる. 0番はA, 20番はB, 30番はCに対応する. BからCはx座標が変わらないから, 歯車は回転せず, 台車だけが下がるのが分かる. また歯車はちょうど1回転した場所なので, Psは左端に, Pcは上端にあることも理解できる.



Psのx座標, y座標, Pcのx座標, y座標は次のようであったので, その閉曲線の面積を計算すると, Fls=-5729.58301, Flc=0.00170898438. よってB1=-5729.58301/90=-63.6620331(ミリ). Maderの例の例の理論値と一致した. やれやれ. もちろんB7までもすべて一致している. 機械を使った実測値も2,3桁の精度なのはすごい. (上の右の図で面積が0なのを見るのは難しい.)

(-28.65 -27.25 -23.18 -16.84 -8.85 0. 8.85 16.84 23.18 27.25
28.65 27.25 23.18 16.84 8.85 0. -8.85 -16.84 -23.18 -27.25
-28.65 -28.65 -28.65 -28.65 -28.65 -28.65 -28.65 -28.65
-28.65 -28.65 -28.65 -27.25 -23.18 -16.84 -8.85 0. 8.85
16.84 23.18 27.25 28.65 27.25 23.18 16.84 8.85 0. -8.85
-16.84 -23.18 -27.25)

(311.77 340.34 366.78 390.41 410.66 427.22 439.97 449.1
455.03 458.4 460. 460.7 461.36 462.75 465.48 469.92 476.17
484.05 493.11 502.64 511.77 491.77 471.77 451.77 431.77
411.77 391.77 371.77 351.77 331.77 311.77 312.64 313.11
314.05 316.17 319.92 325.48 332.75 341.36 350.7 360. 368.4
375.03 379.1 379.97 377.22 370.66 360.41 346.78 330.34)

(0. 8.85 16.84 23.18 27.25 28.65 27.25 23.18 16.84 8.85
0. -8.85 -16.84 -23.18 -27.25 -28.65 -27.25 -23.18 -16.84
-8.85 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -8.85 -16.84 -23.18
-27.25 -28.65 -27.25 -23.18 -16.84 -8.85 0. 8.85 16.84 23.18
27.25 28.65 27.25 23.18 16.84 8.85)

(340.42 358.74 373.12 384.07 392.27 398.57 403.87 409.09
415.02 422.3 431.35 442.3 455.02 469.09 483.87 498.57 512.27
524.07 533.12 538.74 540.42 520.42 500.42 480.42 460.42
440.42 420.42 400.42 380.42 360.42 340.42 348.74 353.12
354.07 352.27 348.57 343.87 339.09 335.02 332.3 331.35 332.3
335.02 339.09 343.87 348.57 352.27 354.07 353.12 348.74)
ところでこれで調和解析ができることが直観的に分かる方法はないだろうか.

2013年6月13日木曜日

十六進算盤

このブログに八進法算盤のことを2回書いた. 2008年9月14日と2011年4月15日だ. そのとき十六進の算盤も考えられるとも書いた.

当時私の考えた十六進算盤は, 天(梁の上)に4珠が3つ, 地(梁の下)に1珠が3つのものだ. ところが最近Wikipediaを見ていたら, 昔の中国の算盤に5珠が2つ, 1珠が5つのものがあり, これは5X2+1X5=15だから16進法のものだという記述があって驚いた.

十六進算盤はどちらが使い易いか. ちょっと検討してみたい.

まず私が小学生の頃習った5珠1つ, 1珠4つの標準的十進算盤で, 加算の仕方をおさらいしよう.

0から9が置いてあるある桁に, 0から9を足すわけだが, 加算のアルゴリズムとしては(A,B)が置いてあるある桁に, (C,D)を足す; ただし

A,C={0,5} (天の値)
B,D={0,1,2,3,4} (地の値)
である.

まず地のCとDを足す.

0<=C+D<5なら地に1珠をD個上げてC+Dにする
9>C+D>=5なら地から(5-D)を下げ, C-(5-D)=C+D-5にし天に5を足す.

次に天でAとBを足す.

0<=A+B<10なら天をA+Bにする. A+B=10なら天を0にして1桁上に1を足す.

足すアルゴリズムの中に「足す」が現れるから, 再帰定義になっている.

私は小学生のころ, 学校で算盤を習ったが, 大体はこのアルゴリズムに従って考え考えやったように思う. しかし普通の人は反射的にやれるように練習しているに違いない.

ところで十六進の算盤の出番だ. 前のブログの図を再掲する.


上が十六進の算盤で, 左から右へ0からfの表わし方を示す. これをとにかく使うには, 0からfについて, 上の絵のように4で割った剰余(地の値B)と, 剰余を引いた差(天の値A)が反射的に思い出せる必要がある. 6なら2と4, 9なら1と8, cなら0と12という具合だ. (n=(A,B) ただしB = n mod 4, A = n - B)

足し算のアルゴリズムは下のようだ. 左端の見出しはいま算盤に置いてある数<1610で, 上端の見出しは足す数<1610である. 横が長いので, 足す数の途中で上下に分割してある.

多す数0の場合は何もしないから, 表は空白だ. 0,1,2に1を足すには単に1珠を1個上げればよい. それが「1↑」の表示だ.
上の算盤の絵で, 1を足すのはある桁と右隣の桁の図を見比べればよい. 従って3に1を足す時は, 下の3を払って(3↓), 上の4珠を1個下げる(4↓).

この調子で14+1まではよい. 15+1は図の右端の桁が左端の桁になるのだから, 3を払い, 上の4珠を3個払うことになる(4珠3個を払うのをc↑と示す.

それと同時に繰上げが生じたことを*で示す. つまり*は16, c↑は-12, 3↓は-3で, 都合+1になったわけだ.

他の場所もだいたいそういう方針で出来ている.



次の表は引き算のアルゴリズムである.



この表は目がしばしばするようなものだが, schemeでlatexのtabularのデータを発生させて書いている.

天と地に分解した表を作ると, 要するに下の表になってしまう.


これらの表も実は本質的に同じもので, 足される方は0,1,2,3, 足す方も0,1,2,3. 素直に足せれば足すが, 足せないときは4との差を引き, 繰り上げ処理をする. ある桁の地での足し算の繰り上げは天で処理し, 天での繰り上げは上の桁の地で処理するというだけのことである.

これだけ理解出来れば, 十進の算盤で1から10まで足して遊んだように, 十六進の算盤で0に1から1016まで足して8816にしてみよう.

次の絵がその様子である. 上の段が左から順に1から8まで足したところ(足した数は下に表示してある. 上は部分和を示す.) 下の段が9から1016まで足したところである.



こういう算盤を作ってみたくなった.

ところで中国の算盤を十六進で使うにはどうするか. まず0からfの表現法だが



上は天を6珠, 地を1珠と思う方法で, これだと十六進ではなく十八進の算盤になっている. その最後の16と17を使わない方法である.

下は天を5珠とみるもので, 天の値が0と5の時は地は0から4までとして使い, 0の時だけ0から5までを使うものである.

中国の人たちがどちらを使ったかは分からないが, 多分下の方であろう. しかし, これで加減算をするのは辛そうである. この図を描いただけでモティベーションが消え, まだアルゴリズムを考える気もしない.

2013年5月26日日曜日

サイコロの問題

数式処理Vol.19,No.2を眺めていたら, 算数オリンピックの話題があった. Mathematicaで解くという話だが, ちょっとやってみた.

2004年のジュニア算数オリンピックの問題だそうだ.

「すべて同じサイコロを図のように積み上げました. このサイコロのA, B, Cの面と反対側の面のサイコロのアルファベットをそれぞれ求めなさい.」



一番上のサイコロから, A,B,Cの面で作る頂点の方向から下の図のpのように見える;
また中段の右のサイコロから, D,E,Fの面で作る頂点の方向からはqのように見える;
さらに中段左から面C,Dは隣接する;
下のサイコロから面A,Eは隣接することが分る.



サイコロはどれかの面が上にあるとしてよいから, pのように面Aが上だとすると, B,Cが側面になり, D,E,Fの2つが側面に, 残ったのが底になるわけだ.

Fが底だとすると, qをFが底になるように回転してpにつけるとrのようになる. Eが底だとするとsになり, Dが底だとするとtになる.

これらの図で隣接関係を満すものを探すと, rではA, Eは隣接するがC,Dは裏の関係だからだめ. sではC,Dは隣接だがA,Eは裏でだめ. tはA,EもC,Dも隣接していてよさそうである.

というわけでAの裏はD, Bの裏はF, Cの裏はEであった.

おそらくMathematicaで解く方が面倒くさい.



サイコロの上のような展開図よりも, Rubicキューブで使ったこういう図の方が隣接関係はよく分るように思う.

ところで右の展開図は一筆描きができ, PostScriptで描くのも楽だった.

2013年5月13日月曜日

閏月

前回のブログの切っ掛けはCalendrical Calculationsに9,10,11,1月の閏は稀で, 12月の閏はほとんどないと書いてあったことだ.

これは前回のブログのような理由によるが, Aslaksenが1654年から2644年まで1000年の中国の太陰太陽暦を計算した結果もある. 日本と中国とは標準時が違うので, 節気の時刻が1時間違い, 多少は異なるかもしれないが, とにかくよく計算したものだ.

この閏月の頻度に関連した図がCalendrical Calculationsに載っている(p.249). たいした情報はないと思われるが, 気になるのでその説明をしたい. 下の図は私流に書き直したが, 実質的には同じである.



Calendrical Calculationsの図には上半分しかなく, 右上向きの斜めの線と上の12→1→...の右向きの線は破線. 更に閏月の頻度に従い, →閏11→, →閏1→, →閏9→, →閏10→は灰色, →閏12→はもっと淡い灰色にしてある. またその斜めの破線にその閏月の起きる確率が併記してある.

閏11月 0.005
閏12月 0.000
閏 1月 0.006
閏 2月 0.023
閏 3月 0.047
閏 4月 0.061
閏 5月 0.074
閏 6月 0.059
閏 7月 0.051
閏 8月 0.026
閏 9月 0.008
閏10月 0.009
横向きに1年の時間軸がある. 朔の字の下の縦線が朔の時刻で, それで区切られる短冊型がある暦月だ.

左端の暦月に11と書いてあるから, この間に冬至があるはず. 短冊の下の斜め線は冬至が11暦月の最初にあるか(左上)最後にあるか(右下)を示す.

冬至が11月の最初にあれば, その後の中気はその位置から右に辿った線上にあるはずで, 大寒は次の暦日の始まる翌日くらいにあり, 雨水はさらに次の暦月の2日目くらい, 冬至はだいぶ先の(右から2番目の)暦月の12日目くらいになる.

この場合にはどの暦月の下にも斜め線が存在するから, 閏月はなく, 月名は短冊の下の矢印を右に辿るようにすすむ.

冬至が斜め線の右下, つまり11暦月の最後の日だったらどうか. 大寒は「次の次」の暦月の最初になり, 「次」の暦月は閏月になる. 中気のない部分をこの図では灰色の帯で示す.

これが上の短冊で11から閏11へ向かう斜め矢印になり, その年は11月, 閏11月, 12月と進むことになる.

ある年はこの下の斜め線のある高さを左から右へすすみ, 途中で灰色の帯に遭遇すると, そこが閏月になり, 上の矢印は斜め右上方向へ進み, 13ヶ月の年になる. 灰色に出会わなければ12ヶ月だ.

図から分かるように, 閏月があるのは斜め線の領域の下の方だけで, その幅は7/19に対応するわけだ.

2013年5月10日金曜日

閏月

太陰太陽暦(lunisolar calendar)では「閏四月」というようなのがある. 12朔望月が1太陽年より短いから臨時に挿入する月だが, Calendrical Calculationsに9,10,11,1月の閏は稀で, 12月の閏はほとんどないと書いてあった. まぁ近日点に近いからとは思うが調べてみることにした. (英語では閏月をintercalary monthというらしい.)

インターネットで探したら旧暦を計算してくれるページがあった.

そこで早速1995年から2013年までの旧暦の各月の1日が太陽暦の何月何日かを書き出した.

表の上の方, 旧暦1995年1月1日は新暦の1月31日であると読む. 8月の所に2段あるのは, 旧暦8月1日は新暦8月26日. 旧暦閏8月1日は新暦9月25日であるということだ.



この期間に閏月は2月に1回, 3月に1回, 4月に1回, 5月に2回, 7月に1回, 8月に1回あり, 確かに9月から1月までには, 稀などころかまったくなかった.

また, 1995年から2013年までは19年あり, いわゆる19年7閏法のとおりに閏月が7回あることも確認できた.

(こんな苦労をせずとも私の書棚の「新こよみ便利帳」には1870年から2020年までの対照表があった. それを見ると明治6年(1873年)は閏6月があり, 13ヶ月なので, 政府が13ヶ月分の月給を払いたくないから新暦に改正した理由も分かる. 2033年の閏月は例外的といわれているが, そこまでは表がない.)

そもそも閏月はどこに入れるかをおさらいしよう.

A)1年の時間軸上に, 太陽と月の黄経が一致する時刻(朔)を決める.
B)平均太陽の南中時刻の12時間前から12時間後までを1暦日とする.
C)朔の時刻を含む暦日から次の朔の時刻を含む暦日の前の暦日までを1暦月とする.
一方,
D)太陽の黄経が30度の倍数になる時刻を中気という.
0 春分 30 穀雨 60 小満 90 夏至 120 大暑 150 処暑 180 秋分 210 霜降 240 小雪 270 冬至 300 大寒 330 雨水
春分を含む暦月を2月, 夏至を含む暦月を5月, 秋分を含む暦月を8月, 冬至を含む暦月を11月とし, その前後に中気を含む月の名前を決める.

1朔望月は29.5日, 中気の平均間隔は30.5日なので, 中気の割り当てられない暦月ができることがある. それを閏月として前の月の名前の上に閏をつける.

こういうわけだから, 閏月はどこにでも入り得るが, そうならないのは中気と中気の間隔(solar monthと書いてあったりする)が一定ではないからである. 地球が一定の角速度で公転するなら, 黄経が30度増える時間は一定であるが, Keplerの法則で地球が近日点付近にいる時は角度の増え方が大きく, 従って中気の間隔が短かく, solar monthの中に暦月がすっぽり入る確率が小さい. (正確に計算すると2033年のようなことも起きるが.)

この話は森口繁一先生の「数理つれづれ」にも書いてあるが, 自分で計算したのは私のブログ2008年6月25日の「夏至の日に」に書いた.

なお, 春分の暦月を2月にするということは, 春分が2月の後半にあるならその45日前の立春は1月にある, つまり新年立春のわけだし, 前半にあるなら立春は12月にある, つまり年内立春のことになり, 古今和歌集の「年のうちに春は来にけり」も確率1/2なことが理解できる.

2013年4月26日金曜日

割り算の九九

「二一天作の五」, 「二進の一十(にっちんのいんじゅう と発音する.)」などの割り算の九九というものがある. 私が計数工学科にいた頃, 森口繁一先生が時々その話をされた. 森口先生が小学生の頃は, 学校でこれを教えていた.

吉田光由の塵劫記にも, そろばんの絵とともとこの九九の説明がある. 塵劫記では割りのこえという. そういえば森口先生の「数理つれづれ」にもこの話があったことを思い出した.

要するに1桁の除数2≤d, 被除数n≤dに対する, 10n/dの商と剰余の表である.

使い方は後まわしにし, 2≤d≤9, 1≤n≤dについて商qと剰余rの表を作ろう.
(do ((d 2 (+ d 1))) ((> d 9))
(do ((n 1 (+ n 1))) ((> n d))
(let* ((nn (* n 10)) (q (quotient nn d)) (r (modulo nn d)))
(display (list d n q r)))) (newline))
実行すると(最上段のnと左端のdは後から私が書き込んだ. n=7から先は下に移動した.) 色分けについては下の変換規則A, B, ...などで説明する.
d\n  1        2        3        4        5        6  
2(2 1 5 0)(2 2 10 0)
3(3 1 3 1)(3 2 6 2)(3 3 10 0)
4(4 1 2 2)(4 2 5 0)(4 3 7 2)(4 4 10 0)
5(5 1 2 0)(5 2 4 0)(5 3 6 0)(5 4 8 0)(5 5 10 0)
6(6 1 1 4)(6 2 3 2)(6 3 5 0)(6 4 6 4)(6 5 8 2)(6 6 10 0)
7(7 1 1 3)(7 2 2 6)(7 3 4 2)(7 4 5 5)(7 5 7 1)(7 6 8 4)
8(8 1 1 2)(8 2 2 4)(8 3 3 6)(8 4 5 0)(8 5 6 2)(8 6 7 4)
9(9 1 1 1)(9 2 2 2)(9 3 3 3)(9 4 4 4)(9 5 5 5)(9 6 6 6)


     7        8        9 
(7 7 10 0)
(8 7 8 6)(8 8 10 0)
(9 7 7 7)(9 8 8 8)(9 9 10 0)
が得られる. この4つ組を(d n q r)ということにする. このパターンにより読み方が変わるのが分かりにくい理由である.

A. 各行の最後の(d d 10 0)は「d進の一十」と変換する.
B. (d n 5 0)のときは「d n天作の五」と変換する.
C. d=5, q=2n, r=0のときは, 塵劫記では「5 n倍作の2n」と変換するが, 掛け算より 森口先生の本のように「5 n加n」と変換しよう.
D. n=qのときは「d n下加r」と変換する.
E. それ以外は「d n q十r」と変換する.

これらの変換と算用数字を漢数字になおすと, 上の表は割り算の九九になる.

2(2 1天作の五)(2進の一十)
3(3 1 3十1)(3 2 6十2)(3進の一十)
4(4 1 2十2)(4 2天作の五)(4 3 7十2)(4進の一十)
5(5 1加1)(5 2加2)(5 3加3)(5 4加4)(5進の一十)
6(6 1下加4)(6 2 3十2)(6 3天作の五)(6 4 6十4)(6 5 8十2)
  (6進の一十)
7(7 1下加3)(7 2下加6)(7 3 4十2)(7 4 5十5)(7 5 7十1)
  (7 6 8十4)(7進の一十)
8(8 1下加2)(8 2下加4)(8 3下加6)(8 4天作の五)(8 5 6十2)
  (8 6 7十4)(8 7 8十6)(8進の一十)
9(9 1下加1)(9 2下加2)(9 3下加3)(9 4下加4)(9 5下加5)
  (9 6下加6)(9 7下加7)(9 8下加8)(9進の一十)
しかし割り算をシミュレートするプログラムを書くには変換前の表が必要である.

さてこの九九の使い方だが, 下の図の左端は我々のやり方である.


この136割る3を塵劫記のアルゴリズムを使いそろばんではこう計算するらしい. (図の右 a,b,...)

3の段の九九を使うと覚えておいて, そろばんに136と置く.

矢印の百の桁に注目. 1なので三一三十一を使う. 1を商の3に変え, 下の桁の3に剰余の1を足して4にし(b), 注目の桁を下へ移す(c). 注目する桁の数がd未満のときはいつもこうする.

注目の桁は4なので三進一十を使う. 4から3を引き上の桁に1を足す(d). 「d進一十」を使う場合は注目の桁は下へ移さない.

注目の桁が1なのでまた三一三十一を使う. つまりこの桁を3にし下の桁の6に1を足し, 注目の桁を下へ移す(e).

一番下の桁は7だから三進一十が2回使えて十の桁は3から5になり, 一の桁に剰余の1が残る.

注目の桁が0の場合は下の桁へ移る. 九九の表にn=0の列(d 0 0 0)があると思えばよい.

このアルゴリズムをSchemeで書くとこうなる.

sはそろばんに対応する配列で, 塵劫記には一二三...九を割る例があるのでそれをやってみるつもりだ. 九九の表からd行目, sx列目を取り出す添字の計算が面倒だが, 九九の三角の表を一列のリストにしたからである.
(define kuku '(
(2 1 5 0)(2 2 10 0)
...
(9 1 1 1)(9 2 2 2)...(9 9 10 0)))
(define (div d x)
 (display x) (display s) (newline)
 (if (>= x (vector-length s)) s
  (let ((sx (vector-ref s x)))
   (cond ((= sx 0) (div d (+ x 1)))
    ((>= sx d) (vector-set! s x (- sx d))
     (vector-set! s (- x 1) (+ (vector-ref s (- x 1)) 1))
     (div d x))
   (else (let* (
       (ku (list-ref kuku (+ (/ (* d (- d 1)) 2) sx -2)))
       (q (caddr ku)) (r (cadddr ku)))
    (vector-set! s x q)
    (vector-set! s (+ x 1) (+ (vector-ref s (+ x 1)) r))
    (div d (+ x 1))))))))

(define s (list->vector '(0 1 2 3 4 5 6 7 8 9 0)))

(div 2 0)
塵劫記流に米十二万三千四百五十六石七斗八升九合を2で割ると, 下のようになる. 右端の九九やゼロスキップなどコメントは私が追加した.
0#(0 1 2 3 4 5 6 7 8 9 0) ゼロスキップ
1#(0 1 2 3 4 5 6 7 8 9 0) 二一天作の五
2#(0 5 2 3 4 5 6 7 8 9 0) 二進の一十
2#(0 6 0 3 4 5 6 7 8 9 0) ゼロスキップ
3#(0 6 0 3 4 5 6 7 8 9 0) 二進の一十
3#(0 6 1 1 4 5 6 7 8 9 0) 二一天作の五
4#(0 6 1 5 4 5 6 7 8 9 0) 二進の一十
4#(0 6 1 6 2 5 6 7 8 9 0) 二進の一十
4#(0 6 1 7 0 5 6 7 8 9 0) ゼロスキップ
5#(0 6 1 7 0 5 6 7 8 9 0) 二進の一十
5#(0 6 1 7 1 3 6 7 8 9 0) 二進の一十
5#(0 6 1 7 2 1 6 7 8 9 0) 二一天作の五
6#(0 6 1 7 2 5 6 7 8 9 0) 二進の一十
6#(0 6 1 7 2 6 4 7 8 9 0) 二進の一十
6#(0 6 1 7 2 7 2 7 8 9 0) 二進の一十
6#(0 6 1 7 2 8 0 7 8 9 0) ゼロスキップ
7#(0 6 1 7 2 8 0 7 8 9 0) 二進の一十
7#(0 6 1 7 2 8 1 5 8 9 0) 二進の一十
7#(0 6 1 7 2 8 2 3 8 9 0) 二進の一十
7#(0 6 1 7 2 8 3 1 8 9 0) 二進の一十
8#(0 6 1 7 2 8 3 5 8 9 0) 二一天作の五
8#(0 6 1 7 2 8 3 6 6 9 0) 二進の一十
8#(0 6 1 7 2 8 3 7 4 9 0) 二進の一十
8#(0 6 1 7 2 8 3 8 2 9 0) 二進の一十
8#(0 6 1 7 2 8 3 9 0 9 0) ゼロスキップ
9#(0 6 1 7 2 8 3 9 0 9 0) 二進の一十
9#(0 6 1 7 2 8 3 9 1 7 0) 二進の一十
9#(0 6 1 7 2 8 3 9 2 5 0) 二進の一十
9#(0 6 1 7 2 8 3 9 3 3 0) 二進の一十
9#(0 6 1 7 2 8 3 9 4 1 0) 二一天作の五
10#(0 6 1 7 2 8 3 9 4 5 0) ゼロスキップ
11#(0 6 1 7 2 8 3 9 4 5 0) ゼロスキップ
nが大きい場合, 二進の一十をn/2回使うので大変である. 実用的な割り算の九九には2の段には四進の二十, 六進の三十, 八進の四十が, 3の段には六進の二十, 九進の三十, 4の段には八進の二十があったらしい.

でもこれらは異端であって, 初めに書いた九九は最初の数が除数を示しているのに対し, ここで追加したのはそうなっていない.

そもそも八進に一十, 二十, 四十があり, 六進に一十, 二十, 三十があり, 四進や九進も複数できて使いにくかったのではないかと想像する.

この割り算にはもう一つ問題がある. 9の段など特にそうだが下加5とか下の桁に剰余を足すと繰り上げが生じることがある. それは注目している桁, 商のところに繰り上がる. その対策はあまり見たことがない.

上のようなシミュレーションを見ると, 12などという数が現れて驚く. しかしこういう場合の次は「d進一十」になるし, 繰り上げは高々1なので, すぐにdずつの引き算が始まり, 繰り上げは1回でなくなるから, さしたる問題にならないのかもしれない.

除数が2桁以上, つまり多重精度除数による除算については, またの機会に.

2013年4月12日金曜日

微分解析機

前回のブログの微分解析機の絵は, 機素だけをつなげた, いわば回路図であった. 昔の微分解析機の論文を見ても, 解法にはこういう図しか示していない. 前回の接続図を佐々木達治郎著「計算機械」にある微分解析機の図解の絵のように描くと下のようになる. (情報処理学会誌Vol.52,No.3「微分解析機」の解説で使ったPostScriptのプログラムに手を加えた.)


トルクアンプのベルトは怪しいが, 佐々木の原図に従っている.

図の右半分にあるのが縦横にシャフトの並ぶ連結装置である. 解くべき微分方程式に従って, この部分は大々的に組み替える. 積分機や出力卓などと接続するクロスシャフトを持つ大小の ベイがバスボックスを介して配置されている. (ここでの説明に必要なベイまでしか描いていない.) 各ベイではクロスシャフトの上にバスシャフトを直角に置く. そのうちのある縦横軸対は, はすば歯車で連結される.

ベイとベイの間のバスボックスには数本のスタブシャフトが貫通しており, 隣り同士のベイのバスシャフトを連結する役目をはたす. スタブシャフトは普段は外さないらしい.

バスシャフトは両端にカプラーがはめてあり, カプラーの2本のネジを締めることでカプラーを介してバスシャフトとスタブシャフトは同時回転する. このネジを緩めるとカプラーをベイの内側の方へ滑らせて, バスシャフトを外すことが出来る.

積分機は今回の図では, 円板台に乗っている円板の方の位置が変るように描いた. 円板台を動かすのは平歯車n0を介している被積分関数軸である. この軸を回転すると, 円板の4時の方向に描いてあるナットにより, 円板台が左右に動く.

積分機0では回転子は円板の中心から左へ大きく離れているが(y'=1を示しているつもり), 右ネジに切れている被積分関数軸を連結装置側から見て反時計回りに回転すると, 円板は左へ移動し, 回転子と円板の中心の距離は減少する.

さて3本のバスシャフトは左下の出力卓でxに対するyとy'の曲線を描くように接続されている.

出力卓の上の方の3本の横軸は, 中央の右ネジのが水平移動軸で, xのバスシャフトで駆動され, 軸の左寄りにある滑り体が左右に動く. 滑り体にはyとy'の回転で垂直移動軸に取り付けた鉛筆を前後に動かすようにギアが付いている. 垂直移動軸は離れて2本あるが, 鉛筆は内側に設置され, 2つの鉛筆の先端のx座標が揃うようになっている.

ただ出力卓で描かれる2本の曲線が交差するとき, なにが起きるか, 衝突をどう避けるかがまだよく分らない. 例えば鉛筆を下の図のように取り付ければ, かすり傷かニアミスですれ違えるかも知れないがどうかな.


両移動軸のピッチが1ミリ, 描画可能領域が縦横それぞれ500ミリだとすると, 2πまで移動で, 前回の計算ではxは1600回転, yとy'は正負方向に460回転くらいするので, h0, h1, h2の平歯車はそれぞれ1/4,1/2,1/2に設定するのがよさそうである.

2013年4月6日土曜日

微分解析機

微分解析機の原理は分かっているし, Processingでシミュレータも書いたりしたが, 定量的にはどうか.

佐々木達治郎「計算機械」に解説はあるが, かつて東大理工研や阪大理学部で活躍し, 今は理科大で静態保存されている昭和航空計器製の微分解析機の数値をもとに自分で計算してみた.

下の図は積分機2台で正弦波d2y/dx2=-yを描く場合の機素の接続図である.



I0とI1が積分機で, Axは独立変数の縦軸(バスシャフト), xのA倍で回転する. By,-Dy''はI0の出力軸, I1の被積分関数軸でyのB倍で回転する. Cy'はyの1回微分の軸でI0の被積分関数軸, I1の出力軸だ. ダッシュ'が使えないプログラムではyをy0と, y'をy1 と書くことがある.

図でAx軸のすぐ下の斜め線はその縦軸の回転を, 1:1で積分機I0の変数の横軸(クロスシャフト)へ伝えるはすば歯車. その先の平歯車でm0倍になり円板を回転する. 平歯車の先の軸の回転角と円板の回転角は1:1で固定である.

Cy軸から積分機I0へ向かう横軸, 被積分関数は, n0倍になってから1回転が1ミリのピッチで円板と回転子の相対位置を変える. 図ではI0とI1の円板の位置が上下に揃い, 回転子がずれているが, トルクアンプの関係で, 通常は回転子が上下揃うように構成されている.

回転子は半径rで, j0の平歯車を経て出力軸と結ばれる.

下の積分機I1でも接続は同様だが, Byの軸の斜線の方向が違うのは, 被積分関数軸が縦軸に対して1:-1になっているからだ.

出力卓については次回のブログで説明しよう.

上の積分機I0で考えると, 回転子の円板の中心からの距離lはCy'n0ミリ. その状態でAx軸がdAxだけ回転すると, 円板がdAxm0だけ回転するから, 回転子の位置で円板は円周方向にdAxm0Cy'n0動き, 出力の横軸はdAxm0Cy'n0/rj0(=dBy)だけ回転することになる.

つまり, この値がByに足される.

この微分解析機にはまだ制約がある. rは32ミリ. 円板の半径は約120ミリ. 平歯車のギア比は4, 2, 1, 0.5, 0.25しかない.

そこでA, B, Cやギア比をいろいろ試行錯誤で計算してみると, m0=m1=2, n0=n1=0.25, j0=j1=4, A=256, B=C=448という使えそうな解があることが分かった.

もともとyやy'は-1≤y≤1だから, n0とn1が0.25だとBy,Cy'は, -112ミリから112ミリの間になり, 円板の範囲に収まる.

xを0から2πまで変えるとAxは0から512πまで変わる. yとy'の初期値はそれぞれ0と1だから, ByとCy'の初期値は0と448だ.

この1周の積分を1024分割で実行すると, dAxはπ/2である

これでシミュレートしてみる.
(define m 1) (define n 0.25) (define j 2)
(define a 256) (define b 448) (define c 448)
(define ax 0) (define by0 0) (define cy1 448) 
(define pi (* 4 (atan 1))) (define dax (/ pi 2))
(define r 32)
(define (integrator dx y) (/ (* m dx n y) r j))
(define (step i)
 (if (= (modulo i 16) 0)
  (begin (display (list i ax by0 cy1
    (+ (* (/ by0 b) (/ by0 b)) (* (/ cy1 c) (/ cy1 c)))))
   (newline)))
 (let ((dcy1 (integrator dax (- by0)))
       (dby0 (integrator dax cy1)))
  (set! ax (+ ax dax))
  (set! by0 (+ by0 dby0))
  (set! cy1 (+ cy1 dcy1))))
(do ((i 0 (+ i 1))) ((> i 1024)) (step i))
1024ステップ分出力するのは大変だから1/16つまり65行プリントした. 左からステップ番号, Ax, By, Cy', y2+y'2である.
(   0    0.00000    0.00000  448.00000 1.00000)
(  16   25.13274   43.92436  445.97712 1.00060)
(  32   50.26548   87.45205  439.65678 1.00121)
(  48   75.39822  130.16351  429.09730 1.00181)
(  64  100.53096  171.64681  414.39783 1.00241)
(  80  125.66371  211.50158  395.69748 1.00302)
(  96  150.79645  249.34290  373.17400 1.00362)
( 112  175.92919  284.80504  347.04204 1.00423)
( 128  201.06193  317.54493  317.55119 1.00483)
( 144  226.19467  347.24554  284.98349 1.00544)
( 160  251.32741  373.61893  249.65085 1.00604)
( 176  276.46015  396.40903  211.89195 1.00665)
( 192  301.59289  415.39414  172.06908 1.00725)
( 208  326.72564  430.38907  130.56463 1.00786)
( 224  351.85838  441.24697   87.77739 1.00847)
( 240  376.99112  447.86075   44.11879 1.00908)
( 256  402.12386  450.16415     .00887 1.00968)
( 272  427.25660  448.13236  -44.12771 1.01029)
( 288  452.38934  441.78236  -87.86579 1.01090)
( 304  477.52208  431.17271 -130.78379 1.01151)
( 320  502.65482  416.40305 -172.46777 1.01212)
( 336  527.78757  397.61316 -212.51543 1.01273)
( 352  552.92031  374.98162 -250.54001 1.01334)
( 368  578.05305  348.72414 -286.17396 1.01395)
( 384  603.18579  319.09147 -319.07260 1.01456)
( 400  628.31853  286.36704 -348.91733 1.01517)
( 416  653.45127  250.86423 -375.41882 1.01579)
( 432  678.58401  212.92338 -398.31976 1.01640)
( 448  703.71675  172.90852 -417.39737 1.01701)
( 464  728.84950  131.20387 -432.46556 1.01762)
( 480  753.98224   88.21016 -443.37676 1.01824)
( 496  779.11498   44.34079 -450.02335 1.01885)
( 512  804.24772     .01783 -452.33875 1.01946)
( 528  829.38046  -44.33200 -450.29802 1.02008)
( 544  854.51320  -88.28149 -443.91821 1.02069)
( 560  879.64594 -131.40702 -433.25816 1.02131)
( 576  904.77868 -173.29266 -418.41798 1.02192)
( 592  929.91143 -213.53415 -399.53811 1.02254)
( 608  955.04417 -251.74286 -376.79800 1.02315)
( 624  980.17691 -287.54947 -350.41438 1.02377)
( 640 1005.30965 -320.60761 -320.63922 1.02439)
( 656 1030.44239 -350.59717 -287.75730 1.02501)
( 672 1055.57513 -377.22738 -252.08351 1.02562)
( 688 1080.70787 -400.23970 -213.95984 1.02624)
( 704 1105.84061 -419.41026 -173.75205 1.02686)
( 720 1130.97336 -434.55206 -131.84624 1.02748)
( 736 1156.10610 -445.51682  -88.64505 1.02810)
( 752 1181.23884 -452.19640  -44.56390 1.02872)
( 768 1206.37158 -454.52385    -.02688 1.02934)
( 784 1231.50432 -452.47414   44.53723 1.02996)
( 800 1256.63706 -446.06438   88.69916 1.03058)
( 816 1281.76980 -435.35369  132.03323 1.03120)
( 832 1306.90254 -420.44265  174.12149 1.03182)
( 848 1332.03529 -401.47238  214.55776 1.03244)
( 864 1357.16803 -378.62318  252.95148 1.03306)
( 880 1382.30077 -352.11282  288.93159 1.03369)
( 896 1407.43351 -322.19447  322.15002 1.03431)
( 912 1432.56625 -289.15431  352.28509 1.03493)
( 928 1457.69899 -253.30872  379.04465 1.03556)
( 944 1482.83173 -215.00134  402.16889 1.03618)
( 960 1507.96447 -174.59970  421.43285 1.03680)
( 976 1533.09721 -132.49175  436.64863 1.03743)
( 992 1558.22996  -89.08210  447.66722 1.03805)
(1008 1583.36270  -44.78813  454.37993 1.03868)
(1024 1608.49544    -.03601  456.71951 1.03931)
axの最後の値は (* pi 512) => 1608.495438637974だ.

この出力をもとにサークルテストをしたのが下の図だ.



開始点と終了点である右端がちょっとずれている. 448と456だから2%くらいの誤差である.

2013年3月20日水曜日

三角錐の体積

今月初めにあった埼玉県公立高校入試問題が新聞に出ていた.

その中にこういう数学の問題がある.

右の図は, 三角錐の展開図です. △ABCはAB=16cm, BC=8cm, ∠ABC=90°の直角三角形です. また, 点D, Eは, それぞれ辺AB, ACの中点であり, 点Fは, 線分DBの中点です. このとき, 線分DE, EF, FCを折り曲げてできる三角錐の体積を求めなさい.



面白そうだ. これはちょっとやってみようと思いつつものびのびになっていた.

まずは正攻法.



Dの辺りを原点とし, DとBが合わさった頂点をP, その座標をx,y,zとする. 三角錐の底面の頂点E,F,Cからの稜の長さも図に示す通りだから, 式が書ける.

x2+(4-y)2+z2=42
(4-x)2+y2+z2=42
(8-x)2+(8-y)2+z2=82

x=yなのは初めから分る. 解くとx=y=z=8/3が得られ, 底面の△EFC=24(cm2)だから体積は64/3(cm3).

思いのほか簡単な数値であった. ところで待てよと思う. なるほど∠Bと∠ADEが直角なので, EFPの面に対して稜PCは垂線になるわけだ. こちらの三角錐で計算すれば, 底面積は8(cm2), 高さは8(cm)なので暗算で計算できる問題であった.

2013年3月9日土曜日

スマートエイジング

3月6日から8日まで, 東北大学川内キャンパスで, 情報処理学会の第75回全国大会が開催された.

招待講演のひとつが, 東北大学加齢医学研究所の川島隆太教授の「スマート・エイジング〜脳機能解析学が拓く新しい超高齢社会〜」であった.

高齢者のひとりとして聞かざるを得ない話題である.

高齢期には
□脳を積極的に使う
□身体を積極的に動かす
□バランスのとれた栄養を摂る
□社会や他者と積極的に関わる
が肝心であるとのご高説であった.

脳を積極的に使うには, いわゆる脳トレが不可欠で, ふたつのデモがあった.

ひとつは「スパン課題」で, 次のような画面で9個の黒い丸がランダムに一瞬だけ色が変る. その順番を覚えておくものである.



もうひとつは「Nバック課題」で, 簡単な加算の問題が次々と現れる. それをN個あとに答えるものである. つまり計算の結果を長さNのシフトレジスタにいれ, 出口側から出る数を答えるのである.

これらの脳トレ課題は, 任天堂から発売されているそうだが, 私としてはProcessingでプログラムを書きたい衝動に駆られた.

仙台からの帰途は, 金曜の晩だけあって, 自由席を指定席に変える列が長大で, 並んでいる内に指定席がなくなるような状態であった.

止むをえず, やまびこの自由席で帰ることにしたが, 私の好きなE2系の1000番台で仙台始発. 缶ビールを飲みながら, 早速「スパン課題」のプログラムに 着手した. 最初のバージョンは列車が大宮に着くずっと前に完成した. 簡単であった.

今日はそれに機能を追加した. 画面の下方の6から15は1回の試行に色の変る回数である. 6をクリックすると6回色が変わり, 変わった位地が画面の下に表示される. 位置は最上段の左から0,1,2. 中段は左から3,4,5. 最下段は左から6,7,8である.



画面の右に見えるスライドバーは, 下端にあるときは, 色の替わる周期が2分ごとで, 上端では30秒ごとになる.

脳トレに興味があれば, ここに置いてあるので, 使ってみてほしい.

「Nバック課題」はそのうち書いて見よう.

2013年2月28日木曜日

継子立て

Concrete MathematicsのJosephusの問題のところにこういう演習問題があった.

最初に善玉n人を並べ, 次に悪玉n人を並べて2n人の円形にし, m人ごとに除外すると, 最初に悪玉が全部除外されるような, nに依存したmが存在することを示せ.

例えばn=3ならm=5, n=4ならm=30とできる.

Suppose there are 2n people in a circle; the first n are "good guys" and the last n are "bad guys." Show that there is always an integer m (depending on n) such taht, if we go around the circle executing every mth people, all the bad guys are first to go? (For example, when n=3, we can take m=5; when n=4, we can take m=30.)

n=3,m=5でやってみると,

位置の0,1,2が善玉で, 3,4,5が悪玉なので, 位置4,3,5の順に除外され, たしかにそうなる.

計算機の威を借りてさらにmを探すと
n=3 (5 52 60 65 112 120 125 172 180 185 232 240 245 292 ...)
n=4 (30 71 101 175 205 216 246 320 350 391 421 ...)
n=5 (169 217 330 378 979 ...)

Answerを見ると

m be the least (or any) common multiple of 2n, 2n-1, ..., n+1.

n=3のとき, 上の図では, 2nごとに位置5のところへ回ってくるから, 2nの倍数で5が除外され, 人数が1人減るから次は2n-1の倍数で4が除外され, 最後n+1の倍数で3が除外される.

n=3のときはlcm(6,5,4)=60だからm=60,120,180,... m=5はその方法とは違って偶然うまくいく場合だったらしい. 65,125,185もその流れの解である.

n=4だとlcm(8,7,6,5)=840. n=5だとlcmは2520になる.

nから計算でmは得られるのが分かったが, 問題の例に偶然見つかる方を書くのはどうかなぁ.

ところでConcrete Mathematicsの訳書で, この問題は「2n人が環状に並んでいるとする. 半分のn人は「よい人間」で, 残りのn人は「わるい人間」だとする.」と始まる.

この記述では善悪それぞれn人ずつが固まっているとは読めない. 継子立てのように, 30人が環状に入り交って並び, その半分の15人が先妻の子で, 残りの15人が後妻の子であるのと同様に読めてしまう.

私は原書の方を読んでやっと問題の意味が理解できた. もう少し慎重に翻訳して欲しい. (ここを訳したのは誰だ.)

2013年2月26日火曜日

継子立て

わが家で塵劫記が, それも継子立てが話題になった. 私が小学生の頃持っていた本に塵劫記のことが書いてあり, 私は読んでいたから, 継子立ては70年来のお馴染みである.

Concrete MathematicsやMathWorldによると外国ではJosephusの問題というらしい.

概要はこうだ. 先妻の子15人, 後妻の子15人がいる. 後妻はその中から1人選ぶのを自分の子にしたい.

15人ずつの先妻の子, 後妻の子をある順に円形に並ばせ, 出発の位置から順に数えて10人目ごとに除外する.

次の図の外側の円がそれだ. 円周の内側の黒の数字が位置で, 真上の0から時計回り1,2,...,29まで続く. 円周上の小円内の数字が赤のところには後妻の子を, 黒のところには先妻の子を置く.

小円の中の数字は除外される順を示す. 真上の0から時計廻りに数えると10人目は位置9だ. これが最初なので小円には0を記す. 位置9の子を除外し, その次から数えると10人目は位置19 になる. 小円には1を記し, それも除外し,... を続けると, 小円内の14までが黒なので, (不思議にではないが)先妻の子ばかり14人除外される.



先妻の子の15人目が除外されそうになると, その子は, 次は自分から始めてほしいと訴え, 内側の円のように再開する. 今度は不思議に後妻の子ばかり除外され, 最後に真上の先妻の子0 が残るという話だ.

後妻が子の配置を決めるのは簡単である. 外側のような図を書いてシミュレーションし, 15番目に除外される場所までを先妻の子にすればよい.

しかしコンピュータ屋のわれわれが手でシミュレーションするのは業腹だ. やはりプログラムを書こう.

lispでは円形のリストが作れる.

たとえば
(define foo '(0 1 2 3))
(set-cdr! (list-tail foo 3) foo)

位置 0 1 2 3 4 5 6 7 8 9 a b ..
要素(0 1 2 3 0 1 2 3 0 1 2 3 ..)
の無限リストfooができる.

このようにして30人のリストhead (0 1 2 ... 29)を作り, set-cdr!で円形, 無限にする. 10番目, つまり0から数えれば9番目の要素を除外するには
(set-cdr! (list-tail head 8) (list-ref head 10))
で9の両側のリンクが繋がる. これでheadは9の抜けた
(0 1 2 3 4 5 6 7 8 10 11 12 ...)
になる. そうしてheadを9進める.
(set! head (list-tail head 9))
この様子は次の図の左のようだ. headはnextの位置に進む.



しかし円を構成するメンバーの数が少なくなると, これでは失敗する. 上の図の右を見て欲しい. headから0,1,...,と進むと9番はremoveと書いてある1になる. 従ってその両端を繋ぐ. headをnextのところに進めたいが, 9先は0,2,3,4,5,6,7,0,2,3と進んでnext'と書いた3になる.

この失敗をしないためには, 除外する前に10進んでnextの位置を先に確保することだ. そうして書いたプログラムが次だ.

(define (exclude)
(let ((next (list-tail foo 10)))
(set! bar (cons (list-ref foo 9) bar))
(set-cdr! (list-tail foo 8) next)
(set! foo next))
(display bar)
(display (map (lambda (n) (list-ref foo n)) (a2b 0 16))) (newline)
'ok)
後妻も16人の場合をシミュレートすべきであった.

2013年2月18日月曜日

Rubicキューブのシミュレータ

Winning WaysのRubicキューブの直し方(curing)が分ってきた. 実はその理解のために今回のシミュレータを書いたような次第だ.

島内本では6つの面をtop, north, east, south, west. bottomというが, 今回はWinning Ways流(というかDavid Singmaster流)にUp, Back, Right, Front, Left, Downとする.

各面の回転は時計まわりにU,B,R,F,L,D; 反時計まわりにはU',B',R',F',L',D'. またWinning Waysの説明には, 真ん中の段の回転もあり, それにはギリシア文字α,β,δ,γ,ε,ωを使う. εとωはeastとwestのようで覚えやすい.



このブログでは, 各小体の面を下の図の左のように表す.



右の図は操作の順で, まず下段の辺の2面体Aから始めEへとステージを進める.

A) Aの2面体を向きも合わせて定位置に置く.
B) Bの3面体を向きも合わせて定位置に置く.
C) Cの2面体を向きも合わせて定位置に置く.
D) Dの2面体を定位置に置く. 向きは気にしない.
E) Eの3面体を定位置に置く. 向きは気にしない.
F) Dの2面体とEの3面体の向きを合わせる.

Winning Waysではこれらの各ステージを

A) Aloft, Around (Adjust) and About.
B) Bottom Layer Corner Cubelets.
C) Central Layer Edge Cubelets.
D) Domiciling the Top Edge Cublets.
E) Exchanging Pairs of Top Corners.
F) Finishing Flips and Fiddles.

という.

ステージA) 下段は操作しにくいから, 上段のDの位置に集めてから各側面を回転してAに置く. Aのステージは簡単で, Aloft, Aroundなどはあまり関係ない.

ステージB) すでに置き終えたAの2面体に影響しないようにBの3面体を置く. 例えば面の図のpsの3面体(矢印の先)を置くには, それが上段にあればcnq の位置に来るように上段を回転し, 底の面の色がqにあればB1, nにあればB2, cにあればB3を実行する. 3面体が下段にあればいづれかの操作で上段に上げてから今の操作を行う.



図A) B1:F'U'F 底の面の色がqにあるとき
図B) B2:RUR'底の面の色がnにあるとき
図C) B3:F'UFRU2R' 底の面の色がcにあるとき

3面体が下段に移る他, orの2面体や上段のキューブにも影響があるが, そこは後で処理するので問題ではない.

ステージC) etかblにあるDの2面体をorのC(矢印の先)に移動する. 下段が揃っているから, 単純な面の回転は出来ない.



図D) C1:URU'R'U'F'UF oの色がUの面にあるとき
図E) C2:U'F'UFURU'R' rの色がUの面にあるとき

ここまでで下段と中段が揃う. この後は上段の2面体Dと3面体Eの入れ替えになる.

ステージD) 上段の2面体etとgを交換する.



図F) D1:UFRUR'U'F' 2面体bl, 3面体ai, cnq, hvも替わったように見えるが, これらは向きが替わっただけで, 位置は前のままである. 上段の2面体, 3面体の向きは最後のステージFで修正するから, 今は気にしない.

ステージE) 3面体cnqとhvを交換する(mono swapという).



図G) 1月22日のブログのG)巡回 で3面体aiがcnq, cnqがhv, hvがaiに移動している. これを3面体の交換2回でもとに戻す.
図H) 上段のT'の回転.
図I) E:FDF2D2F2D'F'=Msを実行. Hの赤丸と青丸の3面体を交換した. 中段と下段にも変化があるが, 上段はこの2つが変わっただけなのに注意.
図J) 上段のT'の回転.
図K) E=Msを実行. Jの赤丸と青丸の3面体を交換した. 上段の3つの3面体は希望の位置に移った. Ms2=1なので, 中段と下段がもとに戻ることに注意.

ステージF) 上段の3面体と2面体の向きを合わせる.

3面体の場合


図L) 1月22日のブログのH)ねじり で3面体cnqが時計回りに, hvが反時計回りにねじれている.
図M) F1=(F'RFR')2=Maで3面体cnqを反時計回りに回転する. 中段と下段にも変化があるが, 上段はこの3面体が回転しただけなのに注意.
図N) 上段のT'の回転.
図O) F2=(RF'R'F)2=Mcで時計回りに回転する. 中段と下段はMaMc=1で元に戻る. 後は上段をT回転する.

2面体の場合


図P) F)隣辺向き替え で2面体blとetの向きが反転している.
図Q) F2=(εR)4=Meでblを反転する.
図R) 上段をT'回転する.
図S) F2=Me を実行する. 中段と下段はMe2=1で元に戻る. 後は上段をT回転する.

この方法はMs2=1, MaMc=1, Me2=1をうまく利用していて面白い.

2013年2月4日月曜日

Rubicキューブのシミュレータ

シミュレータが出来たので島内先生の本の6面体完成術を実装することにした.

今回はその大体のやり方を書くことにしよう.

島内流にいうとRubicキューブには26個の小体があり, その内訳は隅にある3面体8個, 辺にある2面体12個, 面の中央にあり動かない1面体6個とである.

島内流の完成術は
A. 3面体を向きは無視して正規の位置に移す.
B. 2面体を向きは無視して正規の位置に移す.
C. 2面体の向きを揃える.
D. 3面体の向きを揃える.
である.

この各手順で主に使われるのが, 1月22日のブログにあった
E) 単純3角形 9a(Bで使う)
F) 隣辺向き替え24a(Cで使う)
G) 巡回 28a(Aで使う)
H) ねじり 33a(Dで使う)
の操作である.

まず6面の色を番号で表す.
上t 0(白)
北n 1(緑)
東e 2(橙)
南s 3(黄)
西w 4(赤)
下b 5(青)

隅(赤字で示す)と辺(青字で示す)も番号で表す.

3面体は正規の位置での3つの色を昇順に並べて表す. 隅番号0,1,2,...にある3面体は順に
[012],[023],[034],[014],[125],[235],[345],[145]

2面体は正規の位置での2つの色を昇順に並べて表す. 辺番号0,1,2,...にある2面体は順に
[01],[02],[03],[04],[12],[23],[34],[14],[15],[25],[35],[45]

現在 隅0にある3面体を知るには, 54の面の色のリストから0,15,18の位置の色をとりだし, 昇順にソートする. こうして8つの隅にどの3面体がいるかが分る.


[ランダムにした出発点]

Aの手順では, [012]を隅0に戻すことから始めて, 隅7まで戻すことである. 最初の頃は間単だ. たとえば[012]が隅番号n=0,1,2,3の何れかにあれば, 面tをnだけ左回転する. n=4,5,6,7なら面bをn-4だけ右回転し隅4に置き, 面eを左回転する.

[023]は隅2にいれば, 面sを右回転. 3にいればw右, s右だ. n=4,5,6,7なら(n+3)%4だけ右回転して隅5へ移し, sを左回転する.

このようにして上の4個はなんとかなる. 下の隅はまとめて対処する.

4つの3面体の相対位置により, 巡回 28aを1回か2回使うと, 正規の位置に収まる.


[3面体が元の位置にもどる]

次はBの手順で, 今度は隅に影響を与えないよう, 単純3角形 9aだけで対応する. これも辺0に[01]を移すのがもっとも楽で, 7になるとなかなか苦しい. 下の面の8,9,10,11は, Aの手順の最後のようにまとめて対応する. これも巡回を2回使うことで無事に完了する.


[2面体も元の位置に戻る]

Cの手順は隣の辺と一蓮託生なので, たとえば辺0の向きを変えると, 辺1,4,3,7のどれかも向きが変わる. それでせっかく揃っていた向きも壊れるかもしれない.

まず一筆描きの辺の順を決める. たとえば0,1,2,3,7,8,4,9,5,10,6,11とする. 辺0の2面体の向きが違っていたら, 0と1の辺の隣辺向き替えを行う. 次に辺1の2面体の向きを調べ, 違っていたら1と2の辺の隣辺向き替えを行う. このように順々に向き替えを続けるとこの手順は完了する.


[2面体の向きを揃える]

Dはねじりで, これは向きが揃うまで同じ操作を2回やらなければならないかもしれない. これも一筆描き, 0,1,2,3,7,6,5,4のように実行すれば 終わる.


[3面体も向きを揃える]

手順AやBの小体の移動に較べると, 向き替えははるかに簡単であった.

問題はそれぞれの操作がある位置についてしか書いてないことで, 人手で揃えるときはキューブをその向きになるよう持ち替えるのだが, 私はx,y,z軸について操作を回転した操作を生成するSchemeのプログラムを書き, あらゆる方向の操作を用意した上で好きな位置で処理が出来ようにした.

乱数を発生させてキューブをこねたあと, 元へ戻すプログラムを走らせると, 一瞬にして完成するのが感動的である. もちろん途中経過を出力しているから, 手順を追って完成させていることも確認できる.

2013年1月25日金曜日

Rubicキューブのシミュレータ

Rubicキューブのどの面をどう回せといわれても, 手でやっていれば, 手順が複雑になってくるとつい間違えてしまうのは誰にも経験があろう.

その点, シミュレータはプログラムをしっかり書けば何回でも正しく処理できる.

Winning WaysにはElena Conwayが好んだパターンを作る手順が書いてある. どんなパターンか興味があったが, シミュレータを使えばそのパターンを作るのは簡単であった.

Aは4 Windows. Bは6 Windows. CはChequers. DはHarlequin. EはStripey. FはZigzag. Gは4 Crosses. Hは6 Crossesというらしい.

手順は省略する. Winning Ways, Vol. 4, p. 876参照.

2013年1月22日火曜日

Rubicキューブのシミュレータ

情報処理学会誌にRubicキューブの話を書いたのは2005年7月号であった. その話題は島内先生の名著(島内剛一:ルービック・キューブ免許皆伝, 日本評論社(1981))の置換をHaskellで記述することであった. (この本はルービック・キューブと数学パズル, 日本評論社(2008)として再刊された.)

最近Winning Waysを見ていたら, Rubik's Hungarian Cubeなる話題があり, いろいろ書いてあるが結構ややこしい. 我が家にもRubicキューブはあったが, 最近は見当たらないので, 私のことだからPostScript(とProcessing)でシミュレータを書いた.

それを使ってWinnig Waysの方法を試みているのだが, なかなかうまくはいかない.

今回はとりあえずシミュレータの話だ.

学会誌にあった絵だが, 下の図でAは魔方体を横からみた図で, この方向を南(S)とする. 上はトップ(T)とする. 右にわずかにみえる面は東(E)である. 当然みえない右の方に西(W), 向こう側に北(N), 下にボトム(B)の面がある.

BはTの面を回転しているところで, 時計廻りをt(あるいはt+), 反時計廻りをt-という. CはSの面の回転で, s(s+)とs-になる. DはE面の回転e(e+)とe-である. (どの図もちょうど45度回ったところで, 回転方向は分からない.)

今回のシミュレータはカラーで, 上の図では見えなかった3面が右につなげて描いてある.

Bの図では左半分の上段の色が右回転した(t)後であることを示している. C, Dはsとeの結果である.

この下の図も(図の記号は違うが)学会誌の図で, 島内本にある基本操作
E) 単純3角形 9a
F) 隣辺向き替え24a
G) 巡回 28a
H) ねじり 33a
を示す. そのシミュレータ版が次だ. 54の面には次のように0から53の番号が振ってあり, リストに各面の色が入っている. 回転によりその色を入れ替えている. 例えば回転tについては
((0 2 8 6) (1 5 7 3) (15 24 33 42) (12 21 30 39) (9 18 27 36))
のような置換のリストがある. つまり面0は2へ行き, 2は8へ行き, 8は6へ行き, 6は0へ行くを示す. t-はこのリストを逆に読む. まだ前途遼遠だが, 出来たところまでをまずアップロードしよう.