2022年4月12日火曜日

mex関数

mex関数についてブログに全回書いたのは, もう2年も前だった. その頃, Schemeでの実装も 出来ていたが, なんとなく気に入らなかった.
最近, もう少し気持ちよくしたいと思い, やった書き直したので, ブログに記録して おく. まずinsertion. (ins n t)は, nを, 木tに挿入する関数である.
(define (ins n t)  ;inserts n in a tree t, returns a new tree
  (if (null? t) (list '() (list n n) '())
      (let ((lt (car t)) (l (caadr t)) (r (cadadr t))
            (rt (caddr t)))
  (define (insr t)
    (let ((lt (car t)) (m (caadr t)) (s (cadadr t))
          (rt (caddr t)))
      (if (null? rt)
        (if (= s (- n 1)) (begin (set! l m) lt)
            (begin (set! l n) t))
        (list lt (list m s) (insr rt)))))
  (define (insl t)
    (let ((lt (car t)) (m (caadr t)) (s (cadadr t))
          (rt (caddr t)))
      (if (null? lt)
        (if (= m (+ n 1)) (begin (set! r s) rt)
            (begin (set! r n) t))
        (list (insl lt) (list m s) rt))))
  (cond ((<= l n r) t)
        ((= n (- l 1))
          (if (null? lt) (list lt (list n r) rt)
              (let ((lt (insr lt)))
                (list lt (list l r) rt))))
        ((< n (- l 1)) (list (ins n lt) (list l r) rt))
        ((= n (+ r 1))
          (if (null? rt) (list lt (list l n) rt)
              (let ((rt (insl rt)))
                (list lt (list l r) rt))))
        ((> n (+ r 1)) (list lt (list l r) (ins n rt)))))))
一方, これを使う(mex ms)は次のようだ.
(define (mex ns) (let ((t '()))
 (for-each (lambda (n) (set! t (ins n t))) ns)
  (if (null? t) 0
      (begin (while (not (null? (car t))) (set! t (car t)))
             (if (< 0 (caadr t)) 0 (+ (cadadr t) 1))))))
全回と似たような図だが, 8 4 12 2 6 10 14 7 3 5 1 11 13 15 9 の順に挿入した時の木構造の移り代りの図を示す.
左上の A [8,8] 8 は, 左のAが図の記号, 右の8が今回挿入したn, 中央が(最初は空だった)木にnが挿入された, 新しい木 のように見る.

2021年10月31日日曜日

クイーン支配問題

前回のクイーン支配問題のブログは, 8×8の大きさのチェス盤に, 5個の クイーンを置き, 盤面のすべての場所が支配できる配置の数を計算するものであった. 4860という解が得られたが, どう考えても, 90度の回転や, 左右の反転で同じになる 仲間が沢山含まれている筈である.

そこで次の課題は, 本質的に異なる配置の数を知ることである. 前回のすべての配置の それぞれに, 回転や反転を施し, それで一緒になる配置を仲間とするプログラムを書き, 実行すると, 638という答が出た.

クイーン支配の配置は, かなりランダムらしいから, 回転で仲間が4倍になり, 反転で また2倍になって, 仲間は8人がほとんどになりそうである. しかし, 反転するとどうなるか. 左右反転では, 偶数の列に奇数のクイーンを置くから, クイーンの奇数個を中央の列に置き, 残りの偶数個を, 左右対称の位置に置けばいいわけだが, 中央の列が存在しないから, 左右や上下の反転で一緒になるものはない. 一方, 対角線で反転することを考えると, 対角線には中央があるから, そこに奇数個を置くと, 対称になり得る. この対角線の 中央を通る, もう一方の対角線で, もう一回反転できるかというと, 対角線の長さ が8という偶数なので, そういう反転はあり得ない.

従って, 4人の仲間のグループと, 8人の仲間のグループとで, グループは638あり, メンバーの総数は, 4860人ということになる.

4人と8人のグループはそれぞれ何個あるかは, いわゆる鶴亀算で得られる. a+b=638, 4a+8b=4860 だから a=61, b=577である.

そこで, その638の配置をすべて図示したのが次である. 1行に16個, 縦に20行あって, 最初のページに320個, 次のページに318並んでいる. 638個の配置には, 0番から637番 の番号があり, 最初のページの左上が0, そこから右へ1,2,3,...と続く. 次のページの最下行の右端が, 637番である.

これらの図は, 小さ過ぎるので, それぞれの図で, 各クイーンがどのように全盤を支配する かを示す図を描いた. 下の20の図は, 上の図の左端の配置の一つ置きに描いたもので, 塗りつぶしの丸がクイーン, それから四方八方に出る破線が支配領域である. 8×8 のすべての盤面が破線で蔽われているのが分る.

最後に4人仲間の配置はどれかという話題である. 合格発表のような下の番号が, 2回対称軸を持つ61の4人組である.

(20 25 29 30 59 128 129 137 141 149 163 165 171 180 182 183 204 215 219 229 232 238 242 258 259 261 267 276 283 288 289 365 368 372 391 393 441 468 470 473 474 511 516 517 518 524 526 530 551 568 569 574 582 584 592 594 605 617 620 623 633)

この表の0, 15, 30, 45番目を選び(20, 183, 289, 524番になる), 上と同じように 示したのが, 下の図である. 対角線で反転するのが分る. 中のふたつは面白い.

2021年10月27日水曜日

クイーン支配問題

クイーン支配問題 The Art of Computer Programmingに, クイーン支配問題(queens domination probmel)という話題があった.

8×8の盤に, 5個のクイーンを, 盤のすべての場所を支配するように置く という問題である. 1863年の本にあるそうだから, 古典的問題である.

この解は, 4860通りあるらしいので, 自分でも計算することにした.

私のことだから, Schemeを使う. またSchemeを使うなら, bit-stringが役立つ.

プログラムを書いたら, 思ったより簡単で, 1時間はかからなかった. こんな具合だ.

(define (r n) (map (lambda (i) (+ (* n 8) i)) (range 0 8)))
(define (c n) (map (lambda (i) (+ (* i 8) n)) (range 0 8)))
(define (d n) (if (< n 0)  ;n=i-j -7<=n<8
  (map (lambda (i) (+ (* i 9) (- 0 n))) (range 0 (+ n 8)))
  (map (lambda (i) (+ (* i 9) (* n 8))) (range 0 (- 8 n)))))
(define (u n) (if (< n 7)  ;n=i+j 0<=n<15
  (map (lambda (i) (+ (* i -7) (* n 8))) (range 0 (+ n 1)))
  (map (lambda (i) (+ (* i -7) (+ n 49))) (range 0 (- 15 n)))))
(define (bs n)
  (let ((i (quotient n 8)) (j (modulo n 8))
	(b (make-bit-string 64 #t)))
    (for-each (lambda (m) (bit-string-clear! b m))
      (append (r i) (c j) (d (- i j)) (u (+ i j))))
    b))
(define bss (map bs (range 0 64)))
(do ((i 0 (+ i 1)) (count 0) (bs (make-bit-string 64 #t)))
    ((= i 60) count) (display (list 'i i))
(let ((bs (bit-string-copy
       (bit-string-and bs (list-ref bss i)))))
  (do ((j (+ i 1) (+ j 1))) ((= j 61))
   (let ((bs (bit-string-copy
          (bit-string-and bs (list-ref bss j)))))
    (do ((k (+ j 1) (+ k 1))) ((= k 62))
     (let ((bs (bit-string-copy
            (bit-string-and bs (list-ref bss k)))))
      (do ((l (+ k 1) (+ l 1))) ((= l 63))
       (let ((bs (bit-string-copy
              (bit-string-and bs (list-ref bss l)))))
	(do ((m (+ l 1) (+ m 1))) ((= m 64))
         (if (bit-string-zero?
              (bit-string-and bs (list-ref bss m)))
          (set! count (+ count 1))))))))))))
盤のますには, 左上から右上にかけ, そして上から下へ, 0から63まで番号を付ける. 行rは上が0, 下が7(その番号をiとする), 列cは左が0, 右が7(その番号をjとする), 左上から右下への対角線dは, ます7を通る(i=0,j=7)のを-7, ます0を通るのを0, ます56を通る(i=7, j=0)のを7, 左下から右上への対角線uは, ます0を通る(i=0,j=0) のを0, ます56を通る(i=7,j=0)のを7, ます63を通る(i=7,j=7)のを14とする.

先頭の関数, r,c,d,uは, 行n, 列n, 下り対角線n, 上り対角線nにいるクイーンが支配する ますの番号を返す.

bssは, ますnのクイーンが支配するますをビット0, それ以外をビット1にした, 64ビット のビット列の配列である。<br>
これだけ用意し, 後は5個のクイーンを, 全面に置きながら, 支配場所を集め, 5個置いた 時に, すべてのます(に対応するビット)が0なら, countを1増やす.

やってみたら, 18秒くらいで, 4860が得られた.

2021年6月7日月曜日

Piの1000桁

前回のブログのPiの1000桁をもう少し詳しく説明したい. Piの値の計算に Machinの式

π/4 = 4 tan-1(1/5)- tan-1(1/239)
tan-1(1/5)=1/5-1/(3·53)+1/(5·55) -1/(7·57)+...
tan-1(1/239)=1/239-1/(3·2393)+1/(5·239 5)-1/(7·2397)+...

を使った. これは絵を描いてみるとこういうことだ.
円弧の中に底辺:高さが5:1の直角三角形を描く. その中心角をαとする. つまり α=tan-1(1/5). αは11.3°くらいだ.

αの斜辺の上にまた今の直角三角形を描く. それを3回繰り返すと4αの 角度が得られ, 上の値を使えば4α=45.2°だ. 一方π/4は45°だから, 今度は底辺:高さ=239:1の直角三角形(赤線)を角度を引く方法に描くとその斜辺が 丁度45°になるというのがMachinの式であった. 計算すると次のようになる.
一番上はtanの加法公式である. これを使いαから2α, 4α を 計算すると4α=120/119が得られる. この120と119の差が1なのが 重要だ. 4αとπとの差の角度βのtanはまた加法公式を使うと1/239 になるわけだ. 一番下のようにMachinの式になる.

Machinが分ったから計算に進む. 最初のPiの1000桁のブログでやったように, 40桁計算する様子を見るには, 小数点を40桁右に移動して長い整数(bignum)で計算する. つまり1/5とせず, ss=1040/5 から始める.

それを1で割ってtsに足し, ss=ss/(5×5)にしてから ss/3をtsから引き, ss=ss/(5×5)にしてから ss/5をtsに足し, ... を計算する. しかし最初だけ 5で割り, 後は5×5で割るなら, 最初の分子を5倍しておき, いつも5×5 で割るのがプログラムは簡単だ.

さらにMachinの式はπ/4を計算するが, 我々はπが欲しいから, 最初のss は 4×4×5×1040とする.

従って
(set! ts 0)
(set! ss (* 4 4 5 (expt 10 40)))
続けて
(set! ss (quotient ss 25))
(set! ts (+ ts (quotient ss 1)))
(set! ss (quotient ss 25))
(set! ts (- ts (quotient ss 3)))
...
これをssが0になるまで実行する. 初期化の後 doループを使い
(set! ts 0)
(set! ss (* 4 4 5 (expt 10 40)))
(do ((n 1 (+ n 2)) (ff #t (not ff))) ((= ss 0) n)
  (set! ss (quotient ss 25))
  (set! ts ((if ff + -) ts (quotient ss n))))
doループは, nを1, ffを#tで始め, ループを回る度にnはn+2にし, ffは ¬ffにしてssが0になるまで, 本体を実行すると読む. つまりnは1,3,5,7,... ffは#t,#f,#t,#f,...となる. 終了条件s=0の後のnはdoループの返す値である. nをどこまで進めるとsが0になったかを知るためだ.

((if ff + -) ts (quotient ss n))は, ffが#tなら (+ ts (quotient...)), #fなら(- ts (quotient ...))のことである. つまり1回置きに加算と減算を繰り返す.

初期化したssとループでのssとtsを出力すると
ss 800000000000000000000000000000000000000000(
(n 1)
ss 32000000000000000000000000000000000000000
ts 32000000000000000000000000000000000000000
(n 3)
ss  1280000000000000000000000000000000000000
ts 31573333333333333333333333333333333333334
(n 5)
ss    51200000000000000000000000000000000000
ts 31583573333333333333333333333333333333334
(n 7)
ss     2048000000000000000000000000000000000
ts 31583280761904761904761904761904761904763

...

(n 57)
ss                                        23
ts 31583289575980921339207962431166446951613
(n 59)
ss                                         0
ts 31583289575980921339207962431166446951613
61
最後の61はdoループが返したnの値である. 下の方の23や0は元々は右端に 揃っていたのだが, ブログはそういう風には出してくれない.

同じようにtan-11/239も計算すると
(n 17)
ss                                         1
ts 31415926535897932384626433832795028841973
(n 19)
ss                                         0
ts 31415926535897932384626433832795028841973
21
と円周率が見えてくる.

余談だが, 私の恩師 故高橋秀俊先生は中学生の時円周率を筆算で50桁計算されたそうだ. 40桁でもn=57まで計算するのだから, 50桁はまだまだ先までで大変だっただろう.

次はこの計算を有限語長ので分割しながら実行するわけだが, それは次のブログにしよう.

2021年5月25日火曜日

Piの1000桁

1970年頃のことだ. HITAC 10, FACOM R, MACC 7など, 1語が16ビットの4K 語のミニコンが多く使われていた. これらのミニコンは, 正面に押しボタ ンのようなスイッチが並んでいて, 電源を入れた後, スイッチを操作し, 数語の「初期入力プログラム」を手動で入力する. プログラム全体が入力 し終わると初期入力プログラムが動き出し, 紙テープのプログラムを読み 込めた.

ところがこのメーカー製プログラムの手の入力が面倒で, 誰もがもっと短 いプログラムに出来ないか工夫を凝らし, 徐々に短いプログラムが出現し た. これが一番短いと思っていると, 他から更に短いのが出来たという話 が来る. するとまた知恵を絞り, その語数のプログラムが作れるのであっ た.

つまり「何語で出来る」という情報が重要で, それを聞くとその語数のプ ログラムに向って注力できる.

最近ツイッターを眺めていたら, Automatic Computer@nt_a_c さんの

Automatic Computer@nt_a_c·3月18日
PC-1で円周率、400桁は出せた。もう数日早くできればちょうど良かったのだが。後はどのくらいコードを削れるか。

Automatic Computer@nt_a_c·3月24日
円周率、できた。1000桁、印刷と合わせて3回読み込み、合わせて40分くらい。https://automaticcomputer.github.io/PC1onUnity/pi....

というツイートがあった. 60年も前の計算機 PC-1でプログラムを書く奇 特な人がまだいるかと驚く一方, どうすればあの容量のメモリーでπが1000桁計算 できるか私も考え始めた.

πの1000桁については, このブログでも2012年11月5日に書いている. 簡単に計算するには やはりMachinの式でやることになろう.

π/4 = 4 tan-1(1/5)- tan-1(1/239)

tan-1(1/5)=1/5-1/(3·53)+1/(5·55) -1/(7·57)+...

tan-1(1/239)=1/239-1/(3·2393)+1/(5·239 5)-1/(7·2397)+...

PC-1は1語が18ビット. 計算は2語をつなげた長語36ビットを使う. 十進10 桁で1長語を使うから, 1000桁計算には100長語が必要で, eの1000桁では 300番地から500番地の100長語を作業場所に使っている. 512短語のPC-1で は100長語はどうしても2組しかとれないから,それで上の計算をしなけれ ばならない.

普通に考えると16×tan-1(1/5)を計算するには, atanを置く場所tsを0にし, 多倍長の16/5を作る. これをssとする. それを1で割り(その商をqsとする), qsをtsに足す. ssを25で割りssにし, そのssを3でわりqsとし, tsからqsを引く. ...

これをやるとtan-1(1/5)を計算するだけで100長語が3組必要になりすでに 破綻する.

tsはtan-1(1/5)とtan-1(1/239)で共用出来る. ssは 1/5が済んだ後1/239用に使えばよい.

残る手段はqsを作ってから, 下の桁からtsに足したりtsから引いたりするのを なんとかすることだ.

我々が小学校以来習った加減算は下の桁から始めるものであった. それに対して 上級になって習った算盤では上の桁から始めた. 小学生のころ, 上の桁から加減算 が出来るのが不思議であったが, そのうち繰上げをローカルに処理しながら下へ進んで いるからだと納得した. さらに一度繰上げが通過した桁に, その後から再び繰上げ が来ることはないことも証明できた. 基数が10から2^35と大きくなると, 繰上げが何桁も伝搬する確率は小さくなるが, 無視できない. い

要するにqsの結果はいらない. 各々の桁の商はqsに足すか引くかし, 剰余は次の 桁の除算に使うと作業場所は要らなくなり, なんとかなりそうだ.

そこでSchemeを使ってシミュレートしたのが次のプログラムである.

lnは長語数. b35, d10は二進35桁, 十進10桁の定数. tsはatanの値を蓄積する 場所, ssはatan(1/z)の計算で1/z, 1/z^3, 1/z^5,...を保持する場所である. (zは5と239) ldivはssをd = z^2で次々割る. ldiv+- はssをn(1,3,5,..)で割りながら, 部分商をtsに足すか引くかする.
(define ln 100)
(define b35 (expt 2 35))
(define d10 (expt 10 10))
(define (ldiv ss d)
  (let ((r 0) (n 0))
    (do ((i 0 (+ i 1))) ((= i ln))
      (set! n (+ (fix:lsh r 35) (list-ref ss i)))
      (list-set! ss i (quotient n d))
      (set! r (modulo n d)))))
(define (ldiv+- ts ss ff d)
  (let ((r 0) (n 0) (s 0) (c 0))
    (do ((i 0 (+ i 1))) ((= i ln))
      (set! n (+ (fix:lsh r 35) (list-ref ss i)))
      (set! r (modulo n d))
      (set! s ((if ff + -) (list-ref ts i)
        (quotient n d)))
      (list-set! ts i (fix:and s (- b35 1)))
      (set! c (fix:and s b35))
      (do ((j (- i 1) (- j 1))) ((= c 0))
        (set! s ((if ff + -) (list-ref ts j) 1))
        (list-set! ts j (fix:and s (- b35 1)))
        (set! c (fix:and s b35))))))
(define (lmul ns m)
  (let ((e 0) (p 0))
    (do ((i (- ln 1) (- i 1))) ((= i 0) e)
      (set! p (+ (* (list-ref ns i) m) e))
      (list-set! ns i (modulo p b35))
      (set! e (quotient p b35)))))
(define ts (make-list ln 0))
(define ss (make-list ln 0)) (list-set! ss 0 80)
(do ((n 1 (+ n 2)) (ff #t (not ff))) ((> n 1433)) 
  (ldiv ss (* 5 5))
  (ldiv+- ts ss ff n))
(define ss (make-list ln 0))(list-set! ss 0 (* 4 239))
(do ((n 1 (+ n 2)) (ff #f (not ff))) ((> n 423))
  (ldiv ss (* 239 239))
  (ldiv+- ts ss ff n))
(do ((n 0 (+ n 1))) ((= n 100))
  (display (list (lmul ts d10))))
この出力の始めの方は次のようだ.

1415926535 8979323846 2643383279 5028841971
6939937510 5820974944 5923078164 628620899
8628034825 3421170679 8214808651 3282306647
938446095 5058223172 5359408128 4811174502
8410270193 8521105559 6446229489 5493038196
...
10^10進法の各桁の後に空白があるのは, 十進の一番上が0なのが分るためだ.

これで一応は計算出来るが, 考えてみるに, ssは徐々に小くなっていく. その途中にずーっと先頭から除算を始めるのはやはり無駄である. つまり 除算を開始する位置を制御すれば計算は早く終わる. 下はそのように修正 したプログラムで, 変数bが除算開始位置を覚えている.

(define b 0)
(define (ldiv ss d)
  (let ((r 0) (n 0) (q 0))
    (do ((i b (+ i 1))) ((= i ln))
      (set! n (+ (fix:lsh r 35) (list-ref ss i)))
      (set! q (quotient n d))
      (list-set! ss i q)
      (if (= q 0) (set! b (+ i 1)))
      (set! r (modulo n d)))))
(define (ldiv+- ts ss ff d)
  (let ((r 0) (n 0) (s 0) (c 0))
    (do ((i b (+ i 1))) ((= i ln))
      (set! n (+ (fix:lsh r 35) (list-ref ss i)))
      (set! r (modulo n d))
      (set! s ((if ff + -)
        (list-ref ts i) (quotient n d)))
      (list-set! ts i (fix:and s (- b35 1)))
      (set! c (fix:and s b35))
      (do ((j (- i 1) (- j 1))) ((= c 0))
        (set! s ((if ff + -) (list-ref ts j) 1))
        (list-set! ts j (fix:and s (- b35 1)))
        (set! c (fix:and s b35))))))
(define (lmul ns m)
  (let ((e 0) (p 0))
    (do ((i (- ln 1) (- i 1))) ((= i 0) e)
      (set! p (+ (* (list-ref ns i) m) e))
      (list-set! ns i (modulo p b35))
      (set! e (quotient p b35)))))
(define start (real-time-clock))
(define ts (make-list ln 0))
(define ss (make-list ln 0)) (list-set! ss 0 80)
(set! b 0)
(do ((n 1 (+ n 2)) (ff #t (not ff))) ((> n 1433))
   (ldiv ss (* 5 5)) (ldiv+- ts ss ff n))
(define ss (make-list ln 0)) (list-set! ss 0 (* 4 239))
(set! b 0)
(do ((n 1 (+ n 2)) (ff #f (not ff))) ((> n 423))
   (ldiv ss (* 239 239)) (ldiv+- ts ss ff n))
(do ((n 0 (+ n 1))) ((= n 100))
  (display (lmul ts d10)) (display " "))
上のプログラムをPC-1のプログラムにしたのが次だ. 作業番地tsはeの1000桁と同様, 300番地以降を使う. プログラムの最初にあるインタールードでクリアしておく.

ssは0番地から200番地までを使う. つまりイニシアルオーダーR0は壊す. 計算だけで記憶場所を使うので, 結果tsの出力は別のプログラムにする.

まだコメントが少ないが, 徐々に追加するつもりだ.

      0p=200,
  0   0p=196,0p:p2, インタールード 0をaccへ
  1   tl(300),  作業領域ts(300から498)をクリア
  2   p1r,      
  3   a40,      2を足す
  4   x1r,
  5   s8r,      終了か
  6   z40,      R0へ戻る 
  7   jlr,      ループする
  8   500,jlp.

  0   0p=196,0n=276,0m=496,0p:a7n,
  1   x16r,     atan サブルーチン 帰り番地を入れる
  2   p4n,      作業場所ssをクリア
  3   x7r,
  4   s19n,   終了?
  5   z11r,
  6   p4n,      0をaccに
  7   tl(   ),  作業場所へ
  8   p7r,      アドレスを増やす
  9   a7n,
 10   jl3r,
 11   pl10n,    ssの最初の語に数を入れる
 12   tl0,
 13   pl4n,     atan 計算ループの回数
 14   tl2n,
 15   sl8n,     その回数が0になる迄
 16   zl(   ),   帰り番地
 17   q4n,      z2の除算ループ. 最初の剰余をRregに
 18   p22n,
 19   x23r,     ssの語のアドレスを順に作る
 20   x25r,
 21   s19n,     終り?
 22   z34r,
 23   ql(   ), 被除数をaccとRregに
 24   dl12n,     z^2による除算
 25   tl(   ),  商を格納
 26   sl4n,     商が零かテスト 1を引き零なら負にする
 27   kl31r,
 28   p22n,    除算開始の位置b
 29   a7n,   bを2増やす
 30   x22n,
 31   p23r,
 32   a7n,      ssのアドレスを増やす
 33   jl19r,
 34   q4n,   
 35   p22n,    ssをnで除し, tsに加減する
 36   x43r,
 37   a20n,
 38   x46r,
 39   x49r,
 40   x52r,
 41   s21n,    ループの終か
 42   z70r,
 43   ql(   ),  被除数を作り
 44   dl2n,     nで割る.
 45   tln,      商を一旦しまう
 46   pl(   ),  ts
 47   aln,   aln<->sln 商を足すか引く
 48   k51r,  sign bitが立っていたら繰上げ処理へ
 49   tl(   ), 商をしまってこの語は終り
 50   jl67r,
 51   cl14n,   sign bitを消して
 52   tl(   ), しまう
 53   p46r,  繰上げ処理
 54   s7n,
 55   x58r,  繰上げ先の語のアドレス
 56   x61r,
 57   x64r,
 58   pl(   ), 
 59   al4n,  aln<->sln 繰上げを足すか引く
 60   k63r,  sign bitがあれば繰上げ
 61   tl(   ), 繰上げがなければ処理は終り
 62   jl67r,
 63   cl14n,  sign bitを消す
 64   tl(   ), しまう
 65   p58r,   繰上げ先のアドレスを進める
 66   jl54r,
 67   p43r,   ssの語のアドレスを進める
 68   a7n,
 69   jl36r,
 70   p79r,  バイナリスイッチ切替え
 71   b47r,
 72   t47r,
 73   p79r,
 74   b59r,
 75   t59r,
 76   pl2n,   nを増やす
 77   al6n,
 78   jl14r,
 79   i,      a xor s = i

  0   0n:,,   データ領域
  2   ,,
  4   ,1,
  6   ,2,
  8   ,1435,
 10   ,80,    =5*16 atan(1/5)で使う
 12   ,25,
 14   131071,262143, マスクパターン
 16   425,956,    956=239*4
 18   57121,194,  57121=239*239 
 20   300,494,
 22   0,

  0   0m:it,  主ルーチン
  1   jlp,   atan(1/5)
  2   p16n,  atan(1/239)用に
  3   t9n,   データを入替え
  4   p17n,
  5   t11n,
  6   p18n,
  7   t13n,
  8   p4n,
  9   t22n,
 10   p4n,
 11   t22n,
 12   it,
 13   jlp, atan(1/239)
 14   0,     計算終り 停止する
      jlm.
このプログラムは私のPC-1シミュレータに追加してある. シミュレータを開き, R0 を読み込み, Place Tapeをクリック. メニューの下方にあるpi1000cmpを クリックすると計算のプログラムが設定される. そこでFree Runをクリック; さらにInit Startをクリックするとプログラムが走り出す.

計算は私のMacBook Pro環境では約3分半で終り, シミュレータは停止する.

結果の出力には作業場所 を同じにしてあるので, eの1000桁の出力プログラムを多少修正して使う.

計算プログラムがR0領域を破壊したから, 出力プログラムの読込み前にもう一度 R0をロードする. 後はPlace TapeからPi1000prtを選び, そのプログラムを 読込み起動すればπの1000桁が出力される.

久し振りにPC-1のプログラムを書いたが, インタプリタで必要な情報を取り出し たり, Schemeのプログラムで予備テストが出来たり, 便利な世の中である.

2021年3月29日月曜日

PC-1シミュレータ

このブログにパラメトロン計算機 PC-1のシミュレータの話を書いたのはそろそろ10年前の ことだ. そのシミュレータはProcessingで書いたものだが, その後シミュレータの 置いてある計算機が引退したりで, 長いあいだ動かなかった. ブログを見た方には申し訳 ないことをした.

ずーっとシミュレータの書き直しを考えていたが果せず, 最近やっとJavaScriptで書き, どうやら無事に動くようなので, 公開することにした.

シミュレータは https://www.iijlab.net/~ew/pc1sim2/pc1sim.html にあるので試してみて欲しい.

シミュレータの見掛けは以前のままである. もともとPC-1が出来たころ, 近隣の 計算機と競争していた, 自然対数の底eの値を1000桁計算するプログラムを シミュレートするのが目的だったので, そのプログラムに特化してはあるが, 他の プログラムも少数ながら動く.

シミュレータを開くと次のような画面になる. 右手に3x3の押しボタンがある.



最初にやることは初期入力ルーチンR0の読込みで, まず左下のPlace R0をクリックする. R0のテープが現れたら, 次にその上のInit Loadをクリック. これでR0が読み込まれる.

今度は中下のPlace e1000をクリック. e1000桁のテープが現れる. そこで 中央のFree Runを押し, 続けて中央上のInit Startをクリック. これでe1000桁の テープがR0で読み込まれ, 計算が始まる. 程なくして結果の印字が始まり, 計算が 終了する.

他に用意してあるテープは,

factorize32 231より少し小さい整数の素因数分解
eratosthenes エラトステネスの篩
v2test フーリエ変換
va1test 素数を法とする行列式の計算

で, これらを走らせるには右下のPlace Tapeをクリックする. そうすると右上に プログラム名が表示されるから, それをクリックする. テープが現れたら中央の Free Run, その上のInit Startをクリック. なかにはテープ読込み中に停止する ものもあるが, その時は中央のFree Runと, 右上のRestartをクリックする.

2020年12月10日木曜日

無限クイーン

TAOCPの演習問題7221-38に, 無限クイーン列<qn>の効率的計算法を考案せよというのがある.

面白い解答なので, その説明をしたい.

二進の配列 a,b,cを使う. cは正と負の添字がある.
G1. [初期化.] r←0, s←1, t←0, n←0とする. ( 1≤k≤nについて, qkの計算が済んだ.)
G2. [qn≤nを調べる.] (この時点で1≤k<sについて ak=1でa2=0; また-r<k≤tについてck=1 でc-r=ct+1=0; 各ベクターにはn個の1がある.) n←n+1, k←sとする.
G3. [見つけた?] ak=bk+n=ck-n=0なら G5へ進む. そうでないなら, k←k+1とし, k≤n-rならこのステップを繰り返す.
G4. [qn>nを作る.] t←t+1; qn←n+t; ak←b2n+t←ct←1とし, G2へ戻る.
G5. [qn≤nを作る.] qn←k; ak←bk+n←ck-n←1とする. k=sなら, as=0になるまでs←s+1を繰り返す. k=n-rなら c-r=0になるまでr←r+1を繰り返す. G2へ戻る.

このプログラムを忠実に辿ると, 大体こういうことをやっていると判明する.

TAOCPのプログラムでは, チェス盤の行も列も1から始めるが, 私はやはり0から 始めたい. 前回の無限クイーンの最後にあるように, 解の列は

0 2 4 1 3 8 10 12 14 5 7 18 6 21 9 24 26 28 30 11 13 34 36 38 40 15 17 44 16 47

のように始まる. その始めの方の置きかたを見ると(下の図 前回の図では, 行は下に 延びていたのに, 今回は逆に上に延びるように描いて申し訳けない.)

まず列n=0は行0にクイーンが置ける. それが左下の黒丸だ. そうするとどの行にはもう置けない ことを示す配列aの0を1に, 行-列の対角線に置けないことを示す配列cの0-0=0を1にする. 図ではそれを黒丸のすぐ右の青小丸とすぐ右上の赤小丸で示す. 配列bも0+0=0が 1になるがそれは図には示さない.

続いてn=1にクイーンを置くのだが, 青線が示すようにこの行には置けず, 赤の対角線が 示すようにこの斜線にも置けない. よって赤斜線の上のq=2に置くことになる. それが (1,2)の黒丸で, その右の青小丸と右上の赤小丸で配列の設定も分る.

n=2の列は, 行0は青線でだめ, 行2は赤線でだめ, その間の行1は, すぐ左上にクイーン がいてだめ(これは配列bで分る). 行4も(1,2)からの斜線のためだめで, 結局 行4に置く.

このようにして置いていくのだが, 青線はこれより下はすべて詰っている; 2本の赤線は この間はすべて詰っていることを示す. したがって行0の青線は, 列3の行1に置くと 行2まで詰ったので, ここからは行2から横に進む.

赤の斜線も同様に出来ている. 従ってある列でクイーンが置けるか調べるには, 青線路 の上から始め, 下の赤線の下までである. ここまでで配列a, b, cを調べ, 置けない 時は, 上の赤線の上に置くことになる.

この方針で私流に書いたプログラムが次だ.
(define qs (make-list 30))
(define as (make-list 50 0))
(define bs (make-list 80 0))
(define cs (make-list 40 0))
(define coff 20)
(define (a k) (list-ref as k))
(define (b k) (list-ref bs k))
(define (c k) (list-ref cs (+ k coff)))
(define (a! k) (list-set! as k 1))
(define (b! k) (list-set! bs k 1))
(define (c! k) (list-set! cs (+ k coff) 1))
(define (place k) (list-set! qs n k)
 (a! k) (b! (+ k n)) (c! (- k n)))
(define n 0) (define s -1) (define t 0) (define r 0)
(define k)
(while (< n 30) (set! k (+ s 1))
 (while 
  (and (< k (+ n r))
   (or (= (a k) 1) (= (b (+ k n)) 1) (= (c (- k n)) 1)))
  (set! k (+ k 1)))
 (if (and (= (a k) 0) (= (b (+ k n)) 0)
  (= (c (- k n)) 0))
  (begin (place k)
   (while (= (a (+ s 1)) 1) (set! s (+ s 1)))  
   (while (= (c (- r 1)) 1) (set! r (- r 1))))
  (begin (place (+ n t 1)) (set! t (+ t 1))))
 (set! n (+ n 1)))
qs

=>(0 2 4 1 3 8 10 12 14 5 7 18 6 21 9 24 26 28 30 11 13
  34 36 38 40 15 17 44 16 47)	  
変数sは青線を, rは下の赤線を, tは右の赤線を示す. これは最初のTAOCPのプログラムと 同じである. 変数kの使い方も同様.

下の図は, 上のと同じだが, nの範囲を広げた. また配列を調べた位置を白丸で示した. 白丸の数が少ないことに注目して欲しい. TAOCPのプログラムは判定の方法が違うので, 白丸の数はこれよりかなり多い.