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