面積計を使う調和解析器

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

一般のnについては,

nA_n=-1/π∫sin(nθ)dy,
nB_n=1/π∫cos(nθ)dy,

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

A_n=1/π∫₀^2πy cos(nθ)d&theta,
B_n=1/π∫₀^2πy sin(nθ)d&theta,

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


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

面積計を使う調和解析器

厳密にいうとこれは面積計は使わないが, 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月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

など. 便利である.

面積計を使う調和解析器

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

下の図は適当に描いた周期関数である. 関数は赤, 緑, 青, 橙の各点を通り, 赤に戻る. 左端の赤い縦線はθ=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だった.)

もっと先の方の係数A_n, B_nを求めるには, 元の関数の図を横に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/π)²×π=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年頃)に載っているらしい.

高速フーリエ変換

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

データの個数の少い例なので, プログラムは簡単に書ける. 書きなが, 学生のころ(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のデータX₀,X₁,...,X_n-1(nは2の巾乗)のFourier変換. A_m=∑_i=0^n-1X_icos(2πmi/n) (m=0,1,...,n/2), B_m=∑_i=0^n-1X_i sin(2πmi/n) (m=1,2,...,n/2-1)をすべてのmについて計算する. cosineの最後の変換はプログラム上 sineの最初の変換の場所に入っている. すなわちB₀は常にゼロのために, B₀の場所には求めたA_n/2が入ってくる.」
とある. このような記述も懐しい. なおプログラムの作成は1966年8月22日だ.

sとxの初期化の後のdoループで, n2をn,n/2,n/4,...とlog₂n回まわし, その次の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行目は係数のA₀, 1行目はA₁, 2行目はA₂などである.

A₀は(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だからA₀,

2は

x0 x2 x4  x6 x8  xa xc xe
 1 s2  0 -s2 -1 -s2  0 s2

だからA₁,

4は

x0 x2 x4 x6 x8 xa xc xe
 1  0 -1  0  1  0 -1  0

だからA₂,

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

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

xにsin(4πi/16)つまり2瘤のsin関数を置いてn=16で計算すると, B₂がご覧のように8で, 他が0になる. このB₂を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カーブや鋸歯状波でやってみるのが一番である.

面積計を使う調和解析器

コンピュータや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が求まるとB_n, A_nが得られる.

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)

ところでこれで調和解析ができることが直観的に分かる方法はないだろうか.

十六進算盤

このブログに八進法算盤のことを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)

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

多す数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から10₁₆まで足して88₁₆にしてみよう.

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



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

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



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

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

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