Merge branch 'unstable' of https://github.com/wintersteiger/z3 into unstable

src/util/bit_util.cpp

@@ -351,7 +351,7 @@ bool has_one_at_first_k_bits(unsigned sz, unsigned const * data, unsigned k) {
     }
     if (word_sz < sz) {
         unsigned bit_sz  = k % (8 * sizeof(unsigned));
-        unsigned mask    = (1 << bit_sz) - 1;
+        unsigned mask    = (1u << bit_sz) - 1;
         return (data[word_sz] & mask) != 0;
     }
     return false;

src/util/bit_vector.cpp

@@ -206,6 +206,34 @@ void bit_vector::display(std::ostream & out) const {
 #endif
 }
 
+bool bit_vector::contains(bit_vector const& other) const {
+    unsigned n = num_words();
+    if (n == 0)
+        return true;
+    
+    for (unsigned i = 0; i < n - 1; ++i) {
+        if ((m_data[i] & other.m_data[i]) != other.m_data[i])
+            return false;
+    }
+    unsigned bit_rest = m_num_bits % 32;
+    unsigned mask = (1U << bit_rest) - 1;
+    if (mask == 0) mask = UINT_MAX;
+    unsigned other_data = other.m_data[n-1] & mask;
+    return (m_data[n-1] & other_data) == other_data;
+}
+
+unsigned bit_vector::get_hash() const {
+    return string_hash(reinterpret_cast<char const* const>(m_data), size()/8,  0);
+}
+
+bit_vector& bit_vector::neg() {
+    unsigned n = num_words();
+    for (unsigned i = 0; i < n; ++i) {
+        m_data[i] = ~m_data[i];
+    }
+    return *this;
+}
+
 void fr_bit_vector::reset() {
     unsigned sz = size();
     unsigned_vector::const_iterator it  = m_one_idxs.begin();

src/util/bit_vector.h

@@ -75,8 +75,11 @@ public:
     bit_vector(bit_vector const & source):
         m_num_bits(source.m_num_bits),
         m_capacity(source.m_capacity),
-        m_data(alloc_svect(unsigned, source.m_capacity)) {
-        memcpy(m_data, source.m_data, source.m_capacity * sizeof(unsigned));
+        m_data(0) {
+        if (source.m_data) {
+            m_data = alloc_svect(unsigned, m_capacity);
+            memcpy(m_data, source.m_data, m_capacity * sizeof(unsigned));
+        }
     }
     
     bit_vector(unsigned const * source, int num_bits):
@@ -91,7 +94,8 @@ public:
     }
     
     void reset() {
-        memset(m_data, 0, m_capacity * sizeof(unsigned));
+        if (m_data)
+            memset(m_data, 0, m_capacity * sizeof(unsigned));
         m_num_bits = 0;
     }
 
@@ -123,6 +127,8 @@ public:
     unsigned get_word(unsigned word_idx) const {
         return m_data[word_idx];
     }
+
+    unsigned get_hash() const;
     
     bool get(unsigned bit_idx) const {
         SASSERT(bit_idx < size());
@@ -184,6 +190,9 @@ public:
 
     bit_vector & operator=(bit_vector const & source) {
         m_num_bits = source.m_num_bits;
+        if (!source.m_data)
+            return *this;
+
         if (m_capacity < source.m_capacity) {
             dealloc_svect(m_data);
             m_data     = alloc_svect(unsigned, source.m_capacity);
@@ -196,8 +205,13 @@ public:
     bit_vector & operator|=(bit_vector const & source);
 
     bit_vector & operator&=(bit_vector const & source);
+
+    bit_vector & neg();
     
     void display(std::ostream & out) const;
+
+    bool contains(const bit_vector & other) const;
+
 };
 
 inline std::ostream & operator<<(std::ostream & out, bit_vector const & b) {

src/util/dependency.h

@@ -171,7 +171,7 @@ public:
             m_todo.push_back(d);
             unsigned qhead = 0;
             while (qhead < m_todo.size()) {
-                dependency * d = m_todo[qhead];
+                d = m_todo[qhead];
                 qhead++;
                 if (d->is_leaf()) {
                     if (to_leaf(d)->m_value == v) {
@@ -201,7 +201,7 @@ public:
             m_todo.push_back(d);
             unsigned qhead = 0;
             while (qhead < m_todo.size()) {
-                dependency * d = m_todo[qhead];
+	        d = m_todo[qhead];
                 qhead++;
                 if (d->is_leaf()) {
                     vs.push_back(to_leaf(d)->m_value);

src/util/fixed_bit_vector.cpp

@@ -0,0 +1,168 @@
+/*++
+Copyright (c) 2014 Microsoft Corporation
+
+Module Name:
+
+    fixed_bit_vector.cpp
+
+Abstract:
+
+    Simple bitvector implementation for fixed size bit-vectors.
+
+Author:
+
+    Nikolaj Bjorner (nbjorner) 2014-9-15.
+    Leonardo de Moura (leonardo) 2006-10-03.
+
+Revision History:
+    Based on bit_vector.cpp
+
+--*/
+
+#include<limits.h>
+#include"fixed_bit_vector.h"
+#include"trace.h"
+#include"hash.h"
+
+void fixed_bit_vector::set(fixed_bit_vector const& other, unsigned hi, unsigned lo) {
+    if ((lo % 32) == 0) {
+        unsigned sz32 = (hi-lo+1)/32;
+        unsigned lo32 = lo/32;
+        for (unsigned i = 0; i < sz32; ++i) {
+            m_data[lo32 + i] = other.m_data[i];
+        }
+        for (unsigned i = sz32*32; i < hi - lo + 1; ++i) {
+            set(lo + i, other.get(i));
+        }
+        return;
+    }
+    for (unsigned i = 0; i < hi - lo + 1; ++i) {
+        set(lo + i, other.get(i));
+    }
+}
+
+fixed_bit_vector_manager::fixed_bit_vector_manager(unsigned num_bits):
+    m_alloc("fixed_bit_vector") {
+    m_num_bits = num_bits;
+    m_num_words = num_words(num_bits);
+    m_num_bytes = m_num_words * sizeof(unsigned);
+    unsigned bit_rest = m_num_bits % 32;
+    m_mask = (1U << bit_rest) - 1;
+    if (m_mask == 0) m_mask = UINT_MAX;
+}
+
+
+fixed_bit_vector* fixed_bit_vector_manager::allocate() {
+    if (m_num_bytes == 0) return &m_0;
+    return static_cast<fixed_bit_vector*>(m_alloc.allocate(m_num_bytes));
+}
+
+fixed_bit_vector* fixed_bit_vector_manager::allocate0() {
+    fixed_bit_vector* result = allocate();
+    fill0(*result);
+    return result;
+}
+
+fixed_bit_vector* fixed_bit_vector_manager::allocate1() {
+    fixed_bit_vector* result = allocate();
+    fill1(*result);
+    return result;
+}
+
+fixed_bit_vector* fixed_bit_vector_manager::allocate(fixed_bit_vector const& bv) {
+    fixed_bit_vector* result = allocate();
+    copy(*result, bv);
+    return result;
+}
+
+void fixed_bit_vector_manager::deallocate(fixed_bit_vector* bv) {
+    if (m_num_bytes > 0) m_alloc.deallocate(m_num_bytes, bv);
+}
+
+
+void fixed_bit_vector_manager::copy(fixed_bit_vector& dst, fixed_bit_vector const& src) const {
+    memcpy(dst.m_data, src.m_data, num_bytes());
+}
+
+
+fixed_bit_vector& 
+fixed_bit_vector_manager::fill0(fixed_bit_vector& bv) const {
+    memset(bv.m_data, 0, num_bytes());
+    return bv;
+}
+
+fixed_bit_vector& 
+fixed_bit_vector_manager::fill1(fixed_bit_vector& bv) const {
+    memset(bv.m_data, 0xFF, num_bytes());
+    return bv;
+}
+
+fixed_bit_vector& 
+fixed_bit_vector_manager::set_and(fixed_bit_vector& dst, fixed_bit_vector const& src) const {
+    for (unsigned i = 0; i < m_num_words; i++) 
+        dst.m_data[i] &= src.m_data[i];
+    return dst;
+}
+
+fixed_bit_vector& 
+fixed_bit_vector_manager::set_or(fixed_bit_vector& dst,  fixed_bit_vector const& src) const {
+    for (unsigned i = 0; i < m_num_words; i++) 
+        dst.m_data[i] |= src.m_data[i];
+    return dst;
+}
+
+fixed_bit_vector& 
+fixed_bit_vector_manager::set_neg(fixed_bit_vector& dst) const {
+    for (unsigned i = 0; i < m_num_words; i++) 
+        dst.m_data[i] = ~dst.m_data[i];
+    return dst;
+}
+
+unsigned fixed_bit_vector_manager::last_word(fixed_bit_vector const& bv) const {
+    unsigned n = num_words();
+    if (n == 0) return 0;
+    return bv.m_data[n-1] & m_mask;
+}
+
+bool fixed_bit_vector_manager::equals(fixed_bit_vector const& a, fixed_bit_vector const& b) const {
+    if (&a == &b) return true;
+    unsigned n = num_words();
+    if (n == 0)
+        return true;
+    for (unsigned i = 0; i < n - 1; i++) {
+        if (a.m_data[i] != b.m_data[i])
+            return false;
+    }
+    return last_word(a) == last_word(b);
+}
+unsigned fixed_bit_vector_manager::hash(fixed_bit_vector const& src) const {
+    return string_hash(reinterpret_cast<char const* const>(src.m_data), num_bits()/8, num_bits());
+}
+
+bool fixed_bit_vector_manager::contains(fixed_bit_vector const& a, fixed_bit_vector const& b) const {
+    unsigned n = num_words();
+    if (n == 0)
+        return true;
+    
+    for (unsigned i = 0; i < n - 1; ++i) {
+        if ((a.m_data[i] & b.m_data[i]) != b.m_data[i])
+            return false;
+    }
+    unsigned b_data = last_word(b);
+    return (last_word(a) & b_data) == b_data;
+}
+
+std::ostream& fixed_bit_vector_manager::display(std::ostream& out, fixed_bit_vector const& b) const {
+    unsigned i = num_bits();
+    while (i > 0) {
+        --i;
+        if (b.get(i))
+            out << "1";
+        else
+            out << "0";
+    }
+    return out;
+}
+
+
+

src/util/fixed_bit_vector.h

@@ -0,0 +1,146 @@
+/*++
+Copyright (c) 2014 Microsoft Corporation
+
+Module Name:
+
+    fixed_bit_vector.h
+
+Abstract:
+
+    Simple bitvector implementation for fixed size bit-vectors.
+
+Author:
+
+    Nikolaj Bjorner (nbjorner) 2014-9-15.
+
+Revision History:
+
+    Related to bit_vector, but is based on a manager.
+
+--*/
+#ifndef _FIXED_BIT_VECTOR_H_
+#define _FIXED_BIT_VECTOR_H_
+
+#include<string.h>
+#include"debug.h"
+#include"small_object_allocator.h"
+
+class fixed_bit_vector {
+    friend class fixed_bit_vector_manager;
+    friend class tbv_manager;
+    unsigned    m_data[1];
+
+    static unsigned get_pos_mask(unsigned bit_idx) {
+        return 1 << (bit_idx % 32);
+    }
+
+
+    unsigned get_bit_word(unsigned bit_idx) const {
+        return m_data[bit_idx / 32];
+    }
+
+    unsigned & get_bit_word(unsigned bit_idx) {
+        return m_data[bit_idx / 32];
+    }
+
+public:
+
+    fixed_bit_vector() {}
+    
+    ~fixed_bit_vector() {}
+    
+    unsigned get_word(unsigned word_idx) const { return m_data[word_idx]; }
+
+    bool operator[](unsigned bit_idx) const {
+        return get(bit_idx);
+    }
+    
+    bool get(unsigned bit_idx) const {
+        return (get_bit_word(bit_idx) & get_pos_mask(bit_idx)) != 0;
+    }
+    
+private:
+    void set(unsigned bit_idx) {
+        get_bit_word(bit_idx) |= get_pos_mask(bit_idx);
+    }
+
+    void unset(unsigned bit_idx) {
+        get_bit_word(bit_idx) &= ~get_pos_mask(bit_idx);
+    }
+    
+    void set(unsigned bit_idx, bool val) {
+        int _val = static_cast<int>(val);
+        get_bit_word(bit_idx) ^= (-_val ^ get_bit_word(bit_idx)) & get_pos_mask(bit_idx);
+    }
+
+    // assign bits this[lo:hi] := other[0:hi-lo+1]
+    void set(fixed_bit_vector const& other, unsigned hi, unsigned lo);
+
+};
+
+class fixed_bit_vector_manager {
+    friend class fixed_bit_vector;
+    small_object_allocator m_alloc;
+    unsigned               m_num_bits;
+    unsigned               m_num_bytes;
+    unsigned               m_num_words;
+    unsigned               m_mask;
+    fixed_bit_vector       m_0;
+
+    static unsigned num_words(unsigned num_bits) { 
+        return (num_bits + 31) / 32;
+    }    
+
+public:
+    fixed_bit_vector_manager(unsigned num_bits);
+
+    void reset() { m_alloc.reset(); }
+    fixed_bit_vector* allocate();
+    fixed_bit_vector* allocate1();
+    fixed_bit_vector* allocate0();
+    fixed_bit_vector* allocate(fixed_bit_vector const& bv);
+    void deallocate(fixed_bit_vector* bv);
+
+    void copy(fixed_bit_vector& dst, fixed_bit_vector const& src) const;
+    unsigned num_words() const { return m_num_words; }
+    unsigned num_bytes() const { return m_num_bytes; }
+    unsigned num_bits() const { return m_num_bits; }
+    fixed_bit_vector& reset(fixed_bit_vector& bv) const { return fill0(bv); }
+    fixed_bit_vector& fill0(fixed_bit_vector& bv) const;
+    fixed_bit_vector& fill1(fixed_bit_vector& bv) const;
+    fixed_bit_vector& set_and(fixed_bit_vector& dst, fixed_bit_vector const& src) const;
+    fixed_bit_vector& set_or(fixed_bit_vector& dst,  fixed_bit_vector const& src) const;
+    fixed_bit_vector& set_neg(fixed_bit_vector& dst) const;
+    unsigned last_word(fixed_bit_vector const& bv) const;
+    unsigned get_mask() const { return m_mask; }
+    bool equals(fixed_bit_vector const& a, fixed_bit_vector const& b) const;
+    unsigned hash(fixed_bit_vector const& src) const;
+    bool contains(fixed_bit_vector const& a, fixed_bit_vector const& b) const;
+    std::ostream& display(std::ostream& out, fixed_bit_vector const& b) const;    
+    void set(fixed_bit_vector& dst, unsigned bit_idx) {
+        SASSERT(bit_idx < num_bits());
+        dst.set(bit_idx);
+    }
+    void unset(fixed_bit_vector& dst, unsigned bit_idx) {
+        SASSERT(bit_idx < num_bits());
+        dst.unset(bit_idx);
+    }
+    
+    void set(fixed_bit_vector& dst, unsigned bit_idx, bool val) {
+        SASSERT(bit_idx < num_bits());
+        dst.set(bit_idx, val);
+    }
+
+    // assign bits this[lo:hi] := other[0:hi-lo+1]
+    void set(fixed_bit_vector& dst, fixed_bit_vector const& other, unsigned hi, unsigned lo) {
+        SASSERT(lo <= hi && hi < num_bits());
+        dst.set(other, hi, lo);
+    }
+
+};
+
+
+
+
+#endif /* _FIXED_BIT_VECTOR_H_ */
+

src/util/gparams.cpp

@@ -266,6 +266,7 @@ public:
         }
     }
 
+
     void set(param_descrs const & d, symbol const & param_name, char const * value, symbol const & mod_name) {
         param_kind k = d.get_kind(param_name);
         params_ref & ps = get_params(mod_name);

src/util/hash.cpp

@@ -19,11 +19,18 @@ Revision History:
 
 #include"debug.h"
 #include"hash.h"
+#include <string.h>
+
+static unsigned read_unsigned(const char *s) {
+  unsigned n;
+  memcpy(&n, s, sizeof(unsigned));
+  return n;
+}
 
 // I'm using Bob Jenkin's hash function.
 // http://burtleburtle.net/bob/hash/doobs.html
 unsigned string_hash(const char * str, unsigned length, unsigned init_value) {
-    register unsigned a, b, c, len;
+    unsigned a, b, c, len;
 
     /* Set up the internal state */
     len = length;
@@ -31,10 +38,11 @@ unsigned string_hash(const char * str, unsigned length, unsigned init_value) {
     c = init_value;      /* the previous hash value */
 
     /*---------------------------------------- handle most of the key */
+    SASSERT(sizeof(unsigned) == 4);
     while (len >= 12) {
-        a += reinterpret_cast<const unsigned *>(str)[0];
-        b += reinterpret_cast<const unsigned *>(str)[1];
-        c += reinterpret_cast<const unsigned *>(str)[2];
+        a += read_unsigned(str);
+        b += read_unsigned(str+4);
+        c += read_unsigned(str+8);
         mix(a,b,c);
         str += 12; len -= 12;
     }

src/util/hwf.cpp

@@ -16,7 +16,9 @@ Author:
 Revision History:
 
 --*/
+#include<math.h>
 #include<float.h>
+#include<sstream>
 
 #ifdef _WINDOWS
 #pragma float_control( except, on )   // exception semantics; this does _not_ mean that exceptions are enabled (we want them off!)
@@ -24,19 +26,13 @@ Revision History:
 #pragma fp_contract(off)              // contractions off (`contraction' means x*y+z is turned into a fused-mul-add).
 #pragma fenv_access(on)               // fpu environment sensitivity (needed to be allowed to make FPU mode changes).
 #else
-#ifdef _WINDOWS
-#pragma STDC FENV_ACCESS ON
-#endif
-#include <math.h>
-#include <fenv.h>
+#include<fenv.h>
 #endif
 
 #ifndef _M_IA64
 #define USE_INTRINSICS
 #endif
 
-#include<sstream>
-
 #include"hwf.h"
 
 // Note:
@@ -57,7 +53,7 @@ Revision History:
 
 hwf_manager::hwf_manager() :
     m_mpz_manager(m_mpq_manager)
-{
+{    
 #ifdef _WINDOWS
 #if defined(_AMD64_) || defined(_M_IA64)
     // Precision control is not supported on x64. 
@@ -306,20 +302,31 @@ void hwf_manager::div(mpf_rounding_mode rm, hwf const & x, hwf const & y, hwf &
 #pragma fp_contract(on)
 #endif
 
-void hwf_manager::fused_mul_add(mpf_rounding_mode rm, hwf const & x, hwf const & y, hwf const &z, hwf & o) {
+void hwf_manager::fma(mpf_rounding_mode rm, hwf const & x, hwf const & y, hwf const &z, hwf & o) {
     // CMW: fused_mul_add is not available on most CPUs. As of 2012, only Itanium, 
     // Intel Sandybridge and AMD Bulldozers support that (via AVX).
 
-#ifdef _M_IA64
-    // IA64 (Itanium) will do it, if contractions are on.
     set_rounding_mode(rm);
+
+#ifdef _M_IA64
+    // IA64 (Itanium) will do it, if contractions are on.    
     o.value = x.value * y.value + z.value;
 #else
-    // NOT_IMPLEMENTED_YET();
-    // Just a dummy for now:
-    hwf t;
-    mul(rm, x, y, t);
-    add(rm, t, z, o);
+#if defined(_WINDOWS)    
+#if _MSC_VER >= 1800
+    o.value = ::fma(x.value, y.value, z.value);    
+#else // Windows, older than VS 2013
+  #ifdef USE_INTRINSICS
+      _mm_store_sd(&o.value, _mm_fmadd_sd(_mm_set_sd(x.value), _mm_set_sd(y.value), _mm_set_sd(z.value)));
+  #else
+      // If all else fails, we are imprecise.
+      o.value = (x.value * y.value) + z;
+  #endif
+#endif
+#else
+    // Linux, OSX
+    o.value = ::fma(x.value, y.value, z.value);
+#endif
 #endif
 }
 
@@ -338,8 +345,38 @@ void hwf_manager::sqrt(mpf_rounding_mode rm, hwf const & x, hwf & o) {
 
 void hwf_manager::round_to_integral(mpf_rounding_mode rm, hwf const & x, hwf & o) {
     set_rounding_mode(rm);
-    modf(x.value, &o.value);
-    // Note: on x64, this sometimes produces an SNAN instead of a QNAN?
+    // CMW: modf is not the right function here.
+    // modf(x.value, &o.value);
+
+    // According to the Intel Architecture manual, the x87-instrunction FRNDINT is the 
+    // same in 32-bit and 64-bit mode. The _mm_round_* intrinsics are SSE4 extensions.
+#ifdef _WINDOWS
+#ifdef USE_INTRINSICS
+    switch (rm) {
+    case 0: _mm_store_sd(&o.value, _mm_round_pd(_mm_set_sd(x.value), _MM_FROUND_TO_NEAREST_INT)); break;
+    case 2: _mm_store_sd(&o.value, _mm_round_pd(_mm_set_sd(x.value), _MM_FROUND_TO_POS_INF)); break;
+    case 3: _mm_store_sd(&o.value, _mm_round_pd(_mm_set_sd(x.value), _MM_FROUND_TO_NEG_INF)); break;
+    case 4: _mm_store_sd(&o.value, _mm_round_pd(_mm_set_sd(x.value), _MM_FROUND_TO_ZERO)); break;
+    case 1:
+        UNREACHABLE(); // Note: MPF_ROUND_NEAREST_TAWAY is not supported by the hardware!
+        break;
+    default:
+        UNREACHABLE(); // Unknown rounding mode.
+    }
+#else
+    double xv = x.value;
+    double & ov = o.value;
+
+    __asm {
+        fld     xv
+        frndint
+        fstp    ov // Store result away.
+    }
+#endif
+#else
+    // Linux, OSX.
+    o.value = nearbyint(x.value);
+#endif
 }
 
 void hwf_manager::rem(hwf const & x, hwf const & y, hwf & o) {

src/util/hwf.h

@@ -105,7 +105,7 @@ public:
     void mul(mpf_rounding_mode rm, hwf const & x, hwf const & y, hwf & o);
     void div(mpf_rounding_mode rm, hwf const & x, hwf const & y, hwf & o);    
 
-    void fused_mul_add(mpf_rounding_mode rm, hwf const & x, hwf const & y, hwf const &z, hwf & o);
+    void fma(mpf_rounding_mode rm, hwf const & x, hwf const & y, hwf const &z, hwf & o);
 
     void sqrt(mpf_rounding_mode rm, hwf const & x, hwf & o);
 

src/util/inf_eps_rational.h

@@ -127,6 +127,10 @@ class inf_eps_rational {
         return m_r.get_uint64();
     }
 
+    Numeral const& get_numeral() const {
+        return m_r;
+    }
+
     rational const& get_rational() const {
         return m_r.get_rational();
     }
@@ -139,16 +143,37 @@ class inf_eps_rational {
         return m_infty;
     }
 
+    bool is_finite() const {
+        return m_infty.is_zero();
+    }
+
+    static inf_eps_rational zero() {
+        return inf_eps_rational(Numeral::zero());
+    }
+
+    static inf_eps_rational one() {
+        return inf_eps_rational(Numeral::one());
+    }
+
+    static inf_eps_rational minus_one() {
+        return inf_eps_rational(Numeral::minus_one());
+    }
+
+    static inf_eps_rational infinity() {
+        return inf_eps_rational(rational::one(), Numeral::zero());
+    }
+
+
     inf_eps_rational & operator=(const inf_eps_rational & r) {
         m_infty = r.m_infty;
         m_r = r.m_r;
 	return *this;
     }
 
-    inf_eps_rational & operator=(const rational & r) {
+    inf_eps_rational & operator=(const Numeral & r) {
         m_infty.reset();
         m_r = r;
-        return *this;
+	return *this;
     }
 
     inf_eps_rational & operator+=(const inf_eps_rational & r) { 
@@ -163,6 +188,16 @@ class inf_eps_rational {
 	return *this; 
     }
 
+    inf_eps_rational & operator-=(const inf_rational & r) { 
+        m_r      -= r;
+	return *this; 
+    }
+
+    inf_eps_rational & operator+=(const inf_rational & r) { 
+        m_r      += r;
+	return *this; 
+    }
+
     inf_eps_rational & operator+=(const rational & r) { 
         m_r  += r;
 	return *this; 
@@ -272,12 +307,12 @@ class inf_eps_rational {
     }
 
     friend inline rational floor(const inf_eps_rational & r) {
-        SASSERT(r.m_infty.is_zero());
+        // SASSERT(r.m_infty.is_zero());
         return floor(r.m_r);
     }
 
     friend inline rational ceil(const inf_eps_rational & r) {
-        SASSERT(r.m_infty.is_zero());
+        // SASSERT(r.m_infty.is_zero());
         return ceil(r.m_r);
     }
 

src/util/inf_int_rational.cpp

@@ -19,9 +19,9 @@ Revision History:
 #include<sstream>
 #include"inf_int_rational.h"
 
-inf_int_rational inf_int_rational::m_zero(0);
-inf_int_rational inf_int_rational::m_one(1);
-inf_int_rational inf_int_rational::m_minus_one(-1);
+inf_int_rational inf_int_rational::m_zero;
+inf_int_rational inf_int_rational::m_one;
+inf_int_rational inf_int_rational::m_minus_one;
 
 std::string inf_int_rational::to_string() const {
     if (m_second == 0) {
@@ -39,3 +39,22 @@ std::string inf_int_rational::to_string() const {
     return s.str();
 }
 
+void initialize_inf_int_rational() {
+  inf_int_rational::init();
+}
+
+void inf_int_rational::init() {
+  m_zero.m_first = rational::zero();
+  m_one.m_first = rational::one();
+  m_minus_one.m_first = rational::minus_one();
+}
+
+void finalize_inf_int_rational() {
+  inf_int_rational::finalize();
+}
+
+void inf_int_rational::finalize() {
+  m_zero.~inf_int_rational();
+  m_one.~inf_int_rational();
+  m_minus_one.~inf_int_rational();
+}

src/util/inf_int_rational.h

@@ -33,6 +33,8 @@ class inf_int_rational {
     rational m_first;
     int      m_second;
  public:
+    static void init(); // called from rational::initialize() only
+    static void finalize();  // called from rational::finalize() only
 
     unsigned hash() const { 
         return m_first.hash() ^ (static_cast<unsigned>(m_second) + 1);
@@ -155,6 +157,17 @@ class inf_int_rational {
 	return *this; 
     }
 
+    inf_int_rational & operator*=(const rational & r) { 
+        if (!r.is_int32()) {
+            throw default_exception("multiplication with large rational is not possible");
+        }
+        m_first  *= r;
+        m_second *= r.get_int32();
+	return *this; 
+    }
+
+
+
     inf_int_rational & operator-=(const inf_int_rational & r) { 
         m_first  -= r.m_first;
         m_second -= r.m_second;
@@ -261,7 +274,7 @@ class inf_int_rational {
             if (r.m_second >= 0) {
                 return r.m_first;
             }
-            return r.m_first - rational(1);
+            return r.m_first - rational::one();
         }
         
         return floor(r.m_first);
@@ -272,7 +285,7 @@ class inf_int_rational {
             if (r.m_second <= 0) {
                 return r.m_first;
             }
-            return r.m_first + rational(1);
+            return r.m_first + rational::one();
         }
         
         return ceil(r.m_first);
@@ -344,6 +357,10 @@ inline inf_int_rational operator+(const inf_int_rational & r1, const inf_int_rat
     return inf_int_rational(r1) += r2; 
 }
 
+inline inf_int_rational operator*(const rational & r1, const inf_int_rational & r2) { 
+    return inf_int_rational(r2) *= r1;
+}
+
 inline inf_int_rational operator-(const inf_int_rational & r1, const inf_int_rational & r2) { 
     return inf_int_rational(r1) -= r2; 
 }

src/util/inf_rational.cpp

@@ -18,9 +18,9 @@ Revision History:
 --*/
 #include"inf_rational.h"
 
-inf_rational inf_rational::m_zero(0);
-inf_rational inf_rational::m_one(1);
-inf_rational inf_rational::m_minus_one(-1);
+inf_rational inf_rational::m_zero;
+inf_rational inf_rational::m_one;
+inf_rational inf_rational::m_minus_one;
 
 inf_rational inf_mult(inf_rational const& r1, inf_rational const& r2) 
 {
@@ -128,7 +128,7 @@ inf_rational inf_power(inf_rational const& r, unsigned n)
         // 0 will work.
     }
     else if (r.m_first.is_zero()) {
-        result.m_first = rational(-1);        
+        result.m_first = rational::minus_one();
     }
     else if (r.m_first.is_pos()) {
         result.m_first = rational(r.m_first - r.m_first/rational(2)).expt(n);
@@ -152,7 +152,7 @@ inf_rational sup_power(inf_rational const& r, unsigned n)
         result.m_first = r.m_first.expt(n);
     }
     else if (r.m_first.is_zero() || (n == 0)) {
-        result.m_first = rational(1);        
+        result.m_first = rational::one();
     }
     else if (r.m_first.is_pos() || is_even) {
         result.m_first = rational(r.m_first + r.m_first/rational(2)).expt(n);
@@ -177,3 +177,23 @@ inf_rational sup_root(inf_rational const& r, unsigned n)
     // use r.
     return r;
 }
+
+void initialize_inf_rational() {
+    inf_rational::init();
+}
+
+void inf_rational::init() {
+    m_zero.m_first = rational::zero();
+    m_one.m_first = rational::one();
+    m_minus_one.m_first = rational::minus_one();
+}
+
+void finalize_inf_rational() {
+    inf_rational::finalize();
+}
+
+void inf_rational::finalize() {
+    m_zero.~inf_rational();
+    m_one.~inf_rational();
+    m_minus_one.~inf_rational();
+}

src/util/inf_rational.h

@@ -33,6 +33,8 @@ class inf_rational {
     rational m_first;
     rational m_second;
  public:
+    static void init(); // called from rational::initialize() only
+    static void finalize();  // called from rational::finalize() only
 
     unsigned hash() const { 
         return m_first.hash() ^ (m_second.hash()+1);
@@ -40,7 +42,7 @@ class inf_rational {
 
     struct hash_proc {  unsigned operator()(inf_rational const& r) const { return r.hash(); }  };
 
-	struct eq_proc { bool operator()(inf_rational const& r1, inf_rational const& r2) const { return r1 == r2; } };
+    struct eq_proc { bool operator()(inf_rational const& r1, inf_rational const& r2) const { return r1 == r2; } };
 
     void swap(inf_rational & n) { 
         m_first.swap(n.m_first);
@@ -63,10 +65,7 @@ class inf_rational {
         return s;
     }
 
-    inf_rational():
-        m_first(rational()),
-        m_second(rational())
-     {}
+    inf_rational() {}
 
     inf_rational(const inf_rational & r): 
         m_first(r.m_first),
@@ -85,10 +84,10 @@ class inf_rational {
 
     explicit inf_rational(rational const& r, bool pos_inf):
         m_first(r),
-        m_second(pos_inf?rational(1):rational(-1))
+        m_second(pos_inf ? rational::one() : rational::minus_one())
     {}
 
-    explicit inf_rational(rational const& r):
+    inf_rational(rational const& r):
         m_first(r)
     {
         m_second.reset();
@@ -316,7 +315,7 @@ class inf_rational {
             if (r.m_second.is_nonneg()) {
                 return r.m_first;
             }
-            return r.m_first - rational(1);
+            return r.m_first - rational::one();
         }
         
         return floor(r.m_first);
@@ -327,7 +326,7 @@ class inf_rational {
             if (r.m_second.is_nonpos()) {
                 return r.m_first;
             }
-            return r.m_first + rational(1);
+            return r.m_first + rational::one();
         }
         
         return ceil(r.m_first);

src/util/map.h

@@ -108,6 +108,10 @@ public:
         m_table.insert(key_data(k, v));
     }
   
+    bool insert_if_not_there_core(key const & k, value const & v, entry *& et) {
+        return m_table.insert_if_not_there_core(key_data(k,v), et);
+    }
+
     key_data const & insert_if_not_there(key const & k, value const & v) {
         return m_table.insert_if_not_there(key_data(k, v));
     }

src/util/memory_manager.cpp

@@ -1,3 +1,9 @@
+
+/*++
+Copyright (c) 2015 Microsoft Corporation
+
+--*/
+
 #include<iostream>
 #include<stdlib.h>
 #include<limits.h>

src/util/mpf.cpp

@@ -210,7 +210,7 @@ void mpf_manager::set(mpf & o, unsigned ebits, unsigned sbits, mpf_rounding_mode
 
     std::string f, e;
 
-	f = (e_pos != std::string::npos) ? v.substr(0, e_pos) : v;
+    f = (e_pos != std::string::npos) ? v.substr(0, e_pos) : v;
     e = (e_pos != std::string::npos) ? v.substr(e_pos+1) : "0";    
 
     TRACE("mpf_dbg", tout << " f = " << f << " e = " << e << std::endl;);   
@@ -232,18 +232,18 @@ void mpf_manager::set(mpf & o, unsigned ebits, unsigned sbits, mpf_rounding_mode
 
     o.ebits = ebits;
     o.sbits = sbits;
-    o.sign = m_mpq_manager.is_neg(significand);	
+    o.sign = m_mpq_manager.is_neg(significand);    
     
     if (m_mpq_manager.is_zero(significand))
         mk_zero(ebits, sbits, o.sign, o);
-    else {	
+    else {    
         scoped_mpq sig(m_mpq_manager);
         scoped_mpz exp(m_mpq_manager);
 
         m_mpq_manager.set(sig, significand);
         m_mpq_manager.abs(sig);
-        m_mpz_manager.set(exp, exponent);	
-	
+        m_mpz_manager.set(exp, exponent);    
+    
         // Normalize
         while (m_mpq_manager.ge(sig, 2)) {
             m_mpq_manager.div(sig, mpq(2), sig);
@@ -277,7 +277,7 @@ void mpf_manager::set(mpf & o, unsigned ebits, unsigned sbits, mpf_rounding_mode
 
         TRACE("mpf_dbg", tout << "sig = " << m_mpz_manager.to_string(o.significand) << 
                                  " exp = " << o.exponent << std::endl;);
-	
+    
         if (m_mpz_manager.is_small(exp)) {
             o.exponent = m_mpz_manager.get_int64(exp);
             round(rm, o);
@@ -348,7 +348,7 @@ void mpf_manager::set(mpf & o, unsigned ebits, unsigned sbits, mpf_rounding_mode
             bool sticky = false;
             while (ds < 0)
             {
-                if (!m_mpz_manager.is_even(o.significand)) sticky = true;
+                sticky |= m_mpz_manager.is_odd(o.significand);
                 m_mpz_manager.machine_div2k(o.significand, 1);            
                 ds++;
             }
@@ -1003,9 +1003,38 @@ void mpf_manager::round_to_integral(mpf_rounding_mode rm, mpf const & x, mpf & o
         mk_nan(x.ebits, x.sbits, o);
     else if (is_inf(x))
         set(o, x);
-    else if (x.exponent < 0) 
-        mk_zero(x.ebits, x.sbits, x.sign, o);
-    else if (x.exponent >= x.sbits-1)
+    else if (is_zero(x))
+        mk_zero(x.ebits, x.sbits, x.sign, o); // -0.0 -> -0.0, says IEEE754, Sec 5.9.
+    else if (x.exponent < 0) {
+        if (rm == MPF_ROUND_TOWARD_ZERO)
+            mk_zero(x.ebits, x.sbits, x.sign, o);
+        else if (rm == MPF_ROUND_TOWARD_NEGATIVE) {
+            if (x.sign)
+                mk_one(x.ebits, x.sbits, true, o);
+            else
+                mk_zero(x.ebits, x.sbits, false, o);
+        }
+        else if (rm == MPF_ROUND_TOWARD_POSITIVE) {
+            if (x.sign)
+                mk_zero(x.ebits, x.sbits, true, o);
+            else
+                mk_one(x.ebits, x.sbits, false, o);
+        }
+        else {
+            SASSERT(rm == MPF_ROUND_NEAREST_TEVEN || rm == MPF_ROUND_NEAREST_TAWAY);
+            bool tie = m_mpz_manager.is_zero(x.significand) && x.exponent == -1;
+            TRACE("mpf_dbg", tout << "tie = " << tie << std::endl;);
+            if (tie && rm == MPF_ROUND_NEAREST_TEVEN)
+                mk_zero(x.ebits, x.sbits, x.sign, o);
+            else if (tie && rm == MPF_ROUND_NEAREST_TAWAY)
+                mk_one(x.ebits, x.sbits, x.sign, o);
+            else if (x.exponent < -1)
+                mk_zero(x.ebits, x.sbits, x.sign, o);
+            else
+                mk_one(x.ebits, x.sbits, x.sign, o);
+        }
+    }
+    else if (x.exponent >= x.sbits - 1)
         set(o, x);
     else {
         SASSERT(x.exponent >= 0 && x.exponent < x.sbits-1);
@@ -1016,21 +1045,74 @@ void mpf_manager::round_to_integral(mpf_rounding_mode rm, mpf const & x, mpf & o
 
         scoped_mpf a(*this);
         set(a, x);
-        unpack(a, true);
+        unpack(a, true); // A includes hidden bit
+
+        TRACE("mpf_dbg", tout << "A = " << to_string_raw(a) << std::endl;);
+
+        SASSERT(m_mpz_manager.lt(a.significand(), m_powers2(x.sbits)));
+        SASSERT(m_mpz_manager.ge(a.significand(), m_powers2(x.sbits - 1)));
 
-        TRACE("mpf_dbg", tout << "A = " << to_string(a) << std::endl;);
-                
         o.exponent = a.exponent();
         m_mpz_manager.set(o.significand, a.significand());
+        
+        unsigned shift = (o.sbits - 1) - ((unsigned)o.exponent);
+        const mpz & shift_p = m_powers2(shift);
+        TRACE("mpf_dbg", tout << "shift=" << shift << std::endl;);
+        
+        scoped_mpz div(m_mpz_manager), rem(m_mpz_manager);
+        m_mpz_manager.machine_div_rem(o.significand, shift_p, div, rem);
+        TRACE("mpf_dbg", tout << "div=" << m_mpz_manager.to_string(div) << " rem=" << m_mpz_manager.to_string(rem) << std::endl;);                
 
-        unsigned q = (unsigned) o.exponent;
-        unsigned shift = o.sbits-q-1;
-        TRACE("mpf_dbg", tout << "Q = " << q << " shift=" << shift << std::endl;);
-        m_mpz_manager.machine_div2k(o.significand, shift);
-        m_mpz_manager.mul2k(o.significand, shift+3);
+        const mpz & shift_p1 = m_powers2(shift-1);
+        TRACE("mpf_dbg", tout << "shift_p1=" << m_mpz_manager.to_string(shift_p1) << std::endl;);
 
-        round(rm, o);
-    }    
+        switch (rm) {
+        case MPF_ROUND_NEAREST_TEVEN:
+        case MPF_ROUND_NEAREST_TAWAY: {            
+            bool tie = m_mpz_manager.eq(rem, shift_p1);
+            bool less_than_tie = m_mpz_manager.lt(rem, shift_p1);
+            bool more_than_tie = m_mpz_manager.gt(rem, shift_p1);
+            TRACE("mpf_dbg", tout << "tie= " << tie << "; <tie = " << less_than_tie << "; >tie = " << more_than_tie << std::endl;);
+            if (tie) {
+                if ((rm == MPF_ROUND_NEAREST_TEVEN && m_mpz_manager.is_odd(div)) ||
+                    (rm == MPF_ROUND_NEAREST_TAWAY && m_mpz_manager.is_even(div))) {
+                    TRACE("mpf_dbg", tout << "div++ (1)" << std::endl;);
+                    m_mpz_manager.inc(div);
+                }
+            }
+            else {
+                SASSERT(less_than_tie || more_than_tie);                
+                if (more_than_tie) {
+                    m_mpz_manager.inc(div);
+                    TRACE("mpf_dbg", tout << "div++ (2)" << std::endl;);
+                }
+            }
+            break;
+        }
+        case MPF_ROUND_TOWARD_POSITIVE:
+            if (!m_mpz_manager.is_zero(rem) && !o.sign)
+                m_mpz_manager.inc(div);
+            break;
+        case MPF_ROUND_TOWARD_NEGATIVE:
+            if (!m_mpz_manager.is_zero(rem) && o.sign)
+                m_mpz_manager.inc(div);
+            break;
+        case MPF_ROUND_TOWARD_ZERO:
+        default:
+            /* nothing */;
+        }
+           
+        m_mpz_manager.mul2k(div, shift, o.significand);
+        SASSERT(m_mpz_manager.ge(o.significand, m_powers2(o.sbits - 1)));
+
+        // re-normalize
+        while (m_mpz_manager.ge(o.significand, m_powers2(o.sbits))) {
+            m_mpz_manager.machine_div2k(o.significand, 1);
+            o.exponent++;
+        }
+
+        m_mpz_manager.sub(o.significand, m_powers2(o.sbits - 1), o.significand); // strip hidden bit
+    }
 
     TRACE("mpf_dbg", tout << "INTEGRAL = " << to_string(o) << std::endl;);
 }
@@ -1099,10 +1181,10 @@ void mpf_manager::rem(mpf const & x, mpf const & y, mpf & o) {
         mk_nan(x.ebits, x.sbits, o);
     else if (is_inf(y))
         set(o, x);
-    else if (is_zero(x))        
-        set(o, x);
     else if (is_zero(y))
         mk_nan(x.ebits, x.sbits, o);
+    else if (is_zero(x))        
+        set(o, x);    
     else {        
         o.ebits = x.ebits;
         o.sbits = x.sbits;
@@ -1147,7 +1229,11 @@ void mpf_manager::rem(mpf const & x, mpf const & y, mpf & o) {
 }
 
 void mpf_manager::maximum(mpf const & x, mpf const & y, mpf & o) {
-    if (is_nan(x) || (is_zero(x) && is_zero(y)))
+    if (is_nan(x))
+        set(o, y);
+    else if (is_zero(x) && is_zero(y) && sgn(x) != sgn(y))
+        mk_pzero(x.ebits, x.sbits, o);
+    else if (is_zero(x) && is_zero(y))
         set(o, y);
     else if (is_nan(y))
         set(o, x);
@@ -1158,7 +1244,11 @@ void mpf_manager::maximum(mpf const & x, mpf const & y, mpf & o) {
 }
 
 void mpf_manager::minimum(mpf const & x, mpf const & y, mpf & o) {
-    if (is_nan(x) || (is_zero(x) && is_zero(y)))
+    if (is_nan(x))
+        set(o, y);
+    else if (is_zero(x) && is_zero(y) && sgn(x) != sgn(y))
+        mk_pzero(x.ebits, x.sbits, o);
+    else if (is_zero(x) && is_zero(y))
         set(o, y);
     else if (is_nan(y))
         set(o, x);
@@ -1350,7 +1440,7 @@ bool mpf_manager::is_ninf(mpf const & x) {
 }
 
 bool mpf_manager::is_normal(mpf const & x) {
-    return !(has_top_exp(x) || is_denormal(x));
+    return !(has_top_exp(x) || is_denormal(x) || is_zero(x));
 }
 
 bool mpf_manager::is_denormal(mpf const & x) {
@@ -1449,6 +1539,14 @@ void mpf_manager::mk_nan(unsigned ebits, unsigned sbits, mpf & o) {
     o.sign = false;
 }
 
+void mpf_manager::mk_one(unsigned ebits, unsigned sbits, bool sign, mpf & o) const {
+    o.sbits = sbits;
+    o.ebits = ebits;
+    o.sign = sign;
+    m_mpz_manager.set(o.significand, 0);
+    o.exponent = 0;
+}
+
 void mpf_manager::mk_max_value(unsigned ebits, unsigned sbits, bool sign, mpf & o) {
     o.sbits = sbits;
     o.ebits = ebits;
@@ -1475,6 +1573,9 @@ void mpf_manager::mk_ninf(unsigned ebits, unsigned sbits, mpf & o) {
 
 void mpf_manager::unpack(mpf & o, bool normalize) {
     // Insert the hidden bit or adjust the exponent of denormal numbers.    
+    if (is_zero(o))
+        return;
+
     if (is_normal(o))
         m_mpz_manager.add(o.significand, m_powers2(o.sbits-1), o.significand);
     else {

src/util/mpf.h

@@ -46,7 +46,7 @@ class mpf {
     mpz significand;
     mpf_exp_t exponent;
     mpf & operator=(mpf const & other) { UNREACHABLE(); return *this; }
-    void set(unsigned ebits, unsigned sbits);
+    void set(unsigned _ebits, unsigned _sbits);
 public:    
     mpf();
     mpf(unsigned ebits, unsigned sbits);
@@ -208,7 +208,9 @@ public:
 
     void to_sbv_mpq(mpf_rounding_mode rm, const mpf & x, scoped_mpq & o);
 
-protected:    
+protected:
+    void mk_one(unsigned ebits, unsigned sbits, bool sign, mpf & o) const;
+
     bool has_bot_exp(mpf const & x);
     bool has_top_exp(mpf const & x);
 

src/util/mpff.cpp

@@ -255,7 +255,7 @@ void mpff_manager::set(mpff & n, int64 v) {
     }
     else {
         if (v < 0) {
-            set(n, static_cast<uint64>(-v));
+            set(n, 1 + static_cast<uint64>(-(1+v)));
             n.m_sign = 1;
         }
         else {
@@ -680,7 +680,7 @@ void mpff_manager::add_sub(bool is_sub, mpff const & a, mpff const & b, mpff & c
     
     // Make sure that a and b have the same exponent.
     if (exp_a > exp_b) {
-        unsigned shift = exp_a - exp_b;
+        unsigned shift = (unsigned)exp_a - (unsigned)exp_b;
         n_sig_b = m_buffers[0].c_ptr();
         shr(m_precision, sig_b, shift, m_precision, n_sig_b);
         if (sgn_b != m_to_plus_inf && has_one_at_first_k_bits(m_precision, sig_b, shift)) {

src/util/mpn.h

@@ -27,7 +27,9 @@ Revision History:
 typedef unsigned int mpn_digit;
 
 class mpn_manager {
+#ifndef _NO_OMP_
     omp_nest_lock_t m_lock;
+#endif
 #define MPN_BEGIN_CRITICAL() omp_set_nest_lock(&m_lock);
 #define MPN_END_CRITICAL() omp_unset_nest_lock(&m_lock);
 

src/util/mpq_inf.cpp

@@ -19,7 +19,7 @@ Revision History:
 #include"mpq_inf.h"
 
 template<bool SYNCH>
-std::string mpq_inf_manager<SYNCH>::to_string(mpq_inf const & a) const {
+std::string mpq_inf_manager<SYNCH>::to_string(mpq_inf const & a) {
     if (m.is_zero(a.second))
         return m.to_string(a.first);
     
@@ -38,5 +38,6 @@ std::string mpq_inf_manager<SYNCH>::to_string(mpq_inf const & a) const {
     return s;
 }
 
+
 template class mpq_inf_manager<true>;
 template class mpq_inf_manager<false>;

src/util/mpq_inf.h

@@ -26,12 +26,12 @@ typedef std::pair<mpq, mpq> mpq_inf;
 
 template<bool SYNCH = true>
 class mpq_inf_manager {
-    mpq_manager<SYNCH> & m;
+    mpq_manager<SYNCH>   m;
     double               m_inf;
 public:
     typedef mpq_inf numeral;
 
-    mpq_inf_manager(mpq_manager<SYNCH> & _m, double inf = 0.0001):m(_m) {
+    mpq_inf_manager(double inf = 0.0001) {
         set_inf(inf);
     }
 
@@ -83,6 +83,14 @@ public:
     }
 
     bool is_int(mpq_inf const & a) const { return m.is_int(a.first) && m.is_zero(a.second); }
+
+    bool is_pos(mpq_inf const & a) const { 
+        return m.is_pos(a.first) || (m.is_zero(a.first) && m.is_pos(a.second)); 
+    }
+
+    bool is_neg(mpq_inf const & a) const { 
+        return m.is_neg(a.first) || (m.is_zero(a.first) && m.is_neg(a.second)); 
+    }
     
     bool is_rational(mpq_inf const & a) const { return m.is_zero(a.second); }
 
@@ -104,15 +112,15 @@ public:
         return m.is_zero(a.first) && m.is_zero(a.second);
     }
 
-    bool eq(mpq_inf const & a, mpq_inf const & b) const {
+    bool eq(mpq_inf const & a, mpq_inf const & b) {
         return m.eq(a.first, b.first) && m.eq(a.second, b.second);
     }
 
-    bool eq(mpq_inf const & a, mpq const & b) const {
+    bool eq(mpq_inf const & a, mpq const & b) {
         return m.eq(a.first, b) && m.is_zero(a.second);
     }
         
-    bool eq(mpq_inf const & a, mpq const & b, inf_kind k) const {
+    bool eq(mpq_inf const & a, mpq const & b, inf_kind k) {
         if (!m.eq(a.first, b))
             return false;
         switch (k) {
@@ -124,15 +132,15 @@ public:
         return false;
     }
     
-    bool lt(mpq_inf const & a, mpq_inf const & b) const {
+    bool lt(mpq_inf const & a, mpq_inf const & b) {
         return m.lt(a.first, b.first) || (m.lt(a.second, b.second) && m.eq(a.first, b.first));
     }
 
-    bool lt(mpq_inf const & a, mpq const & b) const {
+    bool lt(mpq_inf const & a, mpq const & b) {
         return m.lt(a.first, b) || (m.is_neg(a.second) && m.eq(a.first, b));
     }
 
-    bool lt(mpq_inf const & a, mpq const & b, inf_kind k) const {
+    bool lt(mpq_inf const & a, mpq const & b, inf_kind k) {
         if (m.lt(a.first, b))
             return true;
         if (m.eq(a.first, b)) {
@@ -146,13 +154,13 @@ public:
         return false;
     }
 
-    bool gt(mpq_inf const & a, mpq_inf const & b) const { return lt(b, a); }
+    bool gt(mpq_inf const & a, mpq_inf const & b) { return lt(b, a); }
     
-    bool gt(mpq_inf const & a, mpq const & b) const {
+    bool gt(mpq_inf const & a, mpq const & b) {
         return m.gt(a.first, b) || (m.is_pos(a.second) && m.eq(a.first, b));
     }
 
-    bool gt(mpq_inf const & a, mpq const & b, inf_kind k) const {
+    bool gt(mpq_inf const & a, mpq const & b, inf_kind k) {
         if (m.gt(a.first, b))
             return true;
         if (m.eq(a.first, b)) {
@@ -166,17 +174,17 @@ public:
         return false;
     }
 
-    bool le(mpq_inf const & a, mpq_inf const & b) const { return !gt(a, b); }
+    bool le(mpq_inf const & a, mpq_inf const & b) { return !gt(a, b); }
 
-    bool le(mpq_inf const & a, mpq const & b) const { return !gt(a, b); }
+    bool le(mpq_inf const & a, mpq const & b) { return !gt(a, b); }
 
-    bool le(mpq_inf const & a, mpq const & b, inf_kind k) const { return !gt(a, b, k); }
+    bool le(mpq_inf const & a, mpq const & b, inf_kind k) { return !gt(a, b, k); }
 
-    bool ge(mpq_inf const & a, mpq_inf const & b) const { return !lt(a, b); }
+    bool ge(mpq_inf const & a, mpq_inf const & b) { return !lt(a, b); }
 
-    bool ge(mpq_inf const & a, mpq const & b) const { return !lt(a, b); }
+    bool ge(mpq_inf const & a, mpq const & b) { return !lt(a, b); }
 
-    bool ge(mpq_inf const & a, mpq const & b, inf_kind k) const { return !lt(a, b, k); }
+    bool ge(mpq_inf const & a, mpq const & b, inf_kind k) { return !lt(a, b, k); }
 
     void add(mpq_inf const & a, mpq_inf const & b, mpq_inf & c) {
         m.add(a.first, b.first, c.first);
@@ -208,6 +216,16 @@ public:
         m.mul(b, a.second, c.second);
     }
 
+    void div(mpq_inf const & a, mpq const & b, mpq_inf & c) {
+        m.div(a.first, b, c.first);
+        m.div(a.second, b, c.second);
+    }
+
+    void div(mpq_inf const & a, mpz const & b, mpq_inf & c) {
+        m.div(a.first, b,  c.first);
+        m.div(a.second, b, c.second);
+    }
+
     void inc(mpq_inf & a) {
         m.inc(a.first);
     }
@@ -221,10 +239,16 @@ public:
         m.neg(a.second);
     }
 
+    void abs(mpq_inf & a) {
+        if (is_neg(a)) {
+            neg(a);
+        }
+    }
+
     void ceil(mpq_inf const & a, mpq & b) {
         if (m.is_int(a.first)) {
             // special cases for  k - delta*epsilon where k is an integer
-            if (m.is_pos(a.first))
+            if (m.is_pos(a.second))
                 m.add(a.first, mpq(1), b); // ceil(k + delta*epsilon) --> k+1
             else
                 m.set(b, a.first);
@@ -246,7 +270,13 @@ public:
         }
     }
 
-    std::string to_string(mpq_inf const & a) const;
+    std::string to_string(mpq_inf const & a);
+
+    void display(std::ostream & out, mpq_inf const & a) {
+        out << to_string(a);
+    }
+
+    mpq_manager<SYNCH>& get_mpq_manager() { return m; }
 };
 
 typedef mpq_inf_manager<true>  synch_mpq_inf_manager;

src/util/mpz.cpp

@@ -336,17 +336,17 @@ void mpz_manager<SYNCH>::set(mpz & target, unsigned sz, digit_t const * digits)
             }
         }
 #else
-	mk_big(target);
-	// reset
-	mpz_set_ui(*target.m_ptr, digits[sz - 1]);
-	SASSERT(sz > 0);
-	unsigned i = sz - 1;
-	while (i > 0) {
-	  --i;
-	  mpz_mul_2exp(*target.m_ptr, *target.m_ptr, 32);
-	  mpz_set_ui(m_tmp, digits[i]);
-	  mpz_add(*target.m_ptr, *target.m_ptr, m_tmp);
-	}
+    mk_big(target);
+    // reset
+    mpz_set_ui(*target.m_ptr, digits[sz - 1]);
+    SASSERT(sz > 0);
+    unsigned i = sz - 1;
+    while (i > 0) {
+      --i;
+      mpz_mul_2exp(*target.m_ptr, *target.m_ptr, 32);
+      mpz_set_ui(m_tmp, digits[i]);
+      mpz_add(*target.m_ptr, *target.m_ptr, m_tmp);
+    }
 #endif        
     }
 }
@@ -2037,16 +2037,16 @@ bool mpz_manager<SYNCH>::decompose(mpz const & a, svector<digit_t> & digits) {
         }
         return a.m_val < 0;
 #else
-	bool r = is_neg(a);
-	mpz_set(m_tmp, *a.m_ptr);
-	mpz_abs(m_tmp, m_tmp);
-	while (mpz_sgn(m_tmp) != 0) {
-	  mpz_tdiv_r_2exp(m_tmp2, m_tmp, 32);
-	  unsigned v = mpz_get_ui(m_tmp2);
-	  digits.push_back(v);
-	  mpz_tdiv_q_2exp(m_tmp, m_tmp, 32);
-	}
-	return r;
+    bool r = is_neg(a);
+    mpz_set(m_tmp, *a.m_ptr);
+    mpz_abs(m_tmp, m_tmp);
+    while (mpz_sgn(m_tmp) != 0) {
+      mpz_tdiv_r_2exp(m_tmp2, m_tmp, 32);
+      unsigned v = mpz_get_ui(m_tmp2);
+      digits.push_back(v);
+      mpz_tdiv_q_2exp(m_tmp, m_tmp, 32);
+    }
+    return r;
 #endif
     }
 }

src/util/rational.cpp

@@ -24,14 +24,14 @@ Revision History:
 #endif
 
 synch_mpq_manager *  rational::g_mpq_manager = 0;
-rational             rational::m_zero(0);
-rational             rational::m_one(1);
-rational             rational::m_minus_one(-1);
+rational             rational::m_zero;
+rational             rational::m_one;
+rational             rational::m_minus_one;
 vector<rational>     rational::m_powers_of_two;
 
-void mk_power_up_to(vector<rational> & pws, unsigned n) {
+static void mk_power_up_to(vector<rational> & pws, unsigned n) {
     if (pws.empty()) {
-        pws.push_back(rational(1));
+        pws.push_back(rational::one());
     }
     unsigned sz = pws.size();
     rational curr = pws[sz - 1];
@@ -53,14 +53,32 @@ rational rational::power_of_two(unsigned k) {
     return result;
 }
 
+// in inf_rational.cpp
+void initialize_inf_rational();
+void finalize_inf_rational();
+
+// in inf_int_rational.cpp
+void initialize_inf_int_rational();
+void finalize_inf_int_rational();
+
 void rational::initialize() {
     if (!g_mpq_manager) {
         g_mpq_manager = alloc(synch_mpq_manager);
+        m().set(m_zero.m_val, 0);
+        m().set(m_one.m_val, 1);
+        m().set(m_minus_one.m_val, -1);
+        initialize_inf_rational();
+        initialize_inf_int_rational();
     }
 }
 
 void rational::finalize() {
+    finalize_inf_rational();
+    finalize_inf_int_rational();
     m_powers_of_two.finalize();
+    m_zero.~rational();
+    m_one.~rational();
+    m_minus_one.~rational();
     dealloc(g_mpq_manager);
     g_mpq_manager = 0;
 }

src/util/rational.h

@@ -111,6 +111,11 @@ public:
         return INT_MIN <= v && v <= INT_MAX;
     }
 
+    int get_int32() const {
+        SASSERT(is_int32());
+        return (int)get_int64();
+    }
+
     double get_double() const { return m().get_double(m_val); }
 
     rational const & get_rational() const { return *this; }

src/util/region.h

@@ -53,6 +53,7 @@ public:
         m_scopes.push_back(m_chuncks.size());
     }
 
+
     void pop_scope() {
         unsigned old_size = m_scopes.back();
         m_scopes.pop_back();

src/util/s_integer.h

@@ -67,7 +67,7 @@ public:
     s_integer const& get_s_integer() const { return *this; }
     s_integer const& get_infinitesimal() const { return zero(); }
     static bool is_rational() { return true; }
-    s_integer const& get_rational() const { return *this; }
+    s_integer const& get_rational() const { return *this; } 
     s_integer & operator=(const s_integer & r) { m_val = r.m_val; return *this; }
     friend inline s_integer numerator(const s_integer & r) { return r; }
     friend inline s_integer denominator(const s_integer & r) { return one(); }

src/util/scoped_numeral.h

@@ -39,7 +39,7 @@ public:
     numeral const & get() const { return m_num; }
     numeral & get() { return m_num; }
     
-    _scoped_numeral & operator=(_scoped_numeral & n) {
+    _scoped_numeral & operator=(_scoped_numeral const & n) {
         if (this == &n)
             return *this;
         m().set(m_num, n.m_num);
@@ -97,6 +97,10 @@ public:
         return a.m().eq(a, b);
     }
 
+    friend bool operator!=(_scoped_numeral const & a, numeral const & b) {
+        return !a.m().eq(a, b);
+    }
+
     friend bool operator<(_scoped_numeral const & a, numeral const & b) {
         return a.m().lt(a, b);
     }
@@ -113,6 +117,26 @@ public:
         return a.m().ge(a, b);
     }
 
+    bool is_zero() const {
+        return m().is_zero(*this);
+    }
+
+    bool is_pos() const {
+        return m().is_pos(*this);
+    }
+
+    bool is_neg() const {
+        return m().is_neg(*this);
+    }
+
+    bool is_nonpos() const {
+        return m().is_nonpos(*this);
+    }
+
+    bool is_nonneg() const {
+        return m().is_nonneg(*this);
+    }
+
     friend bool is_zero(_scoped_numeral const & a) {
         return a.m().is_zero(a);
     }
@@ -133,6 +157,12 @@ public:
         return a.m().is_nonpos(a);
     }
 
+    friend _scoped_numeral abs(_scoped_numeral const& a) {
+        _scoped_numeral res(a);
+        a.m().abs(res);
+        return res;
+    }
+
     void neg() {
         m().neg(m_num);
     }

src/util/scoped_numeral_vector.h

@@ -46,6 +46,10 @@ public:
         m_manager.set(this->back(), v);
     }
 
+    void pop_back() {
+        shrink(this->size()-1);
+    }
+
     void shrink(unsigned sz) {
         unsigned old_sz = this->size();
         if (old_sz == sz)

src/util/sorting_network.h

@@ -0,0 +1,767 @@
+/*++
+Copyright (c) 2013 Microsoft Corporation
+
+Module Name:
+
+    sorting_network.h
+
+Abstract:
+
+    Utility for creating a sorting network.
+
+Author:
+
+    Nikolaj Bjorner (nbjorner) 2013-11-07
+
+Notes:
+
+    Same routine is used in Formula.
+
+--*/
+
+#include "vector.h"
+
+#ifndef _SORTING_NETWORK_H_
+#define _SORTING_NETWORK_H_
+
+
+    template <typename Ext>
+    class sorting_network {
+        typedef typename Ext::vector vect;
+        Ext&               m_ext;
+        svector<unsigned>  m_currentv;
+        svector<unsigned>  m_nextv;
+        svector<unsigned>* m_current;
+        svector<unsigned>* m_next;
+
+        unsigned& current(unsigned i) { return (*m_current)[i]; }
+        unsigned& next(unsigned i) { return (*m_next)[i]; }
+       
+        void exchange(unsigned i, unsigned j, vect& out) {
+            SASSERT(i <= j);
+            if (i < j) {
+                typename Ext::T ei = out.get(i);
+                typename Ext::T ej = out.get(j);
+                out.set(i, m_ext.mk_ite(m_ext.mk_le(ei, ej), ei, ej));
+                out.set(j, m_ext.mk_ite(m_ext.mk_le(ej, ei), ei, ej));
+            }
+        }
+        
+        void sort(unsigned k, vect& out) {
+            SASSERT(is_power_of2(k) && k > 0);
+            if (k == 2) {
+                for (unsigned i = 0; i < out.size()/2; ++i) {
+                    exchange(current(2*i), current(2*i+1), out);
+                    next(2*i) = current(2*i);
+                    next(2*i+1) = current(2*i+1);
+                }
+                std::swap(m_current, m_next);
+            }
+            else {
+                
+                for (unsigned i = 0; i < out.size()/k; ++i) {
+                    unsigned ki = k * i;
+                    for (unsigned j = 0; j < k / 2; ++j) {
+                        next(ki + j) = current(ki + (2 * j));
+                        next(ki + (k / 2) + j) = current(ki + (2 * j) + 1);
+                    }
+                }
+                
+                std::swap(m_current, m_next);
+                sort(k / 2, out);
+                for (unsigned i = 0; i < out.size() / k; ++i) {
+                    unsigned ki = k * i;
+                    for (unsigned j = 0; j < k / 2; ++j) {
+                        next(ki + (2 * j)) = current(ki + j);
+                        next(ki + (2 * j) + 1) = current(ki + (k / 2) + j);
+                    }
+                    
+                    for (unsigned j = 0; j < (k / 2) - 1; ++j) {
+                        exchange(next(ki + (2 * j) + 1), next(ki + (2 * (j + 1))), out);
+                    }
+                }
+                std::swap(m_current, m_next);
+            }
+        }        
+
+        bool is_power_of2(unsigned n) const {
+            return n != 0 && ((n-1) & n) == 0;
+        }
+        
+    public:
+        sorting_network(Ext& ext):
+            m_ext(ext),
+            m_current(&m_currentv),
+            m_next(&m_nextv)
+        {}
+        
+        void operator()(vect const& in, vect& out) {
+            out.reset();
+            out.append(in);
+            if (in.size() <= 1) {
+                return;
+            }
+            while (!is_power_of2(out.size())) {
+                out.push_back(m_ext.mk_default());
+            }
+            for (unsigned i = 0; i < out.size(); ++i) {
+                m_currentv.push_back(i);
+                m_nextv.push_back(i);
+            }
+            unsigned k = 2;
+            while (k <= out.size()) {
+                sort(k, out);
+                k *= 2;
+            }
+        }        
+    };
+
+    // parametric sorting network
+    // Described in Abio et.al. CP 2013.
+    template<class psort_expr>
+    class psort_nw {
+        typedef typename psort_expr::literal literal;
+        typedef typename psort_expr::literal_vector literal_vector;
+
+        class vc {
+            unsigned v; // number of vertices
+            unsigned c; // number of clauses
+            static const unsigned lambda = 5;
+        public:
+            vc(unsigned v, unsigned c):v(v), c(c) {}
+
+            bool operator<(vc const& other) const {
+                return to_int() < other.to_int();
+            }
+            vc operator+(vc const& other) const {
+                return vc(v + other.v, c + other.c);
+            }
+            unsigned to_int() const {
+                return lambda*v + c;
+            }
+            vc operator*(unsigned n) const {
+                return vc(n*v, n*c);
+            }
+        };
+
+        static vc mk_min(vc const& v1, vc const& v2) {
+            return (v1.to_int() < v2.to_int())?v1:v2;
+        }
+
+
+        enum cmp_t { LE, GE, EQ, GE_FULL, LE_FULL };
+        psort_expr&  ctx;
+        cmp_t        m_t;
+
+        // for testing
+        static const bool m_disable_dcard    = false;
+        static const bool m_disable_dsorting = false;
+        static const bool m_disable_dsmerge  = false;
+        static const bool m_force_dcard      = false;
+        static const bool m_force_dsorting   = false;
+        static const bool m_force_dsmerge    = false;
+
+    public:
+        struct stats {
+            unsigned m_num_compiled_vars;
+            unsigned m_num_compiled_clauses;
+            void reset() { memset(this, 0, sizeof(*this)); }
+            stats() { reset(); }
+        };
+        stats        m_stats;
+
+        psort_nw(psort_expr& c): ctx(c) {}
+
+        literal ge(bool full, unsigned k, unsigned n, literal const* xs) {
+            if (k > n) {
+                return ctx.mk_false();
+            }
+            if (k == 0) {
+                return ctx.mk_true();
+            }
+            SASSERT(0 < k && k <= n);
+            literal_vector in, out;
+            if (dualize(k, n, xs, in)) {
+                return le(full, k, in.size(), in.c_ptr());
+            }
+            else {
+                SASSERT(2*k <= n);
+                m_t = full?GE_FULL:GE;
+                psort_nw<psort_expr>::card(k, n, xs, out);                
+                return out[k-1]; 
+            }
+        }
+        
+        literal le(bool full, unsigned k, unsigned n, literal const* xs) {
+            if (k >= n) {
+                return ctx.mk_true();
+            }
+            SASSERT(k < n);
+            literal_vector in, out;
+            if (dualize(k, n, xs, in)) {
+                return ge(full, k, n, in.c_ptr());
+            }
+            else {
+                SASSERT(2*k <= n);
+                m_t = full?LE_FULL:LE;
+                card(k + 1, n, xs, out);
+                return ctx.mk_not(out[k]); 
+            }
+        }
+
+        literal eq(unsigned k, unsigned n, literal const* xs) {
+            if (k > n) {
+                return ctx.mk_false();
+            }
+            SASSERT(k <= n);
+            literal_vector in, out;
+            if (dualize(k, n, xs, in)) {
+                return eq(k, n, in.c_ptr());
+            }
+            else {
+                SASSERT(2*k <= n);
+                m_t = EQ;
+                card(k+1, n, xs, out);
+                SASSERT(out.size() >= k+1);
+                return ctx.mk_min(out[k-1], ctx.mk_not(out[k]));
+            }
+        }
+
+        
+    private:
+
+        std::ostream& pp(std::ostream& out, unsigned n, literal const* lits) {
+            for (unsigned i = 0; i < n; ++i) ctx.pp(out, lits[i]) << " ";
+            return out;
+        }
+
+        std::ostream& pp(std::ostream& out, literal_vector const& lits) {
+            for (unsigned i = 0; i < lits.size(); ++i) ctx.pp(out, lits[i]) << " ";
+            return out;
+        }
+
+        // 0 <= k <= N
+        //     SUM x_i >= k 
+        // <=>
+        //     SUM ~x_i <= N - k
+        // suppose k > N/2, then it is better to solve dual.
+
+        bool dualize(unsigned& k, unsigned N, literal const* xs, literal_vector& in) {
+            SASSERT(0 <= k && k <= N);
+            if (2*k <= N) {
+                return false;
+            }
+            k = N - k;
+            for (unsigned i = 0; i < N; ++i) {
+                in.push_back(ctx.mk_not(xs[i]));
+            }
+            TRACE("pb", 
+                  pp(tout << N << ": ", in);                  
+                  tout << " ~ " << k << "\n";);
+            return true;
+        }
+        
+
+        bool even(unsigned n) const { return (0 == (n & 0x1)); }
+        bool odd(unsigned n) const { return !even(n); }
+        unsigned ceil2(unsigned n) const { return n/2 + odd(n); }
+        unsigned floor2(unsigned n) const { return n/2; }
+        unsigned power2(unsigned n) const { SASSERT(n < 10); return 1 << n; }
+
+
+        literal mk_max(literal a, literal b) {
+            if (a == b) return a;
+            m_stats.m_num_compiled_vars++;
+            return ctx.mk_max(a, b);
+        }
+
+        literal mk_min(literal a, literal b) {
+            if (a == b) return a;
+            m_stats.m_num_compiled_vars++;
+            return ctx.mk_min(a, b);
+        }
+
+        literal fresh() {
+            m_stats.m_num_compiled_vars++;
+            return ctx.fresh();
+        }
+        void add_clause(literal l1, literal l2, literal l3) {
+            literal lits[3] = { l1, l2, l3 };
+            add_clause(3, lits);
+        }
+        void add_clause(literal l1, literal l2) {
+            literal lits[2] = { l1, l2 };
+            add_clause(2, lits);
+        }
+        void add_clause(unsigned n, literal const* ls) {
+            m_stats.m_num_compiled_clauses++;
+            literal_vector tmp(n, ls);
+            ctx.mk_clause(n, tmp.c_ptr());
+        }
+
+        // y1 <= mk_max(x1,x2)
+        // y2 <= mk_min(x1,x2)
+        void cmp_ge(literal x1, literal x2, literal y1, literal y2) {
+            add_clause(ctx.mk_not(y2), x1);
+            add_clause(ctx.mk_not(y2), x2);
+            add_clause(ctx.mk_not(y1), x1, x2);            
+        }
+
+        // mk_max(x1,x2) <= y1
+        // mk_min(x1,x2) <= y2
+        void cmp_le(literal x1, literal x2, literal y1, literal y2) {
+            add_clause(ctx.mk_not(x1), y1);
+            add_clause(ctx.mk_not(x2), y1);
+            add_clause(ctx.mk_not(x1), ctx.mk_not(x2), y2);
+        }
+
+        void cmp_eq(literal x1, literal x2, literal y1, literal y2) {
+            cmp_ge(x1, x2, y1, y2);
+            cmp_le(x1, x2, y1, y2);
+        }
+
+        void cmp(literal x1, literal x2, literal y1, literal y2) {
+            switch(m_t) {
+            case LE: case LE_FULL: cmp_le(x1, x2, y1, y2); break;
+            case GE: case GE_FULL: cmp_ge(x1, x2, y1, y2); break;
+            case EQ: cmp_eq(x1, x2, y1, y2); break;
+            }
+        }
+        vc vc_cmp() {
+            return vc(2, (m_t==EQ)?6:3);
+        }
+
+        void card(unsigned k, unsigned n, literal const* xs, literal_vector& out) {
+            TRACE("pb", tout << "card k:" << k << " n: " << n << "\n";);
+            if (n <= k) {
+                psort_nw<psort_expr>::sorting(n, xs, out);
+            }
+            else if (use_dcard(k, n)) {
+                dsorting(k, n, xs, out);
+            }
+            else {
+                literal_vector out1, out2;
+                unsigned l = n/2; // TBD
+                card(k, l, xs, out1);
+                card(k, n-l, xs + l, out2);
+                smerge(k, out1.size(), out1.c_ptr(), out2.size(), out2.c_ptr(), out);
+            }
+            TRACE("pb", tout << "card k:" << k << " n: " << n << "\n";
+                  pp(tout << "in:", n, xs) << "\n";
+                  pp(tout << "out:", out) << "\n";);
+
+        }        
+        vc vc_card(unsigned k, unsigned n) {
+            if (n <= k) {
+                return vc_sorting(n);
+            }
+            else if (use_dcard(k, n)) {
+                return vc_dsorting(k, n);
+            }
+            else {
+                return vc_card_rec(k, n);
+            }
+        }
+        vc vc_card_rec(unsigned k, unsigned n) {
+            unsigned l = n/2;
+            return vc_card(k, l) + vc_card(k, n-l) + vc_smerge(k, l, n-l);
+        }
+        bool use_dcard(unsigned k, unsigned n) {
+            return m_force_dcard || (!m_disable_dcard && n < 10 && vc_dsorting(k, n) < vc_card_rec(k, n));
+        }
+
+
+        void merge(unsigned a, literal const* as, 
+                   unsigned b, literal const* bs,                    
+                   literal_vector& out) {
+            TRACE("pb", tout << "merge a: " << a << " b: " << b << "\n";);
+            if (a == 1 && b == 1) {
+                literal y1 = mk_max(as[0], bs[0]);
+                literal y2 = mk_min(as[0], bs[0]);
+                out.push_back(y1);
+                out.push_back(y2);
+                psort_nw<psort_expr>::cmp(as[0], bs[0], y1, y2);
+            }
+            else if (a == 0) {
+                out.append(b, bs);
+            }
+            else if (b == 0) {
+                out.append(a, as);
+            }
+            else if (use_dsmerge(a, b, a + b)) {
+                dsmerge(a + b, a, as, b, bs, out);
+            }
+            else if (even(a) && odd(b)) {
+                merge(b, bs, a, as, out);                
+            }
+            else {
+                literal_vector even_a, odd_a;
+                literal_vector even_b, odd_b;
+                literal_vector out1, out2;
+                SASSERT(a > 1 || b > 1);
+                split(a, as, even_a, odd_a);
+                split(b, bs, even_b, odd_b);
+                SASSERT(!even_a.empty());
+                SASSERT(!even_b.empty());
+                merge(even_a.size(), even_a.c_ptr(),
+                      even_b.size(), even_b.c_ptr(), out1);
+                merge(odd_a.size(), odd_a.c_ptr(),
+                      odd_b.size(), odd_b.c_ptr(), out2);
+                interleave(out1, out2, out); 
+            }
+            TRACE("pb", tout << "merge a: " << a << " b: " << b << "\n";
+                  pp(tout << "a:", a, as) << "\n";
+                  pp(tout << "b:", b, bs) << "\n";
+                  pp(tout << "out:", out) << "\n";);
+        }
+        vc vc_merge(unsigned a, unsigned b) {
+            if (a == 1 && b == 1) {
+                return vc_cmp();
+            }
+            else if (a == 0 || b == 0) {
+                return vc(0, 0);
+            }
+            else if (use_dsmerge(a, b, a + b)) {
+                return vc_dsmerge(a, b, a + b);
+            }
+            else {
+                return vc_merge_rec(a, b);
+            }
+        }
+        vc vc_merge_rec(unsigned a, unsigned b) {
+            return 
+                vc_merge(ceil2(a),  ceil2(b)) +
+                vc_merge(floor2(a), floor2(b)) +
+                vc_interleave(ceil2(a) + ceil2(b), floor2(a) + floor2(b));
+        }
+        void split(unsigned n, literal const* ls, literal_vector& even, literal_vector& odd) {
+            for (unsigned i = 0; i < n; i += 2) {
+                even.push_back(ls[i]);
+            }
+            for (unsigned i = 1; i < n; i += 2) {
+                odd.push_back(ls[i]);
+            }
+        }
+
+        void interleave(literal_vector const& as, 
+                        literal_vector const& bs,
+                        literal_vector& out) {
+            TRACE("pb", tout << "interleave: " << as.size() << " " << bs.size() << "\n";);
+            SASSERT(as.size() >= bs.size());
+            SASSERT(as.size() <= bs.size() + 2);
+            SASSERT(!as.empty());
+            out.push_back(as[0]);
+            unsigned sz = std::min(as.size()-1, bs.size());
+            for (unsigned i = 0; i < sz; ++i) {
+                literal y1 = mk_max(as[i+1],bs[i]);
+                literal y2 = mk_min(as[i+1],bs[i]);
+                psort_nw<psort_expr>::cmp(as[i+1], bs[i], y1, y2);
+                out.push_back(y1);
+                out.push_back(y2);
+            }            
+            if (as.size() == bs.size()) {
+                out.push_back(bs[sz]);
+            }
+            else if (as.size() == bs.size() + 2) {
+                out.push_back(as[sz+1]);
+            }
+            SASSERT(out.size() == as.size() + bs.size());
+            TRACE("pb", tout << "interleave: " << as.size() << " " << bs.size() << "\n";
+                  pp(tout << "a: ", as) << "\n";
+                  pp(tout << "b: ", bs) << "\n";
+                  pp(tout << "out: ", out) << "\n";);
+
+        }
+        vc vc_interleave(unsigned a, unsigned b) {
+            return vc_cmp()*std::min(a-1,b);
+        }
+        
+        void sorting(unsigned n, literal const* xs, literal_vector& out) {
+            TRACE("pb", tout << "sorting: " << n << "\n";);
+            switch(n) {
+            case 0: 
+                break;
+            case 1: 
+                out.push_back(xs[0]); 
+                break;
+            case 2: 
+                psort_nw<psort_expr>::merge(1, xs, 1, xs+1, out); 
+                break;
+            default:
+                if (use_dsorting(n)) {
+                    dsorting(n, n, xs, out);
+                }
+                else {
+                    literal_vector out1, out2;
+                    unsigned l = n/2;  // TBD
+                    sorting(l, xs, out1);
+                    sorting(n-l, xs+l, out2);
+                    merge(out1.size(), out1.c_ptr(), 
+                          out2.size(), out2.c_ptr(),
+                          out);                
+                }
+                break;
+            }
+            TRACE("pb", tout << "sorting: " << n << "\n";
+                  pp(tout << "in:", n, xs) << "\n"; 
+                  pp(tout << "out:", out) << "\n";);
+
+        }
+        vc vc_sorting(unsigned n) {
+            switch(n) {
+            case 0: return vc(0,0);
+            case 1: return vc(0,0);
+            case 2: return vc_merge(1,1);
+            default:
+                if (use_dsorting(n)) {
+                    return vc_dsorting(n, n);
+                }
+                else {
+                    return vc_sorting_rec(n);
+                }
+            }
+        }
+        vc vc_sorting_rec(unsigned n) {
+            SASSERT(n > 2);
+            unsigned l = n/2;
+            return vc_sorting(l) + vc_sorting(n-l) + vc_merge(l, n-l);
+        }
+
+        bool use_dsorting(unsigned n) {
+            SASSERT(n > 2);            
+            return m_force_dsorting || 
+                (!m_disable_dsorting && n < 10 && vc_dsorting(n, n) < vc_sorting_rec(n));
+        }
+
+        void smerge(unsigned c,
+                    unsigned a, literal const* as,
+                    unsigned b, literal const* bs, 
+                    literal_vector& out) {
+            TRACE("pb", tout << "smerge: c:" << c << " a:" << a << " b:" << b << "\n";);
+            if (a == 1 && b == 1 && c == 1) {
+                literal y = mk_max(as[0], bs[0]);
+                if (m_t != GE) {
+                    // x1 <= mk_max(x1,x2)
+                    // x2 <= mk_max(x1,x2)
+                    add_clause(ctx.mk_not(as[0]), y);
+                    add_clause(ctx.mk_not(bs[0]), y);
+                }
+                if (m_t != LE) {
+                    // mk_max(x1,x2) <= x1, x2
+                    add_clause(ctx.mk_not(y), as[0], bs[0]);
+                }
+                out.push_back(y);
+            }
+            else if (a == 0) {
+                out.append(std::min(c, b), bs);
+            }
+            else if (b == 0) {
+                out.append(std::min(c, a), as);
+            }
+            else if (a > c) {
+                smerge(c, c, as, b, bs, out);
+            }
+            else if (b > c) {
+                smerge(c, a, as, c, bs, out);
+            }
+            else if (a + b <= c) {
+                merge(a, as, b, bs, out);
+            }
+            else if (use_dsmerge(a, b, c)) {
+                dsmerge(c, a, as, b, bs, out);
+            }
+            else {
+                literal_vector even_a, odd_a;
+                literal_vector even_b, odd_b;
+                literal_vector out1, out2;               
+                split(a, as, even_a, odd_a);
+                split(b, bs, even_b, odd_b);
+                SASSERT(!even_a.empty());
+                SASSERT(!even_b.empty());
+                unsigned c1, c2;
+                if (even(c)) {
+                    c1 = 1 + c/2; c2 = c/2;
+                }
+                else {
+                    c1 = (c + 1)/2; c2 = (c - 1)/2;
+                }
+                smerge(c1, even_a.size(), even_a.c_ptr(),
+                       even_b.size(), even_b.c_ptr(), out1);
+                smerge(c2, odd_a.size(), odd_a.c_ptr(),
+                       odd_b.size(), odd_b.c_ptr(), out2);
+                SASSERT(out1.size() == std::min(even_a.size()+even_b.size(), c1));
+                SASSERT(out2.size() == std::min(odd_a.size()+odd_b.size(), c2));
+                literal y;
+                if (even(c)) {
+                    literal z1 = out1.back();
+                    literal z2 = out2.back();
+                    out1.pop_back();
+                    out2.pop_back();
+                    y = mk_max(z1, z2);
+                    if (m_t != GE) {
+                        add_clause(ctx.mk_not(z1), y);
+                        add_clause(ctx.mk_not(z2), y);
+                    }
+                    if (m_t != LE) {
+                        add_clause(ctx.mk_not(y), z1, z2);
+                    }
+                }
+                interleave(out1, out2, out);                 
+                if (even(c)) {
+                    out.push_back(y);
+                }
+            }
+            TRACE("pb", tout << "smerge: c:" << c << " a:" << a << " b:" << b << "\n";
+                  pp(tout << "a:", a, as) << "\n";
+                  pp(tout << "b:", b, bs) << "\n";
+                  pp(tout << "out:", out) << "\n";
+                  );
+            SASSERT(out.size() == std::min(a + b, c));
+        }
+
+        vc vc_smerge(unsigned a, unsigned b, unsigned c) {
+            if (a == 1 && b == 1 && c == 1) {
+                vc v(1,0);
+                if (m_t != GE) v = v + vc(0, 2);
+                if (m_t != LE) v = v + vc(0, 1);
+                return v;
+            }
+            if (a == 0 || b == 0) return vc(0, 0);
+            if (a > c) return vc_smerge(c, b, c);
+            if (b > c) return vc_smerge(a, c, c);
+            if (a + b <= c) return vc_merge(a, b);
+            if (use_dsmerge(a, b, c)) return vc_dsmerge(a, b, c);
+            return vc_smerge_rec(a, b, c);
+        }
+        vc vc_smerge_rec(unsigned a, unsigned b, unsigned c) {
+            return 
+                vc_smerge(ceil2(a), ceil2(b), even(c)?(1+c/2):((c+1)/2)) +
+                vc_smerge(floor2(a), floor2(b), even(c)?(c/2):((c-1)/2)) +
+                vc_interleave(ceil2(a)+ceil2(b),floor2(a)+floor2(b)) +
+                vc(1, 0) +
+                ((m_t != GE)?vc(0, 2):vc(0, 0)) +
+                ((m_t != LE)?vc(0, 1):vc(0, 0));
+        }
+        bool use_dsmerge(unsigned a, unsigned b, unsigned c) {
+            return 
+                m_force_dsmerge ||
+                (!m_disable_dsmerge && 
+                 a < (1 << 15) && b < (1 << 15) && 
+                 vc_dsmerge(a, b, a + b) < vc_smerge_rec(a, b, c));
+        }
+
+        void dsmerge(            
+            unsigned c,
+            unsigned a, literal const* as, 
+            unsigned b, literal const* bs,                    
+            literal_vector& out) {
+            TRACE("pb", tout << "dsmerge: c:" << c << " a:" << a << " b:" << b << "\n";);
+            SASSERT(a <= c);
+            SASSERT(b <= c);
+            SASSERT(a + b >= c);
+            for (unsigned i = 0; i < c; ++i) {
+                out.push_back(fresh());
+            }
+            if (m_t != GE) {
+                for (unsigned i = 0; i < a; ++i) {
+                    add_clause(ctx.mk_not(as[i]), out[i]);
+                }
+                for (unsigned i = 0; i < b; ++i) {
+                    add_clause(ctx.mk_not(bs[i]), out[i]);
+                }
+                for (unsigned i = 1; i <= a; ++i) {
+                    for (unsigned j = 1; j <= b && i + j <= c; ++j) {
+                        add_clause(ctx.mk_not(as[i-1]),ctx.mk_not(bs[j-1]),out[i+j-1]);
+                    }
+                }
+            }
+            if (m_t != LE) {
+                literal_vector ls;
+                for (unsigned k = 0; k < c; ++k) {
+                    ls.reset();
+                    ls.push_back(ctx.mk_not(out[k]));
+                    if (a <= k) {
+                        add_clause(ctx.mk_not(out[k]), bs[k-a]);
+                    }
+                    if (b <= k) {
+                        add_clause(ctx.mk_not(out[k]), as[k-b]);
+                    }
+                    for (unsigned i = 0; i < std::min(a,k + 1); ++i) {
+                        unsigned j = k - i;
+                        SASSERT(i + j == k);
+                        if (j < b) {
+                            ls.push_back(as[i]);
+                            ls.push_back(bs[j]);
+                            add_clause(ls.size(), ls.c_ptr());
+                            ls.pop_back();
+                            ls.pop_back();
+                        }
+                    }
+                }
+            }
+        }
+        vc vc_dsmerge(unsigned a, unsigned b, unsigned c) {
+            vc v(c, 0);
+            if (m_t != GE) {
+                v = v + vc(0, a + b + std::min(a, c)*std::min(b, c)/2);
+            }
+            if (m_t != LE) {
+                v = v + vc(0, std::min(a, c)*std::min(b, c)/2);
+            }
+            return v;
+        }
+
+
+        void dsorting(unsigned m, unsigned n, literal const* xs, 
+                      literal_vector& out) {
+            TRACE("pb", tout << "dsorting m: " << m << " n: " << n << "\n";);
+            SASSERT(m <= n);
+            literal_vector lits;
+            for (unsigned i = 0; i < m; ++i) {
+                out.push_back(fresh());
+            }
+            if (m_t != GE) {
+                for (unsigned k = 1; k <= m; ++k) {
+                    lits.push_back(out[k-1]);
+                    add_subset(true, k, 0, lits, n, xs);
+                    lits.pop_back();
+                }
+            }
+            if (m_t != LE) {                
+                for (unsigned k = 1; k <= m; ++k) {
+                    lits.push_back(ctx.mk_not(out[k-1]));
+                    add_subset(false, n-k+1, 0, lits, n, xs);
+                    lits.pop_back();
+                }
+            }
+        }
+        vc vc_dsorting(unsigned m, unsigned n) {
+            SASSERT(m <= n && n < 10);            
+            vc v(m, 0);
+            if (m_t != GE) {
+                v = v + vc(0, power2(n-1));
+            }
+            if (m_t != LE) {
+                v = v + vc(0, power2(n-1));
+            }
+            return v;
+        }
+
+        void add_subset(bool polarity, unsigned k, unsigned offset, literal_vector& lits, 
+                        unsigned n, literal const* xs) {
+            TRACE("pb", tout << "k:" << k << " offset: " << offset << " n: " << n << " ";
+                  pp(tout, lits) << "\n";);
+            SASSERT(k + offset <= n);
+            if (k == 0) {
+                add_clause(lits.size(), lits.c_ptr());
+                return;
+            }
+            for (unsigned i = offset; i < n - k + 1; ++i) {
+                lits.push_back(polarity?ctx.mk_not(xs[i]):xs[i]);
+                add_subset(polarity, k-1, i+1, lits, n, xs);
+                lits.pop_back();
+            }
+        }
+    };
+
+#endif

src/util/stopwatch.h

@@ -69,6 +69,7 @@ public:
 #undef ARRAYSIZE
 #define ARRAYSIZE ARRAYSIZE_TEMP
 #undef max
+#undef min
 
 
 #elif defined(__APPLE__) && defined (__MACH__) // Mac OS X

src/util/trace.cpp

@@ -23,8 +23,8 @@ Revision History:
 std::ofstream tout(".z3-trace"); 
 #endif
 
-bool g_enable_all_trace_tags = false;
-str_hashtable* g_enabled_trace_tags = 0;
+static bool g_enable_all_trace_tags = false;
+static str_hashtable* g_enabled_trace_tags = 0;
 
 static str_hashtable& get_enabled_trace_tags() {
     if (!g_enabled_trace_tags) {

src/util/trace.h

@@ -24,6 +24,10 @@ Revision History:
 #undef max
 #undef min
 #endif
+#ifdef __APPLE__
+#undef max
+#undef min
+#endif
 #include<fstream>
 
 #ifdef _TRACE

src/util/trail.h

@@ -0,0 +1,360 @@
+/*++
+Copyright (c) 2006 Microsoft Corporation
+
+Module Name:
+
+    trail.h
+
+Abstract:
+
+    <abstract>
+
+Author:
+
+    Leonardo de Moura (leonardo) 2008-06-02.
+
+Revision History:
+
+--*/
+#ifndef _TRAIL_H_
+#define _TRAIL_H_
+
+#include"obj_hashtable.h"
+#include"region.h"
+#include"obj_ref.h"
+
+template<typename Ctx>
+class trail {
+public:
+    virtual ~trail() {
+    }
+    virtual void undo(Ctx & ctx) = 0;
+};
+
+template<typename Ctx, typename T>
+class value_trail : public trail<Ctx> {
+    T & m_value;
+    T   m_old_value;
+    
+public:
+    value_trail(T & value):
+        m_value(value),
+        m_old_value(value) {
+    }
+    
+    virtual ~value_trail() {
+    }
+    
+    virtual void undo(Ctx & ctx) {
+        m_value = m_old_value;
+    }
+};
+
+template<typename Ctx>
+class reset_flag_trail : public trail<Ctx> {
+    bool & m_value;
+public:
+    reset_flag_trail(bool & value):
+        m_value(value) {
+    }
+    
+    virtual ~reset_flag_trail() {
+    }
+    
+    virtual void undo(Ctx & ctx) {
+        m_value = false;
+    }
+};
+
+template<typename Ctx, typename T>
+class set_ptr_trail : public trail<Ctx> {
+    T * & m_ptr;
+public:
+    set_ptr_trail(T * & ptr):
+        m_ptr(ptr) {
+        SASSERT(m_ptr == 0);
+    }
+    
+    virtual void undo(Ctx & ctx) {
+        m_ptr = 0;
+    }
+};
+
+template<typename Ctx, typename T, bool CallDestructors=true>
+class restore_size_trail : public trail<Ctx> {
+    vector<T, CallDestructors> & m_vector;
+    unsigned                     m_old_size;
+public:
+    restore_size_trail(vector<T, CallDestructors> & v, unsigned sz):
+        m_vector(v),
+        m_old_size(sz) {
+    }
+    restore_size_trail(vector<T, CallDestructors> & v):
+        m_vector(v),
+        m_old_size(v.size()) {
+    }
+    virtual ~restore_size_trail() {
+    }
+    virtual void undo(Ctx & ctx) {
+        m_vector.shrink(m_old_size);
+    }
+};    
+        
+template<typename Ctx, typename T, bool CallDestructors=true>
+class vector_value_trail : public trail<Ctx> {
+    vector<T, CallDestructors> & m_vector;
+    unsigned                     m_idx;
+    T                            m_old_value;
+public:
+    vector_value_trail(vector<T, CallDestructors> & v, unsigned idx):
+        m_vector(v),
+        m_idx(idx),
+        m_old_value(v[idx]) {
+    }
+    
+    virtual ~vector_value_trail() {
+    }
+    
+    virtual void undo(Ctx & ctx) {
+        m_vector[m_idx] = m_old_value;
+    }
+};
+
+
+template<typename Ctx, typename D, typename R>
+class insert_obj_map : public trail<Ctx> {
+    obj_map<D,R>&     m_map;
+    D*                m_obj;
+public:
+    insert_obj_map(obj_map<D,R>& t, D* o) : m_map(t), m_obj(o) {}
+    virtual ~insert_obj_map() {}
+    virtual void undo(Ctx & ctx) { m_map.remove(m_obj); }
+};
+
+template<typename Ctx, typename M, typename D>
+class insert_map : public trail<Ctx> {
+    M&            m_map;
+    D             m_obj;
+public:
+    insert_map(M& t, D o) : m_map(t), m_obj(o) {}
+    virtual ~insert_map() {}
+    virtual void undo(Ctx & ctx) { m_map.remove(m_obj); }
+};
+
+
+
+template<typename Ctx, typename V>
+class push_back_vector : public trail<Ctx> {
+    V & m_vector;
+public:
+    push_back_vector(V & v):
+        m_vector(v) {
+    }
+    
+    virtual void undo(Ctx & ctx) {
+        m_vector.pop_back();
+    }
+};
+
+template<typename Ctx, typename T>
+class set_vector_idx_trail : public trail<Ctx> {
+    ptr_vector<T> & m_vector;
+    unsigned                         m_idx;
+public:
+    set_vector_idx_trail(ptr_vector<T> & v, unsigned idx):
+        m_vector(v),
+        m_idx(idx) {
+    }
+    
+    virtual ~set_vector_idx_trail() {
+    }
+    
+    virtual void undo(Ctx & ctx) {
+        m_vector[m_idx] = 0;
+    }
+};
+    
+template<typename Ctx, typename T, bool CallDestructors=true>
+class pop_back_trail : public trail<Ctx> {
+    vector<T, CallDestructors> & m_vector;
+    T m_value;
+public:
+    pop_back_trail(vector<T, CallDestructors> & v):
+    m_vector(v),
+    m_value(m_vector.back()) {
+    }
+    
+    virtual void undo(Ctx & ctx) {
+        m_vector.push_back(m_value);
+    }
+};
+
+template<typename Ctx, typename T, bool CallDestructors=true>
+class pop_back2_trail : public trail<Ctx> {
+    vector<T, CallDestructors> & m_vector;
+    typedef vector<vector<T, CallDestructors>, true> vector_t;
+    unsigned m_index;
+    T m_value;
+public:
+    pop_back2_trail(vector<T, CallDestructors> & v, unsigned index):
+    m_vector(v),
+    m_index(index),
+    m_value(m_vector[index].back()) {
+    }
+    
+    virtual void undo(Ctx & ctx) {
+        m_vector[m_index].push_back(m_value);
+    }
+};
+
+
+
+template<typename Ctx, typename T, bool CallDestructors=true>
+class push_back_trail : public trail<Ctx> {
+    vector<T, CallDestructors> & m_vector;
+public:
+    push_back_trail(vector<T, CallDestructors> & v):
+        m_vector(v) {
+    }
+    
+    virtual void undo(Ctx & ctx) {
+        m_vector.pop_back();
+    }
+};
+
+template<typename Ctx, typename T, bool CallDestructors=true>
+class push_back2_trail : public trail<Ctx> {
+    typedef vector<vector<T, CallDestructors>, true> vector_t;
+    vector_t & m_vector;
+    unsigned   m_index;
+public:
+    push_back2_trail(vector_t & v, unsigned index) : 
+    m_vector(v),
+    m_index(index) {
+    }
+    
+    virtual void undo(Ctx & ctx) {
+        m_vector[m_index].pop_back();
+    }
+};
+
+template<typename Ctx>
+class set_bitvector_trail : public trail<Ctx> {
+    svector<bool> & m_vector;
+    unsigned        m_idx;
+public:
+    set_bitvector_trail(svector<bool> & v, unsigned idx):
+        m_vector(v),
+        m_idx(idx) {
+        SASSERT(m_vector[m_idx] == false);
+        m_vector[m_idx] = true;
+    }
+    
+    virtual void undo(Ctx & ctx) {
+        m_vector[m_idx] = false;
+    }
+};
+    
+template<typename Ctx, typename T>
+class new_obj_trail : public trail<Ctx> {
+    T * m_obj;
+public:
+    new_obj_trail(T * obj):
+        m_obj(obj) {
+    }
+    
+    virtual void undo(Ctx & ctx) {
+        dealloc(m_obj);
+    }
+};
+
+template<typename Ctx, typename M, typename T>
+class obj_ref_trail : public trail<Ctx> {
+    obj_ref<T,M> m_obj;
+public:
+    obj_ref_trail(obj_ref<T,M>& obj):
+        m_obj(obj) {
+    }
+    
+    virtual void undo(Ctx & ctx) {
+        m_obj.reset();
+    }
+};
+
+template<typename Ctx, typename T>
+class insert_obj_trail : public trail<Ctx> {
+    obj_hashtable<T>& m_table;
+    T*                m_obj;
+public:
+    insert_obj_trail(obj_hashtable<T>& t, T* o) : m_table(t), m_obj(o) {}
+    virtual ~insert_obj_trail() {}
+    virtual void undo(Ctx & ctx) { m_table.remove(m_obj); }
+};
+
+
+
+template<typename Ctx, typename T>
+class remove_obj_trail : public trail<Ctx> {
+    obj_hashtable<T>& m_table;
+    T*                m_obj;
+public:
+    remove_obj_trail(obj_hashtable<T>& t, T* o) : m_table(t), m_obj(o) {}
+    virtual ~remove_obj_trail() {}
+    virtual void undo(Ctx & ctx) { m_table.insert(m_obj); } 
+};
+
+
+template<typename Ctx>
+void undo_trail_stack(Ctx & ctx, ptr_vector<trail<Ctx> > & s, unsigned old_size) {
+    SASSERT(old_size <= s.size());
+    typename ptr_vector<trail<Ctx> >::iterator begin = s.begin() + old_size;
+    typename ptr_vector<trail<Ctx> >::iterator it    = s.end();
+    while (it != begin) {
+        --it;
+        (*it)->undo(ctx);
+    }
+    s.shrink(old_size);
+}
+
+template<typename Ctx>
+class trail_stack {
+    Ctx &                   m_ctx;
+    ptr_vector<trail<Ctx> > m_trail_stack;
+    unsigned_vector         m_scopes;
+    region                  m_region;
+public:
+    trail_stack(Ctx & c):m_ctx(c) {}
+
+    ~trail_stack() {}
+    
+    region & get_region() { return m_region; }
+    
+    void reset() { 
+        pop_scope(m_scopes.size()); 
+        // Undo trail objects stored at lvl 0 (avoid memory leaks if lvl 0 contains new_obj_trail objects).
+        undo_trail_stack(m_ctx, m_trail_stack, 0); 
+    }
+
+    void push_ptr(trail<Ctx> * t) { m_trail_stack.push_back(t); }
+
+    template<typename TrailObject>
+    void push(TrailObject const & obj) { m_trail_stack.push_back(new (m_region) TrailObject(obj)); }
+
+    unsigned get_num_scopes() const { return m_scopes.size(); }
+
+    void push_scope() { m_region.push_scope(); m_scopes.push_back(m_trail_stack.size()); }
+
+    void pop_scope(unsigned num_scopes) {
+        if (num_scopes == 0) return;
+        unsigned lvl      = m_scopes.size();
+        SASSERT(num_scopes <= lvl);
+        unsigned new_lvl  = lvl - num_scopes;
+        unsigned old_size = m_scopes[new_lvl];
+        undo_trail_stack(m_ctx, m_trail_stack, old_size);
+        m_scopes.shrink(new_lvl);
+        m_region.pop_scope(num_scopes);
+    }
+};
+
+#endif /* _TRAIL_H_ */
+

src/util/union_find.h

@@ -0,0 +1,243 @@
+/*++
+Copyright (c) 2006 Microsoft Corporation
+
+Module Name:
+
+    union_find.h
+
+Abstract:
+
+    <abstract>
+
+Author:
+
+    Leonardo de Moura (leonardo) 2008-05-31.
+
+Revision History:
+
+--*/
+#ifndef _UNION_FIND_H_
+#define _UNION_FIND_H_
+
+#include "trail.h"
+#include "trace.h"
+
+class union_find_default_ctx {
+public:
+    typedef trail_stack<union_find_default_ctx> _trail_stack;
+    union_find_default_ctx() : m_stack(*this) {}
+       
+    void unmerge_eh(unsigned, unsigned) {}
+    void merge_eh(unsigned, unsigned, unsigned, unsigned) {}
+    void after_merge_eh(unsigned, unsigned, unsigned, unsigned) {}
+
+    _trail_stack& get_trail_stack() { return m_stack; }
+
+private:
+    _trail_stack m_stack;
+};
+
+template<typename Ctx = union_find_default_ctx>
+class union_find {
+    Ctx &                         m_ctx;
+    trail_stack<Ctx> &            m_trail_stack;
+    svector<unsigned>             m_find;
+    svector<unsigned>             m_size;
+    svector<unsigned>             m_next;
+    
+    class mk_var_trail;
+    friend class mk_var_trail;
+
+    class mk_var_trail : public trail<Ctx> {
+        union_find & m_owner;
+    public:
+        mk_var_trail(union_find & o):m_owner(o) {}
+        virtual ~mk_var_trail() {}
+        virtual void undo(Ctx & ctx) {
+            m_owner.m_find.pop_back();
+            m_owner.m_size.pop_back();
+            m_owner.m_next.pop_back();
+        }
+    };
+
+    mk_var_trail                  m_mk_var_trail;
+
+    class merge_trail;
+    friend class merge_trail;
+
+    class merge_trail : public trail<Ctx> {
+        union_find & m_owner;
+        unsigned     m_r1;
+    public:
+        merge_trail(union_find & o, unsigned r1):m_owner(o), m_r1(r1) {}
+        virtual ~merge_trail() {}
+        virtual void undo(Ctx & ctx) { m_owner.unmerge(m_r1); }
+    };
+
+    void unmerge(unsigned r1) {
+        unsigned r2 = m_find[r1];
+        TRACE("union_find", tout << "unmerging " << r1 << " " << r2 << "\n";);
+        SASSERT(find(r2) == r2);
+        m_size[r2] -= m_size[r1];
+        m_find[r1]  = r1;
+        std::swap(m_next[r1], m_next[r2]);
+        m_ctx.unmerge_eh(r2, r1);
+        CASSERT("union_find", check_invariant());
+    }
+
+public:    
+    union_find(Ctx & ctx):m_ctx(ctx), m_trail_stack(ctx.get_trail_stack()), m_mk_var_trail(*this) {}
+
+    unsigned mk_var() {
+        unsigned r = m_find.size();
+        m_find.push_back(r);
+        m_size.push_back(1);
+        m_next.push_back(r);
+        m_trail_stack.push_ptr(&m_mk_var_trail);
+        return r;
+    }
+
+    unsigned get_num_vars() const { return m_find.size(); }
+
+
+    unsigned find(unsigned v) const {
+        while (true) {
+            unsigned new_v = m_find[v];
+            if (new_v == v)
+                return v;
+            v = new_v;
+        }
+    }
+
+    unsigned next(unsigned v) const { return m_next[v]; }
+
+    unsigned size(unsigned v) const { return m_size[find(v)]; }
+
+    bool is_root(unsigned v) const { return m_find[v] == v; }
+
+    void merge(unsigned v1, unsigned v2) {
+        unsigned r1 = find(v1);
+        unsigned r2 = find(v2);
+        TRACE("union_find", tout << "merging " << r1 << " " << r2 << "\n";);
+        if (r1 == r2)
+            return;
+        if (m_size[r1] > m_size[r2])
+            std::swap(r1, r2);
+        m_ctx.merge_eh(r2, r1, v2, v1);
+        m_find[r1] = r2;
+        m_size[r2] += m_size[r1];
+        std::swap(m_next[r1], m_next[r2]);
+        m_trail_stack.push(merge_trail(*this, r1));
+        m_ctx.after_merge_eh(r2, r1, v2, v1);
+        CASSERT("union_find", check_invariant());
+    }
+
+    // dissolve equivalence class of v
+    // this method cannot be used with backtracking.
+    void dissolve(unsigned v) {
+        unsigned w;
+        do {
+            w = next(v);                        
+            m_size[v] = 1;
+            m_find[v] = v;
+            m_next[v] = v;            
+        }
+        while (w != v);
+    }
+
+    void display(std::ostream & out) const {
+        unsigned num = get_num_vars(); 
+        for (unsigned v = 0; v < num; v++) {
+            out << "v" << v << " --> v" << m_find[v] << " (" << size(v) << ")\n";
+        }
+    }
+
+#ifdef Z3DEBUG
+    bool check_invariant() const {
+        unsigned num = get_num_vars(); 
+        for (unsigned v = 0; v < num; v++) {
+            if (is_root(v)) {
+                unsigned curr = v;
+                unsigned sz   = 0;
+                do {
+                    SASSERT(find(curr) == v);
+                    sz++;
+                    curr = next(curr);
+                }
+                while (curr != v);
+                SASSERT(m_size[v] == sz);
+            }
+        }
+        return true;
+    }
+#endif
+};
+
+
+class basic_union_find {
+    unsigned_vector   m_find;
+    unsigned_vector   m_size;
+    unsigned_vector   m_next;
+    
+    void ensure_size(unsigned v) {
+        while (v >= get_num_vars()) {
+            mk_var();
+        }
+    }
+ public:
+    unsigned mk_var() {
+        unsigned r = m_find.size();
+        m_find.push_back(r);
+        m_size.push_back(1);
+        m_next.push_back(r);
+        return r;
+    }
+    unsigned get_num_vars() const { return m_find.size(); }
+    
+    unsigned find(unsigned v) const {
+        if (v >= get_num_vars()) {
+            return v;
+        }
+        while (true) {
+            unsigned new_v = m_find[v];
+            if (new_v == v)
+                return v;
+            v = new_v;
+        }
+    }
+    
+    unsigned next(unsigned v) const { 
+        if (v >= get_num_vars()) {
+            return v;
+        }
+        return m_next[v]; 
+    }
+    
+    bool is_root(unsigned v) const { 
+        return v >= get_num_vars() || m_find[v] == v; 
+    }
+    
+    void merge(unsigned v1, unsigned v2) {
+        unsigned r1 = find(v1);
+        unsigned r2 = find(v2);
+        if (r1 == r2)
+            return;
+        ensure_size(v1);
+        ensure_size(v2);
+        if (m_size[r1] > m_size[r2])
+            std::swap(r1, r2);
+        m_find[r1] = r2;
+        m_size[r2] += m_size[r1];
+        std::swap(m_next[r1], m_next[r2]);
+    }
+    
+    void reset() {
+        m_find.reset();
+        m_next.reset();
+        m_size.reset();
+    }
+};
+
+
+#endif /* _UNION_FIND_H_ */
+