閏月

前回のブログの切っ掛けは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に対応するわけだ.

閏月

太陰太陽暦(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なことが理解できる.

割り算の九九

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

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

要するに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桁以上, つまり多重精度除数による除算については, またの機会に.

微分解析機

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


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

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

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

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

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

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

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

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

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


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

微分解析機

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

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

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



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

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

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

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

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

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

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

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

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

そこでA, B, Cやギア比をいろいろ試行錯誤で計算してみると, m₀=m₁=2, n₀=n₁=0.25, j₀=j₁=4, A=256, B=C=448という使えそうな解があることが分かった.

もともとyやy'は-1≤y≤1だから, n₀とn₁が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', y²+y'²である.

(   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%くらいの誤差である.

三角錐の体積

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

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

右の図は, 三角錐の展開図です. △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からの稜の長さも図に示す通りだから, 式が書ける.

x²+(4-y)²+z²=4²
(4-x)²+y²+z²=4²
(8-x)²+(8-y)²+z²=8²

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

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

スマートエイジング

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バック課題」はそのうち書いて見よう.