z3py

Z3はMicrosoft Researchの開発した，SAT・SMTソルバである． ブール代数の充足問題に限らず，整数や実数の制約充足問題，さらには整数計画問題などにも適用が可能である．

z3pyはZ3のPythonバインディングであり，Z3に付属している．

ドキュメントはPythonモジュールのhelpドキュメント，あるいはそれをHTMLにまとめたものが良い．rise4funにチュートリアルもあるらしいのだが接続できたことがない．

インストール

Z3のGithubページからgit cloneしてビルド・インストールする． z3pyは /usr/lib/python2.7/dist-packages 以下に同時にインストールされる． Python3用はないの？

基本的な流れ

   1 # モジュールをimport
   2 from z3 import *
   3 
   4 # 変数を作成．引数は人間が見てわかりやすい変数名．
   5 x, y = Bools(["x", "y"])
   6 
   7 # ソルバのインスタンスを生成して
   8 s = Solver()
   9 
  10 # 制約を追加
  11 s.add(x == True)
  12 s.add(x != y)
  13 
  14 # 解を探索，モデルを取得
  15 s.check()
  16 m = s.model()
  17 
  18 # 解を取得
  19 solution_x = is_true(m[x])
  20 solution_y = is_true(m[y])
  21 
  22 print solution_x, solution_y

変数の作成

ブール変数を作る

一つ作成

   1 x = Bool("x")

複数同時作成

変数名をスペースで区切って渡すか

   1 x, y = Bools("x y")

あるいは次のようにする．

   1 x, y = Bools(["x", "y"])

配列で同時作成

変数を128個含むリストを生成するには次のようにする．

   1 xs = BoolVector("x", 128)

変数名には，与えた変数名に続けて__0, __1, ...と番号が追加される．

あるいはPythonのリストの内包記法を用いて次のようにすればよい．

   1 xs = [Bool("x%d" % i) for i in range(128)]

整数変数を作る

一つ作成

   1 x = Int("x")

複数同時作成

変数名をスペースで区切って渡すか

   1 x, y = Ints("x y")

あるいは次のようにする．

   1 x, y = Ints(["x", "y"])

配列で同時作成

変数を128個含むリストを生成するには次のようにする．

   1 xs = IntVector("x", 128)

変数名には，与えた変数名に続けて__0, __1, ...と番号が追加される．

あるいはPythonのリストの内包記法を用いて次のようにすればよい．

   1 xs = [Int("x%d" % i) for i in range(128)]

実数変数を作る

一つ作成

   1 x = Real("x")

複数同時作成

変数名をスペースで区切って渡すか

   1 x, y = Reals("x y")

あるいは次のようにする．

   1 x, y = Reals(["x", "y"])

配列で同時作成

変数を128個含むリストを生成するには次のようにする．

   1 xs = RealVector("x", 128)

変数名には，与えた変数名に続けて__0, __1, ...と番号が追加される．

あるいはPythonのリストの内包記法を用いて次のようにすればよい．

   1 xs = [Real("x%d" % i) for i in range(128)]

制約の設定

制約の追加

次節以降で紹介する論理式やブール変数により表される命題$P_1, P_2, \dots$があるとき， それらが真になるという制約を加えるには次のようにする．

   1 s.add(P1, P2, ...)

追加した制約を全て満たすかどうか判定される．つまり，全てANDで接続した論理式が評価される．

他にappend, assert_exprs, insertというメソッドもあるが単なる別名である．

論理AND

$P_1 \wedge P_2 \wedge \dots$ を表すには次のようにする．

   1 And(P1, P2, ...)

あるいは次のようにする．

   1 And([P1, P2, ...])

注意: Pythonのandでは求める動作をしない．

論理OR

$P_1 \vee P_2 \vee \dots$ を表すには次のようにする．

   1 Or(P1, P2, ...)

あるいは次のようにする．

   1 Or([P1, P2, ...])

注意: Pythonのorでは求める動作をしない．

論理XOR

$P \oplus Q$ を表すには次のようにする．

   1 Xor(P, Q)

引数はちょうど2個のみで配列も受け付けないので注意する．

注意: Pythonの^演算子では求める動作をしない． ^演算子はビットXORの項を参照せよ．

論理NOT

$\overline{P}$ を表すには次のようにする．

   1 Not(P)

注意: Pythonのnotでは求める動作をしない．

Implies (=>, ならば)

$P \Rightarrow Q$ を表すには次のようにする．

   1 Implies(P, Q)

同値(Equals)

$P \Leftrightarrow Q$ を表すには次のようにする．

   1 P == Q

または数値に対し$x = y$も同様である．

   1 x == y

不等式

普通のPythonの不等式と同じ． ただし不等式のチェーンはできない．

駄目な例(不等式のチェーン)

>>> 0 < x < 1
x < 1

算術演算

普通のPythonの算術演算と同じ．

• 加算 +

• 減算 -

• 乗算 *

• 除算 /

• 冪乗 **

• 剰余 %

なお整数型に対して除算を行うと，少数以下は切り下げられる．

>>> x = Int("x")
>>> s = Solver()
>>> s.add(x / 10 == -1)
>>> s.add(x % 10 == 1)
>>> s.check()
sat
>>> s.model()
[x = -9]

全て異なる値を持つ(Distinct)

式 $x_1, x_2, \dots$ が異なる値を持つことを表すには次のようにする．

   1 Distinct(x1, x2, ...)

あるいは次のようにする．

   1 Distinct([x1, x2, ...])

条件分岐(If)

if $P$ then $x$ else $y$ を表すには次のようにする．

   1 If(P, x, y)

型の変換(Int <-> Real)

Int型とReal型の間はToIntとToRealで変換できる．

ToIntでは小数点以下は切り下げられる．

>>> x = Real("x")
>>> s = Solver()
>>> s.add(-1 < x, x < 0)
>>> s.add(ToInt(x) == -1)
>>> s.check()
sat
>>> s.model()
[x = -1/2]

充足可能性を判定する方法

制約を追加したソルバオブジェクトに対し，checkメソッドを実行することで，与えた制約に対する充足可能性が判定できる．

checkメソッドはsat, unsat, unknownの三種類の値を返す．これらは定数値として定義されている．

>>> x = Int("x")
>>> s = Solver()
>>> s.add(x < 0, x > 1)
>>> r = s.check(); r
unsat
>>> r == unsat
True

充足解を得る方法

充足解を得るためには，まずcheckメソッドで充足可能性を判定し，その後にmodelメソッドでモデルを取得する． 得られたモデルを文字列化して表示すると充足解が表示されるため，手作業で行うときはこれでも良いが， プログラム上で解を利用するためには，次のようにしてモデルから対応する変数を取り出し， z3pyの型から通常のPythonの型に変換する必要がある．

>>> x, y = Ints("x y")
>>> s = Solver()
>>> s.add(0 <= x, x <= 1)
>>> s.add(x < y)
>>> s.check()
sat
>>> m = s.model(); m
[x = 0, y = 1]

まず，モデルmから変数xを表すオブジェクトに変換するためには辞書のように m[x] とする． ここでxは変数名ではなく変数のオブジェクトであることに注意する． m[x]を表示すると一見，整数や真偽値が得られているように見えるが， 実際にはz3pyのオブジェクトなので変換が必要である．

>>> m[x]
0
>>> type(m[x])
<type 'instance'>
>>> m[x] == 0
0 == 0

以下では，m[x]から各種の型に変換する方法を述べる．

ブール変数の解を得る

xがブール変数であるとする．このとき次のようにしてm[x]の真偽値を得ることが出来る．

   1 is_true(m[x])

例

>>> x = Bool("x")
>>> s = Solver()
>>> s.add(Not(x))
>>> s.check()
sat
>>> m = s.model()
>>> is_true(m[x])
False

整数変数の解を得る

xが整数変数であるとする．このとき次のようにしてm[x]の整数値を得ることが出来る．

   1 m[x].as_long()

例

>>> x = Int("x")
>>> s = Solver()
>>> s.add(x * x - x - 2 == 0)
>>> s.check()
sat
>>> m = s.model()
>>> m[x].as_long()
2

実数変数の解を正確な有理数として得る

xが実数変数であり，充足解が有理数として得られているとする． このとき，次のようにしてm[x]の分数表現を得ることが出来る．

   1 m[x].as_fraction()

分母と分子をそれぞれ整数値として得たい場合には次のようにする．

   1 # 分母
   2 m[x].denominator_as_long()
   3 # 分子
   4 m[x].numerator_as_long()

m[x]が有理数として得られているかどうかは次のようにして判定できる．

   1 is_rational_value(m[x])

例

>>> x = Real("x")
>>> s = Solver()
>>> s.add(x * 2 == 1)
>>> s.check()
sat
>>> m = s.model()
>>> is_rational_value(m[x])
True
>>> m[x].as_fraction()
Fraction(1, 2)
>>> m[x].denominator_as_long()
2
>>> m[x].numerator_as_long()
1

実数変数の解を近似された有理数として得る

xが実数変数であり，充足解が無理数として得られており， 誤差$10^{-p}$以下で近似された有理数を得たいとする． このとき，次のようにしてm[x]の近似された分数表現を得ることが出来る．

   1 m[x].approx(p).as_fraction()

逆に解が有理数として得られている時に，approxメソッドを呼び出すとエラーになるため注意が必要である． つまり，実際には次のようにするのが良い．

   1 (m[x] if is_rational_value(m[x]) else m[x].approx(p)).as_fraction()

例

>>> x = Real("x")
>>> s = Solver()
>>> s.add(x ** 2 == 2)
>>> s.check()
sat
>>> m = s.model()
>>> is_rational_value(m[x])
False
>>> m[x].approx(10)
-3109888511975/2199023255552
>>> m[x].approx(10).as_fraction()
Fraction(-3109888511975, 2199023255552)

実数変数の解を浮動小数点数として得る

十分な精度でFraction型に変換した後にfloat関数で浮動小数点型に変換すれば良い．

例

魔法陣

$n \times n$の魔方陣を生成する．

   1 from z3 import *
   2 import sys
   3 
   4 n = 4
   5 
   6 x = [[Int("x[%d,%d]" % (i,j)) for j in range(n)] for i in range(n)]
   7 s = Solver()
   8 
   9 for i in range(n):
  10     for j in range(n):
  11         s.add(1 <= x[i][j], x[i][j] <= n*n)
  12 
  13 # all numbers are unique
  14 s.add(Distinct(sum(x, [])))
  15 
  16 # all rows have same sum
  17 for i in range(1, n):
  18     s.add(sum(x[0]) == sum(x[i]))
  19 
  20 # all columns have same sum
  21 for j in range(1, n):
  22     s.add(sum(map(lambda row: row[0], x)) == sum(map(lambda row: row[j], x)))
  23 
  24 s.check()
  25 m = s.model()
  26 
  27 for i in range(n):
  28     for j in range(n):
  29         sys.stdout.write(" %2d" % m[ x[i][j] ].as_long())
  30     sys.stdout.write("\n")

この方法では$n=5$を解くのにi7-4771プロセッサで15秒を要した．6以上は現実的な時間では無理のようだ．

なお，各行・各列の和を定数として与えてやれば早くなるかと思いきや，むしろ遅くなった． 余計な判定処理が加わったためだろうか．以下の制約を追加したところ$n=5$で49秒を要した．

   1 # sum is n(n*n+1)/2
   2 s.add(sum(x[0]) == n*(n*n+1)/2)
   3 s.add(sum(map(lambda row: row[0], x)) == n*(n*n+1)/2)

虫食い算

Googleのアレ

  WWWDOT
- GOOGLE
--------
  DOTCOM

値$x = x_n \dots x_1,\, y = y_n \dots y_1,\, z = z_n \dots z_1$に関する減算$x - y = z$は，各$i$桁目について，上の桁からのボロー$b_{i+1}$を用いて$x_i + b_{i+1} - y_i - b_i = z_i$と表すことができる． Z3は整数を扱うことができるため変数を整数とし，各桁について値の範囲が0から9となるよう制約を与え，この関係式を制約として与える． さらに最上位の桁は0にならないように制約を与え，同じアルファベットは同じ数字，異なるアルファベットは別の数字となるように制約を与える．

   1 # -*- coding: utf-8 -*-
   2 # 大学院のレポートに出したやつ:-p
   3 from z3 import *
   4 
   5 W, D, O, T, G, L, E, C, M = Ints("W D O T G L E C M")
   6 
   7 # 筆算の変数
   8 #   x_5 ... x_0
   9 # - y_5 ... y_0
  10 # -------------
  11 #   z_5 ... z_0
  12 x = IntVector("x", 6)
  13 y = IntVector("y", 6)
  14 z = IntVector("z", 6)
  15 b = IntVector("b", 7) # borrow
  16 
  17 s = Solver()
  18 
  19 # アルファベットは全て別の数字
  20 s.add(Distinct(W, D, O, T, G, L, E, C, M))
  21 
  22 # 筆算のロジック
  23 s.add(b[0] == 0)
  24 for i in range(6):
  25     s.add(And(0 <= x[i], x[i] < 10))
  26     s.add(And(0 <= y[i], y[i] < 10))
  27     s.add(And(0 <= z[i], z[i] < 10))
  28     s.add(Or(b[i] == 0, b[i] == 1))
  29     s.add(x[i] + b[i+1]*10 - y[i] - b[i] == z[i])
  30 # 最上位桁は0にならない
  31 s.add(x[5] != 0)
  32 s.add(y[5] != 0)
  33 s.add(z[5] != 0)
  34 s.add(b[6] == 0)
  35 
  36 # アルファベットを筆算に割り当て
  37 for (u, v) in zip([W,W,W,D,O,T], x[::-1]):
  38     s.add(u == v)
  39 for (u, v) in zip([G,O,O,G,L,E], y[::-1]):
  40     s.add(u == v)
  41 for (u, v) in zip([D,O,T,C,O,M], z[::-1]):
  42     s.add(u == v)
  43 
  44 # 判定・出力
  45 s.check()
  46 m = s.model()
  47 for v in [W, D, O, T, G, L, E, C, M]:
  48     print "%s = %d" % (v, m[v].as_long())

なお，記事を書くにあたって，制約をx - y == zの形で記述して，除算と剰余で各桁の値を取り出してアルファベットに割り当てる手法でも実装を試みたが，こちらは短時間では停止しなかった．

   1 # 引き算の変数
   2 xx = Int("x")
   3 yy = Int("y")
   4 zz = Int("z")
   5 s.add(100000 <= xx, xx <= 999999)
   6 s.add(100000 <= yy, yy <= 999999)
   7 s.add(100000 <= zz, zz <= 999999)
   8 s.add(xx - yy == zz)
   9 
  10 # 各桁の値
  11 for i in range(6):
  12     s.add(x[i] == xx / pow(10, i) % 10)
  13     s.add(y[i] == yy / pow(10, i) % 10)
  14     s.add(z[i] == zz / pow(10, i) % 10)

整数計画問題