Z3-Python 中的泰勒展开式三角函数

Question

我需要在 Z3 中设计一个余弦（和正弦）函数，但这通常很难且无法确定（例如，请参阅如何在 Z3 Python 中使用内置三角函数？）。

但是，我对近似精度方法很满意。因此，我设计了以下使用泰勒展开式的函数（对于

cos(a)

，其中

a

以弧度表示）：

from z3 import *

def calculate_cosine(angle):
    solver = Solver()
    result = Real('result')

    #Note the Taylor expansion for cos (e.g., next term is "+ angle**8/40320")
    #Also, note 2!=2, 4!=24 etc.
    solver.add(result == 1 - angle**2/2 + angle**4/24 - angle**6/720)

    if solver.check() == sat:
        model = solver.model()
        cosine = model[result]
        return cosine

    return None

angle = math.pi / 8  # Angle in radians
cosine = calculate_cosine(angle)
print(cosine)

但是，这只是Python代码，而我想执行像

Exists y. 0<=approx_cos(y)<=1

这样的查询，我认为我必须构建一个Z3函数

approx_cos

（具有典型语法

If...

），它将在这个中调用案例，一些 NRA 求解器。

但是我不熟悉如何构建这样的功能，你能帮助我吗？

PS：Z3有内置阶乘运算吗？

Answer 1

您可以使用简单的迭代循环来编写泰勒级数展开式，并随时添加各项。下面是一个实现，它将您想要包含的术语数量作为参数：

from z3 import *

# Approximate cosine using taylor expansion
# n is the number of terms we want to include
# Essentially, we implement the formula
#
#     sum_{k=0}^{n-1} (-1)^k * x^2k / (2k)!
#
# in an iterative way, to reduce calculations.
def taylor_cosine(x, n):
    if n < 2:
        raise Exception("taylor_cosine: n must be at least 2, got: %d" % n)

    s    = 1
    sign = 1
    term = 1

    # In each iteration:
    #   sign flips
    #   term gets multiplied by x^2 / (2k * (2k-1))
    for k in range(1, n):
        sign *= -1
        term *= x*x
        term /= 2*k*(2*k - 1)
        s    += sign * term

    return s

注意这个函数与z3无关，你可以直接使用它：

>> taylor_cosine(math.pi, 5)
-0.9760222126236076

包含的术语越多，结果就越准确。

但是这个函数的好处是我们小心翼翼地以某种方式编写它，以便它也可以象征性地使用。也就是说，参数

x

可以是z3符号变量。 （请注意，

n

必须是一个常量：没有简单/可判定的方法来支持符号

n

。但这完全是另一个讨论。）

作为示例，这里有一个小测试程序，尝试通过询问

pi

来近似

z3

的值，逐次逼近，其中调用余弦的结果是

-1

:

def test(numberOfTerms):
  s  = Solver()
  x  = Real('x')
  s.add(And(0 <= x, x <= 2*math.pi))
  cx = taylor_cosine(x, numberOfTerms)
  s.add(And(cx == -1))

  print("Testing with %d terms:" % numberOfTerms)

  r = s.check()
  if r == sat:
      mx = s.model()[x]
      if(is_algebraic_value(mx)):
          mx = mx.approx()
      vx = float(mx.as_fraction())
      print("  x                    = %s"  % vx)
      print("  taylor_cosine(x, %2d) = %s" % (numberOfTerms, taylor_cosine(vx, numberOfTerms)))
      print("  cosine(x)            = %s"  % math.cos(vx))
  else:
      print("  Solver said: %s" % r)

我们只需创建一个

Real

，并要求求解器确保其余弦为

-1

。我们小心地处理实数和代数实数，因为 z3 可以为涉及实数的查询生成这两种结果。 （前者当解是有理数时，后者是代数解，即作为多项式的根。）

让我们看看实际效果：

for i in range(2, 10):
    test(i)
    print("====")

打印：

Testing with 2 terms:
  x                    = 2.0
  taylor_cosine(x,  2) = -1.0
  cosine(x)            = -0.4161468365471424
====
Testing with 3 terms:
  Solver said: unsat
====
Testing with 4 terms:
  x                    = 2.751711543190595
  taylor_cosine(x,  4) = -1.000000000000076
  cosine(x)            = -0.9249542537694367
====
Testing with 5 terms:
  Solver said: unsat
====
Testing with 6 terms:
  x                    = 3.087083093614183
  taylor_cosine(x,  6) = -1.0000000000000144
  cosine(x)            = -0.9985147217565723
====
Testing with 7 terms:
  Solver said: unsat
====
Testing with 8 terms:
  x                    = 3.138726428355767
  taylor_cosine(x,  8) = -0.9999999999999999
  cosine(x)            = -0.999995892379266
====
Testing with 9 terms:
  Solver said: unsat
====

观察：

1. 求解器可以说

   unsat

   （情况 3、5、7 和 9）：这是因为对于给定运行中我们包含的项数，它可能能够推断出调用所产生的表达式

   Taylor_cosine

   永远不会完全是

   -1

   。如果你不希望这样，你应该小心并断言你的约束，以留出一些回旋的空间。 （即，在用户选择的 epsilon 内。）

2. 随着项数的增加，结果通常会越来越准确。我们可以以一种惊人的方式观察到这一点：有两项，我们得到的值相差甚远：

   2

   。有了 8 个项，我们得到

   3.139

   ，这还不算太糟糕。我们可以看到，随着项数的增加，我们对

   pi

   的近似变得越来越好：

   2

   、

   2.75

   、

   3.08

   、

   3.139

   等。

3. 不用说，包含的术语越多，z3 的问题就会变得越复杂。如果您只坚持实值约束，那么您总是会从 z3 得到答案，因为它确实有一个非线性实数算术的决策过程。 （当然，这并不意味着它会“很快”。只是它能够在给定足够的时间/内存等的情况下决定它。）如果您在那里混合整数约束，那么所有的赌注都会被取消因为你将处于半可判定片段中。我的猜测是，您会比您想要的更频繁地获得

   unknown

   ，但这一切都取决于周围的其他限制以及您的实际应用程序是什么。

无论如何..希望这能让你起步。要记住的关键一点是，您需要更多的术语来提高准确性，但代价是计算复杂性。由于您正在处理近似值，请确保您的约束也允许近似结果；即，不要要求“恰好是这个值”，而是要求“在这个值周围加上您选择的 epsilon”。否则，您将得到

unsat

，如上所示。