Тригонометрические функции разложения Тейлора в Z3-Python

Вы можете запрограммировать разложение ряда Тейлора с помощью простого итерационного цикла, складывая члены по ходу дела. Ниже приведена реализация, которая принимает количество терминов, которые вы хотите включить в качестве параметра:

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. Если вы этого не хотите, вам следует быть осторожным и отстаивать свои ограничения, чтобы оставить некоторую свободу действий. (т. е. в пределах выбранного пользователем эпсилона.)

2. Результат, как правило, будет становиться все более и более точным по мере увеличения количества терминов. Мы можем наблюдать это поразительным образом: с двумя членами значение, которое мы получаем, очень далеко: 2. При 8 слагаемых получаем 3.139, что не так уж и страшно. И мы видим, что наше приближение pi становится все лучше и лучше с увеличением количества терминов: 2, 2.75, 3.08, 3.139 и т. д.

3. Излишне говорить, что чем больше терминов вы включите, тем сложнее станет проблема для z3. Если вы придерживаетесь только вещественнозначных ограничений, вы всегда получите ответ от z3, поскольку у него есть процедура принятия решения для нелинейной вещественной арифметики. (Конечно, это не означает, что он будет «быстрым». Просто он сможет принимать решения, имея достаточно времени/памяти и т. д.) Если вы смешаете туда целочисленные ограничения, тогда все ставки отключены, поскольку вы окажетесь в полуразрешимом фрагменте. Я предполагаю, что вы будете получать unknown чаще, чем хотелось бы, но все зависит от других ограничений и вашего фактического приложения.

В любом случае... Надеюсь, это поможет вам сдвинуться с мертвой точки. Важно помнить, что вам нужно больше условий для точности за счет вычислительной сложности. А поскольку вы имеете дело с приближениями, убедитесь, что ваши ограничения также позволяют получить приблизительные результаты; то есть вместо того, чтобы запрашивать «именно это значение», попросите «около этого значения с эпсилоном по вашему выбору». В противном случае вы получите unsat, как показано выше.