From 25f75b60fefd5c8be1f905bc3aa2f8c844732a3d Mon Sep 17 00:00:00 2001 From: "Christoph M. Wintersteiger" Date: Thu, 10 Sep 2015 11:37:06 +0100 Subject: [PATCH] Bugfix for fp.mul and fp.fma for small numbers of e/s bits. Fixes #215. --- src/ast/rewriter/fpa_rewriter.cpp | 2 +- src/util/mpf.cpp | 19 ++++++++++++++----- src/util/mpf.h | 2 +- 3 files changed, 16 insertions(+), 7 deletions(-) diff --git a/src/ast/rewriter/fpa_rewriter.cpp b/src/ast/rewriter/fpa_rewriter.cpp index 07a39bcae..1d12621cd 100644 --- a/src/ast/rewriter/fpa_rewriter.cpp +++ b/src/ast/rewriter/fpa_rewriter.cpp @@ -468,7 +468,7 @@ br_status fpa_rewriter::mk_fma(expr * arg1, expr * arg2, expr * arg3, expr * arg scoped_mpf v2(m_fm), v3(m_fm), v4(m_fm); if (m_util.is_numeral(arg2, v2) && m_util.is_numeral(arg3, v3) && m_util.is_numeral(arg4, v4)) { scoped_mpf t(m_fm); - m_fm.fused_mul_add(rm, v2, v3, v4, t); + m_fm.fma(rm, v2, v3, v4, t); result = m_util.mk_value(t); return BR_DONE; } diff --git a/src/util/mpf.cpp b/src/util/mpf.cpp index 57e5b56b5..36974519c 100644 --- a/src/util/mpf.cpp +++ b/src/util/mpf.cpp @@ -638,7 +638,11 @@ void mpf_manager::mul(mpf_rounding_mode rm, mpf const & x, mpf const & y, mpf & // Remove the extra bits, keeping a sticky bit. scoped_mpz sticky_rem(m_mpz_manager); - m_mpz_manager.machine_div_rem(o.significand, m_powers2(x.sbits-4), o.significand, sticky_rem); + if (o.sbits >= 4) + m_mpz_manager.machine_div_rem(o.significand, m_powers2(o.sbits - 4), o.significand, sticky_rem); + else + m_mpz_manager.mul2k(o.significand, 4-o.sbits, o.significand); + if (!m_mpz_manager.is_zero(sticky_rem) && m_mpz_manager.is_even(o.significand)) m_mpz_manager.inc(o.significand); @@ -728,7 +732,7 @@ void mpf_manager::div(mpf_rounding_mode rm, mpf const & x, mpf const & y, mpf & } } -void mpf_manager::fused_mul_add(mpf_rounding_mode rm, mpf const & x, mpf const & y, mpf const &z, mpf & o) { +void mpf_manager::fma(mpf_rounding_mode rm, mpf const & x, mpf const & y, mpf const &z, mpf & o) { SASSERT(x.sbits == y.sbits && x.ebits == y.ebits && x.sbits == y.sbits && z.ebits == z.ebits); @@ -861,7 +865,7 @@ void mpf_manager::fused_mul_add(mpf_rounding_mode rm, mpf const & x, mpf const & o.exponent = mul_res.exponent(); - unsigned extra = x.sbits-4; + unsigned extra = 0; // Result could overflow into 4.xxx ... SASSERT(m_mpz_manager.lt(o.significand, m_powers2(2 * x.sbits + 2))); if(m_mpz_manager.ge(o.significand, m_powers2(2 * x.sbits + 1))) @@ -871,8 +875,13 @@ void mpf_manager::fused_mul_add(mpf_rounding_mode rm, mpf const & x, mpf const & TRACE("mpf_dbg", tout << "Addition overflew!" << std::endl;); } - // Get rid of the extra bits. - m_mpz_manager.machine_div_rem(o.significand, m_powers2(extra), o.significand, sticky_rem); + // Remove the extra bits, keeping a sticky bit. + m_mpz_manager.set(sticky_rem, 0); + if (o.sbits >= (4+extra)) + m_mpz_manager.machine_div_rem(o.significand, m_powers2(4+extra), o.significand, sticky_rem); + else + m_mpz_manager.mul2k(o.significand, (4+extra) - o.sbits, o.significand); + if (!m_mpz_manager.is_zero(sticky_rem) && m_mpz_manager.is_even(o.significand)) m_mpz_manager.inc(o.significand); diff --git a/src/util/mpf.h b/src/util/mpf.h index b99f8ab08..059f2ab25 100644 --- a/src/util/mpf.h +++ b/src/util/mpf.h @@ -121,7 +121,7 @@ public: void mul(mpf_rounding_mode rm, mpf const & x, mpf const & y, mpf & o); void div(mpf_rounding_mode rm, mpf const & x, mpf const & y, mpf & o); - void fused_mul_add(mpf_rounding_mode rm, mpf const & x, mpf const & y, mpf const &z, mpf & o); + void fma(mpf_rounding_mode rm, mpf const & x, mpf const & y, mpf const &z, mpf & o); void sqrt(mpf_rounding_mode rm, mpf const & x, mpf & o);