mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-04-23 17:15:31 +00:00
Code Activity

mpf partial remainder draft

Browse source
This commit is contained in:
Christoph M. Wintersteiger 2016-05-03 18:20:18 +01:00
parent 6895cc7cc6
commit a7c66356ae
2 changed files with 182 additions and 42 deletions
Download patch file Download diff file Expand all files Collapse all files

214
src/util/mpf.cpp
View file

@ -1214,6 +1214,154 @@ void mpf_manager::to_ieee_bv_mpz(const mpf & x, scoped_mpz & o) {
}
}
void mpf_manager::prem(scoped_mpf & x, scoped_mpf const & y, mpf_exp_t const & exp_diff, bool partial) {
unsigned ebits = x.get().ebits;
unsigned sbits = x.get().sbits;
bool sign = x.get().sign;
signed int D = (unsigned)(partial ? 0 : exp_diff);
mpf_exp_t N = sbits/2;
TRACE("mpf_dbg_rem", tout << "x=" << to_string(x) << std::endl;
tout << "y=" << to_string(y) << std::endl;
tout << "d=" << D << std::endl;
tout << "partial=" << partial << std::endl;);
// 1. Compute a/b
unsigned extra_bits = sbits + 2;
scoped_mpz x_sig_shifted(m_mpz_manager), x_div_y_sig_lrg(m_mpz_manager);
scoped_mpz x_div_y_sig(m_mpz_manager), x_rem_y_sig(m_mpz_manager);
m_mpz_manager.mul2k(x.significand(), sbits + extra_bits, x_sig_shifted);
m_mpz_manager.machine_div(x_sig_shifted, y.significand(), x_div_y_sig_lrg);
m_mpz_manager.machine_div_rem(x_div_y_sig_lrg, m_powers2(extra_bits-2), x_div_y_sig, x_rem_y_sig);
TRACE("mpf_dbg_rem", tout << "X/Y_sig=" << m_mpz_manager.to_string(x_div_y_sig) << std::endl;
tout << "X%Y_sig=" << m_mpz_manager.to_string(x_rem_y_sig) << std::endl;
tout << "X/Y=" << to_string_hexfloat(
to_packed_mpf(x.sign(), D, x_div_y_sig, ebits, sbits, 3)) << std::endl;);
SASSERT(m_mpz_manager.lt(x_div_y_sig, m_powers2(sbits + 3))); // 3 extra bits in x_div_y_sig.
SASSERT(m_mpz_manager.ge(x_div_y_sig, m_powers2(sbits - 1 + 3)));
// Bug: a_rem_b_sig is not used.
// 2. Round a/b to integer Q/QQ
mpf_exp_t Q_exp = partial ? N : D;
scoped_mpz Q_sig(m_mpz_manager), Q_rem(m_mpz_manager);
if (partial) {
// Round according to MPF_ROUND_TOWARD_ZERO
SASSERT(N < D && D < INT_MAX);
unsigned D_N = (unsigned)(D-N);
unsigned Q_shft = (sbits - 1) - D_N;
m_mpz_manager.machine_div_rem(x_div_y_sig, m_powers2(Q_shft+3), Q_sig, Q_rem);
m_mpz_manager.mul2k(Q_sig, Q_shft);
}
else {
// Round according to MPF_ROUND_NEAREST_TEVEN
unsigned Q_shft = (sbits - 1) - D;
m_mpz_manager.machine_div_rem(x_div_y_sig, m_powers2(Q_shft+3), Q_sig, Q_rem);
const mpz & shiftm1_p = m_powers2(Q_shft-1);
bool tie = m_mpz_manager.eq(Q_rem, shiftm1_p);
bool more_than_tie = m_mpz_manager.gt(Q_rem, shiftm1_p);
TRACE("mpf_dbg_rem", tout << "tie= " << tie << "; >tie= " << more_than_tie << std::endl;);
if ((tie && m_mpz_manager.is_odd(Q_sig)) || more_than_tie)
m_mpz_manager.inc(Q_sig);
m_mpz_manager.mul2k(Q_sig, Q_shft);
SASSERT(m_mpz_manager.ge(Q_sig, m_powers2(sbits-1)));
}
TRACE("mpf_dbg_rem", tout << "Q_sig'=" << m_mpz_manager.to_string(Q_sig) << std::endl;
tout << "Q_rem=" << m_mpz_manager.to_string(Q_rem) << std::endl; );
// re-normalize Q
while (m_mpz_manager.ge(Q_sig, m_powers2(sbits))) {
m_mpz_manager.machine_div2k(Q_sig, 1);
Q_exp++;
}
TRACE("mpf_dbg_rem", tout << "Q_exp=" << Q_exp << std::endl;
tout << "Q_sig=" << m_mpz_manager.to_string(Q_sig) << std::endl;
scoped_mpf & Q = to_packed_mpf(x.sign(), Q_exp, Q_sig, ebits, sbits, 0);
tout << "Q=" << to_string_hexfloat(Q) << std::endl;);
// 3. Compute Y*Q / Y*QQ*2^{D-N}
scoped_mpz YQ_sig(m_mpz_manager);
mpf_exp_t YQ_exp = (Q_exp + y.exponent()) + (partial ? D-N : 0);
m_mpz_manager.mul(y.significand(), Q_sig, YQ_sig);
if (sbits >= 7) {
scoped_mpz sticky_rem(m_mpz_manager);
m_mpz_manager.machine_div_rem(YQ_sig, m_powers2(sbits-7), YQ_sig, sticky_rem);
if (!m_mpz_manager.is_zero(sticky_rem) && m_mpz_manager.is_even(YQ_sig))
m_mpz_manager.inc(YQ_sig);
}
else
m_mpz_manager.mul2k(YQ_sig, 7-sbits, YQ_sig);
TRACE("mpf_dbg_rem", tout << "YQ_exp=" << YQ_exp << std::endl;
tout << "YQ_sig=" << m_mpz_manager.to_string(YQ_sig) << std::endl;
scoped_mpf & YQ = to_packed_mpf(x.sign(), YQ_exp, YQ_sig, ebits, sbits, 6);
tout << "YQ=" << to_string_hexfloat(YQ) << std::endl;);
// 4. Compute X-Y*Q
bool YQ_sgn = !x.sign();
mpf_exp_t X_YQ_exp = x.exponent() - 3;
scoped_mpz X_YQ_sig(m_mpz_manager);
mpf_exp_t exp_delta = exp(x) - YQ_exp;
TRACE("mpf_dbg_rem", tout << "exp_delta=" << exp_delta << std::endl;);
scoped_mpz minuend(m_mpz_manager), subtrahend(m_mpz_manager);
if (exp_delta >= 0) {
m_mpz_manager.set(minuend, x.significand());
m_mpz_manager.set(subtrahend, YQ_sig);
}
else if (exp_delta < 0) {
m_mpz_manager.set(minuend, YQ_sig);
m_mpz_manager.set(subtrahend, x.significand());
YQ_sgn = !YQ_sgn;
exp_delta = -exp_delta;
}
if (exp_delta != 0) {
scoped_mpz sticky_rem(m_mpz_manager);
if (exp_delta > sbits+5) {
m_mpz_manager.set(sticky_rem, 0);
sticky_rem.swap(subtrahend);
}
else
m_mpz_manager.machine_div_rem(subtrahend, m_powers2((unsigned)exp_delta), subtrahend, sticky_rem);
if (!m_mpz_manager.is_zero(sticky_rem) && m_mpz_manager.is_even(subtrahend))
m_mpz_manager.inc(subtrahend);
TRACE("mpf_dbg_rem", tout << "aligned subtrahend=" << m_mpz_manager.to_string(subtrahend) << std::endl;);
}
scoped_mpz minuend_lrg(m_mpz_manager);
m_mpz_manager.mul2k(minuend, 6, minuend_lrg); // + 6 extra bits into X
m_mpz_manager.sub(minuend_lrg, subtrahend, X_YQ_sig);
TRACE("mpf_dbg_rem", tout << "minuend=" << m_mpz_manager.to_string(minuend) << std::endl;
tout << "subtrahend=" << m_mpz_manager.to_string(subtrahend) << std::endl;
tout << "X_YQ_exp=" << X_YQ_exp << std::endl;
tout << "X_YQ_sig=" << m_mpz_manager.to_string(X_YQ_sig) << std::endl;
scoped_mpf & X_YQ = to_packed_mpf(YQ_sgn, X_YQ_exp, X_YQ_sig, ebits, sbits, 3);
tout << "X-YQ=" << to_string_hexfloat(X_YQ) << std::endl;);
SASSERT(m_mpz_manager.lt(X_YQ_sig, m_powers2(sbits+5)));
// 5. Rounding
if (m_mpz_manager.is_zero(X_YQ_sig)) {
TRACE("mpf_dbg_rem", tout << "sig zero" << std::endl;);
mk_zero(ebits, sbits, x.sign(), x);
}
else {
bool neg = m_mpz_manager.is_neg(X_YQ_sig);
if (neg) m_mpz_manager.neg(X_YQ_sig);
bool X_YQ_sgn = ((!x.sign() && YQ_sgn && neg) ||
(x.sign() && !YQ_sgn && !neg) ||
(x.sign() && YQ_sgn));
set(x, ebits, sbits, X_YQ_sgn, X_YQ_exp, X_YQ_sig);
round(MPF_ROUND_NEAREST_TEVEN, x);
}
TRACE("mpf_dbg_rem", tout << "partial remainder = " << to_string_hexfloat(x) << std::endl;);
}
void mpf_manager::rem(mpf const & x, mpf const & y, mpf & o) {
SASSERT(x.sbits == y.sbits && x.ebits == y.ebits);
@ -1231,54 +1379,32 @@ void mpf_manager::rem(mpf const & x, mpf const & y, mpf & o) {
else if (is_zero(x))
set(o, x);
else {
// This is a generalized version of the algorithm for FPREM1 in the Intel
SASSERT(is_regular(x) && is_regular(y));
SASSERT(x.exponent >= y.exponent);
// This is a generalized version of the algorithm for FPREM1 in the `Intel
// 64 and IA-32 Architectures Software Developer's Manual',
// Section 3-402 Vol. 2A FPREM1-Partial Remainder'.
// Section 3-402 Vol. 2A `FPREM1-Partial Remainder'.
scoped_mpf ST0(*this), ST1(*this);
set(ST0, x);
set(ST1, y);
unpack(ST0, true);
unpack(ST1, true);
const mpf_exp_t B = x.sbits-1; // max bits per iteration.
const mpf_exp_t B = x.sbits; // max bits per iteration.
mpf_exp_t D;
do {
D = ST0.exponent() - ST1.exponent();
TRACE("mpf_dbg_rem", tout << "st0=" << to_string_hexfloat(ST0) << std::endl;
tout << "st1=" << to_string_hexfloat(ST1) << std::endl;
SASSERT(0 <= D && D < INT64_MAX);
TRACE("mpf_dbg_rem", tout << "ST0=" << to_string_raw(ST0) << "=" << to_string_hexfloat(ST0) << std::endl;
tout << "ST1=" << to_string_raw(ST1) << "=" << to_string_hexfloat(ST1) << std::endl;
tout << "D=" << D << std::endl;);
if (D < B) {
scoped_mpf ST0_DIV_ST1(*this), Q(*this), ST1_MUL_Q(*this), ST0p(*this);
div(MPF_ROUND_NEAREST_TEVEN, ST0, ST1, ST0_DIV_ST1);
round_to_integral(MPF_ROUND_NEAREST_TEVEN, ST0_DIV_ST1, Q);
mul(MPF_ROUND_NEAREST_TEVEN, ST1, Q, ST1_MUL_Q);
sub(MPF_ROUND_NEAREST_TEVEN, ST0, ST1_MUL_Q, ST0p);
TRACE("mpf_dbg_rem", tout << "ST0/ST1=" << to_string_hexfloat(ST0_DIV_ST1) << std::endl;
tout << "Q=" << to_string_hexfloat(Q) << std::endl;
tout << "ST1*Q=" << to_string_hexfloat(ST1_MUL_Q) << std::endl;
tout << "ST0'=" << to_string_hexfloat(ST0p) << std::endl;);
set(ST0, ST0p);
}
else {
const mpf_exp_t N = B;
scoped_mpf ST0_DIV_ST1(*this), QQ(*this), ST1_MUL_QQ(*this), ST0p(*this);
div(MPF_ROUND_TOWARD_ZERO, ST0, ST1, ST0_DIV_ST1);
ST0_DIV_ST1.get().exponent -= D - N;
round_to_integral(MPF_ROUND_TOWARD_ZERO, ST0_DIV_ST1, QQ);
mul(MPF_ROUND_NEAREST_TEVEN, ST1, QQ, ST1_MUL_QQ);
ST1_MUL_QQ.get().exponent += D - N;
sub(MPF_ROUND_NEAREST_TEVEN, ST0, ST1_MUL_QQ, ST0p);
TRACE("mpf_dbg_rem", tout << "ST0/ST1/2^{D-N}=" << to_string_hexfloat(ST0_DIV_ST1) << std::endl;
tout << "QQ=" << to_string_hexfloat(QQ) << std::endl;
tout << "ST1*QQ*2^{D-N}=" << to_string_hexfloat(ST1_MUL_QQ) << std::endl;
tout << "ST0'=" << to_string_hexfloat(ST0p) << std::endl;);
SASSERT(!eq(ST0, ST0p));
set(ST0, ST0p);
}
SASSERT(ST0.exponent() - ST1.exponent() <= D);
prem(ST0, ST1, D, (D >= B));
} while (D >= B);
set(o, ST0);
set(o, x.ebits, x.sbits, MPF_ROUND_TOWARD_ZERO, ST0);
if (is_zero(o))
o.sign = x.sign;
}
@ -1364,9 +1490,9 @@ std::string mpf_manager::to_string(mpf const & x) {
}
}
//DEBUG_CODE(
// res += " " + to_string_hexfloat(x);
//);
DEBUG_CODE(
res += " " + to_string_raw(x);
);
return res;
}
@ -1407,6 +1533,16 @@ std::string mpf_manager::to_string_raw(mpf const & x) {
return res;
}
scoped_mpf mpf_manager::to_packed_mpf(bool sgn, mpf_exp_t exp, scoped_mpz const & sig, unsigned ebits, unsigned sbits, unsigned rbits) {
scoped_mpf q(*this);
scoped_mpz q_sig(m_mpz_manager);
m_mpz_manager.set(q_sig, sig);
if (rbits != 0) m_mpz_manager.div(q_sig, m_powers2(rbits), q_sig); // restore scale
m_mpz_manager.sub(q_sig, m_powers2(sbits-1), q_sig); // strip hidden bit
set(q, ebits, sbits, sgn, exp, q_sig);
return q;
}
std::string mpf_manager::to_string_hexfloat(mpf const & x) {
std::stringstream ss("");
std::ios::fmtflags ff = ss.setf(std::ios_base::hex | std::ios_base::uppercase |

10
src/util/mpf.h
View file

@ -185,9 +185,6 @@ public:
void mk_pinf(unsigned ebits, unsigned sbits, mpf & o);
void mk_ninf(unsigned ebits, unsigned sbits, mpf & o);
std::string to_string_raw(mpf const & a);
std::string to_string_hexfloat(mpf const & a);
unsynch_mpz_manager & mpz_manager(void) { return m_mpz_manager; }
unsynch_mpq_manager & mpq_manager(void) { return m_mpq_manager; }
@ -226,6 +223,8 @@ protected:
void round(mpf_rounding_mode rm, mpf & o);
void round_sqrt(mpf_rounding_mode rm, mpf & o);
void prem(scoped_mpf & x, scoped_mpf const & y, mpf_exp_t const & exp_diff, bool partial);
void mk_round_inf(mpf_rounding_mode rm, mpf & o);
// Convert x into a mpz numeral. zm is the manager that owns o.
@ -284,6 +283,10 @@ protected:
}
};
std::string to_string_raw(mpf const & a);
std::string to_string_hexfloat(mpf const & a);
scoped_mpf to_packed_mpf(bool sgn, mpf_exp_t exp, scoped_mpz const & sig, unsigned ebits, unsigned sbits, unsigned rbits);
public:
powers2 m_powers2;
};
@ -291,6 +294,7 @@ public:
class scoped_mpf : public _scoped_numeral<mpf_manager> {
friend class mpf_manager;
mpz & significand() { return get().significand; }
const mpz & significand() const { return get().significand; }
bool sign() const { return get().sign; }
mpf_exp_t exponent() const { return get().exponent; }
unsigned sbits() const { return get().sbits; }