数式お絵描き

前回のブログでこうもりの描き方がすっかり分ったかといえば, まだすっきりしない部分がある. 下の線を描くプログラムの11から15行目が思いのほか複雑なのだ.

この式の構造はこうなっている. 11行目の(1/2)*に掛けられる因子は4つの項で出来ている.

3*Sqrt[1-(x/7)^2]            (楕円)
+Sqrt[1-(Abs[Abs[x]-2]-1)^2] (4つの小楕円)
+Abs[x/2]                    (折線)
-((3*Sqrt[33]-7)/112)*x^2    (放物線)
-3

最初はお馴染の楕円で高さは3にしてある. 次が4個の楕円で高さはどれも1だ. その次が(0,0)で折り返す直線で, 最後がx²の抛物線である.



4個の楕円はx=-4,-2,0,2,4の5個所でy=0である. このうち中央のx=0では, 他の項の楕円の3を足し, 最後の3を引き, 折線と放物線はx=0で0だから, 0のままである.

x=4では楕円は3√33/7, 折線は2, 放物線は16*(3√33-7)/112, それから3を引くが, これは計算すると結局0になる.

x=2のところは面倒になったので, Schemeで計算した.

(define (a x)
 (define (sq x) (* x x))
 (+ (* 3 (sqrt (- 1 (sq (/ x 7)))))
    (abs (/ x 2))
    (- (* (/ (- (* 3 (sqrt 33)) 7) 112) (sq x)))
    (- 3)))

(a 2) => .5094556875135123
(a 3) => .3881737849967646
(a 1) => .3778577420858067

という次第で, 0,±4では0, ±2では0.5くらいになる. 1と3では1.37程度になろうか.



この図を見ると高さの値は微妙に違うようだが, まぁいいか. この図では±4の外側では定義されていないということで, 青線がないが, 前回のフィルタを掛けた下の図では

のように外側の方もよくみるとx軸上に青の線が見える. 高さも2倍になっているのが分かる. これから楕円の下半分を引くと, 思ったとおりの下の曲線が出来る.

ところでこの下の図は, 前回のこうもりの絵とはすこし違っている. 最後のAspectRatio->Automaticをつけなかったのである.



AspectRatioの詳細はまだわからないが, Plot関数では, AspectRatio->Automatic がないと縦横比が黄金比のφになるらしい.

ディスプレイ上で測ってみると, 左右は56ミリ, 上下は33ミリであった.

(/ 56. 33) => 1.696969696969697

確かに黄金比に近いといえなくはない. x軸の0から2と, y軸の0から2の長さは明らかに違う. Automaticを指定した前回のこうもりの絵では, 横:縦は正しく7:3になっている.

黄金比にするような見栄えを工夫するのがAutomaticの機能のように思うのだが, Mathematicaの文化はこういうものか. 前回のSinの図でも, Automaticがないので, 横/縦は1.7くらい. Automaticを追加すると, 縦と横の長さは1:1になり, もっと平べったい, 正確になった. (世の中のSinカーブでは, 前回の図も含めてx=0での勾配が1になっていない図がかなり多い.)

数式お絵描き

Raspberry PiでMathematicaが使えるのはうれしい.

そのアナウンスにこうもりのような絵があった. もちろんMathematicaのPlot関数で描いたものだ.



その描画プログラムも添えてあったので, Mathematica屋の絵を描く方法が分るのではないかと期待してプログラムを解読した. とりあえずは目がくらくらするが,

である.

そもそもMathematicaで標準関数の絵を描くのは簡単だ. "Plot["と書き, 描くべき式を書き, ","で区切り, "{x,-Pi,Pi}"のように変数の範囲を書き, "]"で閉じる.

Plot[Sin[x],{x,-Pi,Pi}]

こうもりの絵のようなのは, 上の線と下の線を{,}で囲み区切ればよい.

Plot[{Sin[x],Cos[x]},{x,-Pi,Pi}]

こうもりプログラムでは, 1行目の"With"から10行目の"]]"までが上の曲線, 11行目から15行目の"2]"までが下の曲線である.

上の曲線のプログラムは"With[{"から w(3行目), l(2,3行目), h(4,5,6行目), r(7,8行目)を局所定義し(8行目の"},"で終る), それらを使い9,10行目の本体で描画する.

w(3*Sqrt[1-(x/7)^2])はこういう曲線だ.

これは要するに長軸(半径)7, 短軸1のx²/49+y²=1なる楕円の, 平方根は正をとるからその上半分を描いたわけである. あとから縦軸を3倍した.

次はl.(6/7)*Sqrt[10]+(3+x)/2-(3/7)*Sqrt[10]*Sqrt[4-(x+1)^2] のような形だ. 下に左右反対のrがある. xで変る項は(3+x)2とSqrt[4-(x+1)^2]で, 前者は斜線, 後者は例によって楕円である.

大きさを合せて合成するとこうなる. 上の図の斜線はxの範囲-7から7まであるが, 下の図では引くべき楕円が定義されている範囲が-3から1までなので, その外の 線は消えている.



4行目のhはこうもりの耳のところで, ( (1/2)*(3*(Abs[x-1/2]+Abs[x+1/2]+6)-11*(Abs[x-3/4]+Abs[x+3/4])))といろいろの斜線を組合せているに過ぎない.

ここまで準備したところで, 上の曲線は-7から-3,-3から-1,-1から1,1から3,3から7の区間をUnitStepで接続する. UnitStep[x]はx<0では0, 0≤xでは1になる関数で, UnitStep[x+3]とすると, -3から右が1になる.

w+(l-w)*UnitStep[x+3]+(h-l)*UnitStep[x+1]+(r-h)*UnitStep[x-1]+ (w-r)*UnitStep[x-3]]

x=-7では4個のどのUnitStepも0だから最初の項wで描きだす. x=-3になるとUnitStep[x+3]が1になり, それに(l-w)が掛るから, wが消えlの曲線を描き始める. x=-1 でlが終りhが始まる. x=1でrに, x=3で(w-r)があるので, rも終り, wが復活して最後を描く.


下の線もなかなか凝っているが, 全体的には上と同じ楕円に, 4個の楕円が組み合わさっている構造だ.
ここで一番主要な線はSqrt[1-(Abs[Abs[x]-2]-1)^2]だ.



すぐ上の図は左からAbs[x], Abs[x]-2, Abs[Abs[x]-2], Abs[Abs[x]-2]-1で, ここでそれぞれの斜線部分に対応して例の楕円が4個できる仕掛けである.



UnitStepの代わりにこちらではバンド幅フィルタのような関数でこの4個の楕円を取り出している. それが((x+4)/Abs[x+4]-(x-4)/Abs[x-4])で, 左がx<-4で-1, -4<xで1, 右がx<4で-1, 4<xで1だから, 結局x<-4で0, -4<x<4で2, 4<xで0という形が得られる.



これで4個の楕円を切り出しておき, それから最初と同じ楕円を引いて下の線が描ける.



勿論面倒なのでは区間の切れ目で位置を合わせることで, Sqrt[33]などはそういう計算の結果であろうが, プログラムの骨組みの説明では無視させて貰った. 面白かった.

SATソルバ

Algorithm7222Dも動くようになった. 充足可能なのには可能にする値を, 不可能なのには`unsatisfiable'を出力する.

こんなにgoto文だらけのプログラムがよく動くと感心するが, それがこの本のアルゴリズムの特徴だ. とにかくプログラムが動くようになれば, その機能が理解できるという期待から, 私はいつも実装に心掛けている.

今回お目にかけるこのプログラムも, 中に説明は書いてないが, TAOCPの記述に忠実に実装してあるので, 本文と読み較べてほしい. もっともその方の説明も分り難いが.

プログラムリストについているのは, Waerden(3,3;9)用のデータである.

(define (and1 n)
 (if (even? n) 0 1))
(define (xor1 n)
 ((if (even? n) + -) n 1))
(define (rsh1 n)
 (fix:lsh n -1))
(define (iverson b) (if b 1 0))

;Waerden(3,3;9)用のデータ
(define n 9) (define m 32) ;waerden(3,3;9)
(define ls '(3 11 19 7 13 19 5 11 17 3 9 15 11 15 19 9 13 17
7 11 15 5 9 13 3 7 11 15 17 19 13 15 17 11 13 15 9 11 13 7 9
11 5 7 9 3 5 7 2 10 18 6 12 18 4 10 16 2 8 14 10 14 18 8 12 
16 6 10 14 4 8 12 2 6 10 14 16 18 12 14 16 10 12 14 8 10 12 
6 8 10 4 6 8 2 4 6) )
(define ws '
(0 0 1 17 2 18 3 19 4 20 5 21 6 22 7 23 0 0 0 0))
(define link '
(0 8 9 10 11 12 0 0 13 14 15 0 0 16 0 0 0 24 25 26 27 28 0 0
29 30 31 0 0 32 0 0 0))
(define start '(96 93 90 87 84 81 78 75 72 69 66 63 60 57 54
51 48 45 42 39 36 33 30 27 24 21 18 15 12 9 6 3 0))

;途中の様子を出力する関数
(define (printnext next h t)
 (let ((ring (list)))
  (define (ploop)
   (if (eq? h t) (display (reverse (cons t ring)))
    (begin (set! ring (cons h ring))
     (set! h (list-ref next h)) (ploop))))
  (ploop)))
(define (printcl j)
 (define (pin x)
  (string->symbol (string-append
   (number->string (quotient x 2))
   (list-ref '("+" "-") (and1 x)))))
 (display (map (lambda (x) (pin (list-ref ls x))) 
  (a2b (list-ref start j) (list-ref start (- j 1))))))
(define (printx xs) 
 (display (map (lambda (x) (if (< x 0) '- x)) xs)))

(define d 0) (define h 0) (define t 0) (define f 0) 
(define k 0) (define p 0) (define b 0) (define l 0) 
(define j 0) (define j1 0) (define i 0) (define l1 0) 
(define q 0)
(define xs (map (lambda (x) 0) (a2b 0 (+ n 1))))
(define next (map (lambda (x) 0) (a2b 0 (+ n 1))))
(define hs (map (lambda (x) 0) (a2b 0 (+ n 1))))
(define ms (map (lambda (x) 0) (a2b 0 (+ m 1))))

(call-with-current-continuation
 (lambda (exit)

(define (d1)
 (list-set! ms 0 0)
 (set! d 0)
 (set! h 0)
 (set! t 0)
 (do ((k n (- k 1))) ((= k 0)) 
  (list-set! xs k -1)
  (if (or (not (= (list-ref ws (* 2 k)) 0))
          (not (= (list-ref ws (+ (* 2 k) 1)) 0)))
  (begin
   (list-set! next k h) (set! h k)
   (if (= t 0) (set! t k)))))
 (if (not (= t 0)) (list-set! next t h))
(d2))

(define (d2) 
 (newline) (if (not (= t 0)) (printnext next h t)) 
 (printx (cdr xs)) ;xsを出力
 (if (= t 0) (exit 'satisfiable)
  (begin (set! k t) (d3))))

(define (ex135 l)
 (let ((j (list-ref ws l)))
 (define (un0)
  (set! p (+ (list-ref start j) 1))
   (un1))
 (define (un1)
  (if (= p (list-ref start (- j 1))) 1 (un2)))
 (define (un2)
  (if (= (list-ref xs (rsh1 (list-ref ls p)))
         (and1 (list-ref ls p)))
   (begin (set! p (+ p 1)) (un1)) (un3)))
 (define (un3)
  (set! j (list-ref link j))
  (if (= j 0) 0 (un0)))
 (if (not (= j 0)) (un0) 0)))

(define (d3)
 (set! h (list-ref next k))
 (let ((f0 (ex135 (* 2 h))) (f1 (ex135 (+ (* 2 h) 1))))
 (set! f (+  f0 (* 2 f1))) )
 (cond ((= f 3) (display (list 'un 
   (string-append (number->string h) "+")
   (string-append (number->string h) "-"))) (d7))
       ((or (= f 1) (= f 2))
        (display (list 'un (string-append (number->string h) 
         (list-ref '(0 "+" "-") f)))) 
        (list-set! ms (+ d 1) (+ f 3)) (set! t k) (d5))
       ((not (= h t)) (set! k h) (d3))
       (else (d4))))

(define (d4)
 (set! h (list-ref next t))
 (list-set! ms (+ d 1)
  (iverson (or (= (list-ref ws (* 2 h)) 0)
               (not (= (list-ref ws (+ (* 2 h) 1)) 0)))))
 (d5))

(define (d5)
 (set! d (+ d 1))
 (list-set! hs d h)
 (set! k h)
 (if (= t k) (set! t 0)
   (begin (list-set! next t (list-ref next k))
    (set! h (list-ref next k))))
 (d6))

(define (ex136)
  (define (ud0)
   (set! j1 (list-ref link j))
   (set! i (list-ref start j))
   (set! p (+ i 1))
   (ud1))
  (define (ud1)
   (if (= (list-ref xs (rsh1 (list-ref ls p)))
          (and1 (list-ref ls p)))
    (begin (set! p (+ p 1)) (ud1))
    (ud2)))
  (define (ud2)
   (set! l1 (list-ref ls p))
   (list-set! ls p l)
   (list-set! ls i l1)
   (ud3))  
  (define (ud3)
   (set! p (list-ref ws l1))
   (set! q (list-ref ws (xor1 l1))) 
   (if (or (not (= p 0)) (not (= q 0)) 
       (>= (list-ref xs (rsh1 l1)) 0))
    (ud5) (ud4)))
  (define (ud4)
   (if (= t 0) (begin (set! t (rsh1 l1)) (set! h (rsh1 l1))
      (list-set! next t h))
    (begin (list-set! next (rsh1 l1) h) (set! h (rsh1 l1))
      (list-set! next t h)))
    (ud5))
  (define (ud5)
   (list-set! link j p) 
   (list-set! ws l1 j)
   (ud6))
  (define (ud6)
   (printcl j)
   (set! j j1)
   (if (not (= j 0)) (ud0)))
 (set! l (+ (* 2 k) b))
 (set! j (list-ref ws l))
 (list-set! ws l 0)
 (if (not (= j 0)) (ud0)))

(define (d6)
 (set! b (modulo (+ (list-ref ms d) 1) 2))
 (list-set! xs k b)
 (ex136) 
 (d2))

(define (d7loop)
 (if (>= (list-ref ms d) 2) (begin 
  (set! k (list-ref hs d))
  (list-set! xs k -1)
  (if (or (not (= (list-ref ws (* 2 k)) 0))
          (not (= (list-ref ws (+ (* 2 k) 1)) 0)))
   (begin (list-set! next k h)
    (set! h k)
    (list-set! next t h)))
  (set! d (- d 1)) (d7loop))))
   
(define (d7)
 (set! t k)
 (d7loop)
 (d8))

(define (d8)
 (if (> d 0) (begin
  (list-set! ms d (- 3 (list-ref ms d)))
  (set! k (list-ref hs d)) (d6))
 (exit 'unsatisfiable)))

(d1)))

このアルゴリズムの出力を, TAOCPにあるのと同じ場所まで書いてみるとこうなる. (プログラムの実際の出力はもっと冗長だが, 不要な空白など手でけしてある.) 左からアクティブリング, xの値, unはユニットになったリテラル(これが同じ変数の正負のリテラルだとバックトラックに入る,) 書き換えた 項を示す.

(1234567) - - - - - - - - - 2+1+3+ 3+1+5+ 4+1+7+ 5+1+9+
(234567)  0 - - - - - - - - 3+1+2+ 3+2+4+ 4+2+6+ 5+2+8+
(34567)   0 0 - - - - - - -(un 3+) 4-3-5- 5-3-7- 6-3-9-
(4567)    0 0 1 - - - - - - 6+2+4+ 7+1+4+ 5+4+6+ 6+4+8+
(567)     0 0 1 0 - - - - -(un 6+) 9-3-6- 7-6-8-
(975)     0 0 1 0 - 1 - - -(un 9-)
(75)      0 0 1 0 - 1 - - 0(un 7+) 8-6-7- 8-7-9-
(85)      0 0 1 0 - 1 1 - 0(un 8-)
(5)       0 0 1 0 - 1 1 0 0(un 5+5-) 5-3-4- 5-4-6- 6-4-8-
(69785)   0 0 1 1 - - - - -(un 5-) 4+5+6+ 8+2+5+ 9+1+5+ 
                                   6+5+7+ 7+5+9+
(6978)    0 0 1 1 0 - - - -(un 9+) 6-3-9-

3行目の2番目の項はTAOCPのと違う. もちろん向こうが違っている. クヌース先生には連絡済みだ.

SATソルバ

アルゴリズムBでLangford(3)を実行したのが次だ.

Langford(3)だから, 1,2,3を2個ずつ, 例えば323121のように並べる. 但し2個のnの間には他の数がn個なければならないという制限があるから, 先ほどのは駄目で, 231213がひとつの解である.

このくらいならちょっとやっただけで出来るが, 長くなる(nが増える)と大事だ.

TAOCPでの方法はこうだ. まずこういう表を作る.


1行目 1を1番と3番に置く. 間隔は1だ. 2行目 1を2番と4番に置く. ...

5行目 2を1番と4番に置く. ...

8行目 3を1番と5番に置く.

3を2番と6番に置くというのもありだが, 対称だから8行目だけにする.

表が出来たらd₁の列に1のあるx_jから1つ選ぶ. つまりx₁, x₂, x₃, x₄から1つ選ぶ. d₂からもx₅, x₆, x₇から1つ選ぶ. これをs₆の列まで行う.

1つ選ぶ式をTAOCPではS₁と表記するから, この関係は

の式になり, S₁の展開式を当てはめると, SAT


が得られる. 2-4-, 1-3- が重複しているからそれを除き, リテラルの内部表現にすると

(define sat '((2 4 6 8) (3 5) (3 7) (3 9) (5 7) (5 9) (7 9)
(10 12 14) (11 13) (11 15) (13 15) (16) 
(2 10 16) (3 11) (3 17) (11 17) (4 12) (5 13) (2 6 14) (3 15)
 (7 15) 
(4 8 10) (5 11) (9 11) (6 12 16) (7 13) (7 17) (13 17) (8 14)
 (9 15)))

となる. 従って

(define ls '(9 15 8 14 13 17 7 17 7 13 6 12 16 9 11 5 11 4 8 
10 7 15 3 15 2 6 14 5 13 4 12 11 17 3 17 3 11 2 10 16 16 13
 15 11 15 11 13 10 12 14 7 9 5 9 5 7 3 9 3 7 3 5 2 4 6 8))

前回と同様にsatの先頭だけとると, sが

(0 2 3 3 3 5 5 7 10 11 11 13 16 2 3 3 11 4 5 2 3 7 4 5 9 
 6 7 7 13 8 9)

になり, zも作れ, wやlinkも得られる.

((0) () (19 13 1) (20 15 14 4 3 2) (22 17) (23 18 6 5) 
(25) (27 26 21 7) (29) (30 24) (8) (16 10 9) () (28 11) 
() () (12) ())

ws=
(0 0 1 2 17 5 25 7 29 24 8 9 0 11 0 0 12 0)

link=
(0 13 3 4 14 6 18 21 0 10 16 28 0 19 15 20 0 22 23 0 0 26 0 0
 30 0 27 0 0 0 0)

節の長さがばらばらだから, STARTも作らなければならない.

(define a (map length sat))
a => (4 2 2 2 2 2 2 3 2 2 2 1 3 2 2 2 2 2 3 2 2 3 2 2 3 2
 2 2 2 2)

(define (ac ls)
 (if (null? ls) '(0)
  (let ((a (ac (cdr ls))))
    (cons (+ (car ls) (car a)) a))))

(define start (ac a))
start =>
(66 62 60 58 56 54 52 50 47 45 43 41 40 37 35 33 31 29 27 24
 22 20 17 15 13 10 8 6 4 2 0))

これで準備が完了し, アルゴリズムBを実行すると01000011が得られ, x₂,x₇,x₈を取るのが解と分かり, 泰山鳴動の一幕であった.

SATソルバ

TAOCP 7.2.2.2にあるSATソルバのうち, 今回はAlgorithm 7222Bの話である.

これも結構難しかった. 演習問題128の解答の(7)の式の初期値が頼りであった.

Algorithm Aに較べると, F, B, C, SIZEの配列はなくなり, 代わりにWとLINKが新たに出来た.

SATの問題では, リテラルは節のいろいろな位置に現れるが, このアルゴリズムでは, 節の先頭に現れるものだけについて初期値に設定する. Wはそれぞれのリテラルについて, 最初に現れた節の番号を示す. どの先頭にも現れなければ, 終という意味で0を入れる.

同じリテラルが複数の節の先頭に現れた時は, 続きをLINKで示す. 具体的には下の例を見てほしい.

7.2.2.2の初めの方にSimple Exampleというのがある. 8ビットのx₁...x₈で3個の0の間隔が等しくなく, 3個の1の間隔も等しくないものを探せというのである.

すぐ下に9ビットの場合のSATの式があるから, それの8ビット版を書けばよいわけで, 実装に適した私流の記法では

1+2+3+,2+3+4+,3+4+5+,4+5+6+,5+6+7+,6+7+8+,
1+3+5+,2+4+6+,3+5+7+,4+6+8+,
1+4+7+,2+5+8+,
1-2-3-,2-3-4-,3-4-5-,4-5-6-,5-6-7-,6-7-8-,
1-3-5-,2-4-6-,3-5-7-,4-6-8-,
1-4-7-,2-5-8-,

つまり8変数, 24節になる. (n=8, m=24)

これをリテラルの内部表現にすると

(define n 8)

(define sat '
((2 4 6) (4 6 8) (6 8 10) (8 10 12) (10 12 14) (12 14 16) 
(2 6 10) (4 8 12) (6 10 14) (8 12 16) (2 8 14) (4 10 16) 
(3 5 7) (5 7 9) (7 9 11) (9 11 13) (11 13 15) (13 15 17) 
(3 7 11) (5 9 13) (7 11 15) (9 13 17) (3 9 15) (5 11 17)))

(define m (length sat)) => 24

配列Lはこれをreverseし, appendしたものである.

先頭の節を1番として, 2から17までのリテラルがどの節の先頭に現れるかをみると (listは0から始まるのでダミーの0をconsしてある)

(define s (cons 0 (map car sat)))
=>
(0 2 4 6 8 10 12 2 4 6 8 2 4 3 5 7 9 11 13 3 5 7 9 3 5)

リテラル2は1番の節に, 4は3番の節に現れる. リテラル l が現れる節番号のリストを l 番目にもつリストzを作る.

(define z (map (lambda (x) '()) (a2b 0 (+ (* n 2) 2))))

z => (() () () () () () () () () () () () () () () () () ())

(for-each (lambda (x y)
 (list-set! z x (cons y (list-ref z x)))) s (a2b 0 (+ m 1)))

z=>
((0) () (11 7 1) (23 19 13) (12 8 2) (24 20 14) (9 3) (21 15)
 (10 4) (22 16) (5) (17) (6) (18) () () () ())

リテラル2は1,7,11の節に, リテラル3は13,19,23の節に現れるということだ. これからWとLINKを作る.

(define ws (map (lambda (x) 0) (a2b 0 (+ (* 2 n) 2))))

(for-each (lambda (x y) 
(if (not (null? x)) (list-set! ws y (car (reverse x)))))
 z (a2b 0 (+ (* 2 n) 2)))

ws => (0 0 1 13 2 14 3 15 4 16 5 17 6 18 0 0 0 0)

(define link (map (lambda (x) 0) (a2b 0 (+ m 1))))

(for-each (lambda (x) (let ((t 0))
(map (lambda (y) (list-set! link y t) (set! t y)) x))) z)

link =>
(0 7 8 9 10 0 0 11 12 0 0 0 0 19 20 21 22 0 0 23 24 0 0 0 0)

この解は00110011, 01011010,01100110,10011001,10100101,11001100だが, このアルゴリズムは最初の解(00110011)しか返さない.

(define (algorithm7222b)
(define (and1 n) (if (even? n) 0 1))
;1とのandをとる
(define (xor1 n) ((if (even? n) + -) n 1))
;1とのxorをとる. リテラルの正負を反転する
(define (rsh1 n) (fix:lsh n -1))
;1ビット右シフトする
(define (? b) (if b 1 0))
;iverson notation [b]=b?1:0

;simple example
(define n 8) (define m 24)
(define sat '
((2 4 6) (4 6 8) (6 8 10) (8 10 12) (10 12 14) (12 14 16) 
(2 6 10) (4 8 12) (6 10 14) (8 12 16) (2 8 14) (4 10 16) 
(3 5 7) (5 7 9) (7 9 11) (9 11 13) (11 13 15) (13 15 17) 
(3 7 11) (5 9 13) (7 11 15) (9 13 17) (3 9 15) (5 11 17)))
(define ls (apply append (reverse sat)))
(define ws '(0 0 1 13 2 14 3 15 4 16 5 17 6 18 0 0 0 0))
(define link '(0 7 8 9 10 0 0 11 12 0 0 0 0 19 20 21 22 0 0 
 23 24 0 0 0 0))
(define start '(72 69 66 63 60 57 54 51 48 45 42 39 36 33 30
 27 24 21 18 15 12 9 6 3 0))
(define ms '(0 0 0 0 0 0 0 0 0))

(define d 0) (define j 0) (define i 0) (define i1 0)
(define j1 0) (define k 0) (define l 0) (define l1 0)

(call-with-current-continuation
 (lambda (exit)

(define (b1) (set! d 1) (b2))
(define (b2) 
 (if (> d n) (begin (do ((d 1 (+ d 1))) ((> d n))
  (display (xor1 (and1 (list-ref ms d))))) (exit 'ok))
 (begin (list-set! ms d (? (or (= (list-ref ws (* 2 d)) 0)
  (not (= (list-ref ws (+ (* 2 d) 1)) 0))))) 
 (set! l (+ (* 2 d) (list-ref ms d))) 
 (b3))))

(define (b3kloop)
 (cond ((= k i1) (list-set! ws (xor1 l) j) (b5))
  ((< k i1) 
  (set! l1 (list-ref ls k))
  (if (or (> (rsh1 l1) d)
    (even? (+ l1 (list-ref ms (rsh1 l1)))))
   (begin
    (list-set! ls i l1)
    (list-set! ls k (xor1 l))
    (list-set! link j (list-ref ws l1))
    (list-set! ws l1 j)
    (set! j j1))
  (begin (set! k (+ k 1)) (b3kloop))))))

(define (b3jloop)
 (if (not (= j 0)) (begin 
  (set! i (list-ref start j))
  (set! i1 (list-ref start (- j 1)))
  (set! j1 (list-ref link j))
  (set! k (+ i 1))
  (b3kloop)
  (b3jloop))))

(define (b3)
 (set! j (list-ref ws (xor1 l)))
 (b3jloop)
 (b4))

(define (b4)
 (list-set! ws (xor1 l) 0)
 (set! d (+ d 1))
 (b2))

(define (b5)
 (if (< (list-ref ms d) 2) (begin
  (list-set! ms d (- 3 (list-ref ms d)))
  (set! l (+ (* 2 d) (and1 (list-ref ms d)))) 
  (b3))
 (b6)))

(define (b6)
 (if (= d 1) (exit "unsatisfiable")
  (begin (set! d (- d 1)) (b5))))

(b1)))

(algorithm7222b) => 00110011

他の実行例については次回にでも.

軌跡作図器

JavaScriptによるアニメーションができた.



ここにあります.

PがO'を中心とする円上を移動すると, Qは直線上を移動する.

軌跡作図器

先日のブログの最後にJavaScriptで動かしたいと書いた. まだそこまではいっていないが, Pの位置を変えるとQがどう動くが調べたくなって, PostScriptでプログラムを書いた.

方針としてはこうだ. AB, AO, ADの長さを決める. O, O'の位置を決める.

OO'の長さに等しくO'Pをとるわけだが, つまりO'を中心とする半径OO'の円上をPが動くわけだが, Oの反対側をO''とし, ∠O''O'Pを0から順に増やすことで, Pの位置を決めることにした.

OとPの位置やOBとPBの長さが分っているので, Bの位置が決まる. Bが決まるとAが決まり, BとPが分るからCも決まる. AとCからの長さが分っているから, Dの位置が決まり, Qも得られる.

このようにして, ∠O''O'Pを0°, 30°, 60°, 90°, 120°, 135°と変えて描いたのが次の図である.


0°


30°


60°


90°


120°


135°

この先Pはどこまで進めるかというと, OBとPBの長さが決まっているから, OPの間隔がPB-OBになったときに一直線になって以後破綻する.