initial outline of exponentiation

Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>

src/math/lp/nla_core.h

@@ -465,6 +465,8 @@ public:
     
     lbool  check(vector<lemma>& l_vec);
 
+    void set_lemma_vec(vector<lemma>& l_vec) { m_lemma_vec = &l_vec; }
+
     bool  no_lemmas_hold() const;
     
     lbool  test_check(vector<lemma>& l);

src/math/lp/nla_solver.cpp

@@ -13,6 +13,7 @@
 #include "math/lp/factorization.h"
 #include "math/lp/nla_solver.h"
 #include "math/lp/nla_core.h"
+#include "math/polynomial/algebraic_numbers.h"
 
 namespace nla {
 
@@ -54,15 +55,17 @@ solver::~solver() {
     
     std::ostream& solver::display(std::ostream& out) const {    
         m_core->print_monics(out);
-    if( use_nra_model()) {
+        if (use_nra_model()) 
             m_core->m_nra.display(out);
-    }
         return out;
     }
     
     bool solver::use_nra_model() const { return m_core->use_nra_model(); }
+
     core& solver::get_core() { return *m_core; }
+
     nlsat::anum_manager& solver::am() { return m_core->m_nra.am(); }
+
     nlsat::anum const& solver::am_value(lp::var_index v) const {
         SASSERT(use_nra_model());
         return m_core->m_nra.value(v);
@@ -72,4 +75,106 @@ void solver::collect_statistics(::statistics & st) {
         m_core->collect_statistics(st);
     }
 
+    // ensure r = x^y, add abstraction/refinement lemmas
+    lbool solver::check_power(lpvar r, lpvar x, lpvar y, vector<lemma>& lemmas) {
+        if (x == null_lpvar || y == null_lpvar || r == null_lpvar)
+            return l_undef;
+        
+        if (use_nra_model()) 
+            return l_undef;
+
+        auto xval = m_core->val(x);
+        auto yval = m_core->val(y);
+        auto rval = m_core->val(r);
+
+        core& c = get_core();
+        c.set_lemma_vec(lemmas);
+        lemmas.reset();
+
+        // x >= x0 > 0, y >= y0 > 0 => r >= x0^y0
+        // x >= x0 > 0, y <= y0 => r <= x0^y0
+        // x != 0, y = 0 => r = 1
+        // x = 0, y != 0 => r = 0
+        // 
+        // for x fixed, it is exponentiation
+        // => use tangent lemmas and error tolerance.
+
+        if (xval > 0 && yval.is_unsigned()) {
+            auto r2val = power(xval, yval.get_unsigned());
+            if (rval == r2val)
+                return l_true;
+            if (xval != 0 && yval == 0) {
+                new_lemma lemma(c, "x != 0 => x^0 = 1");
+                lemma |= ineq(x, llc::EQ, rational::zero());
+                lemma |= ineq(y, llc::NE, rational::zero());
+                lemma |= ineq(r, llc::EQ, rational::one());
+                return l_false;
+            }
+            if (xval == 0 && yval > 0) {
+                new_lemma lemma(c, "y != 0 => 0^y = 0");
+                lemma |= ineq(x, llc::NE, rational::zero());
+                lemma |= ineq(y, llc::EQ, rational::zero());
+                lemma |= ineq(r, llc::EQ, rational::zero());
+                return l_false;
+            }
+            if (xval > 0 && r2val < rval) {
+                SASSERT(yval > 0);
+                new_lemma lemma(c, "x >= x0 > 0, y >= y0 > 0 => r >= x0^y0");
+                lemma |= ineq(x, llc::LT, xval);
+                lemma |= ineq(y, llc::LT, yval);
+                lemma |= ineq(r, llc::GE, r2val);
+                return l_false;
+            }
+            if (xval > 0 && r2val < rval) {
+                new_lemma lemma(c, "x >= x0 > 0, y <= y0 => r <= x0^y0");
+                lemma |= ineq(x, llc::LT, xval);
+                lemma |= ineq(y, llc::GT, yval);
+                lemma |= ineq(r, llc::LE, r2val);
+                return l_false;
+            }
+        }
+        if (xval > 0 && yval > 0 && !yval.is_int()) {
+            auto ynum = numerator(yval);
+            auto yden = denominator(yval);
+            if (!ynum.is_unsigned())
+                return l_undef;
+            if (!yden.is_unsigned())
+                return l_undef;
+            //   r = x^{yn/yd}
+            // <=>
+            //   r^yd = x^yn
+            auto ryd = power(rval, yden.get_unsigned());
+            auto xyn = power(xval, ynum.get_unsigned());
+            if (ryd == xyn)
+                return l_true;
+#if 0
+            // try some root approximation instead?
+            if (ryd > xyn) {
+                // todo
+            }
+            if (ryd < xyn) {
+                // todo
+            }
+#endif
+
+        }
+
+
+        return l_undef;
+
+        // anum isn't initialized unless nra_solver is invoked.
+        // there is no proviso for using algebraic numbers outside of the nra solver.
+        // so either we have a rational refinement version _and_ an algebraic numeral refinement
+        // loop or we introduce algebraic numerals outside of the nra_solver
+
+#if 0
+        scoped_anum xval(am()), yval(am()), rval(am());
+
+        am().set(xval, am_value(x));
+        am().set(yval, am_value(y));
+        am().set(rval, am_value(r));
+#endif
+
+    }
+
 }

src/math/lp/nla_solver.h

@@ -32,6 +32,7 @@ public:
         void pop(unsigned scopes);
         bool need_check();
         lbool check(vector<lemma>&);
+        lbool check_power(lpvar r, lpvar x, lpvar y, vector<lemma>&);
         bool is_monic_var(lpvar) const;
         bool influences_nl_var(lpvar) const;
         std::ostream& display(std::ostream& out) const;

src/smt/theory_lra.cpp

@@ -797,8 +797,17 @@ class theory_lra::imp {
         return internalize_linearized_def(term, st);
     }
 
+    lpvar get_lpvar(expr* e) {
+        return get_lpvar(get_enode(e));
+    }
+
+    lpvar get_lpvar(enode* n)  {
+        ensure_column(n);
+        return n ? get_lpvar(n->get_th_var(get_id())) : lp::null_lpvar;
+    }
+
     lpvar get_lpvar(theory_var v) const {
-        return lp().external_to_local(v);
+        return v == null_theory_var ? lp::null_lpvar : lp().external_to_local(v);
     }
 
     lp::tv get_tv(theory_var v) const {
@@ -1424,10 +1433,13 @@ public:
         return v != null_theory_var && lp().external_is_used(v);
     }
 
-    void ensure_column(theory_var v) {
+    void ensure_column(enode* n) {
-        if (!lp().external_is_used(v)) {
+        ensure_column(n->get_th_var(get_id()));
-            register_theory_var_in_lar_solver(v);
     }
+
+    void ensure_column(theory_var v) {
+        if (!lp().external_is_used(v) && v != null_theory_var) 
+            register_theory_var_in_lar_solver(v);
     }
 
     mutable vector<std::pair<lp::tv, rational>> m_todo_terms;
@@ -1570,47 +1582,18 @@ public:
     final_check_status eval_power(expr* e) {
         expr* x, * y;
         VERIFY(a.is_power(e, x, y));
-        enode* n = get_enode(e);
-        enode* xn = get_enode(x);
-        enode* yn = get_enode(y);
-        if (!n || !xn || !yn)
-            return FC_GIVEUP;
-        rational valx, valy, valn;
-        theory_var v = n->get_th_var(get_id());
-        if (!get_value(xn, valx))
-            return FC_GIVEUP;
-        if (!get_value(yn, valy))
-            return FC_GIVEUP;
-        if (!get_value(n, valn))
-            return FC_GIVEUP;
 
-        // TBD - check that values align so return FC_DONE.
+        switch (m_nla->check_power(get_lpvar(e), get_lpvar(x), get_lpvar(y), m_nla_lemma_vector)) {
-
+        case l_true:
-        if (valn < 0 && valx > 0 && valy > 0) {
+            return FC_DONE;
-            mk_axiom(mk_literal(a.mk_le(x, a.mk_numeral(rational(0), x->get_sort()))),
+        case l_false:
-                     ~mk_literal(a.mk_le(e, a.mk_numeral(rational(0), e->get_sort()))));
+            for (const nla::lemma & l : m_nla_lemma_vector) 
+                false_case_of_check_nla(l);
             return FC_CONTINUE;
-        }
+        case l_undef:
-
-        // add tangent lemmas for power:
-        // establish equality using algebraic numerals
-
-        // appears to be useless:
-        if (false && !valy.is_int() && numerator(valy) == 1 && denominator(valy) > 0) {
-            rational d = denominator(valy);
-            if (!d.is_unsigned())
             return FC_GIVEUP;
-            unsigned den = d.get_unsigned();
+        default:
-            // e = x^{1/y}
+            break;
-            // => e^y = x
-
-            ptr_vector<expr> es;
-            for (unsigned i = 0; i < den; ++i)
-                es.push_back(e);
-            expr* em = a.mk_mul(es.size(), es.data());            
-            mk_axiom(~mk_literal(m.mk_eq(y, a.mk_real(valy))), mk_literal(m.mk_eq(em, x)));
-            return FC_CONTINUE;
-            
         }
         return FC_GIVEUP;
     }