2015年7月16日木曜日

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のと違う. もちろん向こうが違っている. クヌース先生には連絡済みだ.

2015年7月2日木曜日

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行目だけにする.

表が出来たらd1の列に1のあるxjから1つ選ぶ. つまりx1, x2, x3, x4から1つ選ぶ. d2からもx5, x6, x7から1つ選ぶ. これをs6の列まで行う.

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

の式になり, S1の展開式を当てはめると, 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が得られ, x2,x7,x8を取るのが解と分かり, 泰山鳴動の一幕であった.

2015年7月1日水曜日

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ビットのx1...x8で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
他の実行例については次回にでも.

2015年6月30日火曜日

軌跡作図器

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



ここにあります.

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

2015年6月27日土曜日

軌跡作図器

先日のブログの最後に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°と変えて描いたのが次の図である.





30°


60°


90°


120°


135°

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

2015年6月26日金曜日

軌跡作図器

學士會会報No.912 (2015-III)の表紙の図に興味をもった. 説明には独国マルチン・シリング(ハーレ社)製「軌跡作図器(Serie XXIV, Modell nr.11)」とあるが, それに続く説明は要領を得ないので改めて説明したい.

同じような図をネットで探して頂いてきたのが次の図だ.(http://goo.gl/Fpfcqt)



リンクABとCDの長さは等しく, ADとBCの長さも等しい. A,B,C,Dはピボットで連結されていて角度は変わる. AB, CB, ADのそれぞれをm : (1-m)に内分した点をO,P,Qとする(0 < m < 1).

Oは固定点で, ABはOを支点に回転する. PにPO'がOO'に等しくなるようなリンクPO'を取付け, O'も固定する. つまりPはO'を中心に円を描き, この円はOを通る.

こうしてPが円(の一部)を描くと, Qは直線の軌跡を描くという道具であった(Hartのinversorという名前である).

上の図を次のように書き直す. 太線がリンク機構である. 三角形ACBと三角形CADは逆向きだが合同だから, ACとBDは平行. ABDCは等脚台形になる. OPQは一直線でこれもACに平行.



BD上にA, CからACに垂線を立ててE,Fとする. EB=DFだから
AC×BDはEF×BD=(ED+EB)×(ED-EB)=ED2-EB2.
ED2+AE2=AD2
EB2+AE2=AB2
だから
AC×BC=AD2-AB2
AO/AB=CP/CB=AQ/AD=mだから
OP=m×AC, OQ=(1-m)×BD
従って
OP×OQ=m×(1-m)×(AD2-AB2)=一定.

OP'×OQ'=OP×OQ=一定 でP'がOPの中点O'を中心とする円を動くと, Q'はQでOQに直交する線上を動く.



上の図で∠OP'P=∠R. ∠PQQ'=∠R だから4点PQQ'P'は円周上にあり, 方冪の定理(定点を過ぎる直線が定円と交わるとき, 定点より2つの交点までの長さの積は一定) によりOP'×OQ'=OP×OQ.

そのうちJavaScriptで動かしたいと思う.

2015年6月22日月曜日

SATソルバ

TAOCP 7.2.2.2にはSATソルバのアルゴリズムがいくつか登場する.

今回はそのうちAlgorithm 7222Aを話題としたい. ただこのアルゴリズムの背景を初めから述べるのは大変なので, TAOCPのその辺を読んである読者を対象に書くことで勘弁して貰いたい.

TAOCPのアルゴリズムは自然言語であるだけでなく, 非常に簡潔な記述なので, 前後の説明をよく理解していないと, なんのことだか分からないことが多い. そういう場合にSchemeで記述したこのブログのプログラムは理解の助けになのではないか. とにかく実際に動くプログラムとして存在するから だ.

TAOCPにも, あるデータを例としての説明はある. 7.2.2.2-6と7の式



で, (6)は充足不能, (7)は充足可能な例である.

まず(7)について, 次のようなデータを用意する.



配列L, F, B, Cは添字が0から30で, Lには(7)の節が置いてある. L(10)から始まる9,7,3は(7)の最後の節の3-,4-,1-を4,3,1の順にし, 正のリテラルは変数番号を2倍し, 負のリテラルは2倍して1足したものなので, 9,7,3になっている. Cの対応する場所の7は, 7番目の節を示す.

L(13)以降も同様である.

FとBは両方向リンクで, F(2)つまりリテラル1+が現れる場所はL(30)であり, F(30)=24は次がL(24)であることを示す. Bは逆向きの情報である.

リンク構造なのは, リテラルや節を削除/挿入するのに便利なためである.

C(2)からC(9)まではそれぞれのリテラルが何個あるかを示す.

それらの下のSTARTは, それぞれの節の先頭の位置; SIZEは節にあるリテラルの個数である.

Schemeで実装したプログラムは以下のようだ.


(define (algorithm7222a)
(define (and1 n) (if (even? n) 0 1))
;1とのandをとる
(define (xor1 n) ((if (even? n) + -) n 1))
;1とのxorをとる. リテラルの正負を反転する
(define (? b) (if b 1 0))
;iverson notation [b]=b?1:0

(define n 4) (define m 7)
(define ls '(0 0 0 0 0 0 0 0 0 0 9 7 3 8 7 5 6 5 3 8 4 3
 8 6 2 9 6 4 7 4 2))
(define fs '(0 0 30 21 29 17 26 28 22 25 9 7 3 8 11 5 6 15 12
 13 4 18 19 16 2 10 23 20 14 27 24))
(define bs '(0 0 24 12 20 15 16 11 13 10 25 14 18 19 28 17 23
 5 21 22 27 3 8 26 30 9 6 29 7 4 2))
(define cs '(0 0 2 3 3 2 3 3 3 2 7 7 7 6 6 6 5 5 5 4 4 4
 3 3 3 2 2 2 1 1 1))
(define size '(0 3 3 3 3 3 3 3))
(define start '(0 28 25 22 19 16 13 10))

(define ms (map (lambda (x) 0) (a2b 0 (+ n 1))))

(define l 0) (define a 0) (define d 0) (define p 0) 
(define i 0) (define j 0) (define q 0) (define r 0) 
(define s 0)

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

(define (a1)
 (set! a m) (set! d 1)
 (a2))
(define (a2)
 (set! l (* 2 d))
 (if (<= (list-ref cs l) (list-ref cs (+ l 1)))
  (set! l (+ l 1)))
 (list-set! ms d (+ (and1 l) 
  (* 4 (? (= (list-ref cs (xor1 l)) 0)))))
 (newline) (display (list 'ms (take (+ d 1) ms)))
 (if (= (list-ref cs l) a) (begin (newline)
  (do ((i 1 (+ i 1))) ((> i n)) 
   (display (and1 (xor1 (list-ref ms i)))))
  (exit 'satisfiable)))
 (a3))

(define (a3back)
 (if (>= p (+ n n 2)) (begin
  (set! j (list-ref cs p))
  (set! i (list-ref size j))
  (list-set! size j (+ i 1))
  (set! p (list-ref bs p))
  (a3back)) (a5)))

(define (a3loop)
 (if (>= p (+ n n 2)) (begin
 (set! j (list-ref cs p))
 (set! i (list-ref size j))
 (if (> i 1) (begin
   (list-set! size j (- i 1)) 
   (set! p (list-ref fs p)) (a3loop))
  (begin (set! p (list-ref bs p)) (a3back))))
(a4)))

(define (a3)
 (set! p (list-ref fs (xor1 l)))
 (a3loop))

(define (a4loop)
 (if (>= p (+ n n 2)) (begin
  (set! j (list-ref cs p))
  (set! i (list-ref start j))
  (set! p (list-ref fs p)) 
  (do ((s i (+ s 1))) ((= s (+ i (list-ref size j) -1)))
   (set! q (list-ref fs s))
   (set! r (list-ref bs s))
   (list-set! bs q r)
   (list-set! fs r q)
   (list-set! cs (list-ref ls s)
    (- (list-ref cs (list-ref ls s)) 1)))
  (a4loop))))

(define (a4)
 (set! p (list-ref fs l))
 (a4loop)
 (set! a (- a (list-ref cs l))) (set! d (+ d 1)) 
 (a2))

(define (a5)
 (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)))) (a3))
 (a6)))

(define (a6)
 (if (= d 1) (exit "unsatisfiable")
  (begin (set! d (- d 1))
   (set! l (+ (* 2 d) (and1 (list-ref ms d)))) (a7))))

(define (a7loop)
 (if (>= p (+ n n 2)) (begin
  (set! j (list-ref cs p))
  (set! i (list-ref start j))
  (set! p (list-ref bs p))
  (do ((s i (+ s 1))) ((= s (+ i (list-ref size j) -1)))
   (set! q (list-ref fs s))
   (set! r (list-ref bs s))
   (list-set! bs q s)
   (list-set! fs r s)
   (list-set! cs (list-ref ls s)
    (+ (list-ref cs (list-ref ls s)) 1)))
 (a7loop))))

(define (a7)
 (set! a (+ a (list-ref cs l)))
 (set! p (list-ref bs l))
 (a7loop)
 (a8))

(define (a8loop)
 (if (>= p (+ n n 2)) (begin
  (set! j (list-ref cs p))
  (set! i (list-ref size j))
  (list-set! size j (+ i 1))
  (set! p (list-ref fs p))
  (a8loop))))

(define (a8)
 (set! p (list-ref fs (xor1 l)))
 (a8loop)
 (a5))

(a1)))
)
(algorithm7222a)
途中, ステップA2で動作状況を示すmsを印字しながらの実行結果は
(ms (0 1))
(ms (0 1 0))
(ms (0 1 0 1))
(ms (0 1 0 1 4))
0101
;... done
;Value: satisfiable
で, 解はx1=0, x2=1, x3=0, x4=1 である.

TAOCPには(6)の実行中のmsの値も掲載してあるから
(define n 4) (define m 8)
(define ls '(0 0 0 0 0 0 0 0 0 0 9 5 2 9 7 3 8 7 5 6 5 3 8 4
 3 8 6 2 9 6 4 7 4 2))
(define fs '(0 0 33 24 32 20 29 31 25 28 9 5 2 10 7 3 8 14 11
 6 18 15 16 4 21 22 19 12 13 26 23 17 30 27))
(define bs '(0 0 12 15 23 11 19 14 16 10 13 18 27 28 17 21 22
 31 20 26 5 24 25 30 3 8 29 33 9 6 32 7 4 2))
(define cs '(0 0 3 3 3 3 3 3 3 3 8 8 8 7 7 7 6 6 6 5 5 5 4 4
 4 3 3 3 2 2 2 1 1 1))
(define size '(0 3 3 3 3 3 3 3 3))
(define start '(0 31 28 25 22 19 16 13 10))
として(6)もやってみると
(ms (0 1))
(ms (0 1 1))
(ms (0 1 1 0))
(ms (0 1 1 3 1))
(ms (0 1 2 1))
(ms (0 1 2 1 1))
(ms (0 1 2 2 1))
(ms (0 2 1))
(ms (0 2 1 1))
(ms (0 2 1 1 1))
(ms (0 2 1 2 1))
(ms (0 2 2 1))
(ms (0 2 2 2 1))
;... done
;Value 11: "unsatisfiable"
のように, TAOCPにあるのと同様な出力が得られた. まぁこの実装は間違ってはいないようだ.