Commits

Leonardo de Moura committed 683687b

more cleanup

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

  • Participants
  • Parent commits c2e95bb

Comments (0)

Files changed (24)

File scripts/mk_project.py

     add_lib('polynomial', ['util'], 'math/polynomial')
     add_lib('sat', ['util'])
     add_lib('nlsat', ['polynomial', 'sat'])
-    add_lib('subpaving', ['util'], 'math/subpaving')
+    add_lib('interval', ['util'], 'math/interval')
+    add_lib('subpaving', ['interval'], 'math/subpaving')
     add_lib('ast', ['util', 'polynomial'])
     add_lib('rewriter', ['ast', 'polynomial'], 'ast/rewriter')
     add_lib('model', ['rewriter'])

File src/ast/ast_dag_pp.cpp

-/*++
-Copyright (c) 2007 Microsoft Corporation
-
-Module Name:
-
-    ast_dag_pp.cpp
-
-Abstract:
- 
-    AST low level pretty printer.
-    
-Author:
-
-    Leonardo de Moura (leonardo) 2006-10-19.
-    Nikolaj Bjorner (nbjorner) 2007-07-17
-
-Revision History:
-
---*/
-
-#include<iostream>
-#include"for_each_ast.h"
-#include"arith_decl_plugin.h"
-#include"bv_decl_plugin.h"
-#include"ast_pp.h"
-
-class dag_printer {
-    ast_manager&    m_manager;
-    std::ostream &  m_out;
-    ast_mark&       m_mark;
-    bool            m_initialized;
-    svector<symbol> m_names;
-    family_id       m_basic_fid;
-    family_id       m_bv_fid;
-    family_id       m_arith_fid;
-    family_id       m_array_fid;
-    arith_util      m_arith;
-    bv_util         m_bv;
-    bool            m_enable_shortcut;
-
-    void process_ast(ast* a) {
-        for_each_ast(*this, m_mark, a);
-    }
-
-    void process_info(decl_info* info) {
-        if (!info) {
-            return;
-        }
-        unsigned num_params = info->get_num_parameters();
-        for (unsigned i = 0; i < num_params; ++i) {
-            parameter const& p = info->get_parameter(i);
-            
-            if (p.is_ast() && !m_mark.is_marked(p.get_ast())) {
-                process_ast(p.get_ast());
-            }
-        }
-    }
-        
-    template<typename T>
-    void display_children(unsigned num_children, T * const * children) {
-        for (unsigned i = 0; i < num_children; i++) {
-            display_node_id(children[i]);
-        }
-    }
-
-    void display_node_id(ast* n) {
-        unsigned id = n->get_id();
-        switch(n->get_kind()) {
-        case AST_FUNC_DECL:
-        case AST_SORT:
-            m_out << "$d" << (id - (1 << 31)) << " ";
-            break;
-        default:
-            m_out << "$" << id << " ";
-            break;
-        }        
-    }
-
-    void display_parameter(parameter const& p) 
-    {
-        if (p.is_int()) {
-            m_out << p.get_int() << " ";
-        }
-        else if (p.is_ast()) {
-            SASSERT(p.is_ast());
-            display_node_id(p.get_ast());
-        }
-        else if (p.is_rational()) {
-            m_out << p.get_rational() << " ";
-        }
-        else if (p.is_symbol()) {
-            display_symbol(p.get_symbol());
-        }
-        else {
-            UNREACHABLE();
-        }
-    }
-
-    // Display: 
-    //   App name [ parameters] arguments
-    //
-    void display_builtin(app* n) {
-        func_decl* d = n->get_decl();
-        unsigned num_params = d->get_num_parameters();
-
-        m_out << "App ";
-        display_node_id(n);        
-        display_symbol(d->get_name());
-        if (num_params > 0) {
-            m_out << "[ ";
-            for (unsigned i = 0; i < num_params; ++i) {
-                display_parameter(d->get_parameter(i));
-            }
-            m_out << "] ";
-        }
-        display_children(n->get_num_args(), n->get_args());
-        m_out << "\n";
-    }
-
-    void display_info(func_decl_info* info) {
-        if (!info) {
-            return;
-        }
-        m_out << "BUILTIN " << get_family_name(info->get_family_id()) << " " << info->get_decl_kind() << " ";
-
-        if (info->is_associative()) {
-            m_out << ":assoc ";
-        }
-        if (info->is_commutative()) {
-            m_out << ":comm ";
-        }
-        if (info->is_injective()) {
-            m_out << ":inj ";
-        }
-        for (unsigned i = 0; i < info->get_num_parameters(); ++i) {
-            display_parameter(info->get_parameter(i));
-        }
-    }
-
-    void display_info(sort_info* info) {
-        if (!info) {
-            return;
-        }
-        m_out << "BUILTIN " << get_family_name(info->get_family_id()) << " " << info->get_decl_kind() << " ";
-        // TODO: remove this if... it doesn't make sense...
-        if (!info->is_infinite() && !info->is_very_big()) {
-            m_out << "Size " << info->get_num_elements().size() << " ";
-        }
-        for (unsigned i = 0; i < info->get_num_parameters(); ++i) {
-            display_parameter(info->get_parameter(i));
-        }
-    }
-
-    symbol get_family_name(family_id id) {
-        if (id == null_family_id) {
-            return symbol("null");
-        }
-        if (!m_initialized) {
-            svector<symbol> names;
-            svector<family_id> range;
-            m_manager.get_dom(names);
-            m_manager.get_range(range);
-            m_names.resize(range.size());
-            for (unsigned i = 0; i < range.size(); ++i) {
-                SASSERT(range[i] < static_cast<family_id>(range.size()));
-                m_names[range[i]] = names[i];
-            }
-            m_initialized = true;
-        }
-        SASSERT(id < static_cast<family_id>(m_names.size()));
-        return m_names[id];
-    }
-
-    bool has_special_char(char const* s) {
-        while (s && *s) {
-            if (*s == ' ') {
-                return true;
-            }
-            ++s;
-        }
-        return false;
-    }
-
-    void display_symbol(symbol const& s) {
-         if (s.is_numerical()) {
-             m_out << s << " ";
-         }
-         else if (!(s.bare_str()[0])) {
-              m_out << "\"null\" ";
-         }
-         else if (!has_special_char(s.bare_str())) {
-             m_out << s << " ";
-         }
-         else {
-             char const* r = s.bare_str();
-             m_out << "\"";
-             while (*r) {
-                 if (*r == ' ' || *r == '\n' ||
-                     *r == '\t' || *r == '\r') {
-                     m_out << "\\" << ((unsigned)(*r));
-                 }
-                 else {
-                     m_out << *r;
-                 }
-                 ++r;
-             }
-             m_out << "\" ";
-         }        
-    }
-
-public:
-
-    dag_printer(ast_manager& mgr, std::ostream & out, ast_mark& mark):
-        m_manager(mgr),
-        m_out(out),
-        m_mark(mark),
-        m_initialized(false),
-        m_basic_fid(mgr.get_basic_family_id()),
-        m_bv_fid(mgr.get_family_id("bv")),
-        m_arith_fid(mgr.get_family_id("arith")),
-        m_array_fid(mgr.get_family_id("array")),
-        m_arith(mgr),
-        m_bv(mgr),
-        m_enable_shortcut(true)
-    {
-    }
-
-    void operator()(sort * n) {
-        process_info(n->get_info());
-        m_out << "Ty ";
-        display_node_id(n);        
-        display_symbol(n->get_name());
-        display_info(n->get_info());
-        m_out << "\n";
-    }
-     
-    void pp_num(app* n, rational const& r) {
-        m_out << "Num ";
-        display_node_id(n);
-        m_out << r << " ";
-        display_node_id(m_manager.get_sort(n));
-        m_out << "\n";
-    }
-
-    void operator()(var * n) {
-        process_ast(n->get_sort());
-        m_out << "Var ";
-        display_node_id(n);        
-        m_out << n->get_idx() << " ";
-        display_node_id(n->get_sort());
-        m_out << "\n";
-    }
-
-    void operator()(func_decl * n) {
-
-        process_info(n->get_info());
-
-        family_id fid = n->get_family_id();
-        if (m_arith_fid == fid &&
-            n->get_info()->get_decl_kind() == OP_NUM) {
-            return;
-        }
-
-        if (m_bv_fid == fid &&
-            n->get_info()->get_decl_kind() == OP_BV_NUM) {
-            return;
-        }
-
-        m_out << "Dec ";
-        display_node_id(n);
-        display_symbol(n->get_name());
-        unsigned dom_size = n->get_arity();
-        for (unsigned i = 0; i < dom_size; ++i) {
-            display_node_id(n->get_domain(i));
-        }
-        display_node_id(n->get_range());
-        display_info(n->get_info());
-        m_out << "\n";
-    }
-
-    void operator()(app * n) {   
-        process_ast(n->get_decl());
-        family_id fid = n->get_family_id();
-        unsigned bv_size;
-        rational val;
-        if (m_arith.is_numeral(n, val)) {
-            pp_num(n, val);
-        }
-        else if (m_bv.is_numeral(n, val, bv_size)) {
-            pp_num(n, val);
-        }
-        else if (m_enable_shortcut &&
-            fid != null_family_id && 
-            (fid == m_basic_fid || fid == m_bv_fid || fid == m_array_fid || fid == m_arith_fid)) {
-            display_builtin(n);
-        }
-        else if (n->get_num_args() == 0 && fid == null_family_id) {
-            func_decl* d = n->get_decl();
-            m_out << "Const ";
-            display_node_id(n);
-            display_symbol(d->get_name());
-            display_node_id(d->get_range());
-            m_out << "\n";
-        }
-        else {
-            m_out << "Fun ";
-            display_node_id(n);
-            display_node_id(n->get_decl());
-            display_children(n->get_num_args(), n->get_args());
-            m_out << "\n";
-        }
-    }
-
-    void operator()(quantifier * n) {
-        m_out << "Qua ";
-        display_node_id(n);
-        m_out << (n->is_forall() ? "FORALL" : "EXISTS") << " ";
-        m_out << n->get_weight() << " ";
-        if (symbol::null != n->get_skid()) {            
-            m_out << "\"" << n->get_skid()   << "\" ";
-        }
-        else {
-            m_out << "\"null\" ";
-        }
-        if (symbol::null != n->get_qid()) {
-            m_out << "\"" << n->get_qid()    << "\" ";
-        }
-        else {
-            m_out << "\"null\" ";
-        }
-        unsigned num_decls = n->get_num_decls();
-        m_out << num_decls << " ";
-        for (unsigned i = 0; i < num_decls; i++) {
-            m_out << n->get_decl_name(i) << " ";
-            display_node_id(n->get_decl_sort(i));
-        }
-        m_out << n->get_num_patterns() << " ";
-        display_children(n->get_num_patterns(), n->get_patterns());
-        display_node_id(n->get_expr());
-        m_out << "\n";
-    }
-};
-
-void ast_dag_pp(std::ostream & out, ast_manager& mgr, ast_mark& mark, ast * n) {    
-    dag_printer p(mgr, out, mark);
-    for_each_ast(p, mark, n, true);
-}
-
-void ast_dag_pp(std::ostream & out, ast_manager& mgr, ast * n) {
-    ast_mark mark;
-    ast_dag_pp(out, mgr, mark, n);
-}
-

File src/ast/ast_dag_pp.h

-/*++
-Copyright (c) 2006 Microsoft Corporation
-
-Module Name:
-
-    ast_dag_pp.h
-
-Abstract:
-
-    AST low level pretty printer.
-
-Author:
-
-    Leonardo de Moura (leonardo) 2006-10-19.
-
-Revision History:
-
---*/
-#ifndef _AST_DAG_PP_H_
-#define _AST_DAG_PP_H_
-
-#include<iostream>
-
-class ast;
-
-void ast_dag_pp(std::ostream & out, ast_manager& mgr, ast * n);
-
-void ast_dag_pp(std::ostream & out, ast_manager& mgr, ast_mark& mark, ast * n);
-
-#endif /* _AST_DAG_PP_H_ */
-

File src/ast/expr2dot.cpp

-/*++
-Copyright (c) 2006 Microsoft Corporation
-
-Module Name:
-
-    expr2dot.cpp
-
-Abstract:
-
-    <abstract>
-
-Author:
-
-    Leonardo de Moura (leonardo) 2008-03-07.
-
-Revision History:
-
---*/
-#include"expr2dot.h"
-#include"for_each_expr.h"
-
-class dot_printer {
-    std::ostream & m_out;
-    ast_manager &  m_manager;
-    bool           m_proofs_only;
-public:
-    dot_printer(std::ostream & out, ast_manager & m, bool proofs):
-        m_out(out),
-        m_manager(m),
-        m_proofs_only(proofs) {
-    }
-
-    char const * get_color(app * n) {
-        if (m_manager.is_unit_resolution(n)) 
-            return "blue";
-        else if (m_manager.is_lemma(n))
-            return "gold";
-        else if (m_manager.is_transitivity(n))
-            return "red";
-        else if (m_manager.is_monotonicity(n))
-            return "green";
-        else 
-            return "black";
-    }
-
-    void operator()(var * n) {
-        if (!m_proofs_only) {
-            m_out << n->get_id() << "[label=\"\",shape=circle];\n";
-        }
-    }
-    
-    void operator()(quantifier * n) {
-        if (!m_proofs_only) {
-            m_out << n->get_id() << "[label=\"\",shape=circle,color=gray];\n";
-            m_out << n->get_expr()->get_id() << " -> " << n->get_id() << ";\n";
-        }
-    }
-
-    void operator()(app * n) {
-        if (!m_proofs_only || m_manager.is_proof(n)) {
-            char const * c = get_color(n);
-            m_out << n->get_id() << "[label=\"\",shape=circle,color=" << c << "];\n";
-            if (m_proofs_only) {
-                unsigned num = m_manager.get_num_parents(n);
-                for (unsigned i = 0; i < num; i++)
-                    m_out << m_manager.get_parent(n, i)->get_id() << " -> " << n->get_id() << " [color=" << c << "];\n";
-            }
-            else {
-                unsigned num = n->get_num_args();
-                for (unsigned i = 0; i < num; i++)
-                    m_out << n->get_arg(i)->get_id() << " -> " << n->get_id() << " [color=" << c << "];\n";
-            }
-        }
-    }
-};
-
-void expr2dot(std::ostream & out, expr * n, ast_manager & m, bool proofs) {
-    out << "digraph \"ast\" {\n";
-    dot_printer p(out, m, proofs);
-    for_each_expr(p, n);
-    out << "}\n";
-}
-
-

File src/ast/expr2dot.h

-/*++
-Copyright (c) 2006 Microsoft Corporation
-
-Module Name:
-
-    expr2dot.h
-
-Abstract:
-
-    Convert expressions into a .DOT file
-
-Author:
-
-    Leonardo de Moura (leonardo) 2008-03-07.
-
-Revision History:
-
---*/
-#ifndef _EXPR2DOT_H_
-#define _EXPR2DOT_H_
-
-#include"ast.h"
-
-void expr2dot(std::ostream & out, expr * a, ast_manager & m, bool proofs = false);
-
-#endif /* _AST2DOT_H_ */
-

File src/math/interval/interval.h

+/*++
+Copyright (c) 2012 Microsoft Corporation
+
+Module Name:
+
+    interval.h
+
+Abstract:
+
+    Goodies/Templates for interval arithmetic
+
+Author:
+
+    Leonardo de Moura (leonardo) 2012-07-19.
+
+Revision History:
+
+--*/
+#ifndef _INTERVAL_H_
+#define _INTERVAL_H_
+
+#include"mpq.h"
+#include"ext_numeral.h"
+
+/**
+   \brief Default configuration for interval manager.
+   It is used for documenting the required interface.
+*/
+class im_default_config {
+    unsynch_mpq_manager &       m_manager;
+public:
+    typedef unsynch_mpq_manager numeral_manager;
+    typedef mpq                 numeral;
+
+    // Every configuration object must provide an interval type.
+    // The actual fields are irrelevant, the interval manager
+    // accesses interval data using the following API.
+    struct interval {
+        numeral   m_lower;
+        numeral   m_upper;
+        unsigned  m_lower_open:1;
+        unsigned  m_upper_open:1;
+        unsigned  m_lower_inf:1;
+        unsigned  m_upper_inf:1;
+    };
+
+    // Should be NOOPs for precise numeral types.
+    // For imprecise types (e.g., floats) it should set the rounding mode.
+    void round_to_minus_inf() {}
+    void round_to_plus_inf() {}
+    void set_rounding(bool to_plus_inf) {}
+
+    // Getters
+    numeral const & lower(interval const & a) const { return a.m_lower; }
+    numeral const & upper(interval const & a) const { return a.m_upper; }
+    numeral & lower(interval & a) { return a.m_lower; }
+    numeral & upper(interval & a) { return a.m_upper; }
+    bool lower_is_open(interval const & a) const { return a.m_lower_open; }
+    bool upper_is_open(interval const & a) const { return a.m_upper_open; }
+    bool lower_is_inf(interval const & a) const { return a.m_lower_inf; }
+    bool upper_is_inf(interval const & a) const { return a.m_upper_inf; }
+
+    // Setters
+    void set_lower(interval & a, numeral const & n) { m_manager.set(a.m_lower, n); }
+    void set_upper(interval & a, numeral const & n) { m_manager.set(a.m_upper, n); }
+    void set_lower_is_open(interval & a, bool v) { a.m_lower_open = v; }
+    void set_upper_is_open(interval & a, bool v) { a.m_upper_open = v; }
+    void set_lower_is_inf(interval & a, bool v) { a.m_lower_inf = v; }
+    void set_upper_is_inf(interval & a, bool v) { a.m_upper_inf = v; }
+    
+    // Reference to numeral manager
+    numeral_manager & m() const { return m_manager; }
+
+    im_default_config(numeral_manager & m):m_manager(m) {}
+};
+
+#define DEP_IN_LOWER1 1
+#define DEP_IN_UPPER1 2
+#define DEP_IN_LOWER2 4
+#define DEP_IN_UPPER2 8
+
+typedef short bound_deps;
+inline bool dep_in_lower1(bound_deps d) { return (d & DEP_IN_LOWER1) != 0; }
+inline bool dep_in_lower2(bound_deps d) { return (d & DEP_IN_LOWER2) != 0; }
+inline bool dep_in_upper1(bound_deps d) { return (d & DEP_IN_UPPER1) != 0; }
+inline bool dep_in_upper2(bound_deps d) { return (d & DEP_IN_UPPER2) != 0; }
+inline bound_deps dep1_to_dep2(bound_deps d) { 
+    SASSERT(!dep_in_lower2(d) && !dep_in_upper2(d));
+    bound_deps r = d << 2; 
+    SASSERT(dep_in_lower1(d) == dep_in_lower2(r));
+    SASSERT(dep_in_upper1(d) == dep_in_upper2(r));
+    SASSERT(!dep_in_lower1(r) && !dep_in_upper1(r));
+    return r;
+}
+
+/**
+   \brief Interval dependencies for unary and binary operations on intervals.
+   It contains the dependencies for the output lower and upper bounds 
+   for the resultant interval.
+*/ 
+struct interval_deps {
+    bound_deps m_lower_deps;
+    bound_deps m_upper_deps;
+};
+
+template<typename C>
+class interval_manager {
+    typedef typename C::numeral_manager       numeral_manager;
+    typedef typename numeral_manager::numeral numeral;
+    typedef typename C::interval              interval;
+    
+    C          m_c;
+    numeral    m_result_lower;
+    numeral    m_result_upper;
+    numeral    m_mul_ad;
+    numeral    m_mul_bc;
+    numeral    m_mul_ac;
+    numeral    m_mul_bd;
+    numeral    m_one;
+    numeral    m_minus_one;
+    numeral    m_inv_k;
+
+    unsigned   m_pi_n;
+    interval   m_pi_div_2;
+    interval   m_pi;
+    interval   m_3_pi_div_2;
+    interval   m_2_pi;
+
+    volatile bool m_cancel;     
+    
+    void round_to_minus_inf() { m_c.round_to_minus_inf(); }
+    void round_to_plus_inf() { m_c.round_to_plus_inf(); }
+    void set_rounding(bool to_plus_inf) { m_c.set_rounding(to_plus_inf); }
+
+    ext_numeral_kind lower_kind(interval const & a) const { return m_c.lower_is_inf(a) ? EN_MINUS_INFINITY : EN_NUMERAL; }
+    ext_numeral_kind upper_kind(interval const & a) const { return m_c.upper_is_inf(a) ? EN_PLUS_INFINITY : EN_NUMERAL; }
+
+    void set_lower(interval & a, numeral const & n) { m_c.set_lower(a, n); }
+    void set_upper(interval & a, numeral const & n) { m_c.set_upper(a, n); }
+    void set_lower_is_open(interval & a, bool v) { m_c.set_lower_is_open(a, v); }
+    void set_upper_is_open(interval & a, bool v) { m_c.set_upper_is_open(a, v); }
+    void set_lower_is_inf(interval & a, bool v) { m_c.set_lower_is_inf(a, v); }
+    void set_upper_is_inf(interval & a, bool v) { m_c.set_upper_is_inf(a, v); }
+
+    void nth_root_slow(numeral const & a, unsigned n, numeral const & p, numeral & lo, numeral & hi);
+    void A_div_x_n(numeral const & A, numeral const & x, unsigned n, bool to_plus_inf, numeral & r);
+    void rough_approx_nth_root(numeral const & a, unsigned n, numeral & o);
+    void approx_nth_root(numeral const & a, unsigned n, numeral const & p, numeral & o);
+    void nth_root_pos(numeral const & A, unsigned n, numeral const & p, numeral & lo, numeral & hi);
+    void nth_root(numeral const & a, unsigned n, numeral const & p, numeral & lo, numeral & hi);
+    
+    void pi_series(int x, numeral & r, bool to_plus_inf);
+    void fact(unsigned n, numeral & o);
+    void sine_series(numeral const & a, unsigned k, bool upper, numeral & o);
+    void cosine_series(numeral const & a, unsigned k, bool upper, numeral & o);
+    void e_series(unsigned k, bool upper, numeral & o);
+
+    void div_mul(numeral const & k, interval const & a, interval & b, bool inv_k);
+
+    void checkpoint();
+
+public:    
+    interval_manager(C const & c);
+    ~interval_manager();
+
+    void set_cancel(bool f) { m_cancel = f; }
+
+    numeral_manager & m() const { return m_c.m(); }    
+
+    void del(interval & a);
+
+    numeral const & lower(interval const & a) const { return m_c.lower(a); }
+    numeral const & upper(interval const & a) const { return m_c.upper(a); }
+    numeral & lower(interval & a) { return m_c.lower(a); }
+    numeral & upper(interval & a) { return m_c.upper(a); }
+    bool lower_is_open(interval const & a) const { return m_c.lower_is_open(a); }
+    bool upper_is_open(interval const & a) const { return m_c.upper_is_open(a); }
+    bool lower_is_inf(interval const & a) const { return m_c.lower_is_inf(a); }
+    bool upper_is_inf(interval const & a) const { return m_c.upper_is_inf(a); }
+    
+    bool lower_is_neg(interval const & a) const { return ::is_neg(m(), lower(a), lower_kind(a)); }
+    bool lower_is_pos(interval const & a) const { return ::is_pos(m(), lower(a), lower_kind(a)); }
+    bool lower_is_zero(interval const & a) const { return ::is_zero(m(), lower(a), lower_kind(a)); }
+
+    bool upper_is_neg(interval const & a) const { return ::is_neg(m(), upper(a), upper_kind(a)); }
+    bool upper_is_pos(interval const & a) const { return ::is_pos(m(), upper(a), upper_kind(a)); }
+    bool upper_is_zero(interval const & a) const { return ::is_zero(m(), upper(a), upper_kind(a)); }
+
+    bool is_P(interval const & n) const { return lower_is_pos(n) || lower_is_zero(n); }
+    bool is_P0(interval const & n) const { return lower_is_zero(n) && !lower_is_open(n); }
+    bool is_P1(interval const & n) const { return lower_is_pos(n) || (lower_is_zero(n) && lower_is_open(n)); }
+    bool is_N(interval const & n) const { return upper_is_neg(n) || upper_is_zero(n); }
+    bool is_N0(interval const & n) const { return upper_is_zero(n) && !upper_is_open(n); }
+    bool is_N1(interval const & n) const { return upper_is_neg(n) || (upper_is_zero(n) && upper_is_open(n)); }
+    bool is_M(interval const & n) const { return lower_is_neg(n) && upper_is_pos(n); }
+    bool is_zero(interval const & n) const { return lower_is_zero(n) && upper_is_zero(n); }
+
+    void set(interval & t, interval const & s);
+
+    bool eq(interval const & a, interval const & b) const;
+
+    /**
+       \brief Set lower bound to -oo.
+    */
+    void reset_lower(interval & a);
+
+    /**
+       \brief Set upper bound to +oo.
+    */
+    void reset_upper(interval & a);
+
+    /**
+       \brief Set interval to (-oo, oo)
+    */
+    void reset(interval & a);
+
+    /**
+       \brief Return true if the given interval contains 0.
+    */
+    bool contains_zero(interval const & n) const;
+    
+    /**
+       \brief Return true if n contains v.
+    */
+    bool contains(interval const & n, numeral const & v) const;
+    
+    void display(std::ostream & out, interval const & n) const;
+
+    bool check_invariant(interval const & n) const;
+
+    /**
+       \brief b <- -a
+    */
+    void neg(interval const & a, interval & b, interval_deps & b_deps);
+    void neg(interval const & a, interval & b);
+    void neg_jst(interval const & a, interval_deps & b_deps);
+
+    /**
+       \brief c <- a + b
+    */
+    void add(interval const & a, interval const & b, interval & c, interval_deps & c_deps);
+    void add(interval const & a, interval const & b, interval & c);
+    void add_jst(interval const & a, interval const & b, interval_deps & c_deps);
+
+    /**
+       \brief c <- a - b
+    */
+    void sub(interval const & a, interval const & b, interval & c, interval_deps & c_deps);
+    void sub(interval const & a, interval const & b, interval & c);
+    void sub_jst(interval const & a, interval const & b, interval_deps & c_deps);
+
+    /**
+       \brief b <- k * a
+    */
+    void mul(numeral const & k, interval const & a, interval & b, interval_deps & b_deps);
+    void mul(numeral const & k, interval const & a, interval & b) { div_mul(k, a, b, false); }
+    void mul_jst(numeral const & k, interval const & a, interval_deps & b_deps);
+    /**
+       \brief b <- (n/d) * a
+    */
+    void mul(int n, int d, interval const & a, interval & b);
+
+    /**
+       \brief b <- a/k
+       
+       \remark For imprecise numerals, this is not equivalent to 
+       m().inv(k)
+       mul(k, a, b)
+    
+       That is, we must invert k rounding towards +oo or -oo depending whether we
+       are computing a lower or upper bound.
+    */
+    void div(interval const & a, numeral const & k, interval & b, interval_deps & b_deps);
+    void div(interval const & a, numeral const & k, interval & b) { div_mul(k, a, b, true); }
+    void div_jst(interval const & a, numeral const & k, interval_deps & b_deps) { mul_jst(k, a, b_deps); }
+
+    /**
+       \brief c <- a * b
+    */
+    void mul(interval const & a, interval const & b, interval & c, interval_deps & c_deps);
+    void mul(interval const & a, interval const & b, interval & c);
+    void mul_jst(interval const & a, interval const & b, interval_deps & c_deps);
+
+    /**
+       \brief b <- a^n
+    */
+    void power(interval const & a, unsigned n, interval & b, interval_deps & b_deps);
+    void power(interval const & a, unsigned n, interval & b);
+    void power_jst(interval const & a, unsigned n, interval_deps & b_deps);
+
+    /**
+       \brief b <- a^(1/n) with precision p.
+
+       \pre if n is even, then a must not contain negative numbers.
+    */
+    void nth_root(interval const & a, unsigned n, numeral const & p, interval & b, interval_deps & b_deps);
+    void nth_root(interval const & a, unsigned n, numeral const & p, interval & b);
+    void nth_root_jst(interval const & a, unsigned n, numeral const & p, interval_deps & b_deps);
+    
+    /**
+       \brief Given an equation x^n = y and an interval for y, compute the solution set for x with precision p.
+       
+       \pre if n is even, then !lower_is_neg(y)
+    */
+    void xn_eq_y(interval const & y, unsigned n, numeral const & p, interval & x, interval_deps & x_deps);
+    void xn_eq_y(interval const & y, unsigned n, numeral const & p, interval & x); 
+    void xn_eq_y_jst(interval const & y, unsigned n, numeral const & p, interval_deps & x_deps);
+   
+    /**
+       \brief b <- 1/a
+       \pre !contains_zero(a)
+    */
+    void inv(interval const & a, interval & b, interval_deps & b_deps);
+    void inv(interval const & a, interval & b);
+    void inv_jst(interval const & a, interval_deps & b_deps);
+
+    /**
+       \brief c <- a/b
+       \pre !contains_zero(b)
+       \pre &a == &c (that is, c should not be an alias for a)
+    */
+    void div(interval const & a, interval const & b, interval & c, interval_deps & c_deps);
+    void div(interval const & a, interval const & b, interval & c);
+    void div_jst(interval const & a, interval const & b, interval_deps & c_deps);
+    
+    /**
+       \brief Store in r an interval that contains the number pi.
+       The size of the interval is (1/15)*(1/16^n)
+    */
+    void pi(unsigned n, interval & r);
+
+    /**
+       \brief Set the precision of the internal interval representing pi.
+    */
+    void set_pi_prec(unsigned n);
+
+    /**
+       \brief Set the precision of the internal interval representing pi to a precision of at least n.
+    */
+    void set_pi_at_least_prec(unsigned n);
+
+    void sine(numeral const & a, unsigned k, numeral & lo, numeral & hi);
+
+    void cosine(numeral const & a, unsigned k, numeral & lo, numeral & hi);
+
+    /**
+       \brief Store in r the Euler's constant e.
+       The size of the interval is 4/(k+1)!
+    */
+    void e(unsigned k, interval & r);
+};
+
+#endif

File src/math/interval/interval_def.h

+/*++
+Copyright (c) 2012 Microsoft Corporation
+
+Module Name:
+
+    interval_def.h
+
+Abstract:
+
+    Goodies/Templates for interval arithmetic
+
+Author:
+
+    Leonardo de Moura (leonardo) 2012-07-19.
+
+Revision History:
+
+--*/
+#ifndef _INTERVAL_DEF_H_
+#define _INTERVAL_DEF_H_
+
+#include"interval.h"
+#include"debug.h"
+#include"trace.h"
+#include"scoped_numeral.h"
+#include"cooperate.h"
+
+#define DEFAULT_PI_PRECISION 2
+
+// #define TRACE_NTH_ROOT
+
+template<typename C>
+interval_manager<C>::interval_manager(C const & c):m_c(c) { 
+    m().set(m_minus_one, -1); 
+    m().set(m_one, 1);
+    m_pi_n = 0;
+    m_cancel = false;
+}
+
+template<typename C>
+interval_manager<C>::~interval_manager() {
+    del(m_pi_div_2);
+    del(m_pi);
+    del(m_3_pi_div_2);
+    del(m_2_pi);
+    m().del(m_result_lower);
+    m().del(m_result_upper);
+    m().del(m_mul_ad);
+    m().del(m_mul_bc);
+    m().del(m_mul_ac);
+    m().del(m_mul_bd);
+    m().del(m_minus_one);
+    m().del(m_one);
+    m().del(m_inv_k);
+}
+
+template<typename C>
+void interval_manager<C>::del(interval & a) {
+    m().del(lower(a));
+    m().del(upper(a));
+}
+
+
+template<typename C>
+void interval_manager<C>::checkpoint() {
+    if (m_cancel)
+        throw default_exception("canceled");
+    cooperate("interval");
+}
+
+/*
+  Compute the n-th root of a with precision p. The result hi - lo <= p
+  lo and hi are lower/upper bounds for the value of the n-th root of a.
+  That is, the n-th root is in the interval [lo, hi]
+
+  If n is even, then a is assumed to be nonnegative.
+
+  If numeral_manager is not precise, the procedure does not guarantee the precision p.
+*/
+template<typename C>
+void interval_manager<C>::nth_root_slow(numeral const & a, unsigned n, numeral const & p, numeral & lo, numeral & hi) {
+#ifdef TRACE_NTH_ROOT
+    static unsigned counter = 0;
+    static unsigned loop_counter = 0;
+    counter++;
+    if (counter % 1000 == 0) 
+        std::cerr << "[nth-root] " << counter << " " << loop_counter << " " << ((double)loop_counter)/((double)counter) << std::endl;
+#endif
+
+    bool n_is_even = (n % 2 == 0);
+    SASSERT(!n_is_even || m().is_nonneg(a));
+    if (m().is_zero(a) || m().is_one(a) || (!n_is_even && m().eq(a, m_minus_one))) {
+        m().set(lo, a);
+        m().set(hi, a);
+        return;
+    }
+    if (m().lt(a, m_minus_one)) {
+        m().set(lo, a);
+        m().set(hi, -1);
+    }
+    else if (m().is_neg(a)) {
+        m().set(lo, -1);
+        m().set(hi, 0);
+    }
+    else if (m().lt(a, m_one)) {
+        m().set(lo, 0);
+        m().set(hi, 1);
+    }
+    else {
+        m().set(lo, 1);
+        m().set(hi, a);
+    }
+    SASSERT(m().le(lo, hi));
+    _scoped_numeral<numeral_manager> c(m()), cn(m());
+    _scoped_numeral<numeral_manager> two(m());
+    m().set(two, 2);
+    while (true) {
+        checkpoint();
+#ifdef TRACE_NTH_ROOT
+        loop_counter++;
+#endif
+        m().add(hi, lo, c);
+        m().div(c, two, c);
+        if (m().precise()) {
+            m().power(c, n, cn);
+            if (m().gt(cn, a)) {
+                m().set(hi, c);
+            }
+            else if (m().eq(cn, a)) {
+                // root is precise
+                m().set(lo, c);
+                m().set(hi, c);
+                return;
+            }
+            else {
+                m().set(lo, c);
+            }
+        }
+        else {
+            round_to_minus_inf();
+            m().power(c, n, cn);
+            if (m().gt(cn, a)) {
+                m().set(hi, c);
+            }
+            else {
+                round_to_plus_inf();
+                m().power(c, n, cn);
+                if (m().lt(cn, a)) {
+                    m().set(lo, c);
+                }
+                else {
+                    // can't improve, numeral_manager is not precise enough,
+                    // a is between round-to-minus-inf(c^n) and round-to-plus-inf(c^n)
+                    return;
+                }
+            }
+        }
+        round_to_plus_inf();
+        m().sub(hi, lo, c);
+        if (m().le(c, p))
+            return; // result is precise enough
+    }
+}
+
+/**
+   \brief Store in o a rough approximation of a^1/n.
+   
+   It uses 2^Floor[Floor(Log2(a))/n]
+
+   \pre is_pos(a)
+*/
+template<typename C>
+void interval_manager<C>::rough_approx_nth_root(numeral const & a, unsigned n, numeral & o) {
+    SASSERT(m().is_pos(a));
+    SASSERT(n > 0);
+    round_to_minus_inf();
+    unsigned k = m().prev_power_of_two(a);
+    m().set(o, 2);
+    m().power(o, k/n, o);
+}
+
+/*
+  Compute the n-th root of \c a with (suggested) precision p. 
+  The only guarantee provided by this method is that a^(1/n) is in [lo, hi].
+
+  If n is even, then a is assumed to be nonnegative.
+*/
+template<typename C>
+void interval_manager<C>::nth_root(numeral const & a, unsigned n, numeral const & p, numeral & lo, numeral & hi) {
+    // nth_root_slow(a, n, p, lo, hi);
+    // return;
+    SASSERT(n > 0);
+    SASSERT(n % 2 != 0 || m().is_nonneg(a));
+    if (n == 1 || m().is_zero(a) || m().is_one(a) || m().is_minus_one(a)) {
+        // easy cases: 1, -1, 0
+        m().set(lo, a);
+        m().set(hi, a);
+        return;
+    }
+    bool is_neg = m().is_neg(a);
+    _scoped_numeral<numeral_manager> A(m());
+    m().set(A, a);
+    m().abs(A);
+
+    nth_root_pos(A, n, p, lo, hi);
+    STRACE("nth_root_trace", 
+           tout << "[nth-root] ("; m().display(tout, A); tout << ")^(1/" << n << ") >= "; m().display(tout, lo); tout << "\n"; 
+           tout << "[nth-root] ("; m().display(tout, A); tout << ")^(1/" << n << ") <= "; m().display(tout, hi); tout << "\n";);
+    if (is_neg) {
+        m().swap(lo, hi);
+        m().neg(lo);
+        m().neg(hi);
+    }
+}
+
+/**
+   r <- A/(x^n)
+   
+   If to_plus_inf,      then r >= A/(x^n)
+   If not to_plus_inf,  then r <= A/(x^n)
+
+*/
+template<typename C>
+void interval_manager<C>::A_div_x_n(numeral const & A, numeral const & x, unsigned n, bool to_plus_inf, numeral & r) {
+    if (n == 1) {
+        if (m().precise()) {
+            m().div(A, x, r);
+        }
+        else {
+            set_rounding(to_plus_inf);
+            m().div(A, x, r);
+        }
+    }
+    else {
+        if (m().precise()) {
+            m().power(x, n, r);
+            m().div(A, r, r);
+        }
+        else {
+            set_rounding(!to_plus_inf);
+            m().power(x, n, r);
+            set_rounding(to_plus_inf);
+            m().div(A, r, r);
+        }
+    }
+}
+
+/**
+   \brief Compute an approximation of A^(1/n) using the sequence
+
+   x' = 1/n((n-1)*x + A/(x^(n-1)))
+
+   The computation stops when the difference between current and new x is less than p.
+   The procedure may not terminate if m() is not precise and p is very small.
+   
+*/
+template<typename C>
+void interval_manager<C>::approx_nth_root(numeral const & A, unsigned n, numeral const & p, numeral & x) {
+    SASSERT(m().is_pos(A));
+    SASSERT(n > 1);
+#ifdef TRACE_NTH_ROOT
+    static unsigned counter = 0;
+    static unsigned loop_counter = 0;
+    counter++;
+    if (counter % 1000 == 0) 
+        std::cerr << "[nth-root] " << counter << " " << loop_counter << " " << ((double)loop_counter)/((double)counter) << std::endl;
+#endif
+    
+    _scoped_numeral<numeral_manager> x_prime(m()), d(m());
+    
+    m().set(d, 1);
+    if (m().lt(A, d)) 
+        m().set(x, A);
+    else
+        rough_approx_nth_root(A, n, x);
+    
+    round_to_minus_inf();
+
+    if (n == 2) {
+        _scoped_numeral<numeral_manager> two(m());
+        m().set(two, 2);
+        while (true) {
+            checkpoint();
+#ifdef TRACE_NTH_ROOT
+            loop_counter++;
+#endif
+            m().div(A, x, x_prime);
+            m().add(x, x_prime, x_prime);
+            m().div(x_prime, two, x_prime);
+            m().sub(x_prime, x, d);
+            m().abs(d);
+            m().swap(x, x_prime);
+            if (m().lt(d, p))
+                return;
+        }
+    }
+    else {
+        _scoped_numeral<numeral_manager> _n(m()), _n_1(m());
+        m().set(_n, n);   // _n contains n
+        m().set(_n_1, n);
+        m().dec(_n_1);    // _n_1 contains n-1 
+        
+        while (true) {
+            checkpoint();
+#ifdef TRACE_NTH_ROOT
+            loop_counter++;
+#endif
+            m().power(x, n-1, x_prime);
+            m().div(A, x_prime, x_prime);
+            m().mul(_n_1, x, d);
+            m().add(d, x_prime, x_prime);
+            m().div(x_prime, _n, x_prime);
+            m().sub(x_prime, x, d);
+            m().abs(d);
+            TRACE("nth_root", 
+                  tout << "A:       "; m().display(tout, A); tout << "\n";
+                  tout << "x:       "; m().display(tout, x); tout << "\n";
+                  tout << "x_prime: "; m().display(tout, x_prime); tout << "\n";
+                  tout << "d:       "; m().display(tout, d); tout << "\n";
+                  );
+            m().swap(x, x_prime);
+            if (m().lt(d, p))
+                return;
+        }
+    }
+}
+
+template<typename C>
+void interval_manager<C>::nth_root_pos(numeral const & A, unsigned n, numeral const & p, numeral & lo, numeral & hi) {
+    approx_nth_root(A, n, p, hi);
+    if (m().precise()) {
+        // Assuming hi has a upper bound for A^(n-1)
+        // Then, A/(x^(n-1)) must be lower bound
+        A_div_x_n(A, hi, n-1, false, lo);
+        // Check if we were wrong
+        if (m().lt(hi, lo)) {
+            // swap if wrong
+            m().swap(lo, hi);
+        }
+    }
+    else {
+        // Check if hi is really a upper bound for A^(n-1)
+        A_div_x_n(A, hi, n-1, true /* lo will be greater than the actual lower bound */, lo);
+        TRACE("nth_root_bug", 
+              tout << "Assuming upper\n";
+              tout << "A: "; m().display(tout, A); tout << "\n";
+              tout << "hi: "; m().display(tout, hi); tout << "\n";
+              tout << "lo: "; m().display(tout, hi); tout << "\n";);
+        if (m().le(lo, hi)) {
+            // hi is really the upper bound
+            // Must compute lo again but approximating to -oo
+            A_div_x_n(A, hi, n-1, false, lo);
+        }
+        else {
+            // hi should be lower bound
+            m().swap(lo, hi);
+            // check if lo is lower bound
+            A_div_x_n(A, lo, n-1, false /* hi will less than the actual upper bound */, hi);
+            if (m().le(lo, hi)) {
+                // lo is really the lower bound
+                // Must compute hi again but approximating to +oo
+                A_div_x_n(A, lo, n-1, true, hi);
+            }
+            else {
+                // we don't have anything due to rounding errors
+                // Be supper conservative
+                // This should not really happen very often.
+                _scoped_numeral<numeral_manager> one(m());
+                if (m().lt(A, one)) {
+                    m().set(lo, 0);
+                    m().set(hi, 1);
+                }
+                else {
+                    m().set(lo, 1);
+                    m().set(hi, A);
+                }
+            }
+        }
+    }
+}
+
+
+/**
+   \brief o <- n!
+*/
+template<typename C>
+void interval_manager<C>::fact(unsigned n, numeral & o) {
+    _scoped_numeral<numeral_manager> aux(m());
+    m().set(o, 1);
+    for (unsigned i = 2; i <= n; i++) {
+        m().set(aux, static_cast<int>(i));
+        m().mul(aux, o, o);
+        TRACE("fact_bug", tout << "i: " << i << ", o: " << m().to_rational_string(o) << "\n";);
+    }
+}
+
+template<typename C>
+void interval_manager<C>::sine_series(numeral const & a, unsigned k, bool upper, numeral & o) {
+    SASSERT(k % 2 == 1);
+    // Compute sine using taylor series up to k
+    // x - x^3/3! + x^5/5! - x^7/7! + ... 
+    // The result should be greater than or equal to the actual value if upper == true
+    // Otherwise it must be less than or equal to the actual value.
+    // The argument upper only matter if the numeral_manager is not precise.
+
+    // Taylor series up to k with rounding to 
+    _scoped_numeral<numeral_manager> f(m());
+    _scoped_numeral<numeral_manager> aux(m());
+    m().set(o, a);
+    bool sign         = true;
+    bool upper_factor = !upper; // since the first sign is negative, we must minimize factor to maximize result
+    for (unsigned i = 3; i <= k; i+=2) {
+        TRACE("sine_bug", tout << "[begin-loop] o: " << m().to_rational_string(o) << "\ni: " << i << "\n";
+              tout << "upper: " << upper << ", upper_factor: " << upper_factor << "\n";
+              tout << "o (default): " << m().to_string(o) << "\n";);
+        set_rounding(upper_factor); 
+        m().power(a, i, f);
+        TRACE("sine_bug", tout << "a^i " << m().to_rational_string(f) << "\n";);
+        set_rounding(!upper_factor);
+        fact(i, aux);
+        TRACE("sine_bug", tout << "i! " << m().to_rational_string(aux) << "\n";);
+        set_rounding(upper_factor);
+        m().div(f, aux, f);
+        TRACE("sine_bug", tout << "a^i/i! " << m().to_rational_string(f) << "\n";);
+        set_rounding(upper);
+        if (sign)
+            m().sub(o, f, o);
+        else
+            m().add(o, f, o);
+        TRACE("sine_bug", tout << "o: " << m().to_rational_string(o) << "\n";);
+        sign         = !sign;
+        upper_factor = !upper_factor;
+    }
+}
+
+template<typename C>
+void interval_manager<C>::sine(numeral const & a, unsigned k, numeral & lo, numeral & hi) {
+    TRACE("sine", tout << "sine(a), a: " << m().to_rational_string(a) << "\na: " << m().to_string(a) << "\n";);
+    SASSERT(&lo != &hi);
+    if (m().is_zero(a)) {
+        m().reset(lo);
+        m().reset(hi);
+        return;
+    }
+    
+    // Compute sine using taylor series
+    // x - x^3/3! + x^5/5! - x^7/7! + ...
+    //
+    // Note that, the coefficient of even terms is 0.
+    // So, we force k to be odd to make sure the error is minimized.
+    if (k % 2 == 0)
+        k++; 
+    
+    // Taylor series error = |x|^(k+1)/(k+1)!
+    _scoped_numeral<numeral_manager> error(m());
+    _scoped_numeral<numeral_manager> aux(m());
+    round_to_plus_inf();
+    m().set(error, a);
+    if (m().is_neg(error))
+        m().neg(error);
+    m().power(error, k+1, error);
+    TRACE("sine", tout << "a^(k+1): " << m().to_rational_string(error) << "\nk : " << k << "\n";);
+    round_to_minus_inf();
+    fact(k+1, aux);
+    TRACE("sine", tout << "(k+1)!: " << m().to_rational_string(aux) << "\n";);
+    round_to_plus_inf();
+    m().div(error, aux, error);
+    TRACE("sine", tout << "error: " << m().to_rational_string(error) << "\n";);
+
+    // Taylor series up to k with rounding to -oo 
+    sine_series(a, k, false, lo);
+
+    if (m().precise()) {
+        m().set(hi, lo);
+        m().sub(lo, error, lo);
+        if (m().lt(lo, m_minus_one)) {
+            m().set(lo, -1);
+            m().set(hi, 1);
+        }
+        else {
+            m().add(hi, error, hi);
+        }
+    }
+    else {
+        // We must recompute the series with rounding to +oo
+        TRACE("sine", tout << "lo before -error: " << m().to_rational_string(lo) << "\n";);
+        round_to_minus_inf();
+        m().sub(lo, error, lo);
+        TRACE("sine", tout << "lo: " << m().to_rational_string(lo) << "\n";);
+        if (m().lt(lo, m_minus_one)) {
+            m().set(lo, -1);
+            m().set(hi, 1);
+            return;
+        }
+        sine_series(a, k, true, hi);
+        round_to_plus_inf();
+        m().add(hi, error, hi);
+        TRACE("sine", tout << "hi: " << m().to_rational_string(hi) << "\n";);
+    }
+}
+
+template<typename C>
+void interval_manager<C>::cosine_series(numeral const & a, unsigned k, bool upper, numeral & o) {
+    SASSERT(k % 2 == 0);
+    // Compute cosine using taylor series up to k
+    // 1 - x^2/2! + x^4/4! - x^6/6! + ...
+    // The result should be greater than or equal to the actual value if upper == true
+    // Otherwise it must be less than or equal to the actual value.
+    // The argument upper only matter if the numeral_manager is not precise.
+
+
+    // Taylor series up to k with rounding to -oo 
+    _scoped_numeral<numeral_manager> f(m());
+    _scoped_numeral<numeral_manager> aux(m());
+    m().set(o, 1);
+    bool sign         = true;
+    bool upper_factor = !upper; // since the first sign is negative, we must minimize factor to maximize result
+    for (unsigned i = 2; i <= k; i+=2) {
+        set_rounding(upper_factor);
+        m().power(a, i, f);
+        set_rounding(!upper_factor);
+        fact(i, aux);
+        set_rounding(upper_factor);
+        m().div(f, aux, f);
+        set_rounding(upper);
+        if (sign)
+            m().sub(o, f, o);
+        else
+            m().add(o, f, o);
+        sign         = !sign;
+        upper_factor = !upper_factor; 
+    }
+}
+
+template<typename C>
+void interval_manager<C>::cosine(numeral const & a, unsigned k, numeral & lo, numeral & hi) {
+    TRACE("cosine", tout << "cosine(a): "; m().display_decimal(tout, a, 32); tout << "\n";);
+    SASSERT(&lo != &hi);
+    if (m().is_zero(a)) {
+        m().set(lo, 1);
+        m().set(hi, 1);
+        return;
+    }
+    
+    // Compute cosine using taylor series
+    // 1 - x^2/2! + x^4/4! - x^6/6! + ...
+    //
+    // Note that, the coefficient of odd terms is 0.
+    // So, we force k to be even to make sure the error is minimized.
+    if (k % 2 == 1)
+        k++; 
+    
+    // Taylor series error = |x|^(k+1)/(k+1)!
+    _scoped_numeral<numeral_manager> error(m());
+    _scoped_numeral<numeral_manager> aux(m());
+    round_to_plus_inf();
+    m().set(error, a);
+    if (m().is_neg(error))
+        m().neg(error);
+    m().power(error, k+1, error);
+    round_to_minus_inf();
+    fact(k+1, aux);
+    round_to_plus_inf();
+    m().div(error, aux, error);
+    TRACE("sine", tout << "error: "; m().display_decimal(tout, error, 32); tout << "\n";);
+
+    // Taylor series up to k with rounding to -oo 
+    cosine_series(a, k, false, lo);
+    
+    if (m().precise()) {
+        m().set(hi, lo);
+        m().sub(lo, error, lo);
+        if (m().lt(lo, m_minus_one)) {
+            m().set(lo, -1);
+            m().set(hi, 1);
+        }
+        else {
+            m().add(hi, error, hi);
+        }
+    }
+    else {
+        // We must recompute the series with rounding to +oo
+        round_to_minus_inf();
+        m().sub(lo, error, lo);
+        if (m().lt(lo, m_minus_one)) {
+            m().set(lo, -1);
+            m().set(hi, 1);
+            return;
+        }
+        cosine_series(a, k, true, hi);
+        round_to_plus_inf();
+        m().add(hi, error, hi);
+    }
+}
+
+template<typename C>
+void interval_manager<C>::reset_lower(interval & a) {
+    m().reset(lower(a));
+    set_lower_is_open(a, true);
+    set_lower_is_inf(a, true);
+}
+
+template<typename C>
+void interval_manager<C>::reset_upper(interval & a) {
+    m().reset(upper(a));
+    set_upper_is_open(a, true);
+    set_upper_is_inf(a, true);
+}
+
+template<typename C>
+void interval_manager<C>::reset(interval & a) {
+    reset_lower(a);
+    reset_upper(a);
+}
+
+template<typename C>
+bool interval_manager<C>::contains_zero(interval const & n) const {
+    return 
+        (lower_is_neg(n) || (lower_is_zero(n) && !lower_is_open(n))) &&
+        (upper_is_pos(n) || (upper_is_zero(n) && !upper_is_open(n)));
+}
+
+    
+template<typename C>
+bool interval_manager<C>::contains(interval const & n, numeral const & v) const {
+    if (!lower_is_inf(n)) {
+        if (m().lt(v, lower(n))) return false;
+        if (m().eq(v, lower(n)) && lower_is_open(n)) return false;
+    }
+    if (!upper_is_inf(n)) {
+        if (m().gt(v, upper(n))) return false;
+        if (m().eq(v, upper(n)) && upper_is_open(n)) return false;
+    }
+    return true;
+}
+
+template<typename C>
+void interval_manager<C>::display(std::ostream & out, interval const & n) const {
+    out << (lower_is_open(n) ? "(" : "[");
+    ::display(out, m(), lower(n), lower_kind(n));
+    out << ", ";
+    ::display(out, m(), upper(n), upper_kind(n));
+    out << (upper_is_open(n) ? ")" : "]");
+}
+
+template<typename C>
+bool interval_manager<C>::check_invariant(interval const & n) const {
+    if (::eq(m(), lower(n), lower_kind(n), upper(n), upper_kind(n))) {
+        SASSERT(!lower_is_open(n));
+        SASSERT(!upper_is_open(n));
+    }
+    else {
+        SASSERT(lt(m(), lower(n), lower_kind(n), upper(n), upper_kind(n)));
+    }
+    return true;
+}
+
+template<typename C>
+void interval_manager<C>::set(interval & t, interval const & s) {
+    if (&t == &const_cast<interval&>(s))
+        return;
+    if (lower_is_inf(s)) {
+        set_lower_is_inf(t, true);
+    }
+    else {
+        m().set(lower(t), lower(s));
+        set_lower_is_inf(t, false);
+    }
+    if (upper_is_inf(s)) {
+        set_upper_is_inf(t, true);
+    }
+    else {
+        m().set(upper(t), upper(s));
+        set_upper_is_inf(t, false);
+    }
+    set_lower_is_open(t, lower_is_open(s));
+    set_upper_is_open(t, upper_is_open(s));
+    SASSERT(check_invariant(t));
+}
+
+template<typename C>
+bool interval_manager<C>::eq(interval const & a, interval const & b) const {
+    return 
+        ::eq(m(), lower(a), lower_kind(a), lower(b), lower_kind(b)) &&
+        ::eq(m(), upper(a), upper_kind(a), upper(b), upper_kind(b)) &&
+        lower_is_open(a) == lower_is_open(b) &&
+        upper_is_open(a) == upper_is_open(b);
+}
+
+template<typename C>
+void interval_manager<C>::neg_jst(interval const & a, interval_deps & b_deps) {
+    if (lower_is_inf(a)) {
+        if (upper_is_inf(a)) {
+            b_deps.m_lower_deps = 0;
+            b_deps.m_upper_deps = 0;
+        }
+        else {
+            b_deps.m_lower_deps = DEP_IN_UPPER1;
+            b_deps.m_upper_deps = 0;
+        }
+    }
+    else {
+        if (upper_is_inf(a)) {
+            b_deps.m_lower_deps = 0;
+            b_deps.m_upper_deps = DEP_IN_LOWER1;
+        }
+        else {
+            b_deps.m_lower_deps = DEP_IN_UPPER1;
+            b_deps.m_upper_deps = DEP_IN_LOWER1;
+        }
+    }
+}
+
+template<typename C>
+void interval_manager<C>::neg(interval const & a, interval & b, interval_deps & b_deps) {
+    neg_jst(a, b_deps);
+    neg(a, b);
+}
+
+template<typename C>
+void interval_manager<C>::neg(interval const & a, interval & b) {
+    if (lower_is_inf(a)) {
+        if (upper_is_inf(a)) {
+            reset(b);
+        }
+        else {
+            m().set(lower(b), upper(a));
+            m().neg(lower(b));
+            set_lower_is_inf(b, false);
+            set_lower_is_open(b, upper_is_open(a));
+
+            m().reset(upper(b));
+            set_upper_is_inf(b, true);
+            set_upper_is_open(b, true);
+        }
+    }
+    else {
+        if (upper_is_inf(a)) {
+            m().set(upper(b), lower(a));
+            m().neg(upper(b));
+            set_upper_is_inf(b, false);
+            set_upper_is_open(b, lower_is_open(a));
+
+            m().reset(lower(b));
+            set_lower_is_inf(b, true);
+            set_lower_is_open(b, true);
+        }
+        else {
+            if (&a == &b) {
+                m().swap(lower(b), upper(b));
+            }
+            else {
+                m().set(lower(b), upper(a));
+                m().set(upper(b), lower(a));
+            }
+            m().neg(lower(b));
+            m().neg(upper(b));
+            set_lower_is_inf(b, false);
+            set_upper_is_inf(b, false);
+            bool l_o = lower_is_open(a);
+            bool u_o = upper_is_open(a);
+            set_lower_is_open(b, u_o);
+            set_upper_is_open(b, l_o);
+        }
+    }
+    SASSERT(check_invariant(b));
+}
+
+template<typename C>
+void interval_manager<C>::add_jst(interval const & a, interval const & b, interval_deps & c_deps) {
+    c_deps.m_lower_deps = DEP_IN_LOWER1 | DEP_IN_LOWER2;
+    c_deps.m_upper_deps = DEP_IN_UPPER1 | DEP_IN_UPPER2;
+}
+
+template<typename C>
+void interval_manager<C>::add(interval const & a, interval const & b, interval & c, interval_deps & c_deps) {
+    add_jst(a, b, c_deps);
+    add(a, b, c);
+}
+    
+template<typename C>
+void interval_manager<C>::add(interval const & a, interval const & b, interval & c) {
+    ext_numeral_kind new_l_kind, new_u_kind;
+    round_to_minus_inf();
+    ::add(m(), lower(a), lower_kind(a), lower(b), lower_kind(b), lower(c), new_l_kind);
+    round_to_plus_inf();
+    ::add(m(), upper(a), upper_kind(a), upper(b), upper_kind(b), upper(c), new_u_kind);
+    set_lower_is_inf(c, new_l_kind == EN_MINUS_INFINITY);
+    set_upper_is_inf(c, new_u_kind == EN_PLUS_INFINITY);
+    set_lower_is_open(c, lower_is_open(a) || lower_is_open(b));
+    set_upper_is_open(c, upper_is_open(a) || upper_is_open(b));
+    SASSERT(check_invariant(c));
+}
+
+template<typename C>
+void interval_manager<C>::sub_jst(interval const & a, interval const & b, interval_deps & c_deps) {
+    c_deps.m_lower_deps = DEP_IN_LOWER1 | DEP_IN_UPPER2;
+    c_deps.m_upper_deps = DEP_IN_UPPER1 | DEP_IN_LOWER2;
+}
+
+template<typename C>
+void interval_manager<C>::sub(interval const & a, interval const & b, interval & c, interval_deps & c_deps) {
+    sub_jst(a, b, c_deps);
+    sub(a, b, c);
+}
+
+template<typename C>
+void interval_manager<C>::sub(interval const & a, interval const & b, interval & c) {
+    ext_numeral_kind new_l_kind, new_u_kind;
+    round_to_minus_inf();
+    ::sub(m(), lower(a), lower_kind(a), upper(b), upper_kind(b), lower(c), new_l_kind);
+    round_to_plus_inf();
+    ::sub(m(), upper(a), upper_kind(a), lower(b), lower_kind(b), upper(c), new_u_kind);
+    set_lower_is_inf(c, new_l_kind == EN_MINUS_INFINITY);
+    set_upper_is_inf(c, new_u_kind == EN_PLUS_INFINITY);
+    set_lower_is_open(c, lower_is_open(a) || upper_is_open(b));
+    set_upper_is_open(c, upper_is_open(a) || lower_is_open(b));
+    SASSERT(check_invariant(c));
+}
+
+template<typename C>
+void interval_manager<C>::mul_jst(numeral const & k, interval const & a, interval_deps & b_deps) {
+    if (m().is_zero(k)) {
+        b_deps.m_lower_deps = 0;
+        b_deps.m_upper_deps = 0;
+   }
+    else if (m().is_neg(k)) {
+        b_deps.m_lower_deps = DEP_IN_UPPER1;
+        b_deps.m_upper_deps = DEP_IN_LOWER1;
+    }
+    else {
+        b_deps.m_lower_deps = DEP_IN_LOWER1;
+        b_deps.m_upper_deps = DEP_IN_UPPER1;
+    }
+}
+
+template<typename C>
+void interval_manager<C>::div_mul(numeral const & k, interval const & a, interval & b, bool inv_k) {
+    if (m().is_zero(k)) {
+        reset(b);
+    }
+    else {
+        numeral const & l = lower(a); ext_numeral_kind l_k = lower_kind(a);
+        numeral const & u = upper(a); ext_numeral_kind u_k = upper_kind(a);
+        numeral & new_l_val = m_result_lower;
+        numeral & new_u_val = m_result_upper;
+        ext_numeral_kind new_l_kind, new_u_kind;
+        bool l_o = lower_is_open(a);
+        bool u_o = upper_is_open(a);
+        if (m().is_pos(k)) {
+            set_lower_is_open(b, l_o);
+            set_upper_is_open(b, u_o);
+            if (inv_k) {
+                round_to_minus_inf();
+                m().inv(k, m_inv_k);
+                ::mul(m(), l, l_k, m_inv_k, EN_NUMERAL, new_l_val, new_l_kind);
+                
+                round_to_plus_inf();
+                m().inv(k, m_inv_k);
+                ::mul(m(), u, u_k, m_inv_k, EN_NUMERAL, new_u_val, new_u_kind);
+            }
+            else {
+                round_to_minus_inf();
+                ::mul(m(), l, l_k, k, EN_NUMERAL, new_l_val, new_l_kind);
+                round_to_plus_inf();
+                ::mul(m(), u, u_k, k, EN_NUMERAL, new_u_val, new_u_kind);
+            }
+        }
+        else {
+            set_lower_is_open(b, u_o);
+            set_upper_is_open(b, l_o);
+            if (inv_k) {
+                round_to_minus_inf();
+                m().inv(k, m_inv_k);
+                ::mul(m(), u, u_k, m_inv_k, EN_NUMERAL, new_l_val, new_l_kind);
+
+                round_to_plus_inf();
+                m().inv(k, m_inv_k);
+                ::mul(m(), l, l_k, m_inv_k, EN_NUMERAL, new_u_val, new_u_kind);
+            }
+            else {
+                round_to_minus_inf();
+                ::mul(m(), u, u_k, k, EN_NUMERAL, new_l_val, new_l_kind);
+                round_to_plus_inf();
+                ::mul(m(), l, l_k, k, EN_NUMERAL, new_u_val, new_u_kind);
+            }
+        }
+        m().swap(lower(b), new_l_val);
+        m().swap(upper(b), new_u_val);
+        set_lower_is_inf(b, new_l_kind == EN_MINUS_INFINITY);
+        set_upper_is_inf(b, new_u_kind == EN_PLUS_INFINITY);
+    }
+}
+
+template<typename C>
+void interval_manager<C>::mul(numeral const & k, interval const & a, interval & b, interval_deps & b_deps) {
+    mul_jst(k, a, b_deps);
+    mul(k, a, b);
+}
+
+template<typename C>
+void interval_manager<C>::mul(int n, int d, interval const & a, interval & b) {
+    _scoped_numeral<numeral_manager> aux(m());
+    m().set(aux, n, d);
+    mul(aux, a, b);
+}
+
+template<typename C>
+void interval_manager<C>::div(interval const & a, numeral const & k, interval & b, interval_deps & b_deps) {
+    div_jst(a, k, b_deps);
+    div(a, k, b);
+}
+
+template<typename C>
+void interval_manager<C>::mul_jst(interval const & i1, interval const & i2, interval_deps & r_deps) {
+    if (is_zero(i1)) {
+        r_deps.m_lower_deps = DEP_IN_LOWER1 | DEP_IN_UPPER1;
+        r_deps.m_upper_deps = DEP_IN_LOWER1 | DEP_IN_UPPER1;
+    }
+    else if (is_zero(i2)) {
+        r_deps.m_lower_deps = DEP_IN_LOWER2 | DEP_IN_UPPER2;
+        r_deps.m_upper_deps = DEP_IN_LOWER2 | DEP_IN_UPPER2;
+    }
+    else if (is_N(i1)) {
+        if (is_N(i2)) {
+            // x <= b <= 0, y <= d <= 0 --> b*d <= x*y
+            // a <= x <= b <= 0, c <= y <= d <= 0 --> x*y <= a*c  (we can use the fact that x or y is always negative (i.e., b is neg or d is neg))
+            r_deps.m_lower_deps = DEP_IN_UPPER1 | DEP_IN_UPPER2;
+            r_deps.m_upper_deps = DEP_IN_LOWER1 | DEP_IN_LOWER2 | DEP_IN_UPPER1; // we can replace DEP_IN_UPPER1 with DEP_IN_UPPER2
+        }
+        else if (is_M(i2)) {
+            // a <= x <= b <= 0,  y <= d, d > 0 --> a*d <= x*y (uses the fact that b is not positive)
+            // a <= x <= b <= 0,  c <= y, c < 0 --> x*y <= a*c (uses the fact that b is not positive)
+            r_deps.m_lower_deps = DEP_IN_LOWER1 | DEP_IN_UPPER2 | DEP_IN_UPPER1;
+            r_deps.m_upper_deps = DEP_IN_LOWER1 | DEP_IN_LOWER2 | DEP_IN_UPPER1;
+        }
+        else {
+            // a <= x <= b <= 0, 0 <= c <= y <= d --> a*d <= x*y (uses the fact that x is neg (b is not positive) or y is pos (c is not negative))
+            // x <= b <= 0,  0 <= c <= y --> x*y <= b*c
+            r_deps.m_lower_deps = DEP_IN_LOWER1 | DEP_IN_UPPER2 | DEP_IN_UPPER1; // we can replace DEP_IN_UPPER1 with DEP_IN_UPPER2
+            r_deps.m_upper_deps = DEP_IN_UPPER1 | DEP_IN_LOWER2;
+        }
+    }
+    else if (is_M(i1)) {
+        if (is_N(i2)) {
+            // b > 0, x <= b,  c <= y <= d <= 0 --> b*c <= x*y (uses the fact that d is not positive)
+            // a < 0, a <= x,  c <= y <= d <= 0 --> x*y <= a*c (uses the fact that d is not positive)
+            r_deps.m_lower_deps = DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2;
+            r_deps.m_upper_deps = DEP_IN_LOWER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2;
+        }
+        else if (is_M(i2)) {
+            r_deps.m_lower_deps = DEP_IN_LOWER1 | DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2;
+            r_deps.m_upper_deps = DEP_IN_LOWER1 | DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2;
+        }
+        else {
+            // a < 0, a <= x, 0 <= c <= y <= d --> a*d <= x*y (uses the fact that c is not negative)
+            // b > 0, x <= b, 0 <= c <= y <= d --> x*y <= b*d (uses the fact that c is not negative)
+            r_deps.m_lower_deps = DEP_IN_LOWER1 | DEP_IN_UPPER2 | DEP_IN_LOWER2;
+            r_deps.m_upper_deps = DEP_IN_UPPER1 | DEP_IN_UPPER2 | DEP_IN_LOWER2;
+        }
+    }
+    else {
+        SASSERT(is_P(i1));        
+        if (is_N(i2)) {
+            // 0 <= a <= x <= b,   c <= y <= d <= 0  -->  x*y <= b*c (uses the fact that x is pos (a is not neg) or y is neg (d is not pos))
+            // 0 <= a <= x,  y <= d <= 0  --> a*d <= x*y
+            r_deps.m_lower_deps = DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_LOWER1; // we can replace DEP_IN_LOWER1 with DEP_IN_UPPER2
+            r_deps.m_upper_deps = DEP_IN_LOWER1 | DEP_IN_UPPER2;
+        }
+        else if (is_M(i2)) {
+            // 0 <= a <= x <= b,  c <= y --> b*c <= x*y (uses the fact that a is not negative)
+            // 0 <= a <= x <= b,  y <= d --> x*y <= b*d (uses the fact that a is not negative)
+            r_deps.m_lower_deps = DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_LOWER1;
+            r_deps.m_upper_deps = DEP_IN_UPPER1 | DEP_IN_UPPER2 | DEP_IN_LOWER1;
+        }
+        else {
+            SASSERT(is_P(i2));
+            // 0 <= a <= x, 0 <= c <= y --> a*c <= x*y 
+            // x <= b, y <= d --> x*y <= b*d (uses the fact that x is pos (a is not negative) or y is pos (c is not negative))
+            r_deps.m_lower_deps = DEP_IN_LOWER1 | DEP_IN_LOWER2;
+            r_deps.m_upper_deps = DEP_IN_UPPER1 | DEP_IN_UPPER2 | DEP_IN_LOWER1; // we can replace DEP_IN_LOWER1 with DEP_IN_LOWER2
+        }
+    }
+}
+
+template<typename C>
+void interval_manager<C>::mul(interval const & i1, interval const & i2, interval & r, interval_deps & r_deps) {
+    mul_jst(i1, i2, r_deps);
+    mul(i1, i2, r);
+}
+
+template<typename C>
+void interval_manager<C>::mul(interval const & i1, interval const & i2, interval & r) {
+#ifdef _TRACE
+    static unsigned call_id = 0;
+#endif
+#if Z3DEBUG || _TRACE
+    bool i1_contains_zero = contains_zero(i1);
+    bool i2_contains_zero = contains_zero(i2);
+#endif
+    if (is_zero(i1)) {
+        set(r, i1);
+        return;
+    }
+    if (is_zero(i2)) {
+        set(r, i2);
+        return;
+    }
+
+    numeral const & a = lower(i1); ext_numeral_kind a_k = lower_kind(i1);
+    numeral const & b = upper(i1); ext_numeral_kind b_k = upper_kind(i1);
+    numeral const & c = lower(i2); ext_numeral_kind c_k = lower_kind(i2);
+    numeral const & d = upper(i2); ext_numeral_kind d_k = upper_kind(i2);
+
+    bool a_o = lower_is_open(i1);
+    bool b_o = upper_is_open(i1);
+    bool c_o = lower_is_open(i2);
+    bool d_o = upper_is_open(i2);
+
+    numeral & new_l_val = m_result_lower;
+    numeral & new_u_val = m_result_upper;
+    ext_numeral_kind new_l_kind, new_u_kind;
+
+    if (is_N(i1)) {
+        if (is_N(i2)) {
+            // x <= b <= 0, y <= d <= 0 --> b*d <= x*y
+            // a <= x <= b <= 0, c <= y <= d <= 0 --> x*y <= a*c  (we can use the fact that x or y is always negative (i.e., b is neg or d is neg))
+            TRACE("interval_bug", tout << "(N, N) #" << call_id << "\n"; display(tout, i1); tout << "\n"; display(tout, i2); tout << "\n";
+                  tout << "a: "; m().display(tout, a); tout << "\n";
+                  tout << "b: "; m().display(tout, b); tout << "\n";
+                  tout << "c: "; m().display(tout, c); tout << "\n";
+                  tout << "d: "; m().display(tout, d); tout << "\n";
+                  tout << "is_N0(i1): " << is_N0(i1) << "\n";
+                  tout << "is_N0(i2): " << is_N0(i2) << "\n";
+                  );
+            set_lower_is_open(r, (is_N0(i1) || is_N0(i2)) ? false : (b_o || d_o));
+            set_upper_is_open(r, a_o || c_o);
+            // if b = 0 (and the interval is closed), then the lower bound is closed
+
+            round_to_minus_inf();
+            ::mul(m(), b, b_k, d, d_k, new_l_val, new_l_kind);
+            round_to_plus_inf();
+            ::mul(m(), a, a_k, c, c_k, new_u_val, new_u_kind);
+        }
+        else if (is_M(i2)) {
+            // a <= x <= b <= 0,  y <= d, d > 0 --> a*d <= x*y (uses the fact that b is not positive)
+            // a <= x <= b <= 0,  c <= y, c < 0 --> x*y <= a*c (uses the fact that b is not positive)
+            TRACE("interval_bug", tout << "(N, M) #" << call_id << "\n";);
+
+            set_lower_is_open(r, a_o || d_o);
+            set_upper_is_open(r, a_o || c_o);
+
+            round_to_minus_inf();
+            ::mul(m(), a, a_k, d, d_k, new_l_val, new_l_kind);
+            round_to_plus_inf();
+            ::mul(m(), a, a_k, c, c_k, new_u_val, new_u_kind);
+        }
+        else {
+            // a <= x <= b <= 0, 0 <= c <= y <= d --> a*d <= x*y (uses the fact that x is neg (b is not positive) or y is pos (c is not negative))
+            // x <= b <= 0,  0 <= c <= y --> x*y <= b*c
+            TRACE("interval_bug", tout << "(N, P) #" << call_id << "\n";);
+            SASSERT(is_P(i2));
+            
+            // must update upper_is_open first, since value of is_N0(i1) and is_P0(i2) may be affected by update 
+            set_upper_is_open(r, (is_N0(i1) || is_P0(i2)) ? false : (b_o || c_o));
+            set_lower_is_open(r, a_o || d_o);
+            
+            round_to_minus_inf();
+            ::mul(m(), a, a_k, d, d_k, new_l_val, new_l_kind);
+            round_to_plus_inf();
+            ::mul(m(), b, b_k, c, c_k, new_u_val, new_u_kind);
+        }
+    }
+    else if (is_M(i1)) {
+        if (is_N(i2)) {
+            // b > 0, x <= b,  c <= y <= d <= 0 --> b*c <= x*y (uses the fact that d is not positive)
+            // a < 0, a <= x,  c <= y <= d <= 0 --> x*y <= a*c (uses the fact that d is not positive)
+            TRACE("interval_bug", tout << "(M, N) #" << call_id << "\n";);
+
+            set_lower_is_open(r, b_o || c_o);
+            set_upper_is_open(r, a_o || c_o); 
+
+            round_to_minus_inf();
+            ::mul(m(), b, b_k, c, c_k, new_l_val, new_l_kind);
+            round_to_plus_inf();
+            ::mul(m(), a, a_k, c, c_k, new_u_val, new_u_kind);
+        }
+        else if (is_M(i2)) {
+            numeral & ad = m_mul_ad; ext_numeral_kind ad_k;
+            numeral & bc = m_mul_bc; ext_numeral_kind bc_k;
+            numeral & ac = m_mul_ac; ext_numeral_kind ac_k;
+            numeral & bd = m_mul_bd; ext_numeral_kind bd_k;
+
+            bool  ad_o = a_o || d_o;
+            bool  bc_o = b_o || c_o;
+            bool  ac_o = a_o || c_o;
+            bool  bd_o = b_o || d_o;
+            
+            round_to_minus_inf();
+            ::mul(m(), a, a_k, d, d_k, ad, ad_k);
+            ::mul(m(), b, b_k, c, c_k, bc, bc_k);
+            round_to_plus_inf();
+            ::mul(m(), a, a_k, c, c_k, ac, ac_k);
+            ::mul(m(), b, b_k, d, d_k, bd, bd_k);
+
+            if (::lt(m(), ad, ad_k, bc, bc_k) || (::eq(m(), ad, ad_k, bc, bc_k) && !ad_o && bc_o)) {
+                m().swap(new_l_val, ad);
+                new_l_kind = ad_k;
+                set_lower_is_open(r, ad_o);
+            }
+            else {
+                m().swap(new_l_val, bc);
+                new_l_kind = bc_k;
+                set_lower_is_open(r, bc_o);
+            }
+
+            
+            if (::gt(m(), ac, ac_k, bd, bd_k) || (::eq(m(), ac, ac_k, bd, bd_k) && !ac_o && bd_o)) {
+                m().swap(new_u_val, ac);
+                new_u_kind = ac_k;
+                set_upper_is_open(r, ac_o);
+            }
+            else {
+                m().swap(new_u_val, bd);
+                new_u_kind = bd_k;
+                set_upper_is_open(r, bd_o);
+            }
+        }
+        else {
+            // a < 0, a <= x, 0 <= c <= y <= d --> a*d <= x*y (uses the fact that c is not negative)
+            // b > 0, x <= b, 0 <= c <= y <= d --> x*y <= b*d (uses the fact that c is not negative)
+            TRACE("interval_bug", tout << "(M, P) #" << call_id << "\n";);
+            SASSERT(is_P(i2));
+
+            set_lower_is_open(r, a_o || d_o);
+            set_upper_is_open(r, b_o || d_o); 
+
+            round_to_minus_inf();
+            ::mul(m(), a, a_k, d, d_k, new_l_val, new_l_kind);
+            round_to_plus_inf();
+            ::mul(m(), b, b_k, d, d_k, new_u_val, new_u_kind);
+        }
+    }
+    else {
+        SASSERT(is_P(i1));        
+        if (is_N(i2)) {
+            // 0 <= a <= x <= b,   c <= y <= d <= 0  -->  x*y <= b*c (uses the fact that x is pos (a is not neg) or y is neg (d is not pos))
+            // 0 <= a <= x,  y <= d <= 0  --> a*d <= x*y
+            TRACE("interval_bug", tout << "(P, N) #" << call_id << "\n";);
+            
+            // must update upper_is_open first, since value of is_P0(i1) and is_N0(i2) may be affected by update 
+            set_upper_is_open(r, (is_P0(i1) || is_N0(i2)) ? false : a_o || d_o);
+            set_lower_is_open(r, b_o || c_o);
+
+            round_to_minus_inf();
+            ::mul(m(), b, b_k, c, c_k, new_l_val, new_l_kind);
+            round_to_plus_inf();
+            ::mul(m(), a, a_k, d, d_k, new_u_val, new_u_kind);
+        }
+        else if (is_M(i2)) {
+            // 0 <= a <= x <= b,  c <= y --> b*c <= x*y (uses the fact that a is not negative)
+            // 0 <= a <= x <= b,  y <= d --> x*y <= b*d (uses the fact that a is not negative)
+            TRACE("interval_bug", tout << "(P, M) #" << call_id << "\n";);
+
+            set_lower_is_open(r, b_o || c_o);
+            set_upper_is_open(r, b_o || d_o);
+
+            round_to_minus_inf();
+            ::mul(m(), b, b_k, c, c_k, new_l_val, new_l_kind);
+            round_to_plus_inf();
+            ::mul(m(), b, b_k, d, d_k, new_u_val, new_u_kind);
+        }
+        else {
+            SASSERT(is_P(i2));
+            // 0 <= a <= x, 0 <= c <= y --> a*c <= x*y 
+            // x <= b, y <= d --> x*y <= b*d (uses the fact that x is pos (a is not negative) or y is pos (c is not negative))
+            TRACE("interval_bug", tout << "(P, P) #" << call_id << "\n";);
+
+            set_lower_is_open(r, (is_P0(i1) || is_P0(i2