3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-04-08 18:31:49 +00:00

Merge branch 'realclosure' into unstable

This commit is contained in:
Leonardo de Moura 2013-01-12 22:03:40 -08:00
commit 93f37bdf9f
27 changed files with 7901 additions and 42 deletions

View file

@ -15,6 +15,7 @@ def init_project_def():
add_lib('sat', ['util'])
add_lib('nlsat', ['polynomial', 'sat'])
add_lib('interval', ['util'], 'math/interval')
add_lib('realclosure', ['interval'], 'math/realclosure')
add_lib('subpaving', ['interval'], 'math/subpaving')
add_lib('ast', ['util', 'polynomial'])
add_lib('rewriter', ['ast', 'polynomial'], 'ast/rewriter')
@ -58,8 +59,8 @@ def init_project_def():
add_lib('ufbv_tactic', ['normal_forms', 'core_tactics', 'macros', 'smt_tactic', 'rewriter'], 'tactic/ufbv')
add_lib('portfolio', ['smtlogic_tactics', 'ufbv_tactic', 'fpa', 'aig_tactic', 'muz_qe', 'sls_tactic', 'subpaving_tactic'], 'tactic/portfolio')
add_lib('smtparser', ['portfolio'], 'parsers/smt')
API_files = ['z3_api.h', 'z3_algebraic.h', 'z3_polynomial.h']
add_lib('api', ['portfolio', 'user_plugin', 'smtparser'],
API_files = ['z3_api.h', 'z3_algebraic.h', 'z3_polynomial.h', 'z3_rcf.h']
add_lib('api', ['portfolio', 'user_plugin', 'smtparser', 'realclosure'],
includes2install=['z3.h', 'z3_v1.h', 'z3_macros.h'] + API_files)
add_exe('shell', ['api', 'sat', 'extra_cmds'], exe_name='z3')
add_exe('test', ['api', 'fuzzing'], exe_name='test-z3', install=False)

View file

@ -1320,11 +1320,13 @@ def mk_config():
check_ar()
CXX = find_cxx_compiler()
CC = find_c_compiler()
SLIBEXTRAFLAGS = ''
if GMP:
test_gmp(CXX)
ARITH = "gmp"
CPPFLAGS = '%s -D_MP_GMP' % CPPFLAGS
LDFLAGS = '%s -lgmp' % LDFLAGS
SLIBEXTRAFLAGS = '%s -lgmp' % SLIBEXTRAFLAGS
else:
CPPFLAGS = '%s -D_MP_INTERNAL' % CPPFLAGS
CXXFLAGS = '%s -c' % CXXFLAGS
@ -1332,10 +1334,9 @@ def mk_config():
if HAS_OMP:
CXXFLAGS = '%s -fopenmp -mfpmath=sse' % CXXFLAGS
LDFLAGS = '%s -fopenmp' % LDFLAGS
SLIBEXTRAFLAGS = '-fopenmp'
SLIBEXTRAFLAGS = '%s -fopenmp' % SLIBEXTRAFLAGS
else:
CXXFLAGS = '%s -D_NO_OMP_' % CXXFLAGS
SLIBEXTRAFLAGS = ''
if DEBUG_MODE:
CXXFLAGS = '%s -g -Wall' % CXXFLAGS
else:

View file

@ -25,6 +25,7 @@ Revision History:
#include"api_log_macros.h"
#include"api_util.h"
#include"reg_decl_plugins.h"
#include"realclosure.h"
// The install_tactics procedure is automatically generated
void install_tactics(tactic_manager & ctx);
@ -138,6 +139,8 @@ namespace api {
if (m_interruptable)
(*m_interruptable)();
m().set_cancel(true);
if (m_rcf_manager.get() == 0)
m_rcf_manager->set_cancel(true);
}
}
@ -391,6 +394,19 @@ namespace api {
m_smtlib_parser_sorts.reset();
}
}
// ------------------------
//
// RCF manager
//
// -----------------------
realclosure::manager & context::rcfm() {
if (m_rcf_manager.get() == 0) {
m_rcf_manager = alloc(realclosure::manager, m_rcf_qm);
}
return *(m_rcf_manager.get());
}
};
@ -476,8 +492,11 @@ extern "C" {
}
void Z3_API Z3_enable_trace(Z3_string tag) {
memory::initialize(UINT_MAX);
LOG_Z3_enable_trace(tag);
enable_trace(tag);
// Tag is a string that was probably not allocated by Z3. Create a copy using symbol.
symbol tag_sym(tag);
enable_trace(tag_sym.bare_str());
}
void Z3_API Z3_disable_trace(Z3_string tag) {

View file

@ -38,6 +38,10 @@ namespace smtlib {
class parser;
};
namespace realclosure {
class manager;
};
namespace api {
Z3_search_failure mk_Z3_search_failure(smt::failure f);
@ -83,7 +87,6 @@ namespace api {
event_handler * m_interruptable; // Reference to an object that can be interrupted by Z3_interrupt
pmanager m_pmanager;
public:
// Scoped obj for setting m_interruptable
class set_interruptable {
@ -175,8 +178,22 @@ namespace api {
// Polynomial manager & caches
//
// -----------------------
private:
pmanager m_pmanager;
public:
polynomial::manager & pm() { return m_pmanager.pm(); }
// ------------------------
//
// RCF manager
//
// -----------------------
private:
unsynch_mpq_manager m_rcf_qm;
scoped_ptr<realclosure::manager> m_rcf_manager;
public:
realclosure::manager & rcfm();
// ------------------------
//
// Solver interface for backward compatibility

293
src/api/api_rcf.cpp Normal file
View file

@ -0,0 +1,293 @@
/*++
Copyright (c) 2013 Microsoft Corporation
Module Name:
api_rcf.cpp
Abstract:
Additional APIs for handling elements of the Z3 real closed field that contains:
- transcendental extensions
- infinitesimal extensions
- algebraic extensions
Author:
Leonardo de Moura (leonardo) 2012-01-05
Notes:
--*/
#include<iostream>
#include"z3.h"
#include"api_log_macros.h"
#include"api_context.h"
#include"realclosure.h"
extern "C" {
static rcmanager & rcfm(Z3_context c) {
return mk_c(c)->rcfm();
}
static void reset_rcf_cancel(Z3_context c) {
rcfm(c).reset_cancel();
}
static rcnumeral to_rcnumeral(Z3_rcf_num a) {
return rcnumeral::mk(a);
}
static Z3_rcf_num from_rcnumeral(rcnumeral a) {
return reinterpret_cast<Z3_rcf_num>(a.c_ptr());
}
void Z3_API Z3_rcf_del(Z3_context c, Z3_rcf_num a) {
Z3_TRY;
LOG_Z3_rcf_del(c, a);
RESET_ERROR_CODE();
rcnumeral _a = to_rcnumeral(a);
rcfm(c).del(_a);
Z3_CATCH;
}
Z3_rcf_num Z3_API Z3_rcf_mk_rational(Z3_context c, Z3_string val) {
Z3_TRY;
LOG_Z3_rcf_mk_rational(c, val);
RESET_ERROR_CODE();
reset_rcf_cancel(c);
scoped_mpq q(rcfm(c).qm());
rcfm(c).qm().set(q, val);
rcnumeral r;
rcfm(c).set(r, q);
RETURN_Z3(from_rcnumeral(r));
Z3_CATCH_RETURN(0);
}
Z3_rcf_num Z3_API Z3_rcf_mk_small_int(Z3_context c, int val) {
Z3_TRY;
LOG_Z3_rcf_mk_small_int(c, val);
RESET_ERROR_CODE();
reset_rcf_cancel(c);
rcnumeral r;
rcfm(c).set(r, val);
RETURN_Z3(from_rcnumeral(r));
Z3_CATCH_RETURN(0);
}
Z3_rcf_num Z3_API Z3_rcf_mk_pi(Z3_context c) {
Z3_TRY;
LOG_Z3_rcf_mk_pi(c);
RESET_ERROR_CODE();
reset_rcf_cancel(c);
rcnumeral r;
rcfm(c).mk_pi(r);
RETURN_Z3(from_rcnumeral(r));
Z3_CATCH_RETURN(0);
}
Z3_rcf_num Z3_API Z3_rcf_mk_e(Z3_context c) {
Z3_TRY;
LOG_Z3_rcf_mk_e(c);
RESET_ERROR_CODE();
reset_rcf_cancel(c);
rcnumeral r;
rcfm(c).mk_e(r);
RETURN_Z3(from_rcnumeral(r));
Z3_CATCH_RETURN(0);
}
Z3_rcf_num Z3_API Z3_rcf_mk_infinitesimal(Z3_context c, Z3_string name) {
Z3_TRY;
LOG_Z3_rcf_mk_infinitesimal(c, name);
RESET_ERROR_CODE();
reset_rcf_cancel(c);
rcnumeral r;
rcfm(c).mk_infinitesimal(name, r);
RETURN_Z3(from_rcnumeral(r));
Z3_CATCH_RETURN(0);
}
unsigned Z3_API Z3_rcf_mk_roots(Z3_context c, unsigned n, Z3_rcf_num const a[], Z3_rcf_num roots[]) {
Z3_TRY;
LOG_Z3_rcf_mk_roots(c, n, a, roots);
RESET_ERROR_CODE();
reset_rcf_cancel(c);
rcnumeral_vector av;
unsigned rz = 0;
for (unsigned i = 0; i < n; i++) {
if (!rcfm(c).is_zero(to_rcnumeral(a[i])))
rz = i + 1;
av.push_back(to_rcnumeral(a[i]));
}
if (rz == 0) {
// it is the zero polynomial
SET_ERROR_CODE(Z3_INVALID_ARG);
return 0;
}
av.shrink(rz);
rcnumeral_vector rs;
rcfm(c).isolate_roots(av.size(), av.c_ptr(), rs);
unsigned num_roots = rs.size();
for (unsigned i = 0; i < num_roots; i++) {
roots[i] = from_rcnumeral(rs[i]);
}
RETURN_Z3_rcf_mk_roots num_roots;
Z3_CATCH_RETURN(0);
}
Z3_rcf_num Z3_API Z3_rcf_add(Z3_context c, Z3_rcf_num a, Z3_rcf_num b) {
Z3_TRY;
LOG_Z3_rcf_add(c, a, b);
RESET_ERROR_CODE();
reset_rcf_cancel(c);
rcnumeral r;
rcfm(c).add(to_rcnumeral(a), to_rcnumeral(b), r);
RETURN_Z3(from_rcnumeral(r));
Z3_CATCH_RETURN(0);
}
Z3_rcf_num Z3_API Z3_rcf_sub(Z3_context c, Z3_rcf_num a, Z3_rcf_num b) {
Z3_TRY;
LOG_Z3_rcf_sub(c, a, b);
RESET_ERROR_CODE();
reset_rcf_cancel(c);
rcnumeral r;
rcfm(c).sub(to_rcnumeral(a), to_rcnumeral(b), r);
RETURN_Z3(from_rcnumeral(r));
Z3_CATCH_RETURN(0);
}
Z3_rcf_num Z3_API Z3_rcf_mul(Z3_context c, Z3_rcf_num a, Z3_rcf_num b) {
Z3_TRY;
LOG_Z3_rcf_mul(c, a, b);
RESET_ERROR_CODE();
reset_rcf_cancel(c);
rcnumeral r;
rcfm(c).mul(to_rcnumeral(a), to_rcnumeral(b), r);
RETURN_Z3(from_rcnumeral(r));
Z3_CATCH_RETURN(0);
}
Z3_rcf_num Z3_API Z3_rcf_div(Z3_context c, Z3_rcf_num a, Z3_rcf_num b) {
Z3_TRY;
LOG_Z3_rcf_div(c, a, b);
RESET_ERROR_CODE();
reset_rcf_cancel(c);
rcnumeral r;
rcfm(c).div(to_rcnumeral(a), to_rcnumeral(b), r);
RETURN_Z3(from_rcnumeral(r));
Z3_CATCH_RETURN(0);
}
Z3_rcf_num Z3_API Z3_rcf_neg(Z3_context c, Z3_rcf_num a) {
Z3_TRY;
LOG_Z3_rcf_neg(c, a);
RESET_ERROR_CODE();
reset_rcf_cancel(c);
rcnumeral r;
rcfm(c).neg(to_rcnumeral(a), r);
RETURN_Z3(from_rcnumeral(r));
Z3_CATCH_RETURN(0);
}
Z3_rcf_num Z3_API Z3_rcf_inv(Z3_context c, Z3_rcf_num a) {
Z3_TRY;
LOG_Z3_rcf_inv(c, a);
RESET_ERROR_CODE();
reset_rcf_cancel(c);
rcnumeral r;
rcfm(c).inv(to_rcnumeral(a), r);
RETURN_Z3(from_rcnumeral(r));
Z3_CATCH_RETURN(0);
}
Z3_rcf_num Z3_API Z3_rcf_power(Z3_context c, Z3_rcf_num a, unsigned k) {
Z3_TRY;
LOG_Z3_rcf_power(c, a, k);
RESET_ERROR_CODE();
reset_rcf_cancel(c);
rcnumeral r;
rcfm(c).power(to_rcnumeral(a), k, r);
RETURN_Z3(from_rcnumeral(r));
Z3_CATCH_RETURN(0);
}
Z3_bool Z3_API Z3_rcf_lt(Z3_context c, Z3_rcf_num a, Z3_rcf_num b) {
Z3_TRY;
LOG_Z3_rcf_lt(c, a, b);
RESET_ERROR_CODE();
reset_rcf_cancel(c);
return rcfm(c).lt(to_rcnumeral(a), to_rcnumeral(b));
Z3_CATCH_RETURN(Z3_FALSE);
}
Z3_bool Z3_API Z3_rcf_gt(Z3_context c, Z3_rcf_num a, Z3_rcf_num b) {
Z3_TRY;
LOG_Z3_rcf_gt(c, a, b);
RESET_ERROR_CODE();
reset_rcf_cancel(c);
return rcfm(c).gt(to_rcnumeral(a), to_rcnumeral(b));
Z3_CATCH_RETURN(Z3_FALSE);
}
Z3_bool Z3_API Z3_rcf_le(Z3_context c, Z3_rcf_num a, Z3_rcf_num b) {
Z3_TRY;
LOG_Z3_rcf_le(c, a, b);
RESET_ERROR_CODE();
reset_rcf_cancel(c);
return rcfm(c).le(to_rcnumeral(a), to_rcnumeral(b));
Z3_CATCH_RETURN(Z3_FALSE);
}
Z3_bool Z3_API Z3_rcf_ge(Z3_context c, Z3_rcf_num a, Z3_rcf_num b) {
Z3_TRY;
LOG_Z3_rcf_ge(c, a, b);
RESET_ERROR_CODE();
reset_rcf_cancel(c);
return rcfm(c).ge(to_rcnumeral(a), to_rcnumeral(b));
Z3_CATCH_RETURN(Z3_FALSE);
}
Z3_bool Z3_API Z3_rcf_eq(Z3_context c, Z3_rcf_num a, Z3_rcf_num b) {
Z3_TRY;
LOG_Z3_rcf_eq(c, a, b);
RESET_ERROR_CODE();
reset_rcf_cancel(c);
return rcfm(c).eq(to_rcnumeral(a), to_rcnumeral(b));
Z3_CATCH_RETURN(Z3_FALSE);
}
Z3_bool Z3_API Z3_rcf_neq(Z3_context c, Z3_rcf_num a, Z3_rcf_num b) {
Z3_TRY;
LOG_Z3_rcf_neq(c, a, b);
RESET_ERROR_CODE();
reset_rcf_cancel(c);
return rcfm(c).eq(to_rcnumeral(a), to_rcnumeral(b));
Z3_CATCH_RETURN(Z3_FALSE);
}
Z3_string Z3_API Z3_rcf_num_to_string(Z3_context c, Z3_rcf_num a) {
Z3_TRY;
LOG_Z3_rcf_num_to_string(c, a);
RESET_ERROR_CODE();
reset_rcf_cancel(c);
std::ostringstream buffer;
rcfm(c).display(buffer, to_rcnumeral(a));
return mk_c(c)->mk_external_string(buffer.str());
Z3_CATCH_RETURN("");
}
Z3_string Z3_API Z3_rcf_num_to_decimal_string(Z3_context c, Z3_rcf_num a, unsigned prec) {
Z3_TRY;
LOG_Z3_rcf_num_to_decimal_string(c, a, prec);
RESET_ERROR_CODE();
reset_rcf_cancel(c);
std::ostringstream buffer;
rcfm(c).display_decimal(buffer, to_rcnumeral(a), prec);
return mk_c(c)->mk_external_string(buffer.str());
Z3_CATCH_RETURN("");
}
};

View file

@ -37,6 +37,7 @@ namespace api {
public:
object():m_ref_count(0) {}
virtual ~object() {}
unsigned ref_count() const { return m_ref_count; }
void inc_ref() { m_ref_count++; }
void dec_ref() { SASSERT(m_ref_count > 0); m_ref_count--; if (m_ref_count == 0) dealloc(this); }
};

153
src/api/python/z3rcf.py Normal file
View file

@ -0,0 +1,153 @@
############################################
# Copyright (c) 2013 Microsoft Corporation
#
# Z3 Python interface for Z3 Real Closed Fields
# that may contain
# - computable transcendentals
# - infinitesimals
# - algebraic extensions
#
# Author: Leonardo de Moura (leonardo)
############################################
from z3 import *
from z3core import *
from z3printer import *
from fractions import Fraction
def _to_rcfnum(num, ctx=None):
if isinstance(num, RCFNum):
return num
else:
return RCFNum(num, ctx)
def Pi(ctx=None):
ctx = z3._get_ctx(ctx)
return RCFNum(Z3_rcf_mk_pi(ctx.ref()), ctx)
def E(ctx=None):
ctx = z3._get_ctx(ctx)
return RCFNum(Z3_rcf_mk_e(ctx.ref()), ctx)
def MkInfinitesimal(name="eps", ctx=None):
ctx = z3._get_ctx(ctx)
return RCFNum(Z3_rcf_mk_infinitesimal(ctx.ref(), name), ctx)
def MkRoots(p, ctx=None):
ctx = z3._get_ctx(ctx)
num = len(p)
_tmp = []
_as = (RCFNumObj * num)()
_rs = (RCFNumObj * num)()
for i in range(num):
_a = _to_rcfnum(p[i], ctx)
_tmp.append(_a) # prevent GC
_as[i] = _a.num
nr = Z3_rcf_mk_roots(ctx.ref(), num, _as, _rs)
r = []
for i in range(nr):
r.append(RCFNum(_rs[i], ctx))
return r
class RCFNum:
def __init__(self, num, ctx=None):
# TODO: add support for converting AST numeral values into RCFNum
if isinstance(num, RCFNumObj):
self.num = num
self.ctx = z3._get_ctx(ctx)
else:
self.ctx = z3._get_ctx(ctx)
self.num = Z3_rcf_mk_rational(self.ctx_ref(), str(num))
def __del__(self):
Z3_rcf_del(self.ctx_ref(), self.num)
def ctx_ref(self):
return self.ctx.ref()
def __repr__(self):
return Z3_rcf_num_to_string(self.ctx_ref(), self.num)
def __add__(self, other):
v = _to_rcfnum(other, self.ctx)
return RCFNum(Z3_rcf_add(self.ctx_ref(), self.num, v.num), self.ctx)
def __radd__(self, other):
v = _to_rcfnum(other, self.ctx)
return RCFNum(Z3_rcf_add(self.ctx_ref(), v.num, self.num), self.ctx)
def __mul__(self, other):
v = _to_rcfnum(other, self.ctx)
return RCFNum(Z3_rcf_mul(self.ctx_ref(), self.num, v.num), self.ctx)
def __rmul__(self, other):
v = _to_rcfnum(other, self.ctx)
return RCFNum(Z3_rcf_mul(self.ctx_ref(), v.num, self.num), self.ctx)
def __sub__(self, other):
v = _to_rcfnum(other, self.ctx)
return RCFNum(Z3_rcf_sub(self.ctx_ref(), self.num, v.num), self.ctx)
def __rsub__(self, other):
v = _to_rcfnum(other, self.ctx)
return RCFNum(Z3_rcf_sub(self.ctx_ref(), v.num, self.num), self.ctx)
def __div__(self, other):
v = _to_rcfnum(other, self.ctx)
return RCFNum(Z3_rcf_div(self.ctx_ref(), self.num, v.num), self.ctx)
def __rdiv__(self, other):
v = _to_rcfnum(other, self.ctx)
return RCFNum(Z3_rcf_div(self.ctx_ref(), v.num, self.num), self.ctx)
def __neg__(self):
return self.__rsub__(0)
def power(self, k):
return RCFNum(Z3_rcf_power(self.ctx_ref(), self.num, k), self.ctx)
def __pow__(self, k):
return self.power(k)
def decimal(self, prec=5):
return Z3_rcf_num_to_decimal_string(self.ctx_ref(), self.num, prec)
def __lt__(self, other):
v = _to_rcfnum(other, self.ctx)
return Z3_rcf_lt(self.ctx_ref(), self.num, v.num)
def __rlt__(self, other):
v = _to_rcfnum(other, self.ctx)
return Z3_rcf_lt(self.ctx_ref(), v.num, self.num)
def __gt__(self, other):
v = _to_rcfnum(other, self.ctx)
return Z3_rcf_gt(self.ctx_ref(), self.num, v.num)
def __rgt__(self, other):
v = _to_rcfnum(other, self.ctx)
return Z3_rcf_gt(self.ctx_ref(), v.num, self.num)
def __le__(self, other):
v = _to_rcfnum(other, self.ctx)
return Z3_rcf_le(self.ctx_ref(), self.num, v.num)
def __rle__(self, other):
v = _to_rcfnum(other, self.ctx)
return Z3_rcf_le(self.ctx_ref(), v.num, self.num)
def __ge__(self, other):
v = _to_rcfnum(other, self.ctx)
return Z3_rcf_ge(self.ctx_ref(), self.num, v.num)
def __rge__(self, other):
v = _to_rcfnum(other, self.ctx)
return Z3_rcf_ge(self.ctx_ref(), v.num, self.num)
def __eq__(self, other):
v = _to_rcfnum(other, self.ctx)
return Z3_rcf_eq(self.ctx_ref(), self.num, v.num)
def __ne__(self, other):
v = _to_rcfnum(other, self.ctx)
return Z3_rcf_neq(self.ctx_ref(), self.num, v.num)

View file

@ -106,3 +106,6 @@ class FuncEntryObj(ctypes.c_void_p):
def __init__(self, e): self._as_parameter_ = e
def from_param(obj): return obj
class RCFNumObj(ctypes.c_void_p):
def __init__(self, e): self._as_parameter_ = e
def from_param(obj): return obj

View file

@ -26,6 +26,7 @@ Notes:
#include"z3_api.h"
#include"z3_algebraic.h"
#include"z3_polynomial.h"
#include"z3_rcf.h"
#undef __in
#undef __out

View file

@ -47,6 +47,7 @@ DEFINE_TYPE(Z3_func_interp);
#define Z3_func_interp_opt Z3_func_interp
DEFINE_TYPE(Z3_func_entry);
DEFINE_TYPE(Z3_fixedpoint);
DEFINE_TYPE(Z3_rcf_num);
DEFINE_VOID(Z3_theory_data);
#endif
@ -1190,6 +1191,7 @@ typedef enum
def_Type('FUNC_ENTRY', 'Z3_func_entry', 'FuncEntryObj')
def_Type('FIXEDPOINT', 'Z3_fixedpoint', 'FixedpointObj')
def_Type('PARAM_DESCRS', 'Z3_param_descrs', 'ParamDescrs')
def_Type('RCF_NUM', 'Z3_rcf_num', 'RCFNumObj')
*/
#ifdef Conly

191
src/api/z3_rcf.h Normal file
View file

@ -0,0 +1,191 @@
/*++
Copyright (c) 2013 Microsoft Corporation
Module Name:
z3_rcf.h
Abstract:
Additional APIs for handling elements of the Z3 real closed field that contains:
- transcendental extensions
- infinitesimal extensions
- algebraic extensions
Author:
Leonardo de Moura (leonardo) 2012-01-05
Notes:
--*/
#ifndef _Z3_RCF_H_
#define _Z3_RCF_H_
#ifdef __cplusplus
extern "C" {
#endif // __cplusplus
/**
\brief Delete a RCF numeral created using the RCF API.
def_API('Z3_rcf_del', VOID, (_in(CONTEXT), _in(RCF_NUM)))
*/
void Z3_API Z3_rcf_del(__in Z3_context c, __in Z3_rcf_num a);
/**
\brief Return a RCF rational using the given string.
def_API('Z3_rcf_mk_rational', RCF_NUM, (_in(CONTEXT), _in(STRING)))
*/
Z3_rcf_num Z3_API Z3_rcf_mk_rational(__in Z3_context c, __in Z3_string val);
/**
\brief Return a RCF small integer.
def_API('Z3_rcf_mk_small_int', RCF_NUM, (_in(CONTEXT), _in(INT)))
*/
Z3_rcf_num Z3_API Z3_rcf_mk_small_int(__in Z3_context c, __in int val);
/**
\brief Return Pi
def_API('Z3_rcf_mk_pi', RCF_NUM, (_in(CONTEXT),))
*/
Z3_rcf_num Z3_API Z3_rcf_mk_pi(__in Z3_context c);
/**
\brief Return e (Euler's constant)
def_API('Z3_rcf_mk_e', RCF_NUM, (_in(CONTEXT),))
*/
Z3_rcf_num Z3_API Z3_rcf_mk_e(__in Z3_context c);
/**
\brief Return a new infinitesimal that is smaller than all elements in the Z3 field.
def_API('Z3_rcf_mk_infinitesimal', RCF_NUM, (_in(CONTEXT), _in(STRING)))
*/
Z3_rcf_num Z3_API Z3_rcf_mk_infinitesimal(__in Z3_context c, __in Z3_string name);
/**
\brief Store in roots the roots of the polynomial <tt>a[n-1]*x^{n-1} + ... + a[0]</tt>.
The output vector \c roots must have size \c n.
It returns the number of roots of the polynomial.
\pre The input polynomial is not the zero polynomial.
def_API('Z3_rcf_mk_roots', UINT, (_in(CONTEXT), _in(UINT), _in_array(1, RCF_NUM), _out_array(1, RCF_NUM)))
*/
unsigned Z3_API Z3_rcf_mk_roots(__in Z3_context c, __in unsigned n, __in_ecount(n) Z3_rcf_num const a[], __out_ecount(n) Z3_rcf_num roots[]);
/**
\brief Return the value a + b.
def_API('Z3_rcf_add', RCF_NUM, (_in(CONTEXT), _in(RCF_NUM), _in(RCF_NUM)))
*/
Z3_rcf_num Z3_API Z3_rcf_add(__in Z3_context c, __in Z3_rcf_num a, __in Z3_rcf_num b);
/**
\brief Return the value a - b.
def_API('Z3_rcf_sub', RCF_NUM, (_in(CONTEXT), _in(RCF_NUM), _in(RCF_NUM)))
*/
Z3_rcf_num Z3_API Z3_rcf_sub(__in Z3_context c, __in Z3_rcf_num a, __in Z3_rcf_num b);
/**
\brief Return the value a * b.
def_API('Z3_rcf_mul', RCF_NUM, (_in(CONTEXT), _in(RCF_NUM), _in(RCF_NUM)))
*/
Z3_rcf_num Z3_API Z3_rcf_mul(__in Z3_context c, __in Z3_rcf_num a, __in Z3_rcf_num b);
/**
\brief Return the value a / b.
def_API('Z3_rcf_div', RCF_NUM, (_in(CONTEXT), _in(RCF_NUM), _in(RCF_NUM)))
*/
Z3_rcf_num Z3_API Z3_rcf_div(__in Z3_context c, __in Z3_rcf_num a, __in Z3_rcf_num b);
/**
\brief Return the value -a
def_API('Z3_rcf_neg', RCF_NUM, (_in(CONTEXT), _in(RCF_NUM)))
*/
Z3_rcf_num Z3_API Z3_rcf_neg(__in Z3_context c, __in Z3_rcf_num a);
/**
\brief Return the value 1/a
def_API('Z3_rcf_inv', RCF_NUM, (_in(CONTEXT), _in(RCF_NUM)))
*/
Z3_rcf_num Z3_API Z3_rcf_inv(__in Z3_context c, __in Z3_rcf_num a);
/**
\brief Return the value a^k
def_API('Z3_rcf_power', RCF_NUM, (_in(CONTEXT), _in(RCF_NUM), _in(UINT)))
*/
Z3_rcf_num Z3_API Z3_rcf_power(__in Z3_context c, __in Z3_rcf_num a, __in unsigned k);
/**
\brief Return Z3_TRUE if a < b
def_API('Z3_rcf_lt', BOOL, (_in(CONTEXT), _in(RCF_NUM), _in(RCF_NUM)))
*/
Z3_bool Z3_API Z3_rcf_lt(__in Z3_context c, __in Z3_rcf_num a, __in Z3_rcf_num b);
/**
\brief Return Z3_TRUE if a > b
def_API('Z3_rcf_gt', BOOL, (_in(CONTEXT), _in(RCF_NUM), _in(RCF_NUM)))
*/
Z3_bool Z3_API Z3_rcf_gt(__in Z3_context c, __in Z3_rcf_num a, __in Z3_rcf_num b);
/**
\brief Return Z3_TRUE if a <= b
def_API('Z3_rcf_le', BOOL, (_in(CONTEXT), _in(RCF_NUM), _in(RCF_NUM)))
*/
Z3_bool Z3_API Z3_rcf_le(__in Z3_context c, __in Z3_rcf_num a, __in Z3_rcf_num b);
/**
\brief Return Z3_TRUE if a >= b
def_API('Z3_rcf_ge', BOOL, (_in(CONTEXT), _in(RCF_NUM), _in(RCF_NUM)))
*/
Z3_bool Z3_API Z3_rcf_ge(__in Z3_context c, __in Z3_rcf_num a, __in Z3_rcf_num b);
/**
\brief Return Z3_TRUE if a == b
def_API('Z3_rcf_eq', BOOL, (_in(CONTEXT), _in(RCF_NUM), _in(RCF_NUM)))
*/
Z3_bool Z3_API Z3_rcf_eq(__in Z3_context c, __in Z3_rcf_num a, __in Z3_rcf_num b);
/**
\brief Return Z3_TRUE if a != b
def_API('Z3_rcf_neq', BOOL, (_in(CONTEXT), _in(RCF_NUM), _in(RCF_NUM)))
*/
Z3_bool Z3_API Z3_rcf_neq(__in Z3_context c, __in Z3_rcf_num a, __in Z3_rcf_num b);
/**
\brief Convert the RCF numeral into a string.
def_API('Z3_rcf_num_to_string', STRING, (_in(CONTEXT), _in(RCF_NUM)))
*/
Z3_string Z3_API Z3_rcf_num_to_string(__in Z3_context c, __in Z3_rcf_num a);
/**
\brief Convert the RCF numeral into a string in decimal notation.
def_API('Z3_rcf_num_to_decimal_string', STRING, (_in(CONTEXT), _in(RCF_NUM), _in(UINT)))
*/
Z3_string Z3_API Z3_rcf_num_to_decimal_string(__in Z3_context c, __in Z3_rcf_num a, __in unsigned prec);
#ifdef __cplusplus
};
#endif // __cplusplus
#endif

View file

@ -105,10 +105,11 @@ struct interval_deps {
template<typename C>
class interval_manager {
public:
typedef typename C::numeral_manager numeral_manager;
typedef typename numeral_manager::numeral numeral;
typedef typename C::interval interval;
private:
C m_c;
numeral m_result_lower;
numeral m_result_upper;
@ -199,6 +200,11 @@ public:
bool eq(interval const & a, interval const & b) const;
/**
\brief Return true if all values in 'a' are less than all values in 'b'.
*/
bool before(interval const & a, interval const & b) const;
/**
\brief Set lower bound to -oo.
*/
@ -350,4 +356,29 @@ public:
void e(unsigned k, interval & r);
};
template<typename Manager>
class _scoped_interval {
public:
typedef typename Manager::interval interval;
private:
Manager & m_manager;
interval m_interval;
public:
_scoped_interval(Manager & m):m_manager(m) {}
~_scoped_interval() { m_manager.del(m_interval); }
Manager & m() const { return m_manager; }
operator interval const &() const { return m_interval; }
operator interval&() { return m_interval; }
interval const & get() const { return m_interval; }
interval & get() { return m_interval; }
interval * operator->() {
return &m_interval;
}
interval const * operator->() const {
return &m_interval;
}
};
#endif

View file

@ -687,6 +687,13 @@ bool interval_manager<C>::eq(interval const & a, interval const & b) const {
upper_is_open(a) == upper_is_open(b);
}
template<typename C>
bool interval_manager<C>::before(interval const & a, interval const & b) const {
if (upper_is_inf(a) || lower_is_inf(b))
return false;
return m().lt(upper(a), lower(b)) || (upper_is_open(a) && m().eq(upper(a), lower(b)));
}
template<typename C>
void interval_manager<C>::neg_jst(interval const & a, interval_deps & b_deps) {
if (lower_is_inf(a)) {

View file

@ -0,0 +1,420 @@
/*++
Copyright (c) 2013 Microsoft Corporation
Module Name:
mpz_matrix.h
Abstract:
Matrix with integer coefficients. This is not a general purpose
module for handling matrices with integer coefficients. Instead,
it is a custom package that only contains operations needed to
implement Sign Determination (Algorithm 10.11) in the Book:
"Algorithms in real algebraic geometry", Basu, Pollack, Roy
Design choices:
- Dense representation. The matrices in Alg 10.11 are small and dense.
- Integer coefficients instead of rational coefficients (it only complicates the solver a little bit).
Remark: in Algorithm 10.11, the coefficients of the input matrices are always in {-1, 0, 1}.
During solving, bigger coefficients are produced, but they are usually very small. It may be
an overkill to use mpz instead of int. We use mpz just to be safe.
Remark: We do not use rational arithmetic. The solver is slightly more complicated with integers, but is saves space.
Author:
Leonardo (leonardo) 2013-01-07
Notes:
--*/
#include"mpz_matrix.h"
#include"buffer.h"
mpz_matrix_manager::mpz_matrix_manager(unsynch_mpz_manager & nm, small_object_allocator & a):
m_nm(nm),
m_allocator(a) {
}
mpz_matrix_manager::~mpz_matrix_manager() {
}
void mpz_matrix_manager::mk(unsigned m, unsigned n, mpz_matrix & A) {
SASSERT(m > 0 && n > 0);
del(A);
A.m = m;
A.n = n;
A.a_ij = new (m_allocator) mpz[m*n];
}
void mpz_matrix_manager::del(mpz_matrix & A) {
if (A.a_ij != 0) {
for (unsigned i = 0; i < A.m; i++)
for (unsigned j = 0; j < A.n; j++)
nm().del(A(i,j));
unsigned sz = sizeof(mpz) * A.m * A.n;
m_allocator.deallocate(sz, A.a_ij);
A.m = 0;
A.n = 0;
A.a_ij = 0;
}
}
void mpz_matrix_manager::set(mpz_matrix & A, mpz_matrix const & B) {
if (&A == &B)
return;
if (A.m != B.m || A.n != B.n) {
del(A);
mk(B.m, B.n, A);
}
SASSERT(A.m == B.m && A.n == B.n);
for (unsigned i = 0; i < B.m; i++)
for (unsigned j = 0; j < B.n; j++)
nm().set(A(i, j), B(i, j));
}
void mpz_matrix_manager::tensor_product(mpz_matrix const & A, mpz_matrix const & B, mpz_matrix & C) {
scoped_mpz_matrix _C(*this);
mk(A.m * B.m, A.n * B.n, _C);
for (unsigned i = 0; i < _C.m(); i++)
for (unsigned j = 0; j < _C.n(); j++)
nm().mul(A(i / B.m, j / B.n),
B(i % B.m, j % B.n),
_C(i, j));
C.swap(_C);
}
void mpz_matrix_manager::swap_rows(mpz_matrix & A, unsigned i, unsigned j) {
if (i != j) {
for (unsigned k = 0; k < A.n; k++)
::swap(A(i, k), A(j, k));
}
}
// If b_i == 0, then method just divides the given row by its GCD
// If b_i != 0
// If the GCD of the row divides *b_i
// divide the row and *b_i by the GCD
// Else
// If int_solver == true ==> return false (the system is unsolvable)
bool mpz_matrix_manager::normalize_row(mpz * A_i, unsigned n, mpz * b_i, bool int_solver) {
scoped_mpz g(nm());
bool first = true;
for (unsigned j = 0; j < n; j++) {
if (nm().is_zero(A_i[j]))
continue;
if (first) {
nm().set(g, A_i[j]);
nm().abs(g);
first = false;
}
else {
nm().gcd(g, A_i[j], g);
}
if (nm().is_one(g))
return true;
}
if (first)
return true; // zero row
if (!nm().is_one(g)) {
if (b_i) {
if (nm().divides(g, *b_i)) {
for (unsigned j = 0; j < n; j++) {
nm().div(A_i[j], g, A_i[j]);
}
nm().div(*b_i, g, *b_i);
}
else {
if (int_solver)
return false; // system does not have an integer solution
}
}
else {
for (unsigned j = 0; j < n; j++) {
nm().div(A_i[j], g, A_i[j]);
}
}
}
return true;
}
/*
Given a matrix of the form
k2
|
V
X X ... X X ... X
0 X ... X X ... X
... ... X X ... X
k1=> 0 0 ... 0 X ... X
0 0 ... 0 X ... X
... ... 0 X ... X
0 0 ... 0 X ... X
It will "zero" the elements a_{k1+1, k2} ... a_{m, k2} by addining multiples of the row k1 to multiples of the
rows k1+1, ..., m
The resultant matrix will look like
k2
|
V
X X ... X X ... X
0 X ... X X ... X
... ... X X ... X
k1=> 0 0 ... 0 X ... X
0 0 ... 0 0 ... X
... ... 0 0 ... X
0 0 ... 0 0 ... X
If b != 0, then the transformations are also applied to b.
If int_solver == true and b != 0, then the method returns false if when
performing the transformations it detected that it is impossible to
solve the integer system of equations A x = b.
*/
bool mpz_matrix_manager::eliminate(mpz_matrix & A, mpz * b, unsigned k1, unsigned k2, bool int_solver) {
// check if first k2-1 positions of row k1 are 0
DEBUG_CODE(for (unsigned j = 0; j < k2; j++) { SASSERT(nm().is_zero(A(k1, j))); });
mpz & a_kk = A(k1, k2);
SASSERT(!nm().is_zero(a_kk));
scoped_mpz t1(nm()), t2(nm());
scoped_mpz a_ik_prime(nm()), a_kk_prime(nm()), lcm_a_kk_a_ik(nm());
// for all rows below pivot
for (unsigned i = k1+1; i < A.m; i++) {
// check if first k-1 positions of row k are 0
DEBUG_CODE(for (unsigned j = 0; j < k2; j++) { SASSERT(nm().is_zero(A(i, j))); });
mpz & a_ik = A(i, k2);
if (!nm().is_zero(a_ik)) {
// a_ik' = lcm(a_kk, a_ik)/a_kk
// a_kk' = lcm(a_kk, a_ik)/a_ik
nm().lcm(a_kk, a_ik, lcm_a_kk_a_ik);
nm().div(lcm_a_kk_a_ik, a_kk, a_ik_prime);
nm().div(lcm_a_kk_a_ik, a_ik, a_kk_prime);
for (unsigned j = k2+1; j < A.n; j++) {
// a_ij <- a_kk' * a_ij - a_ik' * a_kj
nm().mul(a_ik_prime, A(k1, j), t1);
nm().mul(a_kk_prime, A(i, j), t2);
nm().sub(t2, t1, A(i, j));
}
if (b) {
// b_i <- a_kk' * b_i - a_ik' * b_k
nm().mul(a_ik_prime, b[k1], t1);
nm().mul(a_kk_prime, b[i], t2);
nm().sub(t2, t1, b[i]);
}
// a_ik <- 0
nm().set(A(i, k2), 0);
// normalize row i
if (!normalize_row(A.row(i), A.n, b ? &(b[i]) : 0, int_solver))
return false;
}
SASSERT(nm().is_zero(A(i, k2)));
}
return true;
}
bool mpz_matrix_manager::solve_core(mpz_matrix const & _A, mpz * b, bool int_solver) {
SASSERT(_A.n == _A.m);
scoped_mpz_matrix A(*this);
set(A, _A);
for (unsigned k = 0; k < A.m(); k++) {
TRACE("mpz_matrix",
tout << "k: " << k << "\n" << A;
tout << "b:";
for (unsigned i = 0; i < A.m(); i++) {
tout << " ";
nm().display(tout, b[i]);
}
tout << "\n";);
// find pivot
unsigned i = k;
for (; i < A.m(); i++) {
if (!nm().is_zero(A(i, k)))
break;
}
if (i == A.m())
return false; // matrix is singular
// swap rows k and i
swap_rows(A, k, i);
swap(b[k], b[i]);
//
if (!eliminate(A, b, k, k, int_solver))
return false;
}
// Back substitution
unsigned k = A.m();
while (k > 0) {
--k;
DEBUG_CODE(for (unsigned j = 0; j < A.n(); j++) { SASSERT(j == k || nm().is_zero(A(k, j))); });
SASSERT(!nm().is_zero(A(k, k)));
if (nm().divides(A(k, k), b[k])) {
nm().div(b[k], A(k, k), b[k]);
nm().set(A(k, k), 1);
}
else {
if (int_solver)
return false; // no integer solution
if (nm().is_neg(A(k, k))) {
nm().neg(A(k, k));
nm().neg(b[k]);
}
}
if (!int_solver) {
// REMARK:
// For the sign determination algorithm, we only use int_solver == true.
//
// TODO: implement backward substitution when int_solver == false
// In this case, A(k, k) may not be 1.
NOT_IMPLEMENTED_YET();
}
SASSERT(!int_solver || nm().is_one(A(k, k)));
// back substitute
unsigned i = k;
while (i > 0) {
--i;
// Assuming int_solver == true
SASSERT(int_solver); // See comment above
// b_i <- b_i - a_ik * b_k
nm().submul(b[i], A(i, k), b[k], b[i]);
nm().set(A(i, k), 0);
}
}
return true;
}
bool mpz_matrix_manager::solve(mpz_matrix const & A, mpz * b, mpz const * c) {
for (unsigned i = 0; i < A.n; i++)
nm().set(b[i], c[i]);
return solve_core(A, b, true);
}
bool mpz_matrix_manager::solve(mpz_matrix const & A, int * b, int const * c) {
scoped_mpz_matrix _b(*this);
mk(A.n, 1, _b);
for (unsigned i = 0; i < A.n; i++)
nm().set(_b(i,0), c[i]);
bool r = solve_core(A, _b.A.a_ij, true);
if (r) {
for (unsigned i = 0; i < A.n; i++)
b[i] = _b.get_int(i, 0);
}
return r;
}
void mpz_matrix_manager::filter_cols(mpz_matrix const & A, unsigned num_cols, unsigned const * cols, mpz_matrix & B) {
SASSERT(num_cols <= A.n);
// Check pre-condition:
// - All elements in cols are smaller than A.n
// - cols is sorted
// - cols does not contain repeated elements
DEBUG_CODE({
for (unsigned i = 0; i < num_cols; i ++) {
SASSERT(cols[i] < A.n);
SASSERT(i == 0 || cols[i-1] < cols[i]);
}
});
if (num_cols == A.n) {
// keep everything
set(B, A);
}
else {
SASSERT(num_cols < A.n);
scoped_mpz_matrix C(*this);
mk(A.m, num_cols, C);
for (unsigned i = 0; i < A.m; i++)
for (unsigned j = 0; j < num_cols; j++)
nm().set(C(i, j), A(i, cols[j]));
B.swap(C);
}
}
void mpz_matrix_manager::permute_rows(mpz_matrix const & A, unsigned const * p, mpz_matrix & B) {
// Check if p is really a permutation
DEBUG_CODE({
buffer<bool> seen;
seen.resize(A.m, false);
for (unsigned i = 0; i < A.m; i++) {
SASSERT(p[i] < A.m);
SASSERT(!seen[p[i]]);
seen[p[i]] = true;
}
});
scoped_mpz_matrix C(*this);
mk(A.m, A.n, C);
for (unsigned i = 0; i < A.m; i++)
for (unsigned j = 0; j < A.n; j++)
nm().set(C(i, j), A(p[i], j));
B.swap(C);
}
unsigned mpz_matrix_manager::linear_independent_rows(mpz_matrix const & _A, unsigned * r, mpz_matrix & B) {
unsigned r_sz = 0;
scoped_mpz_matrix A(*this);
scoped_mpz g(nm());
scoped_mpz t1(nm()), t2(nm());
scoped_mpz a_ik_prime(nm()), a_kk_prime(nm()), lcm_a_kk_a_ik(nm());
set(A, _A);
sbuffer<unsigned, 128> rows;
rows.resize(A.m(), 0);
for (unsigned i = 0; i < A.m(); i++)
rows[i] = i;
for (unsigned k1 = 0, k2 = 0; k1 < A.m(); k1++) {
TRACE("mpz_matrix", tout << "k1: " << k1 << ", k2: " << k2 << "\n" << A;);
// find pivot
unsigned pivot = UINT_MAX;
for (unsigned i = k1; i < A.m(); i++) {
if (!nm().is_zero(A(i, k2))) {
if (pivot == UINT_MAX) {
pivot = i;
}
else {
if (rows[i] < rows[pivot])
pivot = i;
}
}
}
if (pivot == UINT_MAX)
continue;
// swap rows k and pivot
swap_rows(A, k1, pivot);
std::swap(rows[k1], rows[pivot]);
//
r[r_sz] = rows[k1];
r_sz++;
if (r_sz >= A.n())
break;
eliminate(A, 0, k1, k2, false);
k2++;
}
std::sort(r, r + r_sz);
// Copy linear independent rows to B
mpz_matrix & C = A;
mk(r_sz, _A.n, C);
for (unsigned i = 0; i < r_sz; i++ ) {
for (unsigned j = 0; j < _A.n; j++) {
nm().set(C(i, j), _A(r[i], j));
}
}
B.swap(C);
return r_sz;
}
void mpz_matrix_manager::display(std::ostream & out, mpz_matrix const & A, unsigned cell_width) const {
out << A.m << " x " << A.n << " Matrix\n";
for (unsigned i = 0; i < A.m; i++) {
for (unsigned j = 0; j < A.n; j++) {
if (j > 0)
out << " ";
std::string s = nm().to_string(A(i, j));
if (s.size() < cell_width) {
unsigned space = cell_width - s.size();
for (unsigned k = 0; k < space; k++)
out << " ";
}
out << s;
}
out << "\n";
}
}

View file

@ -0,0 +1,156 @@
/*++
Copyright (c) 2013 Microsoft Corporation
Module Name:
mpz_matrix.h
Abstract:
Matrix with integer coefficients. This is not a general purpose
module for handling matrices with integer coefficients. Instead,
it is a custom package that only contains operations needed to
implement Sign Determination (Algorithm 10.11) in the Book:
"Algorithms in real algebraic geometry", Basu, Pollack, Roy
Design choices:
- Dense representation. The matrices in Alg 10.11 are small and dense.
- Integer coefficients instead of rational coefficients (it only complicates the solver a little bit).
Remark: in Algorithm 10.11, the coefficients of the input matrices are always in {-1, 0, 1}.
During solving, bigger coefficients are produced, but they are usually very small. It may be
an overkill to use mpz instead of int. We use mpz just to be safe.
Remark: We do not use rational arithmetic. The solver is slightly more complicated with integers, but is saves space.
Author:
Leonardo (leonardo) 2013-01-07
Notes:
--*/
#ifndef _MPZ_MATRIX_H_
#define _MPZ_MATRIX_H_
#include"mpz.h"
/**
\brief A mxn matrix.
Remark: Algorithm 10.11 only uses square matrices, but supporting
arbitrary matrices does not increase the complexity of this module.
*/
class mpz_matrix {
friend class mpz_matrix_manager;
friend class scoped_mpz_matrix;
unsigned m;
unsigned n;
mpz * a_ij;
public:
mpz_matrix():m(0), n(0), a_ij(0) {}
mpz const & operator()(unsigned i, unsigned j) const {
SASSERT(i < m);
SASSERT(j < n);
return a_ij[i*n + j]; }
mpz & operator()(unsigned i, unsigned j) { SASSERT(i < m); SASSERT(j < n); return a_ij[i*n + j]; }
void swap(mpz_matrix & B) { std::swap(m, B.m); std::swap(n, B.n); std::swap(a_ij, B.a_ij); }
mpz * row(unsigned i) const { SASSERT(i < m); return a_ij + i*n; }
};
class mpz_matrix_manager {
unsynch_mpz_manager & m_nm;
small_object_allocator & m_allocator;
static void swap_rows(mpz_matrix & A, unsigned i, unsigned j);
bool normalize_row(mpz * A_i, unsigned n, mpz * b_i, bool int_solver);
bool eliminate(mpz_matrix & A, mpz * b, unsigned k1, unsigned k2, bool int_solver);
bool solve_core(mpz_matrix const & A, mpz * b, bool int_solver);
public:
mpz_matrix_manager(unsynch_mpz_manager & nm, small_object_allocator & a);
~mpz_matrix_manager();
unsynch_mpz_manager & nm() const { return m_nm; }
void mk(unsigned m, unsigned n, mpz_matrix & A);
void del(mpz_matrix & r);
void set(mpz_matrix & A, mpz_matrix const & B);
void tensor_product(mpz_matrix const & A, mpz_matrix const & B, mpz_matrix & C);
/**
\brief Solve A*b = c
Return false if the system does not have an integer solution.
\pre A is a square matrix
\pre b and c are vectors of size A.n (== A.m)
*/
bool solve(mpz_matrix const & A, mpz * b, mpz const * c);
/**
\brief Solve A*b = c
Return false if the system does not have an integer solution.
\pre A is a square matrix
\pre b and c are vectors of size A.n (== A.m)
*/
bool solve(mpz_matrix const & A, int * b, int const * c);
/**
\brief Store in B that contains the subset cols of columns of A.
\pre num_cols <= A.n
\pre Forall i < num_cols, cols[i] < A.n
\pre Forall 0 < i < num_cols, cols[i-1] < cols[i]
*/
void filter_cols(mpz_matrix const & A, unsigned num_cols, unsigned const * cols, mpz_matrix & B);
/**
\brief Store in B the matrix obtained after applying the given permutation to the rows of A.
*/
void permute_rows(mpz_matrix const & A, unsigned const * p, mpz_matrix & B);
/**
\brief Store in r the row (ids) of A that are linear independent.
\remark If there is an option between rows i and j,
this method will give preference to the row that occurs first.
\remark The vector r must have at least A.n() capacity
The numer of linear independent rows is returned.
Store the new matrix in B.
*/
unsigned linear_independent_rows(mpz_matrix const & A, unsigned * r, mpz_matrix & B);
// method for debugging purposes
void display(std::ostream & out, mpz_matrix const & A, unsigned cell_width=4) const;
};
class scoped_mpz_matrix {
friend class mpz_matrix_manager;
mpz_matrix_manager & m_manager;
mpz_matrix A;
public:
scoped_mpz_matrix(mpz_matrix_manager & m):m_manager(m) {}
scoped_mpz_matrix(mpz_matrix const & B, mpz_matrix_manager & m):m_manager(m) { m_manager.set(A, B); }
~scoped_mpz_matrix() { m_manager.del(A); }
mpz_matrix_manager & mm() const { return m_manager; }
unsynch_mpz_manager & nm() const { return mm().nm(); }
unsigned m() const { return A.m; }
unsigned n() const { return A.n; }
mpz * row(unsigned i) const { return A.row(i); }
operator mpz_matrix const &() const { return A; }
operator mpz_matrix &() { return A; }
mpz_matrix const & get() const { return A; }
mpz_matrix & get() { return A; }
void swap(mpz_matrix & B) { A.swap(B); }
void set(unsigned i, unsigned j, mpz const & v) { nm().set(A(i, j), v); }
void set(unsigned i, unsigned j, int v) { nm().set(A(i, j), v); }
mpz const & operator()(unsigned i, unsigned j) const { return A(i, j); }
mpz & operator()(unsigned i, unsigned j) { return A(i, j); }
int get_int(unsigned i, unsigned j) const { SASSERT(nm().is_int(A(i, j))); return nm().get_int(A(i, j)); }
};
inline std::ostream & operator<<(std::ostream & out, scoped_mpz_matrix const & m) {
m.mm().display(out, m);
return out;
}
#endif

View file

@ -0,0 +1,8 @@
def_module_params('rcf',
description='real closed fields',
export=True,
params=(('use_prem', BOOL, True, "use pseudo-remainder instead of remainder when computing GCDs and Sturm-Tarski sequences"),
('clean_denominators', BOOL, True, "clean denominators before root isolation"),
('initial_precision', UINT, 24, "a value k that is the initial interval size (as 1/2^k) when creating transcendentals and approximated division"),
('inf_precision', UINT, 24, "a value k that is the initial interval size (i.e., (0, 1/2^l)) used as an approximation for infinitesimal values"),
('max_precision', UINT, 64, "during sign determination we switch from interval arithmetic to complete methods when the interval size is less than 1/2^k, where k is the max_precision")))

File diff suppressed because it is too large Load diff

View file

@ -0,0 +1,406 @@
/*++
Copyright (c) 2013 Microsoft Corporation
Module Name:
realclosure.h
Abstract:
Package for computing with elements of the realclosure of a field containing
- all rationals
- extended with computable transcendental real numbers (e.g., pi and e)
- infinitesimals
Author:
Leonardo (leonardo) 2013-01-02
Notes:
--*/
#ifndef _REALCLOSURE_H_
#define _REALCLOSURE_H_
#include"mpq.h"
#include"params.h"
#include"scoped_numeral.h"
#include"scoped_numeral_vector.h"
#include"interval.h"
#include"z3_exception.h"
namespace realclosure {
class num;
typedef interval_manager<im_default_config> mpqi_manager;
typedef default_exception exception;
class mk_interval {
public:
virtual void operator()(unsigned k, mpqi_manager & im, mpqi_manager::interval & r) = 0;
};
class manager {
public:
struct imp;
private:
friend class save_interval_ctx;
imp * m_imp;
public:
manager(unsynch_mpq_manager & m, params_ref const & p = params_ref(), small_object_allocator * a = 0);
~manager();
typedef num numeral;
typedef svector<numeral> numeral_vector;
typedef _scoped_numeral<manager> scoped_numeral;
typedef _scoped_numeral_vector<manager> scoped_numeral_vector;
static void get_param_descrs(param_descrs & r);
static void collect_param_descrs(param_descrs & r) { get_param_descrs(r); }
void set_cancel(bool f);
void cancel() { set_cancel(true); }
void reset_cancel() { set_cancel(false); }
void updt_params(params_ref const & p);
unsynch_mpq_manager & qm() const;
void del(numeral & a);
/**
\brief Add a new infinitesimal to the current field. The new infinitesimal is smaller than any positive element in the field.
*/
void mk_infinitesimal(char const * name, numeral & r);
void mk_infinitesimal(numeral & r);
/**
\brief Add a new transcendental real value to the field.
The functor \c mk_interval is used to compute approximations of the transcendental value.
This procedure should be used with care, if the value is not really transcendental with respect to the current
field, computations with the new numeral may not terminate.
Example: we extended the field with Pi. Pi is transcendental with respect to a field that contains only algebraic real numbers.
So, this step is fine. Let us call the resultant field F.
Then, we extend the field F with 1 - Pi. 1 - Pi is transcendental with respect to algebraic real numbers, but it is NOT transcendental
with respect to F, since F contains Pi.
*/
void mk_transcendental(char const * name, mk_interval & proc, numeral & r);
void mk_transcendental(mk_interval & proc, numeral & r);
/**
\brief r <- pi
*/
void mk_pi(numeral & r);
/**
\brief r <- e (Euler's constant)
*/
void mk_e(numeral & r);
/**
\brief Isolate the roots of the univariate polynomial as[0] + as[1]*x + ... + as[n-1]*x^{n-1}
The roots are stored in \c roots.
*/
void isolate_roots(unsigned n, numeral const * as, numeral_vector & roots);
/**
\brief a <- 0
*/
void reset(numeral & a);
/**
\brief Return the sign of a.
*/
int sign(numeral const & a);
/**
\brief Return true if a is zero.
*/
bool is_zero(numeral const & a);
/**
\brief Return true if a is positive.
*/
bool is_pos(numeral const & a);
/**
\brief Return true if a is negative.
*/
bool is_neg(numeral const & a);
/**
\brief Return true if a is an integer.
*/
bool is_int(numeral const & a);
/**
\brief Return true if the representation of \c a depends on
infinitesimal extensions.
*/
bool depends_on_infinitesimals(numeral const & a);
/**
\brief a <- n
*/
void set(numeral & a, int n);
void set(numeral & a, mpz const & n);
void set(numeral & a, mpq const & n);
void set(numeral & a, numeral const & n);
void swap(numeral & a, numeral & b);
/**
\brief Return a^{1/k}
Throws an exception if (a is negative and k is even) or (k is zero).
*/
void root(numeral const & a, unsigned k, numeral & b);
/**
\brief Return a^k
Throws an exception if 0^0.
*/
void power(numeral const & a, unsigned k, numeral & b);
/**
\brief c <- a + b
*/
void add(numeral const & a, numeral const & b, numeral & c);
void add(numeral const & a, mpz const & b, numeral & c);
/**
\brief c <- a - b
*/
void sub(numeral const & a, numeral const & b, numeral & c);
/**
\brief c <- a * b
*/
void mul(numeral const & a, numeral const & b, numeral & c);
/**
\brief a <- -a
*/
void neg(numeral & a);
/**
\brief b <- -a
*/
void neg(numeral const & a, numeral & b);
/**
\brief a <- 1/a if a != 0
*/
void inv(numeral & a);
/**
\brief b <- 1/a if a != 0
*/
void inv(numeral const & a, numeral & b);
/**
\brief c <- a/b if b != 0
*/
void div(numeral const & a, numeral const & b, numeral & c);
/**
Return -1 if a < b
Return 0 if a == b
Return 1 if a > b
*/
int compare(numeral const & a, numeral const & b);
/**
\brief a == b
*/
bool eq(numeral const & a, numeral const & b);
bool eq(numeral const & a, mpq const & b);
bool eq(numeral const & a, mpz const & b);
/**
\brief a != b
*/
bool neq(numeral const & a, numeral const & b) { return !eq(a, b); }
bool neq(numeral const & a, mpq const & b) { return !eq(a, b); }
bool neq(numeral const & a, mpz const & b) { return !eq(a, b); }
/**
\brief a < b
*/
bool lt(numeral const & a, numeral const & b);
bool lt(numeral const & a, mpq const & b);
bool lt(numeral const & a, mpz const & b);
/**
\brief a > b
*/
bool gt(numeral const & a, numeral const & b) { return lt(b, a); }
bool gt(numeral const & a, mpq const & b);
bool gt(numeral const & a, mpz const & b);
/**
\brief a <= b
*/
bool le(numeral const & a, numeral const & b) { return !gt(a, b); }
bool le(numeral const & a, mpq const & b) { return !gt(a, b); }
bool le(numeral const & a, mpz const & b) { return !gt(a, b); }
/**
\brief a >= b
*/
bool ge(numeral const & a, numeral const & b) { return !lt(a, b); }
bool ge(numeral const & a, mpq const & b) { return !lt(a, b); }
bool ge(numeral const & a, mpz const & b) { return !lt(a, b); }
void display(std::ostream & out, numeral const & a) const;
/**
\brief Display a real number in decimal notation.
A question mark is added based on the precision requested.
This procedure throws an exception if the \c a is not a real.
*/
void display_decimal(std::ostream & out, numeral const & a, unsigned precision = 10) const;
void display_interval(std::ostream & out, numeral const & a) const;
void clean_denominators(numeral const & a, numeral & p, numeral & q);
};
class value;
class num {
friend class manager;
friend class manager::imp;
value * m_value;
public:
num():m_value(0) {}
// Low level functions for implementing the C API
void * c_ptr() { return m_value; }
static num mk(void * ptr) { num r; r.m_value = reinterpret_cast<value*>(ptr); return r; }
};
};
typedef realclosure::manager rcmanager;
typedef rcmanager::numeral rcnumeral;
typedef rcmanager::numeral_vector rcnumeral_vector;
typedef rcmanager::scoped_numeral scoped_rcnumeral;
typedef rcmanager::scoped_numeral_vector scoped_rcnumeral_vector;
#define RCF_MK_COMPARISON_CORE(EXTERNAL, INTERNAL, TYPE) \
inline bool EXTERNAL(scoped_rcnumeral const & a, TYPE const & b) { \
rcmanager & m = a.m(); \
scoped_rcnumeral _b(m); \
m.set(_b, b); \
return m.INTERNAL(a, _b); \
}
#define RCF_MK_COMPARISON(EXTERNAL, INTERNAL) \
RCF_MK_COMPARISON_CORE(EXTERNAL, INTERNAL, int) \
RCF_MK_COMPARISON_CORE(EXTERNAL, INTERNAL, mpz) \
RCF_MK_COMPARISON_CORE(EXTERNAL, INTERNAL, mpq)
RCF_MK_COMPARISON(operator==, eq);
RCF_MK_COMPARISON(operator!=, neq);
RCF_MK_COMPARISON(operator<, lt);
RCF_MK_COMPARISON(operator<=, le);
RCF_MK_COMPARISON(operator>, gt);
RCF_MK_COMPARISON(operator>=, ge);
#undef RCF_MK_COMPARISON
#undef RCF_MK_COMPARISON_CORE
#define RCF_MK_BINARY_CORE(EXTERNAL, INTERNAL, TYPE) \
inline scoped_rcnumeral EXTERNAL(scoped_rcnumeral const & a, TYPE const & b) { \
rcmanager & m = a.m(); \
scoped_rcnumeral _b(m); \
m.set(_b, b); \
scoped_rcnumeral r(m); \
m.INTERNAL(a, _b, r); \
return r; \
}
#define RCF_MK_BINARY(EXTERNAL, INTERNAL) \
RCF_MK_BINARY_CORE(EXTERNAL, INTERNAL, int) \
RCF_MK_BINARY_CORE(EXTERNAL, INTERNAL, mpz) \
RCF_MK_BINARY_CORE(EXTERNAL, INTERNAL, mpq)
RCF_MK_BINARY(operator+, add)
RCF_MK_BINARY(operator-, sub)
RCF_MK_BINARY(operator*, mul)
RCF_MK_BINARY(operator/, div)
#undef RCF_MK_BINARY
#undef RCF_MK_BINARY_CORE
inline scoped_rcnumeral root(scoped_rcnumeral const & a, unsigned k) {
scoped_rcnumeral r(a.m());
a.m().root(a, k, r);
return r;
}
inline scoped_rcnumeral power(scoped_rcnumeral const & a, unsigned k) {
scoped_rcnumeral r(a.m());
a.m().power(a, k, r);
return r;
}
inline scoped_rcnumeral operator^(scoped_rcnumeral const & a, unsigned k) {
return power(a, k);
}
inline bool is_int(scoped_rcnumeral const & a) {
return a.m().is_int(a);
}
struct rc_sym_pp {
rcmanager & m;
rcnumeral const & n;
rc_sym_pp(rcmanager & _m, rcnumeral const & _n):m(_m), n(_n) {}
rc_sym_pp(scoped_rcnumeral const & _n):m(_n.m()), n(_n.get()) {}
};
inline rc_sym_pp sym_pp(scoped_rcnumeral const & _n) {
return rc_sym_pp(_n);
}
inline std::ostream & operator<<(std::ostream & out, rc_sym_pp const & n) {
n.m.display(out, n.n);
return out;
}
struct rc_decimal_pp {
rcmanager & m;
rcnumeral const & n;
unsigned prec;
rc_decimal_pp(rcmanager & _m, rcnumeral const & _n, unsigned p):m(_m), n(_n), prec(p) {}
rc_decimal_pp(scoped_rcnumeral const & _n, unsigned p):m(_n.m()), n(_n.get()), prec(p) {}
};
inline std::ostream & operator<<(std::ostream & out, rc_decimal_pp const & n) {
n.m.display_decimal(out, n.n, n.prec);
return out;
}
inline rc_decimal_pp decimal_pp(scoped_rcnumeral const & n, unsigned prec = 10) {
return rc_decimal_pp(n, prec);
}
struct rc_interval_pp {
rcmanager & m;
rcnumeral const & n;
rc_interval_pp(rcmanager & _m, rcnumeral const & _n):m(_m), n(_n) {}
rc_interval_pp(scoped_rcnumeral const & _n):m(_n.m()), n(_n.get()) {}
};
inline std::ostream & operator<<(std::ostream & out, rc_interval_pp const & n) {
n.m.display_interval(out, n.n);
return out;
}
inline rc_interval_pp interval_pp(scoped_rcnumeral const & n) {
return rc_interval_pp(n);
}
#endif

View file

@ -176,11 +176,11 @@ static void tst1() {
display_anums(std::cout, rs1);
}
void tst_refine_mpbq() {
void tst_refine_mpbq(int n, int d) {
unsynch_mpq_manager qm;
mpbq_manager bqm(qm);
scoped_mpq q1(qm);
qm.set(q1, 5, 7);
qm.set(q1, n, d);
scoped_mpbq l(bqm);
scoped_mpbq u(bqm);
std::cout << "using refine upper...\n";
@ -207,6 +207,10 @@ void tst_refine_mpbq() {
}
}
void tst_refine_mpbq() {
tst_refine_mpbq(5, 7);
}
void tst_mpbq_root() {
unsynch_mpq_manager qm;
mpbq_manager bqm(qm);

View file

@ -206,6 +206,7 @@ int main(int argc, char ** argv) {
TST(mpff);
TST(horn_subsume_model_converter);
TST(model2expr);
TST(rcf);
}
void initialize_mam() {}

175
src/test/rcf.cpp Normal file
View file

@ -0,0 +1,175 @@
/*++
Copyright (c) 2013 Microsoft Corporation
Module Name:
rcf.cpp
Abstract:
Testing RCF module
Author:
Leonardo (leonardo) 2013-01-04
Notes:
--*/
#include"realclosure.h"
#include"mpz_matrix.h"
static void tst1() {
unsynch_mpq_manager qm;
rcmanager m(qm);
scoped_rcnumeral a(m);
#if 0
a = 10;
std::cout << sym_pp(a) << std::endl;
std::cout << sym_pp(eps) << std::endl;
std::cout << interval_pp(a) << std::endl;
std::cout << interval_pp(eps) << std::endl;
#endif
scoped_rcnumeral eps(m);
m.mk_infinitesimal("eps", eps);
mpq aux;
qm.set(aux, 1, 3);
m.set(a, aux);
#if 0
std::cout << interval_pp(a) << std::endl;
std::cout << decimal_pp(eps, 4) << std::endl;
std::cout << decimal_pp(a) << std::endl;
std::cout << a + eps << std::endl;
std::cout << a * eps << std::endl;
std::cout << (a + eps)*eps - eps << std::endl;
#endif
std::cout << interval_pp(a - eps*2) << std::endl;
std::cout << interval_pp(eps + 1) << std::endl;
scoped_rcnumeral t(m);
t = (a - eps*2) / (eps + 1);
std::cout << t << std::endl;
std::cout << t * (eps + 1) << std::endl;
a = 10;
std::cout << (a + eps > a) << std::endl;
scoped_rcnumeral pi(m);
m.mk_pi(pi);
std::cout << pi + 1 << std::endl;
std::cout << decimal_pp(pi) << std::endl;
std::cout << decimal_pp(pi + 1) << std::endl;
scoped_rcnumeral e(m);
m.mk_e(e);
t = e + (pi + 1)*2;
std::cout << t << std::endl;
std::cout << decimal_pp(t, 10) << std::endl;
std::cout << (eps + 1 > 1) << std::endl;
std::cout << interval_pp((a + eps)/(a - eps)) << std::endl;
}
static void tst2() {
enable_trace("mpz_matrix");
unsynch_mpq_manager nm;
small_object_allocator allocator;
mpz_matrix_manager mm(nm, allocator);
scoped_mpz_matrix A(mm);
mm.mk(3, 3, A);
// Matrix
// 1 1 1
// 0 1 -1
// 0 1 1
A.set(0, 0, 1); A.set(0, 1, 1); A.set(0, 2, 1);
A.set(1, 0, 0); A.set(1, 1, 1); A.set(1, 2, -1);
A.set(2, 0, 0); A.set(2, 1, 1); A.set(2, 2, 1);
std::cout << A;
{
int b[3];
int c[3] = { 10, -2, 8 };
std::cout << "solve: " << mm.solve(A, b, c) << "\n";
for (unsigned i = 0; i < 3; i++) std::cout << b[i] << " ";
std::cout << "\n";
}
scoped_mpz_matrix A2(mm);
mm.tensor_product(A, A, A2);
std::cout << A2;
scoped_mpz_matrix B(mm);
unsigned cols[] = { 1, 3, 7, 8 };
mm.filter_cols(A2, 4, cols, B);
std::cout << B;
scoped_mpz_matrix C(mm);
unsigned perm[] = { 8, 7, 6, 5, 4, 3, 2, 1, 0 };
mm.permute_rows(B, perm, C);
std::cout << C;
}
static void tst_solve(unsigned n, int _A[], int _b[], int _c[], bool solved) {
unsynch_mpq_manager nm;
small_object_allocator allocator;
mpz_matrix_manager mm(nm, allocator);
scoped_mpz_matrix A(mm);
mm.mk(n, n, A);
for (unsigned i = 0; i < n; i++)
for (unsigned j = 0; j < n; j++)
A.set(i, j, _A[i*n + j]);
svector<int> b;
b.resize(n, 0);
if (mm.solve(A, b.c_ptr(), _c)) {
SASSERT(solved);
for (unsigned i = 0; i < n; i++) {
SASSERT(b[i] == _b[i]);
}
}
else {
SASSERT(!solved);
}
}
static void tst_lin_indep(unsigned m, unsigned n, int _A[], unsigned ex_sz, unsigned ex_r[]) {
unsynch_mpq_manager nm;
small_object_allocator allocator;
mpz_matrix_manager mm(nm, allocator);
scoped_mpz_matrix A(mm);
mm.mk(m, n, A);
for (unsigned i = 0; i < m; i++)
for (unsigned j = 0; j < n; j++)
A.set(i, j, _A[i*n + j]);
unsigned_vector r;
r.resize(A.n());
scoped_mpz_matrix B(mm);
mm.linear_independent_rows(A, r.c_ptr(), B);
SASSERT(r.size() == ex_sz);
for (unsigned i = 0; i < ex_sz; i++) {
SASSERT(r[i] == ex_r[i]);
}
}
static void tst_denominators() {
unsynch_mpq_manager qm;
rcmanager m(qm);
scoped_rcnumeral a(m);
scoped_rcnumeral t(m);
scoped_rcnumeral eps(m);
m.mk_pi(a);
m.inv(a);
m.mk_infinitesimal("eps", eps);
t = (a - eps*2) / (a*eps + 1);
// t = t + a * 2;
scoped_rcnumeral n(m), d(m);
std::cout << t << "\n";
m.clean_denominators(t, n, d);
std::cout << "---->\n" << n << "\n" << d << "\n";
}
void tst_rcf() {
enable_trace("rcf_clean");
enable_trace("rcf_clean_bug");
tst_denominators();
return;
tst1();
tst2();
{ int A[] = {0, 1, 1, 1, 0, 1, 1, 1, -1}; int c[] = {10, 4, -4}; int b[] = {-2, 4, 6}; tst_solve(3, A, b, c, true); }
{ int A[] = {1, 1, 1, 0, 1, 1, 0, 1, 1}; int c[] = {3, 2, 2}; int b[] = {1, 1, 1}; tst_solve(3, A, b, c, false); }
{ int A[] = {1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, -1}; unsigned r[] = {0, 1, 4}; tst_lin_indep(5, 3, A, 3, r); }
{ int A[] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, -1}; unsigned r[] = {0, 4}; tst_lin_indep(5, 3, A, 2, r); }
{ int A[] = {1, 1, 1, 1, 1, 1, 1, 0, 1, 2, 1, 2, 3, 1, 3}; unsigned r[] = {0, 2}; tst_lin_indep(5, 3, A, 2, r); }
}

View file

@ -53,11 +53,11 @@ private:
set_data(mem, sz);
}
void init() {
void init(T const & v) {
iterator it = begin();
iterator e = end();
for (; it != e; ++it) {
new (it) T();
new (it) T(v);
}
}
@ -122,7 +122,7 @@ public:
if (m_data) {
if (CallDestructors)
destroy_elements();
a.deallocate(size(), raw_ptr);
a.deallocate(size(), raw_ptr());
m_data = 0;
}
}
@ -137,7 +137,14 @@ public:
void set(Allocator & a, size_t sz, T const * vs) {
SASSERT(m_data == 0);
allocate(a, sz);
init(sz, vs);
init(vs);
}
template<typename Allocator>
void set(Allocator & a, size_t sz, T const & v = T()) {
SASSERT(m_data == 0);
allocate(a, sz);
init(v);
}
size_t size() const {
@ -175,7 +182,13 @@ public:
return m_data + size();
}
T const * c_ptr() const { return m_data; }
T * c_ptr() { return m_data; }
void swap(array & other) {
std::swap(m_data, other.m_data);
}
};
template<typename T>

View file

@ -231,8 +231,13 @@ public:
SASSERT(size() == nsz);
}
private:
buffer& operator=(buffer const&);
buffer & operator=(buffer const & other) {
if (this == &other)
return *this;
reset();
append(other);
return *this;
}
};
template<typename T, unsigned INITIAL_SIZE=16>

View file

@ -855,16 +855,16 @@ void mpbq_manager::approx_div(mpbq const & a, mpbq const & b, mpbq & c, unsigned
unsigned k_prime;
if (m_manager.is_power_of_two(b.m_num, k_prime)) {
// The division is precise, so we ignore k and to_plus_inf
SASSERT(b.m_k == 0); // remark: b.m_num is odd when b.m_k > 0
SASSERT(b.m_k == 0 || k_prime == 0); // remark: b.m_num is odd when b.m_k > 0, since b.m_num is a power of two we have that b.m_k == 0 or b.m_num == 1.
m_manager.set(c.m_num, a.m_num);
if (a.m_k == 0) {
c.m_k = k_prime;
normalize(c);
}
else {
c.m_k = a.m_k + k_prime;
// there is not need to normalize since the least significant bit of a must be 1.
if (b.m_k > 0) {
SASSERT(k_prime == 0);
mpz & pw2 = m_div_tmp1;
m_manager.power(mpz(2), b.m_k, pw2);
m_manager.mul(c.m_num, pw2, c.m_num);
}
c.m_k = a.m_k + k_prime;
normalize(c);
}
else if (m_manager.divides(b.m_num, a.m_num)) {
// result is also precise

View file

@ -251,8 +251,10 @@ class mpz_manager {
}
void mk_big(mpz & a) {
if (a.m_ptr == 0)
if (a.m_ptr == 0) {
a.m_val = 0;
a.m_ptr = allocate();
}
}
#endif
@ -374,7 +376,7 @@ public:
}
// d <- a + b*c
// d <- a - b*c
void submul(mpz const & a, mpz const & b, mpz const & c, mpz & d) {
if (is_one(b)) {
sub(a, c, d);
@ -676,10 +678,18 @@ public:
int64 get_int64(mpz const & a) const;
bool is_uint(mpz const & a) const { return is_uint64(a) && get_uint64(a) < UINT_MAX; }
unsigned get_uint(mpz const & a) const { SASSERT(is_uint(a)); return static_cast<unsigned>(get_uint64(a)); }
bool is_int(mpz const & a) const { return is_int64(a) && INT_MIN < get_int64(a) && get_int64(a) < INT_MAX; }
int get_int(mpz const & a) const { SASSERT(is_int(a)); return static_cast<int>(get_int64(a)); }
double get_double(mpz const & a) const;
std::string to_string(mpz const & a) const;
void display(std::ostream & out, mpz const & a) const;
/**

View file

@ -30,10 +30,10 @@ Revision History:
- void dec_ref(T * obj)
- void inc_ref(T * obj)
*/
template<typename T, typename Ref>
template<typename T, typename Ref, unsigned INITIAL_SIZE=16>
class ref_buffer_core : public Ref {
protected:
ptr_buffer<T> m_buffer;
ptr_buffer<T, INITIAL_SIZE> m_buffer;
void inc_ref(T * o) { Ref::inc_ref(o); }
void dec_ref(T * o) { Ref::dec_ref(o); }
@ -78,11 +78,7 @@ public:
return m_buffer.c_ptr();
}
T const * operator[](unsigned idx) const {
return m_buffer[idx];
}
T * operator[](unsigned idx) {
T * operator[](unsigned idx) const {
return m_buffer[idx];
}
@ -126,6 +122,11 @@ public:
m_buffer.resize(sz, 0);
}
void shrink(unsigned sz) {
SASSERT(sz <= m_buffer.size());
resize(sz);
}
// set pos idx with elem. If idx >= size, then expand.
void setx(unsigned idx, T * elem) {
if (idx >= size()) {
@ -133,19 +134,23 @@ public:
}
set(idx, elem);
}
private:
// prevent abuse:
ref_buffer_core& operator=(ref_buffer_core const & other);
ref_buffer_core & operator=(ref_buffer_core const & other) {
if (this == &other)
return *this;
reset();
append(other);
return *this;
}
};
/**
\brief Buffer of managed references
*/
template<typename T, typename TManager>
class ref_buffer : public ref_buffer_core<T, ref_manager_wrapper<T, TManager> > {
typedef ref_buffer_core<T, ref_manager_wrapper<T, TManager> > super;
template<typename T, typename TManager, unsigned INITIAL_SIZE=16>
class ref_buffer : public ref_buffer_core<T, ref_manager_wrapper<T, TManager>, INITIAL_SIZE> {
typedef ref_buffer_core<T, ref_manager_wrapper<T, TManager>, INITIAL_SIZE> super;
public:
ref_buffer(TManager & m):
super(ref_manager_wrapper<T, TManager>(m)) {
@ -161,8 +166,8 @@ public:
/**
\brief Buffer of unmanaged references
*/
template<typename T>
class sref_buffer : public ref_buffer_core<T, ref_unmanaged_wrapper<T> > {
template<typename T, unsigned INITIAL_SIZE=16>
class sref_buffer : public ref_buffer_core<T, ref_unmanaged_wrapper<T>, INITIAL_SIZE> {
public:
};

View file

@ -243,6 +243,10 @@ public:
m_ptr = 0;
return tmp;
}
void swap(scoped_ptr & p) {
std::swap(m_ptr, p.m_ptr);
}
};
template<typename T1, typename T2>