Add clean_denominators procedure

Signed-off-by: Leonardo de Moura <leonardo@microsoft.com>

src/math/realclosure/realclosure.cpp

@@ -2416,13 +2416,32 @@ namespace realclosure {
             return new (allocator()) rational_value();
         }
 
-        rational_value * mk_rational(mpq & v) {
+        /**
+           \brief Make a rational and swap its value with v
+        */
+        rational_value * mk_rational_and_swap(mpq & v) {
+            SASSERT(!qm().is_zero(v));
             rational_value * r = mk_rational();
             ::swap(r->m_value, v);
             return r;
         }
 
+        rational_value * mk_rational(mpq const & v) {
+            SASSERT(!qm().is_zero(v));
+            rational_value * r = mk_rational();
+            qm().set(r->m_value, v);
+            return r;
+        }
+
+        rational_value * mk_rational(mpz const & v) {
+            SASSERT(!qm().is_zero(v));
+            rational_value * r = mk_rational();
+            qm().set(r->m_value, v);
+            return r;
+        }
+
         rational_value * mk_rational(mpbq const & v) {
+            SASSERT(!bqm().is_zero(v));
             scoped_mpq v_q(qm()); // v as a rational
             ::to_mpq(qm(), v, v_q);
             return mk_rational(v_q);
@@ -2921,7 +2940,7 @@ namespace realclosure {
                 for (unsigned i = 1; i < sz; i++) {
                     mpq i_mpq(i);
                     value_ref a_i(*this);
-                    a_i = mk_rational(i_mpq);
+                    a_i = mk_rational_and_swap(i_mpq);
                     mul(a_i, p[i], a_i);
                     r.push_back(a_i);
                 }
@@ -3086,30 +3105,42 @@ namespace realclosure {
             return true;
         }
 
+        /**
+           \brief r <- p(b)
+        */
+        void mk_polynomial_value(unsigned n, value * const * p, value * b, value_ref & r) {
+            SASSERT(n > 0);
+            if (n == 1 || b == 0) {
+                r = p[0];
+            }
+            else {
+                SASSERT(n >= 2);
+                // We compute the result using the Horner Sequence
+                //  ((a_{n-1}*b + a_{n-2})*b + a_{n-3})*b + a_{n-4} ...
+                // where a_i's are the coefficients of p.
+                mul(p[n - 1], b, r); // r <- a_{n-1} * b
+                unsigned i = n - 1;
+                while (i > 0) {
+                    --i;
+                    if (p[i] != 0)
+                        add(r, p[i], r); // r <- r + a_i
+                    if (i > 0)
+                        mul(r, b, r); // r <- r * b
+                }
+            }
+        }
+
         /**
            \brief Evaluate the sign of p(b) by computing a value object.
         */
         int expensive_eval_sign_at(unsigned n, value * const * p, mpbq const & b) {
             SASSERT(n > 0);
             SASSERT(p[n - 1] != 0);
-            value_ref bv(*this);
+            value_ref _b(*this);
-            bv = mk_rational(b);
+            _b = mk_rational(b);
-            // We compute the result using the Horner Sequence
+            value_ref pb(*this);
-            //  ((a_{n-1}*bv + a_{n-2})*bv + a_{n-3})*bv + a_{n-4} ...
+            mk_polynomial_value(n, p, _b, pb);
-            // where a_i's are the coefficients of p.
+            return sign(pb);
-            value_ref r(*this);
-            // r <- a_{n-1} * bv
-            mul(p[n - 1], bv, r);
-            unsigned i = n - 1;
-            while (i > 0) {
-                checkpoint();
-                --i;
-                if (p[i] != 0)
-                    add(r, p[i], r); // r <- r + a_i
-                if (i > 0)
-                    mul(r, bv, r); // r <- r * bv
-            }
-            return sign(r);
         }
 
         /**
@@ -4290,7 +4321,7 @@ namespace realclosure {
                 if (qm().is_zero(v))
                     r = 0;
                 else 
-                    r = mk_rational(v);
+                    r = mk_rational_and_swap(v);
             }
             else {
                 INC_DEPTH();
@@ -4319,7 +4350,7 @@ namespace realclosure {
                 if (qm().is_zero(v))
                     r = 0;
                 else 
-                    r = mk_rational(v);
+                    r = mk_rational_and_swap(v);
             }
             else {
                 value_ref neg_b(*this);
@@ -4353,7 +4384,7 @@ namespace realclosure {
                 scoped_mpq v(qm());
                 qm().set(v, to_mpq(a));
                 qm().neg(v);
-                r = mk_rational(v);
+                r = mk_rational_and_swap(v);
             }
             else {
                 neg_rf(to_rational_function(a), r);
@@ -4495,7 +4526,7 @@ namespace realclosure {
             else if (is_nz_rational(a) && is_nz_rational(b)) {
                 scoped_mpq v(qm());
                 qm().mul(to_mpq(a), to_mpq(b), v);
-                r = mk_rational(v);
+                r = mk_rational_and_swap(v);
             }
             else {
                 INC_DEPTH();
@@ -4530,7 +4561,7 @@ namespace realclosure {
             else if (is_nz_rational(a) && is_nz_rational(b)) {
                 scoped_mpq v(qm());
                 qm().div(to_mpq(a), to_mpq(b), v);
-                r = mk_rational(v);
+                r = mk_rational_and_swap(v);
             }
             else {
                 value_ref inv_b(*this);
@@ -4561,7 +4592,7 @@ namespace realclosure {
             if (is_nz_rational(a)) {
                 scoped_mpq v(qm());
                 qm().inv(to_mpq(a), v);
-                r = mk_rational(v);
+                r = mk_rational_and_swap(v);
             }
             else {
                 inv_rf(to_rational_function(a), r);
@@ -4651,6 +4682,280 @@ namespace realclosure {
             return compare(a.m_value, b.m_value);
         }
 
+        // ---------------------------------
+        //
+        // Structural equality
+        //
+        // ---------------------------------
+
+        /**
+           \brief Values a and b are said to be "structurally" equal if:
+           - a and b are 0.
+           - a and b are rationals and compare(a, b) == 0
+           - a and b are rational function values p_a(x)/q_a(x) and p_b(y)/q_b(y) where x and y are field extensions, and
+                  * x == y (pointer equality, i.e., they are the same field extension object).
+                  * Every coefficient of p_a is structurally equal to every coefficient of p_b
+                  * Every coefficient of q_a is structurally equal to every coefficient of q_b
+           Clearly structural equality implies equality, but the reverse is not true.
+        */
+        bool struct_eq(value * a, value * b) const {
+            if (a == b)
+                return true;
+            else if (a == 0 || b == 0)
+                return false;
+            else if (is_nz_rational(a) && is_nz_rational(b))
+                return qm().eq(to_mpq(a), to_mpq(b));
+            else if (is_nz_rational(a) || is_nz_rational(b))
+                return false;
+            else {
+                SASSERT(is_rational_function(a));
+                SASSERT(is_rational_function(b));
+                rational_function_value * rf_a = to_rational_function(a);
+                rational_function_value * rf_b = to_rational_function(b);
+                if (rf_a->ext() != rf_b->ext())
+                    return false;
+                return struct_eq(rf_a->num(), rf_b->num()) && struct_eq(rf_a->den(), rf_b->den());
+            }
+        }
+
+        /**
+           Auxiliary method for 
+                bool struct_eq(value * a, value * b)
+        */
+        bool struct_eq(unsigned sz_a, value * const * p_a, unsigned sz_b, value * const * p_b) const {
+            if (sz_a != sz_b)
+                return false;
+            for (unsigned i = 0; i < sz_a; i++) {
+                if (!struct_eq(p_a[i], p_b[i]))
+                    return false;
+            }
+            return true;
+        }
+
+        /**
+           Auxiliary method for 
+                bool struct_eq(value * a, value * b)
+        */
+        bool struct_eq(polynomial const & p_a, polynomial const & p_b) const {
+            return struct_eq(p_a.size(), p_a.c_ptr(), p_b.size(), p_b.c_ptr());
+        }
+
+        // ---------------------------------
+        //
+        // Clean denominators
+        //
+        // ---------------------------------
+
+        /**
+           \brief We say 'a' has "clean" denominators if
+           - a is 0
+           - a is a rational_value that is an integer
+           - a is a rational_function_value of the form p_a(x)/1 where the coefficients of p_a also have clean denominators.
+        */
+        bool has_clean_denominators(value * a) const {
+            if (a == 0)
+                return true;
+            else if (is_nz_rational(a))
+                return qm().is_int(to_mpq(a));
+            else {
+                rational_function_value * rf_a = to_rational_function(a);
+                return is_rational_one(rf_a->den()) && has_clean_denominators(rf_a->num());
+            }
+        }
+
+        /**
+           \brief See comment at has_clean_denominators(value * a)
+         */
+        bool has_clean_denominators(polynomial const & p) const {
+            unsigned sz = p.size();
+            for (unsigned i = 0; i < sz; i++) {
+                if (!has_clean_denominators(p[i]))
+                    return false;
+            }
+            return true;
+        }
+
+        /**
+           \brief "Clean" the denominators of 'a'. That is, return p and q s.t.
+                     a == p/q
+                  and 
+                     has_clean_denominators(p) and has_clean_denominators(q)
+        */
+        void clean_denominators_core(value * a, value_ref & p, value_ref & q) {
+            INC_DEPTH();
+            TRACE("rcf_clean", tout << "clean_denominators_core [" << m_exec_depth << "]\na: "; display(tout, a, false); tout << "\n";);
+            p.reset(); q.reset();
+            if (a == 0) {
+                p = a; 
+                q = one();
+            }
+            else if (is_nz_rational(a)) {
+                p = mk_rational(to_mpq(a).numerator());
+                q = mk_rational(to_mpq(a).denominator());
+            }
+            else {
+                rational_function_value * rf_a = to_rational_function(a);
+                value_ref_buffer p_num(*this), p_den(*this);
+                value_ref d_num(*this), d_den(*this);
+                clean_denominators_core(rf_a->num(), p_num, d_num);
+                clean_denominators_core(rf_a->den(), p_den, d_den);
+                value_ref x(*this);
+                x = mk_rational_function_value(rf_a->ext());
+                mk_polynomial_value(p_num.size(), p_num.c_ptr(), x, p);
+                mk_polynomial_value(p_den.size(), p_den.c_ptr(), x, q);
+                if (!struct_eq(d_den, d_num)) {
+                    mul(p, d_den, p);
+                    mul(q, d_num, q);
+                }
+            }
+        }
+
+        /**
+           \brief Clean the denominators of the polynomial p, it returns clean_p and d s.t.
+                      p = clean_p/d
+                  and has_clean_denominators(clean_p) && has_clean_denominators(d)
+        */
+        void clean_denominators_core(polynomial const & p, value_ref_buffer & norm_p, value_ref & d) {
+            value_ref_buffer nums(*this), dens(*this);
+            value_ref a_n(*this), a_d(*this);
+            bool all_one = true;
+            for (unsigned i = 0; i < p.size(); i++) {
+                if (p[i]) {
+                    clean_denominators_core(p[i], a_n, a_d);
+                    nums.push_back(a_n);
+                    if (!is_rational_one(a_d))
+                        all_one = false;
+                    dens.push_back(a_d);
+                }
+                else {
+                    nums.push_back(0);
+                    dens.push_back(0);
+                }
+            }
+            if (all_one) {
+                norm_p = nums;
+                d      = one();
+            }
+            else {
+                // Compute lcm of the integer elements in dens.
+                // This is a little trick to control the coefficient growth.
+                // We don't compute lcm of the other elements of dens because it is too expensive.
+                scoped_mpq lcm_z(qm()); 
+                bool found_z = false;
+                SASSERT(nums.size() == p.size());
+                SASSERT(dens.size() == p.size());
+                for (unsigned i = 0; i < p.size(); i++) {
+                    if (!dens[i])
+                        continue;
+                    if (is_nz_rational(dens[i])) {
+                        mpq const & _d = to_mpq(dens[i]);
+                        SASSERT(qm().is_int(_d));
+                        if (!found_z) {
+                            found_z = true;
+                            qm().set(lcm_z, _d);
+                        }
+                        else {
+                            qm().lcm(lcm_z, _d, lcm_z);
+                        }
+                    }
+                }
+            
+                value_ref lcm(*this);
+                if (found_z) {
+                    lcm = mk_rational(lcm_z);
+                }
+                else {
+                    lcm = one();
+                }
+            
+                // Compute the multipliers for nums.
+                // Compute norm_p and d 
+                // 
+                // We do NOT use GCD to compute the LCM of the denominators of non-rational values.
+                // However, we detect structurally equivalent denominators.
+                //
+                //   Thus   a/(b+1) + c/(b+1)  is converted into a*c/(b+1) instead of (a*(b+1) + c*(b+1))/(b+1)^2
+                norm_p.reset();
+                d = lcm;
+                value_ref_buffer multipliers(*this);
+                value_ref m(*this);
+                for (unsigned i = 0; i < p.size(); i++) {
+                    if (!nums[i]) {
+                        norm_p.push_back(0);
+                    }
+                    else {
+                        SASSERT(dens[i]);
+                        bool is_z;
+                        if (!is_nz_rational(dens[i])) {
+                            m = lcm;
+                            is_z = false;
+                        }
+                        else {
+                            scoped_mpq num_z(qm());
+                            qm().div(lcm_z, to_mpq(dens[i]), num_z);
+                            SASSERT(qm().is_int(num_z));
+                            m = mk_rational_and_swap(num_z);
+                            is_z = true;
+                        }
+                        bool found_lt_eq = false;
+                        for (unsigned j = 0; j < p.size(); j++) {
+                            TRACE("rcf_clean_bug", tout << "j: " << j << " "; display(tout, m, false); tout << "\n";);
+                            if (!dens[j])
+                                continue;
+                            if (i != j && !is_nz_rational(dens[j])) {
+                                if (struct_eq(dens[i], dens[j])) {
+                                    if (j < i) 
+                                        found_lt_eq = true;
+                                }
+                                else {
+                                    mul(m, dens[j], m);
+                                }
+                            }
+                        }
+                        if (!is_z && !found_lt_eq) {
+                            mul(dens[i], d, d);
+                        }
+                        mul(m, nums[i], m);
+                        norm_p.push_back(m);
+                    }
+                }
+            }
+            SASSERT(norm_p.size() == p.size());
+        }
+
+        void clean_denominators(value * a, value_ref & p, value_ref & q) {
+            if (has_clean_denominators(a)) {
+                p = a;
+                q = one();
+            }
+            else {
+                clean_denominators_core(a, p, q);
+            }
+        }
+        
+        void clean_denominators(polynomial const & p, value_ref_buffer & norm_p, value_ref & d) {
+            if (has_clean_denominators(p)) {
+                norm_p.append(p.size(), p.c_ptr());
+                d = one();
+            }
+            else {
+                clean_denominators_core(p, norm_p, d);
+            }
+        }
+
+        void clean_denominators(numeral const & a, numeral & p, numeral & q) {
+            value_ref _p(*this), _q(*this);
+            clean_denominators(a.m_value, _p, _q);
+            set(p, _p);
+            set(q, _q);
+        }
+
+        // ---------------------------------
+        //
+        // "Pretty printing"
+        //
+        // ---------------------------------
+
         struct collect_algebraic_refs {
             char_vector            m_visited; // Set of visited algebraic extensions.
             ptr_vector<algebraic>  m_found;   // vector/list of visited algebraic extensions.
@@ -4683,11 +4988,20 @@ namespace realclosure {
             }
         };
 
+        static unsigned num_nz_coeffs(polynomial const & p) {
+            unsigned r = 0;
+            for (unsigned i = 0; i < p.size(); i++) {
+                if (p[i])
+                    r++;
+            }
+            return r;
+        }
+
         bool use_parenthesis(value * v) const {
             if (is_zero(v) || is_nz_rational(v)) 
                 return false;
             rational_function_value * rf = to_rational_function(v);
-            return rf->num().size() > 1 || !is_rational_one(rf->den());
+            return num_nz_coeffs(rf->num()) > 1 || !is_rational_one(rf->den());
         }
 
         template<typename DisplayVar>
@@ -5140,6 +5454,10 @@ namespace realclosure {
         m_imp->display_interval(out, a);
     }
 
+    void manager::clean_denominators(numeral const & a, numeral & p, numeral & q) {
+        save_interval_ctx ctx(this);
+        m_imp->clean_denominators(a, p, q);
+    }
 };
 
 void pp(realclosure::manager::imp * imp, realclosure::polynomial const & p, realclosure::extension * ext) {
@@ -5162,6 +5480,10 @@ void pp(realclosure::manager::imp * imp, realclosure::manager::imp::value_ref_bu
         pp(imp, p[i]);
 }
 
+void pp(realclosure::manager::imp * imp, realclosure::manager::imp::value_ref const & v) {
+    pp(imp, v.get());
+}
+
 void pp(realclosure::manager::imp * imp, realclosure::mpbqi const & i) {
     imp->bqim().display(std::cout, i);
     std::cout << std::endl;

src/math/realclosure/realclosure.h

@@ -263,6 +263,8 @@ namespace realclosure {
 
         
         void display_interval(std::ostream & out, numeral const & a) const;
+    
+        void clean_denominators(numeral const & a, numeral & p, numeral & q);
     };
 
     class value;

src/test/rcf.cpp

@@ -135,15 +135,36 @@ static void tst_lin_indep(unsigned m, unsigned n, int _A[], unsigned ex_sz, unsi
             A.set(i, j, _A[i*n + j]);
     unsigned_vector r;
     r.resize(A.n());
-    mm.linear_independent_rows(A, r.c_ptr());
+    scoped_mpz_matrix B(mm);
+    mm.linear_independent_rows(A, r.c_ptr(), B);
     SASSERT(r.size() == ex_sz);
     for (unsigned i = 0; i < ex_sz; i++) {
         SASSERT(r[i] == ex_r[i]);
     }
 }
 
+static void tst_denominators() {
+    unsynch_mpq_manager qm;
+    rcmanager m(qm);
+    scoped_rcnumeral a(m);
+    scoped_rcnumeral t(m);
+    scoped_rcnumeral eps(m);
+    m.mk_pi(a);
+    m.inv(a);
+    m.mk_infinitesimal("eps", eps);
+    t = (a - eps*2) / (a*eps + 1);
+    // t = t + a * 2;
+    scoped_rcnumeral n(m), d(m);
+    std::cout << t << "\n";
+    m.clean_denominators(t, n, d);
+    std::cout << "---->\n" << n << "\n" << d << "\n";
+}
 
 void tst_rcf() {
+    enable_trace("rcf_clean");
+    enable_trace("rcf_clean_bug");
+    tst_denominators();
+    return;
     tst1();
     tst2();
     { int A[] = {0, 1, 1, 1, 0, 1, 1, 1, -1}; int c[] = {10, 4, -4}; int b[] = {-2, 4, 6}; tst_solve(3, A, b, c, true); }