#include "nlsat/nlsat_transcendentals.h"
#include "nlsat/nlsat_solver.h"
#include "math/polynomial/polynomial.h"
#include "math/polynomial/algebraic_numbers.h"
#include <cmath>
#include <limits>
namespace nlsat {
double transcendentals::eval(transcendental_op_kind op, double x) {
return transcendental_eval::eval(op, x);
}
double transcendentals::error_bound(transcendental_op_kind op, double x, double fx) {
return transcendental_eval::error_bound(op, x, fx);
}
rational transcendentals::to_rational(double d) {
return transcendental_eval::to_rational(d);
}
bool transcendentals::exp_taylor_bracket_at(rational const& xr, rational& lo, rational& hi) {
return transcendental_eval::exp_taylor_bracket_at(xr, lo, hi);
}
bool transcendentals::log_taylor_bracket_at(rational const& xr, rational& lo, rational& hi) {
return transcendental_eval::log_taylor_bracket_at(xr, lo, hi);
}
bool transcendentals::atan_taylor_bracket_at(rational const& xr, rational& lo, rational& hi) {
return transcendental_eval::atan_taylor_bracket_at(xr, lo, hi);
}
bool transcendentals::interval_eval(transcendental_op_kind op, double lo, double hi, double& lo_val, double& hi_val) {
return transcendental_eval::wide_interval_eval(op, lo, hi, lo_val, hi_val);
}
static literal bound_literal(solver& s, var x, atom::kind k, rational const& bound) {
rational den = denominator(bound);
rational coeff = den;
rational c = -(bound * den);
polynomial_ref p(s.pm().mk_linear(1, &coeff, &x, c), s.pm());
poly* pp = p.get();
bool is_even = false;
return s.mk_ineq_literal(k, 1, &pp, &is_even);
}
static literal linear_literal(solver& s, var x1, rational const& c1, var x2, rational const& c2,
atom::kind k, rational const& bound) {
rational den = lcm(denominator(c1), lcm(denominator(c2), denominator(bound)));
rational coeffs[2] = { c1 * den, c2 * den };
var vars[2] = { x1, x2 };
rational c = -(bound * den);
polynomial_ref p(s.pm().mk_linear(2, coeffs, vars, c), s.pm());
poly* pp = p.get();
bool is_even = false;
return s.mk_ineq_literal(k, 1, &pp, &is_even);
}
static constexpr double k_pi_2_ub = 1.5707963267948968; static constexpr double k_pi_ub = 3.1415926535897936;
void transcendentals::add(transcendental_op_kind op, var arg, var val) {
m_apps.push_back({ op, arg, val });
m_retry.push_back(0);
add_global_axioms(op, arg, val);
add_range_axioms(op, val);
find_and_add_identity_axiom(op, arg, val);
}
void transcendentals::add_pi(var val) {
if (val == null_var || m_pi_var != null_var)
return;
m_pi_var = val;
rational lo("3.14159265358979");
rational hi("3.14159265358980");
literal lo_lit = ~bound_literal(s, val, atom::LT, lo); literal hi_lit = ~bound_literal(s, val, atom::GT, hi); s.mk_clause(1, &lo_lit, nullptr);
s.mk_clause(1, &hi_lit, nullptr);
}
void transcendentals::add_atan2(var y, var x, var val) {
if (y == null_var || x == null_var || val == null_var)
return;
m_atan2_apps.push_back({ y, x, val });
m_retry_atan2.push_back(0);
rational bnd = to_rational(k_pi_ub);
literal lit1 = ~bound_literal(s, val, atom::GT, bnd);
literal lit2 = ~bound_literal(s, val, atom::LT, -bnd);
s.mk_clause(1, &lit1, nullptr);
s.mk_clause(1, &lit2, nullptr);
}
void transcendentals::add_global_axioms(transcendental_op_kind op, var arg, var val) {
switch (op) {
case transcendental_op_kind::EXP: {
literal lit = ~linear_literal(s, val, rational(1), arg, rational(-1), atom::LT, rational(1));
s.mk_clause(1, &lit, nullptr);
literal lits[2] = {
~bound_literal(s, arg, atom::LT, rational(0)),
bound_literal(s, val, atom::LT, rational(1))
};
s.mk_clause(2, lits, nullptr);
break;
}
case transcendental_op_kind::LOG: {
literal lits[2] = {
~bound_literal(s, arg, atom::GT, rational(0)),
~linear_literal(s, val, rational(1), arg, rational(-1), atom::GT, rational(-1))
};
s.mk_clause(2, lits, nullptr);
break;
}
default:
break; }
}
void transcendentals::add_range_axioms(transcendental_op_kind op, var val) {
switch (op) {
case transcendental_op_kind::SIN:
case transcendental_op_kind::COS: {
auto lit1 = ~bound_literal(s, val, atom::GT, rational(1));
auto lit2 = ~bound_literal(s, val, atom::LT, rational(-1));
s.mk_clause(1, &lit1, nullptr); s.mk_clause(1, &lit2, nullptr); break;
}
case transcendental_op_kind::TANH: {
literal lo = bound_literal(s, val, atom::GT, rational(-1)); literal hi = bound_literal(s, val, atom::LT, rational(1)); s.mk_clause(1, &lo, nullptr);
s.mk_clause(1, &hi, nullptr);
break;
}
case transcendental_op_kind::COSH: {
literal lit = ~bound_literal(s, val, atom::LT, rational(1)); s.mk_clause(1, &lit, nullptr);
break;
}
case transcendental_op_kind::ACOSH: {
literal lit = ~bound_literal(s, val, atom::LT, rational(0)); s.mk_clause(1, &lit, nullptr);
break;
}
case transcendental_op_kind::ASIN: {
rational bnd = to_rational(k_pi_2_ub);
auto lit1 = ~bound_literal(s, val, atom::GT, bnd);
auto lit2 = ~bound_literal(s, val, atom::LT, -bnd);
s.mk_clause(1, &lit1, nullptr);
s.mk_clause(1, &lit2, nullptr);
break;
}
case transcendental_op_kind::ACOS: {
rational bnd = to_rational(k_pi_ub);
auto lit1 = ~bound_literal(s, val, atom::GT, bnd);
auto lit2 = ~bound_literal(s, val, atom::LT, rational(0));
s.mk_clause(1, &lit1, nullptr);
s.mk_clause(1, &lit2, nullptr);
break;
}
default:
break; }
}
void transcendentals::find_and_add_identity_axiom(transcendental_op_kind op, var arg, var val) {
for (unsigned i = 0; i + 1 < m_apps.size(); ++i) {
app const& other = m_apps[i];
if (other.arg != arg)
continue;
if (op == transcendental_op_kind::SIN && other.op == transcendental_op_kind::COS)
add_identity_axiom({ val, other.val }, true, false);
else if (op == transcendental_op_kind::COS && other.op == transcendental_op_kind::SIN)
add_identity_axiom({ other.val, val }, true, false);
else if (op == transcendental_op_kind::SINH && other.op == transcendental_op_kind::COSH)
add_identity_axiom({ other.val, val }, false, true);
else if (op == transcendental_op_kind::COSH && other.op == transcendental_op_kind::SINH)
add_identity_axiom({ val, other.val }, false, true);
else if (op == transcendental_op_kind::TANH && other.op == transcendental_op_kind::COSH)
add_identity_axiom({ other.val, val }, false, false);
else if (op == transcendental_op_kind::COSH && other.op == transcendental_op_kind::TANH)
add_identity_axiom({ val, other.val }, false, false);
}
}
void transcendentals::add_identity_axiom(identity_pair const& pr, bool is_sin_cos, bool is_cosh_sinh) {
auto& pm = s.pm();
polynomial_ref v1(pm.mk_polynomial(pr.v1), pm);
polynomial_ref v2(pm.mk_polynomial(pr.v2), pm);
polynomial_ref v1sq(pm); v1sq = v1 * v1;
polynomial_ref v2sq(pm); v2sq = v2 * v2;
polynomial_ref eq(pm);
if (is_sin_cos)
eq = v1sq + v2sq - rational(1);
else if (is_cosh_sinh)
eq = v1sq - v2sq - rational(1);
else
eq = v1sq - (v1sq * v2sq) - rational(1);
poly* pp = eq.get();
bool is_even = false;
literal lit = s.mk_ineq_literal(atom::EQ, 1, &pp, &is_even);
s.mk_clause(1, &lit, nullptr);
}
bool transcendentals::refine_app(unsigned idx) {
app const& a = m_apps[idx];
rational l, u;
s.am().get_interval(s.value(a.arg), l, u, 64);
auto enclose = [&](rational& xlo, rational& xhi, rational& rlo, rational& rhi) {
bool have_exact = false;
rational lo1, hi1, lo2, hi2;
switch (a.op) {
case transcendental_op_kind::EXP:
have_exact = exp_taylor_bracket_at(xlo, lo1, hi1) && exp_taylor_bracket_at(xhi, lo2, hi2);
break;
case transcendental_op_kind::LOG:
have_exact = log_taylor_bracket_at(xlo, lo1, hi1) && log_taylor_bracket_at(xhi, lo2, hi2);
break;
case transcendental_op_kind::ATAN:
have_exact = atan_taylor_bracket_at(xlo, lo1, hi1) && atan_taylor_bracket_at(xhi, lo2, hi2);
break;
case transcendental_op_kind::SIN:
have_exact = transcendental_eval::sin_cos_taylor_bracket(transcendental_eval::op_kind::SIN, xlo, xhi, rlo, rhi);
break;
case transcendental_op_kind::COS:
have_exact = transcendental_eval::sin_cos_taylor_bracket(transcendental_eval::op_kind::COS, xlo, xhi, rlo, rhi);
break;
default:
break; }
if (have_exact && a.op != transcendental_op_kind::SIN && a.op != transcendental_op_kind::COS) {
rlo = lo1;
rhi = hi2;
}
if (have_exact)
return true;
double lo_d = xlo.get_double(), hi_d = xhi.get_double();
if (!std::isfinite(lo_d) || !std::isfinite(hi_d))
return false;
if (to_rational(lo_d) > xlo)
lo_d = std::nextafter(lo_d, -std::numeric_limits<double>::infinity());
if (to_rational(hi_d) < xhi)
hi_d = std::nextafter(hi_d, std::numeric_limits<double>::infinity());
double lo_val, hi_val;
if (!interval_eval(a.op, lo_d, hi_d, lo_val, hi_val))
return false;
xlo = to_rational(lo_d);
xhi = to_rational(hi_d);
rlo = to_rational(lo_val);
rhi = to_rational(hi_val);
return true;
};
anum const& yv = s.value(a.val);
auto contains_value = [&](rational const& lo, rational const& hi) {
return s.am().ge(yv, lo.to_mpq()) && s.am().le(yv, hi.to_mpq());
};
rational xlo = l, xhi = u, rlo, rhi;
if (!enclose(xlo, xhi, rlo, rhi))
return false;
if (contains_value(rlo, rhi)) {
m_retry[idx] = 0;
return false;
}
rational amid = abs((l + u) / rational(2));
rational scale = std::max(amid, rational(1));
rational growth = rational::power_of_two(std::min(m_retry[idx], 30u));
rational rdelta = rational(1, 1000000) * scale * growth;
for (unsigned attempt = 0; attempt < 64; ++attempt, rdelta /= 2) {
rational lo = l - rdelta, hi = u + rdelta, flo, fhi;
if (enclose(lo, hi, flo, fhi) && !contains_value(flo, fhi)) {
xlo = lo;
xhi = hi;
rlo = flo;
rhi = fhi;
break;
}
}
literal_vector lemma;
lemma.push_back(bound_literal(s, a.arg, atom::LT, xlo));
lemma.push_back(bound_literal(s, a.arg, atom::GT, xhi));
if (s.am().lt(yv, rlo.to_mpq()))
lemma.push_back(~bound_literal(s, a.val, atom::LT, rlo)); else
lemma.push_back(~bound_literal(s, a.val, atom::GT, rhi)); s.mk_clause(lemma.size(), lemma.data(), nullptr);
++m_retry[idx];
return true;
}
bool transcendentals::refine_atan2(unsigned idx) {
atan2_app const& a = m_atan2_apps[idx];
anum const& yv = s.value(a.y);
anum const& xv = s.value(a.x);
anum const& vv = s.value(a.val);
if (s.am().is_pos(yv) && !s.am().is_pos(vv)) {
literal lits[2] = {
~bound_literal(s, a.y, atom::GT, rational(0)), bound_literal(s, a.val, atom::GT, rational(0)) };
s.mk_clause(2, lits, nullptr);
return true;
}
if (s.am().is_neg(yv) && !s.am().is_neg(vv)) {
literal lits[2] = {
~bound_literal(s, a.y, atom::LT, rational(0)), bound_literal(s, a.val, atom::LT, rational(0)) };
s.mk_clause(2, lits, nullptr);
return true;
}
rational yl, yu, xl, xu;
s.am().get_interval(yv, yl, yu, 64);
s.am().get_interval(xv, xl, xu, 64);
rational ymid = (yl + yu) / rational(2); if (ymid.is_neg()) ymid = -ymid;
rational xmid = (xl + xu) / rational(2); if (xmid.is_neg()) xmid = -xmid;
rational scale = std::max(std::max(ymid, xmid), rational(1));
rational growth = rational::power_of_two(std::min(m_retry_atan2[idx], 30u));
rational rdelta = rational(1, 1000000) * scale * growth;
rational ylo = yl - rdelta, yhi = yu + rdelta;
rational xlo = xl - rdelta, xhi = xu + rdelta;
rational pi_lo("3.14159265358979"), pi_hi("3.14159265358980");
rational lower, upper;
if (ylo.is_pos() && xlo.is_pos()) {
lower = rational(0);
upper = pi_hi / rational(2);
}
else if (ylo.is_pos() && xhi.is_neg()) {
lower = pi_lo / rational(2);
upper = pi_hi;
}
else if (yhi.is_neg() && xhi.is_neg()) {
lower = -pi_hi;
upper = -pi_lo / rational(2);
}
else if (yhi.is_neg() && xlo.is_pos()) {
lower = -pi_hi / rational(2);
upper = rational(0);
}
else
return false;
bool below = s.am().lt(vv, lower.to_mpq());
bool above = s.am().gt(vv, upper.to_mpq());
if (!below && !above)
return false;
literal_vector lemma;
lemma.push_back(bound_literal(s, a.y, atom::LT, ylo));
lemma.push_back(bound_literal(s, a.y, atom::GT, yhi));
lemma.push_back(bound_literal(s, a.x, atom::LT, xlo));
lemma.push_back(bound_literal(s, a.x, atom::GT, xhi));
if (below)
lemma.push_back(~bound_literal(s, a.val, atom::LT, lower)); else
lemma.push_back(~bound_literal(s, a.val, atom::GT, upper)); s.mk_clause(lemma.size(), lemma.data(), nullptr);
++m_retry_atan2[idx];
return true;
}
bool transcendentals::refine_monotonicity_pair(app const& a, app const& b) {
anum const& xa = s.value(a.arg);
anum const& xb = s.value(b.arg);
if (a.op == transcendental_op_kind::LOG && (!s.am().is_pos(xa) || !s.am().is_pos(xb)))
return false; anum const& ya = s.value(a.val);
anum const& yb = s.value(b.val);
bool violated = false, a_lt_b = false;
if (s.am().lt(xa, xb) && s.am().ge(ya, yb)) { violated = true; a_lt_b = true; }
else if (s.am().lt(xb, xa) && s.am().ge(yb, ya)) { violated = true; a_lt_b = false; }
if (!violated)
return false;
literal_vector lemma;
if (a.op == transcendental_op_kind::LOG) {
lemma.push_back(~bound_literal(s, a.arg, atom::GT, rational(0))); lemma.push_back(~bound_literal(s, b.arg, atom::GT, rational(0))); }
if (a_lt_b) {
lemma.push_back(~linear_literal(s, a.arg, rational(1), b.arg, rational(-1), atom::LT, rational(0)));
lemma.push_back(linear_literal(s, a.val, rational(1), b.val, rational(-1), atom::LT, rational(0)));
}
else {
lemma.push_back(~linear_literal(s, a.arg, rational(1), b.arg, rational(-1), atom::GT, rational(0)));
lemma.push_back(linear_literal(s, a.val, rational(1), b.val, rational(-1), atom::GT, rational(0)));
}
s.mk_clause(lemma.size(), lemma.data(), nullptr);
return true;
}
bool transcendentals::refine_monotonicity() {
bool added = false;
for (unsigned i = 0; i < m_apps.size(); ++i) {
if (m_apps[i].op != transcendental_op_kind::EXP && m_apps[i].op != transcendental_op_kind::LOG)
continue;
for (unsigned j = i + 1; j < m_apps.size(); ++j) {
if (m_apps[j].op != m_apps[i].op)
continue;
if (refine_monotonicity_pair(m_apps[i], m_apps[j]))
added = true;
}
}
return added;
}
bool transcendentals::refine() {
bool added = false;
for (unsigned i = 0; i < m_apps.size(); ++i)
if (refine_app(i))
added = true;
for (unsigned i = 0; i < m_atan2_apps.size(); ++i)
if (refine_atan2(i))
added = true;
if (refine_monotonicity())
added = true;
return added;
}
}