This commit merges into trunk the branch branches/arithmetic/integers2 from r2650 to r2779.

- This excludes revision 2777.  This revision had some strange performance implications and was delaying the merge.
- This includes the new DioSolver. The DioSolver can discover conflicts, produce substitutions, and produce cuts.
- The DioSolver can be disabled at command line using --disable-dio-solver.
- This includes a number of changes to the arithmetic normal form.
- The Integer class features a number of new number theoretic function.
- This commit includes a few rather loud warning. I will do my best to take care of them today.

5 files changed, 240 insertions, 12 deletions

diff --git a/src/util/bitvector.h b/src/util/bitvector.h
index f05ebaf17..d7f0e13a5 100644
--- a/src/util/bitvector.h
+++ b/src/util/bitvector.h
@@ -91,15 +91,18 @@ public:
   }
 
   BitVector operator ~() const {
+    //is this right? it looks like a no-op?
     return BitVector(d_size, d_value);
   }
 
   BitVector concat (const BitVector& other) const {
-    return BitVector(d_size + other.d_size, (d_value * Integer(2).pow(other.d_size)) + other.d_value);
+    return BitVector(d_size + other.d_size, (d_value.multiplyByPow2(other.d_size)) + other.d_value);
+    //return BitVector(d_size + other.d_size, (d_value * Integer(2).pow(other.d_size)) + other.d_value);
   }
 
   BitVector extract(unsigned high, unsigned low) const {
-    return BitVector(high - low + 1, (d_value % (Integer(2).pow(high + 1))) / Integer(2).pow(low));
+    return BitVector(high - low + 1, d_value.extractBitRange(high - low + 1, low));
+    //return BitVector(high - low + 1, (d_value % (Integer(2).pow(high + 1))) / Integer(2).pow(low));
   }
 
   size_t hash() const {
diff --git a/src/util/integer_cln_imp.h b/src/util/integer_cln_imp.h
index f8ffc0d65..06459e3e1 100644
--- a/src/util/integer_cln_imp.h
+++ b/src/util/integer_cln_imp.h
@@ -167,6 +167,7 @@ public:
     return *this;
   }
 
+  /*
   Integer operator/(const Integer& y) const {
     return Integer( cln::floor1(d_value, y.d_value) );
   }
@@ -182,6 +183,65 @@ public:
     d_value = cln::floor2(d_value, y.d_value).remainder;
     return *this;
   }
+  */
+
+  /**
+   * Return this*(2^pow).
+   */
+  Integer multiplyByPow2(uint32_t pow) const {
+    cln::cl_I ipow(pow);
+    return Integer( d_value << ipow);
+  }
+
+  /** See CLN Documentation. */
+  Integer extractBitRange(uint32_t bitCount, uint32_t low) const {
+    cln::cl_byte range(bitCount, low);
+    return Integer(cln::ldb(d_value, range));
+  }
+
+  /**
+   * Returns the floor(this / y)
+   */
+  Integer floorDivideQuotient(const Integer& y) const {
+    return Integer( cln::floor1(d_value, y.d_value) );
+  }
+
+  /**
+   * Returns r == this - floor(this/y)*y
+   */
+  Integer floorDivideRemainder(const Integer& y) const {
+    return Integer( cln::floor2(d_value, y.d_value).remainder );
+  }
+   /**
+   * Computes a floor quoient and remainder for x divided by y.
+   */
+  static void floorQR(Integer& q, Integer& r, const Integer& x, const Integer& y) {
+    cln::cl_I_div_t res = cln::floor2(x.d_value, y.d_value);
+    q.d_value = res.quotient;
+    r.d_value = res.remainder;
+  }
+
+  /**
+   * Returns the ceil(this / y)
+   */
+  Integer ceilingDivideQuotient(const Integer& y) const {
+    return Integer( cln::ceiling1(d_value, y.d_value) );
+  }
+
+  /**
+   * Returns the ceil(this / y)
+   */
+  Integer ceilingDivideRemainder(const Integer& y) const {
+    return Integer( cln::ceiling2(d_value, y.d_value).remainder );
+  }
+
+  /**
+   * If y divides *this, then exactQuotient returns (this/y)
+   */
+  Integer exactQuotient(const Integer& y) const {
+    Assert(y.divides(*this));
+    return Integer( cln::exquo(d_value, y.d_value) );
+  }
 
   /**
    * Raise this Integer to the power <code>exp</code>.
@@ -208,6 +268,22 @@ public:
   }
 
   /**
+   * Return the least common multiple of this integer with another.
+   */
+  Integer lcm(const Integer& y) const {
+    cln::cl_I result = cln::lcm(d_value, y.d_value);
+    return Integer(result);
+  }
+
+  /**
+   * Return true if *this exactly divides y.
+   */
+  bool divides(const Integer& y) const {
+    cln::cl_I result = cln::rem(y.d_value, d_value);
+    return cln::zerop(result);
+  }
+
+  /**
    * Return the absolute value of this integer.
    */
   Integer abs() const {
@@ -243,6 +319,12 @@ public:
     return output;
   }
 
+  int sgn() const {
+    cln::cl_I sgn = cln::signum(d_value);
+    Assert(sgn == 0 || sgn == -1 || sgn == 1);
+    return cln::cl_I_to_int(sgn);
+  }
+
   //friend std::ostream& operator<<(std::ostream& os, const Integer& n);
 
   long getLong() const {
@@ -281,6 +363,27 @@ public:
     return cln::logbitp(n, d_value);
   }
 
+  /**
+   * If x != 0, returns the unique n s.t. 2^{n-1} <= abs(x) < 2^{n}.
+   * If x == 0, returns 1.
+   */
+  size_t length() const {
+    int s = sgn();
+    if(s == 0){
+      return 1;
+    }else if(s < 0){
+      return cln::integer_length(-d_value);
+    }else{
+      return cln::integer_length(d_value);
+    }
+  }
+
+/*   cl_I xgcd (const cl_I& a, const cl_I& b, cl_I* u, cl_I* v) */
+/* This function ("extended gcd") returns the greatest common divisor g of a and b and at the same time the representation of g as an integral linear combination of a and b: u and v with u*a+v*b = g, g >= 0. u and v will be normalized to be of smallest possible absolute value, in the following sense: If a and b are non-zero, and abs(a) != abs(b), u and v will satisfy the inequalities abs(u) <= abs(b)/(2*g), abs(v) <= abs(a)/(2*g). */
+  static void extendedGcd(Integer& g, Integer& s, Integer& t, const Integer& a, const Integer& b){    
+    g.d_value = cln::xgcd(a.d_value, b.d_value, &s.d_value, &t.d_value);
+  }
+
   friend class CVC4::Rational;
 };/* class Integer */
 
diff --git a/src/util/integer_gmp_imp.h b/src/util/integer_gmp_imp.h
index 16ca8313b..161666df5 100644
--- a/src/util/integer_gmp_imp.h
+++ b/src/util/integer_gmp_imp.h
@@ -135,20 +135,82 @@ public:
     return *this;
   }
 
-  Integer operator/(const Integer& y) const {
-    return Integer( d_value / y.d_value );
+  /**
+   * Return this*(2^pow).
+   */
+  Integer multiplyByPow2(uint32_t pow) const{
+    mpz_class result;
+    mpz_mul_2exp(result.get_mpz_t(), d_value.get_mpz_t(), pow);
+    return Integer( result );
   }
-  Integer& operator/=(const Integer& y) {
-    d_value /= y.d_value;
-    return *this;
+
+  /** See GMP Documentation. */
+  Integer extractBitRange(uint32_t bitCount, uint32_t low) const {
+    // bitCount = high-low+1
+    uint32_t high = low + bitCount-1;
+    //— Function: void mpz_fdiv_r_2exp (mpz_t r, mpz_t n, mp_bitcnt_t b)
+    mpz_class rem, div;
+    mpz_fdiv_r_2exp(rem.get_mpz_t(), d_value.get_mpz_t(), high+1);
+    mpz_fdiv_q_2exp(div.get_mpz_t(), rem.get_mpz_t(), low);
+
+    return Integer(div);
   }
 
-  Integer operator%(const Integer& y) const {
-    return Integer( d_value % y.d_value );
+  /**
+   * Returns the floor(this / y)
+   */
+  Integer floorDivideQuotient(const Integer& y) const {
+    mpz_class q;
+    mpz_fdiv_q(q.get_mpz_t(), d_value.get_mpz_t(), y.d_value.get_mpz_t());
+    return Integer( q );
   }
-  Integer& operator%=(const Integer& y) {
-    d_value %= y.d_value;
-    return *this;
+
+  /**
+   * Returns r == this - floor(this/y)*y
+   */
+  Integer floorDivideRemainder(const Integer& y) const {
+    mpz_class r;
+    mpz_fdiv_r(r.get_mpz_t(), d_value.get_mpz_t(), y.d_value.get_mpz_t());
+    return Integer( r );
+  }
+
+  /**
+   * Computes a floor quoient and remainder for x divided by y.
+   */
+  static void floorQR(Integer& q, Integer& r, const Integer& x, const Integer& y) {
+    mpz_fdiv_qr(q.d_value.get_mpz_t(), r.d_value.get_mpz_t(), x.d_value.get_mpz_t(), y.d_value.get_mpz_t());
+  }
+
+  /**
+   * Returns the ceil(this / y)
+   */
+  Integer ceilingDivideQuotient(const Integer& y) const {
+    mpz_class q;
+    mpz_cdiv_q(q.get_mpz_t(), d_value.get_mpz_t(), y.d_value.get_mpz_t());
+    return Integer( q );
+  }
+
+  /**
+   * Returns the ceil(this / y)
+   */
+  Integer ceilingDivideRemainder(const Integer& y) const {
+    mpz_class r;
+    mpz_cdiv_r(r.get_mpz_t(), d_value.get_mpz_t(), y.d_value.get_mpz_t());
+    return Integer( r );
+  }
+
+  /**
+   * If y divides *this, then exactQuotient returns (this/y)
+   */
+  Integer exactQuotient(const Integer& y) const {
+    Assert(y.divides(*this));
+    mpz_class q;
+    mpz_divexact(q.get_mpz_t(), d_value.get_mpz_t(), y.d_value.get_mpz_t());
+    return Integer( q );
+  }
+
+  int sgn() const {
+    return mpz_sgn(d_value.get_mpz_t());
   }
 
   /**
@@ -172,6 +234,24 @@ public:
   }
 
   /**
+   * Return the least common multiple of this integer with another.
+   */
+  Integer lcm(const Integer& y) const {
+    mpz_class result;
+    mpz_lcm(result.get_mpz_t(), d_value.get_mpz_t(), y.d_value.get_mpz_t());
+    return Integer(result);
+  }
+
+  /**
+   * All non-zero integers z, z.divide(0)
+   * ! zero.divides(zero)
+   */
+  bool divides(const Integer& y) const {
+    int res = mpz_divisible_p(y.d_value.get_mpz_t(), d_value.get_mpz_t());
+    return res != 0;
+  }
+
+  /**
    * Return the absolute value of this integer.
    */
   Integer abs() const {
@@ -217,6 +297,24 @@ public:
     return mpz_tstbit(d_value.get_mpz_t(), n);
   }
 
+  /**
+   * If x != 0, returns the smallest n s.t. 2^{n-1} <= abs(x) < 2^{n}.
+   * If x == 0, returns 1.
+   */
+  size_t length() const {
+    if(sgn() == 0){
+      return 1;
+    }else{
+      return mpz_sizeinbase(d_value.get_mpz_t(),2);
+    }
+  }
+
+  static void extendedGcd(Integer& g, Integer& s, Integer& t, const Integer& a, const Integer& b){
+    //mpz_gcdext (mpz_t g, mpz_t s, mpz_t t, mpz_t a, mpz_t b);
+    mpz_gcdext (g.d_value.get_mpz_t(), s.d_value.get_mpz_t(), t.d_value.get_mpz_t(), a.d_value.get_mpz_t(), b.d_value.get_mpz_t());
+  }
+
+
   friend class CVC4::Rational;
 };/* class Integer */
 
diff --git a/src/util/rational_cln_imp.h b/src/util/rational_cln_imp.h
index 2f2c14ed8..885e6b628 100644
--- a/src/util/rational_cln_imp.h
+++ b/src/util/rational_cln_imp.h
@@ -192,6 +192,18 @@ public:
     }
   }
 
+  Rational abs() const {
+    if(sgn() < 0){
+      return -(*this);
+    }else{
+      return *this;
+    }
+  }
+
+  bool isIntegral() const{
+    return getDenominator() == 1;
+  }
+
   Integer floor() const {
     return Integer(cln::floor1(d_value));
   }
diff --git a/src/util/rational_gmp_imp.h b/src/util/rational_gmp_imp.h
index 37c3c8364..4635ce881 100644
--- a/src/util/rational_gmp_imp.h
+++ b/src/util/rational_gmp_imp.h
@@ -169,6 +169,14 @@ public:
     return mpq_sgn(d_value.get_mpq_t());
   }
 
+  Rational abs() const {
+    if(sgn() < 0){
+      return -(*this);
+    }else{
+      return *this;
+    }
+  }
+
   Integer floor() const {
     mpz_class q;
     mpz_fdiv_q(q.get_mpz_t(), d_value.get_num_mpz_t(), d_value.get_den_mpz_t());
@@ -244,6 +252,10 @@ public:
     return (*this);
   }
 
+  bool isIntegral() const{
+    return getDenominator() == 1;
+  }
+
   /** Returns a string representing the rational in the given base. */
   std::string toString(int base = 10) const {
     return d_value.get_str(base);