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

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

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

で, 解はx₁=0, x₂=1, x₃=0, x₄=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にあるのと同様な出力が得られた. まぁこの実装は間違ってはいないようだ.

SATソルバ

TAOCP 7.2.2.2は充足可能性(Satisfiability)が話題だ. 最初の方に多くの問題が CNF(Conjunctive Normal Form 乗法標準形)で記述できるという話が続く. 私はこの分野にはさほど興味がなく, まったくのど素人で, いつものようにアルゴリズ ムにだけ関心がある. 繰り返すまでもなく, TAOCPのアルゴリズムは自然言語で 記述されているから, その英文解釈が煩わしく, いつもああでもない, こうでもない と思いつつSchemeに変換し, いくつかの例題が無事に解けるのをみてやっと 安心している.

今回もそういう話である.

早速例をTAOCPから拝借すると, 乗法標準形は次のようなものである.


x₁, x₂, ...などが変数(variable). 変数または その上に線を引いたものがリテラル(literal). リテラルを∨でつないだものが 論理和節(clause), そしてCNFはそれらの節を∧でつないだものである.

かっこや∧や∨をいちいち書くのも面倒だということで, TAOCPではこれを


のように書いた.

さてSATソルバはFが真になるように, これらの変数に真(1), 偽(0)を対応させる もので, この例のサイズなら目の子で探せば解が見つかるが, 実用になる問題では 変数の個数がすごく多く, 高速なアルゴリズムが研究対象になっているらしい. ( 4月23日のブログMcGregorグラフを解く場合は440変数の1406論理和節になるとか.)

TAOCPにはまず基本的なアルゴリズムが出ていたので, 練習としてそれを Schemeで実装してみた.

まず上のCNFの例は, 私の実装ではS式で((1+ 2-) (2+ 3+) (1- 3-) (1- 2- 3+)) と書くことにする. 負のリテラルに-がつくだけでなく, 正のリテラルに+が つくのは対称性からきれいではないか. Schemeでは+1は数だが, 1+はシンボルになる. 従って 1+を1-にしたり, その逆にしたりは

(define (neg l)
(let* ((s (symbol->string l)) (w (- (string-length s) 1)))
 (string-set! s w (integer->char 
  (- 88 (char->integer (string-ref s w)))))
 (string->symbol s)))

つまりリテラルlを貰い, その文字列sとその長さ-1のwを 作る. 最後の+か-を取り出し, ASCIIの値とし,

(char->integer #\+) => 43
(char->integer #\-) => 45

だからその和の88から引くと+, -が逆転できる. ついでだが, 変数の記号からリテラルを作る(例えば x からx+, x-, 2 から2+, 2-を作る)手続きは次のようだ.

 (define (lit v s)
  (string->symbol (string-append
   ((if (number? v) number->string symbol->string) v)
   (symbol->string s))))
(lit 'x '+) => x+
(lit 'x '-) => x-
(lit 2 '+) => 2+
(lit 2 '-) => 2-

TAOCPにあるSATソルバの基本的なアルゴリズムB(F)(7222-56)は次のようだ.





考えてみると, このアルゴリズムは当然と思える. Fは節がandで繋がってい るから, 空のandは真だ. (and) => #t.

andで繋がっている節に空のものがあると, 空の節は(or) => #f だから Fは偽になる.

そのどちらでもないなら, まだFに残っている変数の正と偽のリテラル を作り, それで「ただし」以下のF|lを用意してBに 再帰的に作用させ, その結果にリテラルを足したのを返す.

「ただし」の内容はこうだ. lを含む節は真なので, andの並びでは 考慮する必要がないから削除する. 一方lの否定を含む節はorに 偽の項は関係ないから項を削除する. いずれにしろlの変数は 消えることになる.

そうだと思って書いたのが次のSchemeのプログラムである.

(define (b f vs)
 (define (lit v s)
  (string->symbol (string-append
   ((if (number? v) number->string symbol->string) v)
   (symbol->string s))))
 (let ((a '()))
  (define (b1 f ls vs)
   (display (list 'b1 f ls vs )) (newline)
   (cond ((member '() f) #f)
         ((= (length (nub f)) 1)
          (set! a (cons (cons (caar f) ls) a)))
         (else
         (let ((l (lit (car vs) '+)))
           (b1 (apb f l) (cons l ls) (cdr vs)))
         (let ((l (lit (car vs) '-)))
           (b1 (apb f l) (cons l ls) (cdr vs))))))
 (b1 f '() vs) a))

引数のvsは変数のリストである.

だから上の例では

(b '((1+ 2-) (2+ 3+) (1- 3-) (1- 2- 3+)) '(1 2 3))

と呼び出す. 内部手続きb1が途中経過を示すので, 出力は

(b1 ((2+ 3+) (3-) (2- 3+)) (1+) (2 3))
(b1 ((3-) (3+)) (2+ 1+) (3))
(b1 (()) (3+ 2+ 1+) ())
(b1 (()) (3- 2+ 1+) ())
(b1 ((3+) (3-)) (2- 1+) (3))
(b1 (()) (3+ 2- 1+) ())
(b1 (()) (3- 2- 1+) ())
(b1 ((2-) (2+ 3+)) (1-) (2 3))
(b1 (()) (2+ 1-) (3))
(b1 ((3+)) (2- 1-) (3))
=> ((3+ 2- 1-))

トレースの最初の行は, 後ろから見て残りの変数は2と3, リテラルは1+. 従って((1+ 2-) (2+ 3+) (1- 3-) (1- 2- 3+))の左端の項は1+があるから削除. 次は1+, 1-もないからそのまま, 次の2項は1-を削除. そして((2+ 3+) (3-) (2- 3+))を得る.

次の行はリテラルが2+. そこで(2+ 3+)は削除. (2- 3+)の 2-を削除. ((3-) (3+))が残る.

次. 変数3に対して3+, 3-をリテラルとしても()が得られ. εなので 充足不能になる. バックトラックして2-をリテラルにしても失敗.

改めて1-をリテラルにしてみたところ, fが空リストになり成功した.

原理的にはこれでよいが, 実用的には工夫があるらしい. そういう改良 アルゴリズムの話がTAOCPでは続く. その説明はまたいずれ.