// treat src as a homogeneous matrix.
 void dualizeH(matrix& dst, matrix const& src) {
     hilbert_basis hb;
     for (unsigned i = 0; i < src.size(); ++i) {
         vector<rational> v(src.A[i]);
         v.push_back(src.b[i]);
         hb.add_eq(v, rational(0));
     }
     for (unsigned i = 0; i < 1 + src.A[0].size(); ++i) {
         hb.set_is_int(i);
     }
     lbool is_sat = hb.saturate();
     hb.display(std::cout);
     SASSERT(is_sat == l_true);
     dst.reset();
     unsigned basis_size = hb.get_basis_size();
     for (unsigned i = 0; i < basis_size; ++i) {
         bool is_initial;
         vector<rational> soln;
         hb.get_basis_solution(i, soln, is_initial);
         if (!is_initial) {
             dst.b.push_back(soln.back());
             soln.pop_back();
             dst.A.push_back(soln);
         }
     }
 }

//
// Find rationals c, x, such that c + epsilon*x <= r1/r2
//
// let r1 = a + d_1
// let r2 = b + d_2
//
// suppose b != 0:
//
//      r1/b <= r1/r2
// <=> { if b > 0, then r2 > 0, and cross multiplication does not change the sign }
//     { if b < 0, then r2 < 0, and cross multiplication changes sign twice }
//      r1 * r2 <= b * r1
// <=>
//      r1 * (b + d_2) <= r1 * b
// <=>
//      r1 * d_2 <= 0
//
// if r1 * d_2 > 0, then r1/(b + sign_of(r1)*1/2*|b|) <= r1/r2
//
// Not handled here:
// if b = 0, then d_2 != 0
//   if r1 * d_2 = 0 then it's 0.
//   if r1 * d_2 > 0, then result is +oo
//   if r1 * d_2 < 0, then result is -oo
//
inf_rational inf_div(inf_rational const& r1, inf_rational const& r2)
{
    SASSERT(!r2.m_first.is_zero());
    inf_rational result;

    if (r2.m_second.is_neg() && r1.is_neg()) {
        result = r1 / (r2.m_first - (abs(r2.m_first)/rational(2)));
    }
    else if (r2.m_second.is_pos() && r1.is_pos()) {
        result = r1 / (r2.m_first + (abs(r2.m_first)/rational(2)));
    }
    else {
        result = r1 / r2.m_first;
    }
    return result;
}

 app * process_le_ge(func_decl * f, expr * arg1, expr * arg2, bool le) {
     expr * v;
     expr * t;
     if (uncnstr(arg1)) {
         v = arg1;
         t = arg2;
     }
     else if (uncnstr(arg2)) {
         v = arg2;
         t = arg1;
         le = !le;
     }
     else {
         return nullptr;
     }
     app * u;
     if (!mk_fresh_uncnstr_var_for(f, arg1, arg2, u))
         return u;
     if (!m_mc)
         return u;
     // v = ite(u, t, t + 1) if le
     // v = ite(u, t, t - 1) if !le
     add_def(v, m().mk_ite(u, t, m_a_util.mk_add(t, m_a_util.mk_numeral(rational(le ? 1 : -1), m().get_sort(arg1)))));
     return u;
 }

    void transition(
        matrix& dst,
        matrix const& src,
        matrix const& Ab) {
        matrix T;
        // length of rows in Ab are twice as long as
        // length of rows in src.
        SASSERT(2*src.A[0].size() == Ab.A[0].size());
        vector<rational> zeros;
        for (unsigned i = 0; i < src.A[0].size(); ++i) {
            zeros.push_back(rational(0));
        }
        for (unsigned i = 0; i < src.size(); ++i) {
            T.A.push_back(src.A[i]);
            T.A.back().append(zeros);            
        }
        T.b.append(src.b);
        T.append(Ab);

        T.display(std::cout << "T:\n");
        matrix TD;
        dualizeI(TD, T);
        TD.display(std::cout << "TD:\n");
        for (unsigned i = 0; i < TD.size(); ++i) {
            vector<rational> v;            
            v.append(src.size(), TD.A[i].c_ptr() + src.size());
            dst.A.push_back(v);
            dst.b.push_back(TD.b[i]);
        }
        dst.display(std::cout << "dst\n");
    }

 app * process_bv_mul(func_decl * f, unsigned num, expr * const * args) {
     if (num == 0)
         return nullptr;
     if (uncnstr(num, args)) {
         sort * s = m().get_sort(args[0]);
         app * r;
         if (!mk_fresh_uncnstr_var_for(f, num, args, r))
             return r;
         if (m_mc)
             add_defs(num, args, r, m_bv_util.mk_numeral(rational(1), s));
         return r;
     }
     // c * v (c is even) case
     unsigned bv_size;
     rational val;
     rational inv;
     if (num == 2 && 
         uncnstr(args[1]) && 
         m_bv_util.is_numeral(args[0], val, bv_size) &&
         m_bv_util.mult_inverse(val, bv_size, inv)) {
         app * r;
         if (!mk_fresh_uncnstr_var_for(f, num, args, r))
             return r;
         sort * s = m().get_sort(args[1]);
         if (m_mc)
             add_def(args[1], m_bv_util.mk_bv_mul(m_bv_util.mk_numeral(inv, s), r));
         return r;
     }
     return nullptr;
 }

cons_t* proc_add(cons_t *p, environment_t* env)
{
  /*
   * Integers have an IDENTITY, so we can do this,
   * but a more correct approach would be to take
   * the value of the FIRST number we find and
   * return that.
   */
  rational_t sum;
  sum.numerator = 0;
  sum.denominator = 1;
  bool exact = true;

  for ( ; !nullp(p); p = cdr(p) ) {
    cons_t *i = listp(p)? car(p) : p;

    if ( integerp(i) ) {
      if ( !i->number.exact ) exact = false;
      sum += i->number.integer;
    } else if ( rationalp(i) ) {
      if ( !i->number.exact ) exact = false;
      sum += i->number.rational;
    } else if ( realp(i) ) {
      // automatically convert; perform rest of computation in floats
      exact = false;
      return proc_addf(cons(real(sum), p), env);
    } else
      raise(runtime_exception(
        "Cannot add integer with " + to_s(type_of(i)) + ": " + sprint(i)));
  }

  return rational(sum, exact);
}

 void add_ineq() {
     pb_util pb(m);
     expr_ref fml(m), tmp(m);
     th_rewriter rw(m);
     vector<rational> coeffs(vars.size());
     expr_ref_vector args(vars);
     while (true) {
         rational k(rand(6));
         for (unsigned i = 0; i < coeffs.size(); ++i) {
             int v = 3 - rand(5);
             coeffs[i] = rational(v);
             if (coeffs[i].is_neg()) {
                 args[i] = m.mk_not(args[i].get());
                 coeffs[i].neg();
                 k += coeffs[i];
             }
         }       
         fml = pb.mk_ge(args.size(), coeffs.c_ptr(), args.c_ptr(), k);
         rw(fml, tmp);
         rw(tmp, tmp);
         if (pb.is_ge(tmp)) {
             fml = tmp;
             break;
         }
     }
     std::cout << "(assert " << fml << ")\n";
     ctx.assert_expr(fml);
 }

Z3_ast mk_max(Z3_context ctx, Z3_sort bv, bool is_signed) {
    unsigned bvsize = Z3_get_bv_sort_size(ctx, bv);
    unsigned sz = is_signed ? bvsize - 1 : bvsize;
    rational max_bound = power(rational(2), sz);
    --max_bound;
    return Z3_mk_numeral(ctx, max_bound.to_string().c_str(), bv);
}

cons_t* proc_mul(cons_t *p, environment_t *env)
{
  rational_t product;
  product.numerator = 1;
  product.denominator = 1;
  bool exact = true;

  for ( ; !nullp(p); p = cdr(p) ) {
    cons_t *i = listp(p)? car(p) : p;

    if ( integerp(i) ) {
      product *= i->number.integer;
      if ( !i->number.exact ) exact = false;
    } else if ( rationalp(i) ) {
      if ( !i->number.exact ) exact = false;
      product *= i->number.rational;
    } else if ( realp(i) ) {
      // automatically convert; perform rest of computation in floats
      exact = false;
      return proc_mulf(cons(real(product), p), env);
    } else
      raise(runtime_exception("Cannot multiply integer with " + to_s(type_of(i)) + ": " + sprint(i)));
  }

  return rational(product, exact);
}

rational pb_util::to_rational(parameter const& p) const {
    if (p.is_int()) {
        return rational(p.get_int());
    }
    SASSERT(p.is_rational());
    return p.get_rational();
}

 sort * dl_decl_plugin::mk_relation_sort( unsigned num_parameters, parameter const * parameters) {
             bool is_finite = true;
     rational r(1);
     for (unsigned i = 0; is_finite && i < num_parameters; ++i) {
         if (!parameters[i].is_ast() || !is_sort(parameters[i].get_ast())) {
             m_manager->raise_exception("expecting sort parameters");
             return 0;
         }
         sort* s = to_sort(parameters[i].get_ast());
         sort_size sz1 = s->get_num_elements();
         if (sz1.is_finite()) {
             r *= rational(sz1.size(),rational::ui64());
         }
         else {
             is_finite = false;
         }
     }
     sort_size sz;
     if (is_finite && r.is_uint64()) {
         sz = sort_size::mk_finite(r.get_uint64());
     }
     else {
         sz = sort_size::mk_very_big();
     }
     sort_info info(m_family_id, DL_RELATION_SORT, sz, num_parameters, parameters);
     return m_manager->mk_sort(symbol("Table"),info);
 }

rational operator * ( const rational& r1, const rational& r2 )
{
    int a = r1.num * r2.num;
    int b = r1.denum * r2.denum;
    return rational(a,b);

}

cons_t* parse_exact_real(const char* sc, int radix)
{
  if ( radix != 10 )
    raise(parser_exception(
      "Only reals with decimal radix are supported"));

  /*
   * Since the real is already in string form, we can simply turn it into a
   * rational number.
   */
  char *s = strdup(sc);
  char *d = strchr(s, '.');
  *d = '\0';
  const char* left = s;
  const char* right = d+1;

  int decimals = strlen(right);

  /*
   * NOTE: If we overflow here, we're in big trouble.
   * TODO: Throw an error if we overflow.  Or just implement bignums.
   */
  rational_t r;
  r.numerator = to_i(left, radix)*pow10(decimals) + to_i(right, radix);
  r.denominator = pow10(decimals);

  free(s);
  return rational(r, true);
}

        virtual void operator()(model_ref& md) {
            model_ref old_model = alloc(model, m);
            obj_map<func_decl, func_decl*>::iterator it  = m_new2old.begin();
            obj_map<func_decl, func_decl*>::iterator end = m_new2old.end();
            for (; it != end; ++it) {
                func_decl* old_p = it->m_value;
                func_decl* new_p = it->m_key;
                func_interp* old_fi = alloc(func_interp, m, old_p->get_arity());

                if (new_p->get_arity() == 0) {
                    old_fi->set_else(md->get_const_interp(new_p));
                }
                else {
                    func_interp* new_fi = md->get_func_interp(new_p);
                    expr_ref_vector subst(m);
                    var_subst vs(m, false);
                    expr_ref tmp(m);

                    if (!new_fi) {
                        TRACE("dl", tout << new_p->get_name() << " has no value in the current model\n";);
                        dealloc(old_fi);
                        continue;
                    }
                    for (unsigned i = 0; i < old_p->get_arity(); ++i) {
                        subst.push_back(m.mk_var(i, old_p->get_domain(i)));
                    }
                    subst.push_back(a.mk_numeral(rational(1), a.mk_real()));

                    // Hedge that we don't have to handle the general case for models produced
                    // by Horn clause solvers.
                    SASSERT(!new_fi->is_partial() && new_fi->num_entries() == 0);
                    vs(new_fi->get_else(), subst.size(), subst.c_ptr(), tmp);
                    old_fi->set_else(tmp);
                    old_model->register_decl(old_p, old_fi);
                }

 app * process_bv_le(func_decl * f, expr * arg1, expr * arg2, bool is_signed) {
     if (m_produce_proofs) {
         // The result of bv_le is not just introducing a new fresh name,
         // we need a side condition.
         // TODO: the correct proof step
         return nullptr;
     }
     if (uncnstr(arg1)) {
         // v <= t
         expr * v = arg1;
         expr * t = arg2;
         // v <= t --->  (u or t == MAX)   u is fresh
         //     add definition v = ite(u or t == MAX, t, t+1)
         unsigned bv_sz = m_bv_util.get_bv_size(arg1);
         rational MAX;
         if (is_signed)
             MAX = rational::power_of_two(bv_sz - 1) - rational(1);
         else
             MAX = rational::power_of_two(bv_sz) - rational(1);
         app * u;
         bool is_new = mk_fresh_uncnstr_var_for(f, arg1, arg2, u);
         app * r = m().mk_or(u, m().mk_eq(t, m_bv_util.mk_numeral(MAX, bv_sz)));
         if (m_mc && is_new)
             add_def(v, m().mk_ite(r, t, m_bv_util.mk_bv_add(t, m_bv_util.mk_numeral(rational(1), bv_sz))));
         return r;
     }
     if (uncnstr(arg2)) {
         // v >= t
         expr * v = arg2;
         expr * t = arg1;
         // v >= t --->  (u ot t == MIN)  u is fresh
         //    add definition v = ite(u or t == MIN, t, t-1)
         unsigned bv_sz = m_bv_util.get_bv_size(arg1);
         rational MIN;
         if (is_signed)
             MIN = -rational::power_of_two(bv_sz - 1);
         else
             MIN = rational(0);
         app * u;
         bool is_new = mk_fresh_uncnstr_var_for(f, arg1, arg2, u);
         app * r = m().mk_or(u, m().mk_eq(t, m_bv_util.mk_numeral(MIN, bv_sz)));
         if (m_mc && is_new)
             add_def(v, m().mk_ite(r, t, m_bv_util.mk_bv_sub(t, m_bv_util.mk_numeral(rational(1), bv_sz))));
         return r;
     }
     return nullptr;
 }

    bool utvpi_tester::linearize() {

        m_weight.reset();
        m_coeff_map.reset();

        while (!m_terms.empty()) {
            expr* e1, *e2;
            rational num;
            rational mul = m_terms.back().second;
            expr* e = m_terms.back().first;
            m_terms.pop_back();
            if (a.is_add(e)) {
                for (unsigned i = 0; i < to_app(e)->get_num_args(); ++i) {
                    m_terms.push_back(std::make_pair(to_app(e)->get_arg(i), mul));
                }
            }
            else if (a.is_mul(e, e1, e2) && a.is_numeral(e1, num)) {
                m_terms.push_back(std::make_pair(e2, mul*num));
            }
            else if (a.is_mul(e, e2, e1) && a.is_numeral(e1, num)) {
                m_terms.push_back(std::make_pair(e2, mul*num));
            }
            else if (a.is_sub(e, e1, e2)) {
                m_terms.push_back(std::make_pair(e1, mul));
                m_terms.push_back(std::make_pair(e2, -mul));                
            }
            else if (a.is_uminus(e, e1)) {
                m_terms.push_back(std::make_pair(e1, -mul));
            }
            else if (a.is_numeral(e, num)) {
                m_weight += num*mul;
            }
            else if (a.is_to_real(e, e1)) {
                m_terms.push_back(std::make_pair(e1, mul));
            }
            else if (!is_uninterp_const(e)) {
                return false;
            }
            else {
                m_coeff_map.insert_if_not_there2(e, rational(0))->get_data().m_value += mul;
            }
        }
        obj_map<expr, rational>::iterator it  = m_coeff_map.begin();
        obj_map<expr, rational>::iterator end = m_coeff_map.end();
        for (; it != end; ++it) {
            rational r = it->m_value;
            if (r.is_zero()) {
                continue;
            }
            m_terms.push_back(std::make_pair(it->m_key, r));
            if (m_terms.size() > 2) {
                return false;
            }
            if (!r.is_one() && !r.is_minus_one()) {
                return false;
            }
        }
        return true;
    }

Z3_ast mk_min(Z3_context ctx, Z3_sort bv, bool is_signed) {
    unsigned bvsize = Z3_get_bv_sort_size(ctx, bv);
    if (! is_signed) return Z3_mk_numeral(ctx, "0", bv);
    unsigned sz = bvsize - 1;
    rational min_bound = power(rational(2), sz);
    min_bound.neg();
    return Z3_mk_numeral(ctx, min_bound.to_string().c_str(), bv);
}

 Z3_ast Z3_API Z3_mk_pbeq(Z3_context c, unsigned num_args, 
                          Z3_ast const args[], int _coeffs[],
                          int k) {
     Z3_TRY;
     LOG_Z3_mk_pble(c, num_args, args, _coeffs, k);
     RESET_ERROR_CODE();
     pb_util util(mk_c(c)->m());
     vector<rational> coeffs;
     for (unsigned i = 0; i < num_args; ++i) {
         coeffs.push_back(rational(_coeffs[i]));
     }
     ast* a = util.mk_eq(num_args, coeffs.c_ptr(), to_exprs(args), rational(k));
     mk_c(c)->save_ast_trail(a);
     check_sorts(c, a);
     RETURN_Z3(of_ast(a));
     Z3_CATCH_RETURN(0);
 }

br_status fpa_rewriter::mk_to_real_unspecified(expr_ref & result) {
    if (m_hi_fp_unspecified)
        // The "hardware interpretation" is 0.
        result = m_util.au().mk_numeral(rational(0), false);
    else
        result = m_util.mk_internal_to_real_unspecified();

    return BR_DONE;
}

void bv2int_rewriter_ctx::collect_power2(goal const& s) {
    ast_manager& m = m_trail.get_manager();
    arith_util arith(m);
    bv_util bv(m);
    
    for (unsigned j = 0; j < s.size(); ++j) {
        expr* f = s.form(j);
        if (!m.is_or(f)) continue;
        unsigned sz = to_app(f)->get_num_args();
        expr* x, *y, *v = 0;
        rational n;
        vector<rational> bounds;
        bool is_int, ok = true;

        for (unsigned i = 0; ok && i < sz; ++i) {
            expr* e = to_app(f)->get_arg(i);
            if (!m.is_eq(e, x, y)) {
                ok = false;
                break;
            }
            if (arith.is_numeral(y, n, is_int) && is_int &&
                (x == v || v == 0)) {
                v = x;
                bounds.push_back(n);
            }
            else if (arith.is_numeral(x, n, is_int) && is_int &&
                     (y == v || v == 0)) {
                v = y;
                bounds.push_back(n);
            }
            else {
                ok = false;
                break;
            }
        }
        if (!ok || !v) continue;
        SASSERT(!bounds.empty());
        lt_rational lt;
        // lt is a total order on rationals.
        std::sort(bounds.begin(), bounds.end(), lt);
        rational p(1);
        unsigned num_bits = 0;
        for (unsigned i = 0; ok && i < bounds.size(); ++i) {
            ok = (p == bounds[i]); 
            p *= rational(2);
            ++num_bits;
        }
        if (!ok) continue;
        unsigned log2 = 0;
        for (unsigned i = 1; i <= num_bits; i *= 2) ++log2;
        if(log2 == 0) continue;
        expr* logx = m.mk_fresh_const("log2_v", bv.mk_sort(log2));
        logx = bv.mk_zero_extend(num_bits - log2, logx);
        m_trail.push_back(logx);
        TRACE("bv2int_rewriter", tout << mk_pp(v, m) << " |-> " << mk_pp(logx, m) << "\n";);
        m_power2.insert(v, logx);
    }