Verifikation von Algorithmen mit Z3 – Teil 2

Im ersten Artikel dieser Reihe hatten wir den SMT-Solver Z3 kennen gelernt. Wir hatten in der zugehörigen Sprache SMT-LIB-2 eine Spezifikation für die Rasterisierung von Geraden verfasst und dann einen verbesserten Algorithmus entworfen, der mit dem Simplen in seiner Funktionsweise übereinstimmt. Während die Korrespondenz dieser beiden Algorithmen manchen Lesern vielleicht noch offensichtlich erschien, gehen wir in diesem Artikel einen signifikanten Schritt weiter. Wir implementieren den sog. Bresenham-Algorithmus, bei dem auf den ersten Blick ganz und gar nicht mehr klar ist, warum er überhaupt funktioniert.

Der nächste bitte

Der zweite Algorithmus aus unserem ersten Artikel ist zwar nicht mehr ganz so verschwenderisch wie stures Ausmultiplizieren, aber er führt im Algorithmenzustand immer noch eine reelle (rationale) Zahl y mit, welche in jedem Schritt gerundet werden muss. Das ist immer noch unnötig aufwendig. Wir haben uns am Anfang bewusst eingeschränkt auf „flache“ (dx < dy), „aufsteigende“ (dy >= 0) Linien. Bei der Rasterung solcher Linien geht’s für jeden x-Schritt in y-Richtung entweder einen Schritt nach oben, oder wir bleiben auf der gleichen Höhe. Das klingt zum einen plausibel, ist andererseits aber auch beweisbar:¹

(define-fun candidate-next-pixel-top ((st State-1))
  Int

  (+ 1 (state-1-y st)))

(define-fun candidate-next-pixel-bottom ((st State-1))
  Int

  (state-1-y st))

(push)
(declare-const st State-1)

(assert (state-1-valid st))

(assert (not (let ((next-st (step-1 st)))
               (or
                (= (candidate-next-pixel-bottom st)
                   (state-1-y next-st))
                (= (candidate-next-pixel-top st)
                   (state-1-y next-st))))))

(check-sat)
(pop)

Wir wollen im folgenden mehrere Zwischenergebnisse auf unsat prüfen. Damit wir nicht jedes mal Code des vorigen Zwischenschritts wieder auskommentieren müssen, nutzen wir die Ausdrücke (push) und (pop), um einen frischen Kontext einzuführen. Alles, was seit dem letzten (push) deklariert oder assertet wurde, ist nach dem (pop) wieder weg.

Die Entscheidung, welchen dieser beiden Kandidatenpixel wir auswählen, treffen wir anhand der Abstände der jeweiligen Mittelpunkte zur Ideallinie. Dieser Zusammenhang lässt sich am besten anhand einer Grafik illustrieren:

Die Abstände zum oberen Pixel d-top und zum unteren Pixel d-bottom werden wir später als expliziten Teil des Zustands mitführen. Jetzt – in der Phase, wo wir uns noch klar werden müssen über alle Zusammenhänge – können wir die beiden Abstände mit Funktionen ableiten.

(define-fun state-1-d-top ((st State-1))
  Real
  (-
   (candidate-next-pixel-top st)
   (state-1-exact-y (step-1 st))))

(define-fun state-1-d-bottom ((st State-1))
  Real
  (-
   (state-1-exact-y (step-1 st))
   (candidate-next-pixel-bottom st)))

Hierbei ist wichtig, dass sich die beiden Abstände d-bottom und d-top auf den nächsten Zustand beziehen. Deshalb müssen wir zur Berechnung einmal mit step-1 weiterschalten.

Mit diesen Definitionen können wir prüfen, ob das Wissen über die beiden Abstände ausreicht, um die Entscheidung zu treffen, ob wir den oberen Kandidatenpixel auswählen müssen oder den unteren. Eine Entscheidung in SMT-LIB-2 treffen wir mit dem ite-Ausdruck, was für „if-then-else“ steht:

(push)
(declare-const st State-1)
(assert (state-1-valid st))

(assert (not
         (let ((next-y (state-1-y
                        (step-1 st)))
               (dx (line-dx (state-1-line st)))
               (dy (line-dy (state-1-line st)))
               (x (state-1-x st))
               (y (state-1-y st)))

           (ite (< (state-1-d-top st)
                   (state-1-d-bottom st))

                (= next-y
                   (candidate-next-pixel-top st))
                (= next-y
                   (candidate-next-pixel-bottom st))))))

(check-sat)
(pop)

Z3 sagt hierzu sat. Das ist schlecht; unsere eigentliche Aussage gilt nicht. Um herauszufinden, was wir falsch gemacht haben, können wir einen get-value-Aufruf hinter das (check-sat) schreiben. get-value nimmt als Argument eine Liste mit Ausdrücken und druckt deren Werte aus für das gefundene Modell (sprich: die gefundene Variablebelegung). Wir lassen uns ein paar interessante Werte ausdrucken.

(get-value
 (st

  (state-1-y st)

  (state-1-y
   (step-1 st))

  (state-1-d-top st)

  (state-1-d-bottom st)

  (candidate-next-pixel-top st)

  (candidate-next-pixel-bottom st)
  ))

Das Ergebnis davon ist:

sat
((st (mk-state-1 (mk-line 4 3) 1))
 ((state-1-y st) 1)
 ((state-1-y
   (step-1 st)) 2)
 ((state-1-d-top st) (/ 1.0 2.0))
 ((state-1-d-bottom st) (/ 1.0 2.0))
 ((candidate-next-pixel-top st) 2)
 ((candidate-next-pixel-bottom st) 1))

Wir sehen, dass der Abstand zum oberen Kandidaten mit 1/2 genau gleich groß ist wie der Abstand zum unteren Kandidaten. In der Bedingung unseres ite steht ein <, also strikt kleiner, was in diesem Fall nicht gegeben ist. Dann wählen wir den unteren Pixel aus. Richtig gewesen wäre der obere, denn bei 1/2 rundet round-to-nearest auf statt ab. Ein naheliegender Lösungsversuch wäre, das < in ein <= zu ändern. Und siehe da, damit sagt Z3 unsat.

Die Bedingung für das ite ist jetzt (<= (state-1-d-top st) (state-1-d-bottom st)), was immer noch die Berechnung des exakten y-Werts erfordert. Davon wollen wir jetzt weg. Insbesondere wollen wir die reellen Zahlen (Real) in der Implementierung los werden. Die reellen Zahlen sind hervorragend geeignet für die Formulierung der Spezifikation und für die Beschreibung des Zusammenhangs von Implementierung und Spezifikation. In der Implementierung sind sie aber zu rechenaufwendig und auch gar nicht nötig. Zum Zwecke ihrer Loswerdung formulieren wir die ite-Bedingung nach und nach um. Diese Umformungen können wir dabei jeweils von Z3 auf Korrektheit prüfen lassen.

Wir substituieren zunächst von Hand alle Funktionsaufrufe durch deren jeweilige Implementierungen und formen nebenbei die Ungleichung so um, dass links nur noch die Null steht.

(push)

(declare-const st State-1)
(assert (state-1-valid st))

(assert (not
         (=

          ;; Die ursprüngliche Bedingung ...
          (<= (state-1-d-top st)
              (state-1-d-bottom st))

          ;; ... ist äquivalent zu dieser hier
          (<= 0

              (-
               (* 2 (line-y (state-1-line (step-1 st))
                            (state-1-x (step-1 st))))
               (* 2 (state-1-y st))
               1)))))

(check-sat)
(pop)

Ergebnis: unsat.

Wir wissen außerdem, dass die Linie state-1-line sich nicht ändert. Darüberhinaus ist das x nach einem step-1 einfach x + 1. Wenn wir zusätzlich noch die Definition von line-y und dann die von line-slope einsetzen, erhalten wir diese Bedingung:

(push)

(declare-const st State-1)
(assert (state-1-valid st))

(assert (not
         (=

          ;; Die ursprüngliche Bedingung ...
          (<= (state-1-d-top st)
              (state-1-d-bottom st))

          ;; ... ist äquivalent zu dieser hier
          (let ((next-y (state-1-y (step-1 st)))
                (dx (line-dx (state-1-line st)))
                (dy (line-dy (state-1-line st)))
                (x (state-1-x st))
                (y (state-1-y st)))

           (<= 0
               (-
                (* 2 (* (/ dy dx)
                        (+ 1 x)))
                (* 2 y)
                1))))))
(check-sat)
(pop)

Die Ungleichung enthält (/ dy dx) als einzigen Bruch. Den können wir beseitigen, indem wir die Ungleichung mit dx durchmultiplizieren. Das ist zulässig, so lange dx ungleich 0 ist. Anders als auf Papier müssen wir die Zulässigkeit hier gar nicht verargumentieren, sondern wir können die Umformulierung einfach machen und Z3 sagt uns dann, ob das zulässig war. Die Ungleichung sieht jetzt so aus:

(<= 0
    (-
     (* 2 dy (+ 1 x))
     (* 2 y dx)
     dx))

Z3 sagt dazu immer noch unsat. Die Einschränkung (state-1-valid st) enthält die Einschränkung, dass die im Zustand befindliche Linie valide sein soll. In line-valid steht (> (line-dx l) 0). Wir können mal probeweise diese Einschränkung löschen und Z3 fragen, ob dann unsere Aussagen immer noch gelten. Erwartungsgemäß sagt Z3 jetzt sat.

Dem komplizierten Ausdruck (- (* 2 (* dy (+ 1 x))) (* 2 y dx) dx) geben wir den Namen d. Dieses d ist jetzt eine Ganzzahl und zu deren Berechnung müssen wir auch nur Ganzzahlarithmetik betreiben. Nur die Multiplikationen stören uns noch. Aber um eine Multiplikation aufzulösen kennen wir ja schon einen Trick: Wir merken uns den jeweiligen Wert aus dem vorherigen Schritt und addieren für den nächsten Schritt was drauf (oder ziehen was ab). Das sollte hier funktionieren, denn das Produkt (* 2 dy (+ 1 x)) erhöht sich in jedem Schritt um genau 2 dy und das Produkt (* 2 y dx) erhöht sich je nachdem, ob wir den unteren Kandidatenpixel oder den oberen auswählen entweder gar nicht oder um ein konstantes (* 2 dx).

Mit diesen formal verifizierten Überlegungen sind wir gewappnet für den finalen Schlag: Wir bauen einen Algorithmus, der eine Linie zeichnet und sich dabei nur einfacher Ganzzahlarithmetik bedient.

I keep my ends out for the tie that binds

Wir verfolgen wieder das übliche Schema aus dem ersten Artikel. Zuerst definieren wir ein Zustandsobjekt. Das enthält neben dem üblichen Kram jetzt unser d, welches unsere Entscheidung für den oberen oder unteren nächsten Pixel treibt:

(declare-datatypes ()
                   ((State-4 (mk-state-4 (state-4-dx Int)
                                         (state-4-dy Int)
                                         (state-4-x Int)
                                         (state-4-y Int)
                                         (state-4-d Int)
                                         ))))

Ein solcher Zustand ist valide, wenn die üblichen Einschränkungen gelten und zusätzlich unser d dem Wert entspricht, der uns vorschwebt:

(define-fun state-4-m ((st State-4))
  Real

  (/ (state-4-dy st)
     (state-4-dx st)))

(define-fun state-4-valid ((st State-4))
  Bool
  (let ((dx (state-4-dx st))
        (dy (state-4-dy st))
        (m (state-4-m st))
        (x (state-4-x st))
        (y (state-4-y st))
        (d (state-4-d st)))
    (and
     (< 0 dx)
     (<= 0 dy dx)

     (>= x 0)

     ;; y is actually the correct y for m and x
     (= y (round-to-nearest (* m x)))
     
     ;; d is what we want it to describe
     (= d
        (- (* 2 dy (+ 1 x))
           (* 2 y dx)
           dx))
     )))

Fehlt noch die Funktion für den Einstieg:

(define-fun into-state-4 ((st State-1))
  State-4

  (let ((dx (line-dx (state-1-line st)))
        (dy (line-dy (state-1-line st)))
        (x (state-1-x st)))

    (let ((y (round-to-nearest (* (/ dy dx) x))))
      (mk-state-4
       dx
       dy
       x
       y
       (- (* 2 dy (+ 1 x))
          (* 2 y dx)
          dx)))))

Hier findet sich wieder die Bruchrechnung (/ dy dx) und das darauffolgende runden mit round-to-nearest. Im Gegensatz zum langsamen Referenzalgorithmus müssen wir diese aufwendigen Rechenoperationen aber jetzt nur einmal zum Einstieg machen. Wer sowieso immer bei x = 0 anfangen will zu zeichnen, muss noch nicht mal dividieren und runden, sondern kann direkt mit y = 0 loslegen. Auch diese Aussage könnten wir von Z3 verifizieren lassen. We leave this as an exercise to the reader.

Jetzt fehlt noch die Schrittfunktion:

(define-fun step-4 ((st State-4))
  State-4

  (let ((dx (state-4-dx st))
        (dy (state-4-dy st))
        (x (state-4-x st))
        (y (state-4-y st))
        (d (state-4-d st)))

    (ite (<= 0 d)
         ;; Diagonalschritt
         (mk-state-4
          dx
          dy
          (+ x 1)
          (+ y 1)
          (+ d (- (* 2 dy)
                  (* 2 dx))))

         ;; Rechtsschritt
         (mk-state-4
          dx
          dy
          (+ x 1)
          y
          (+ d (* 2 dy))))))

Als Korrektheitsaussage würden wir gern wieder dem oben etablierten Schema entsprechen:

(declare-const st State-1)
(assert (state-1-valid st))

(assert
 (not
  (= (step-4
      (into-state-4 st))
     (into-state-4
      (step-1 st)))))

Das terminiert leider nicht. Das ist leider gar nicht schön und zeigt die Grenzen eines SMT-Solvers als Verifikationsinstrument. Wir können mit etwas Gefrickel aber eine Korrektheitsaussage hinbasteln, die äquivalent ist:

(declare-const before-1 State-1)
(assert (state-1-valid before-1))

(assert (not
         (let ((before-4 (into-state-4 before-1)))
           (let ((after-4 (step-4 before-4))
                 (after-1 (step-1 before-1)))
             (and
              (= (state-4-x after-4)
                 (state-1-x after-1))

              (state-4-valid after-4)

              (= (state-1-y after-1)
                 (state-4-y after-4)))))))

(check-sat)

Die Bedingungen zu x und y bewirken, dass step-4 dasselbe Ergebnis liefert wie step-1 und die Bedingung (step-4-valid after-4) bewirkt, dass wir auch nicht den Pfad der Tugend verlassen.

Because you‘re mine, I walk the line

Unsere finale Formulierung ist eine Instanz des sog. Bresenham-Algorithmus. Wir haben diesen hier nicht nur korrekt implementiert, sondern diese Korrektheit auch sichergestellt und außerdem auf dem Weg dorthin einen Einblick erhalten, warum der Algorithmus funktioniert.

Um die Korrektheit zu zeigen, formulierten wir Bedingungen mit Rückgriff auf die reellen Zahlen – nicht Float, nicht Double, sondern präzise reelle Zahlen. Das ist richtig, denn nur mit reellen Zahlen lässt sich die Rasterisierung einer Geraden akkurat beschreiben. Das bedeutet jedoch nicht, dass eine Implementierung ebenfalls mit reellen Zahlen arbeiten muss. Das wäre für einen Computer sehr aufwendig und mit dem Bresenham-Algorithmus haben wir gezeigt, dass es ja auch ohne geht.

Im nächsten Artikel zeigen wir eine Implementierung des Bresenham-Algorithmus in Haskell. Wir nutzen Liquid Haskell, um sicherzustellen, dass auch diese Haskell-Implementierung korrekt im Sinne unserer obigen Überlegungen ist.