From f1c0ac72e79ce6e16c383fe05d877f0b4c9e7b5b Mon Sep 17 00:00:00 2001 From: "Christoph M. Wintersteiger" Date: Mon, 24 Jul 2017 20:54:29 +0100 Subject: [PATCH 1/5] Fix for fp.fma encoding. Relates to #872. --- src/ast/fpa/fpa2bv_converter.cpp | 117 ++++++++++++++++++------------- src/util/mpf.cpp | 2 +- 2 files changed, 70 insertions(+), 49 deletions(-) diff --git a/src/ast/fpa/fpa2bv_converter.cpp b/src/ast/fpa/fpa2bv_converter.cpp index 99b8c5ac8..5b71b2c2d 100644 --- a/src/ast/fpa/fpa2bv_converter.cpp +++ b/src/ast/fpa/fpa2bv_converter.cpp @@ -1516,7 +1516,8 @@ void fpa2bv_converter::mk_fma(func_decl * f, unsigned num, expr * const * args, mul_sgn = m_bv_util.mk_bv_xor(2, signs); dbg_decouple("fpa2bv_fma_mul_sgn", mul_sgn); - mul_exp = m_bv_util.mk_bv_add(m_bv_util.mk_bv_sub(a_exp_ext, a_lz_ext), + mul_exp = m_bv_util.mk_bv_add( + m_bv_util.mk_bv_sub(a_exp_ext, a_lz_ext), m_bv_util.mk_bv_sub(b_exp_ext, b_lz_ext)); dbg_decouple("fpa2bv_fma_mul_exp", mul_exp); @@ -1529,11 +1530,14 @@ void fpa2bv_converter::mk_fma(func_decl * f, unsigned num, expr * const * args, // The product has the form [-1][0].[2*sbits - 2]. // Extend c - c_sig = m_bv_util.mk_zero_extend(1, m_bv_util.mk_concat(c_sig, m_bv_util.mk_numeral(0, sbits-1))); + expr_ref c_sig_ext(m); + c_sig_ext = m_bv_util.mk_zero_extend(1, m_bv_util.mk_concat(c_sig, m_bv_util.mk_numeral(0, sbits - 1))); c_exp_ext = m_bv_util.mk_bv_sub(c_exp_ext, c_lz_ext); SASSERT(m_bv_util.get_bv_size(mul_sig) == 2 * sbits); - SASSERT(m_bv_util.get_bv_size(c_sig) == 2 * sbits); + SASSERT(m_bv_util.get_bv_size(c_sig_ext) == 2 * sbits); + dbg_decouple("fpa2bv_fma_c_sig_ext", c_sig_ext); + dbg_decouple("fpa2bv_fma_c_exp_ext", c_exp_ext); expr_ref swap_cond(m); swap_cond = m_bv_util.mk_sle(mul_exp, c_exp_ext); @@ -1542,10 +1546,10 @@ void fpa2bv_converter::mk_fma(func_decl * f, unsigned num, expr * const * args, expr_ref e_sgn(m), e_sig(m), e_exp(m), f_sgn(m), f_sig(m), f_exp(m); m_simp.mk_ite(swap_cond, c_sgn, mul_sgn, e_sgn); - m_simp.mk_ite(swap_cond, c_sig, mul_sig, e_sig); // has 2 * sbits + m_simp.mk_ite(swap_cond, c_sig_ext, mul_sig, e_sig); // has 2 * sbits m_simp.mk_ite(swap_cond, c_exp_ext, mul_exp, e_exp); // has ebits + 2 m_simp.mk_ite(swap_cond, mul_sgn, c_sgn, f_sgn); - m_simp.mk_ite(swap_cond, mul_sig, c_sig, f_sig); // has 2 * sbits + m_simp.mk_ite(swap_cond, mul_sig, c_sig_ext, f_sig); // has 2 * sbits m_simp.mk_ite(swap_cond, mul_exp, c_exp_ext, f_exp); // has ebits + 2 SASSERT(is_well_sorted(m, e_sgn)); @@ -1560,6 +1564,11 @@ void fpa2bv_converter::mk_fma(func_decl * f, unsigned num, expr * const * args, SASSERT(m_bv_util.get_bv_size(e_exp) == ebits + 2); SASSERT(m_bv_util.get_bv_size(f_exp) == ebits + 2); + dbg_decouple("fpa2bv_fma_e_sig", e_sig); + dbg_decouple("fpa2bv_fma_e_exp", e_exp); + dbg_decouple("fpa2bv_fma_f_sig", f_sig); + dbg_decouple("fpa2bv_fma_f_exp", f_exp); + expr_ref res_sgn(m), res_sig(m), res_exp(m); expr_ref exp_delta(m); @@ -1582,6 +1591,7 @@ void fpa2bv_converter::mk_fma(func_decl * f, unsigned num, expr * const * args, m_bv_util.mk_zero_extend((3*sbits)-(ebits+2), exp_delta)); shifted_f_sig = m_bv_util.mk_extract(3*sbits-1, sbits, shifted_big); sticky_raw = m_bv_util.mk_extract(sbits-1, 0, shifted_big); + dbg_decouple("fpa2bv_fma_shifted_f_sig", shifted_f_sig); SASSERT(m_bv_util.get_bv_size(shifted_f_sig) == 2 * sbits); SASSERT(is_well_sorted(m, shifted_f_sig)); @@ -1594,14 +1604,16 @@ void fpa2bv_converter::mk_fma(func_decl * f, unsigned num, expr * const * args, expr * or_args[2] = { shifted_f_sig, sticky }; shifted_f_sig = m_bv_util.mk_bv_or(2, or_args); SASSERT(is_well_sorted(m, shifted_f_sig)); + dbg_decouple("fpa2bv_fma_f_sig_or_sticky", shifted_f_sig); // Significant addition. - // Two extra bits for catching the overflow. + // Two extra bits for the sign and for catching overflows. e_sig = m_bv_util.mk_zero_extend(2, e_sig); shifted_f_sig = m_bv_util.mk_zero_extend(2, shifted_f_sig); expr_ref eq_sgn(m); m_simp.mk_eq(e_sgn, f_sgn, eq_sgn); + dbg_decouple("fpa2bv_fma_eq_sgn", eq_sgn); SASSERT(m_bv_util.get_bv_size(e_sig) == 2*sbits + 2); SASSERT(m_bv_util.get_bv_size(shifted_f_sig) == 2*sbits + 2); @@ -1611,27 +1623,24 @@ void fpa2bv_converter::mk_fma(func_decl * f, unsigned num, expr * const * args, expr_ref sum(m), e_plus_f(m), e_minus_f(m); e_plus_f = m_bv_util.mk_bv_add(e_sig, shifted_f_sig); - e_minus_f = m_bv_util.mk_bv_sub(e_sig, shifted_f_sig), - m_simp.mk_ite(eq_sgn, e_plus_f, e_minus_f, sum); - + e_minus_f = m_bv_util.mk_bv_sub(e_sig, shifted_f_sig); + m_simp.mk_ite(eq_sgn, e_plus_f, e_minus_f, sum); + SASSERT(m_bv_util.get_bv_size(sum) == 2*sbits + 2); SASSERT(is_well_sorted(m, sum)); - dbg_decouple("fpa2bv_fma_add_sum", sum); expr_ref sign_bv(m), n_sum(m); sign_bv = m_bv_util.mk_extract(2*sbits+1, 2*sbits+1, sum); n_sum = m_bv_util.mk_bv_neg(sum); + dbg_decouple("fpa2bv_fma_add_sign_bv", sign_bv); + dbg_decouple("fpa2bv_fma_add_n_sum", n_sum); expr_ref res_sig_eq(m), sig_abs(m), one_1(m); one_1 = m_bv_util.mk_numeral(1, 1); m_simp.mk_eq(sign_bv, one_1, res_sig_eq); m_simp.mk_ite(res_sig_eq, n_sum, sum, sig_abs); - dbg_decouple("fpa2bv_fma_add_sign_bv", sign_bv); - dbg_decouple("fpa2bv_fma_add_n_sum", n_sum); dbg_decouple("fpa2bv_fma_add_sig_abs", sig_abs); - res_exp = e_exp; - family_id bvfid = m_bv_util.get_fid(); expr_ref res_sgn_c1(m), res_sgn_c2(m), res_sgn_c3(m); expr_ref not_e_sgn(m), not_f_sgn(m), not_sign_bv(m); @@ -1643,44 +1652,57 @@ void fpa2bv_converter::mk_fma(func_decl * f, unsigned num, expr * const * args, res_sgn_c3 = m.mk_app(bvfid, OP_BAND, e_sgn, f_sgn); expr * res_sgn_or_args[3] = { res_sgn_c1, res_sgn_c2, res_sgn_c3 }; res_sgn = m_bv_util.mk_bv_or(3, res_sgn_or_args); + dbg_decouple("fpa2bv_fma_res_sgn", res_sgn); + + expr_ref is_sig_neg(m); + is_sig_neg = m.mk_eq(one_1, m_bv_util.mk_extract(2 * sbits + 1, 2 * sbits + 1, sig_abs)); + sig_abs = m.mk_ite(is_sig_neg, m_bv_util.mk_bv_neg(sig_abs), sig_abs); + dbg_decouple("fpa2bv_fma_is_sig_neg", is_sig_neg); // Result could have overflown into 4.xxx. SASSERT(m_bv_util.get_bv_size(sig_abs) == 2 * sbits + 2); - expr_ref ovfl_into_4(m); - ovfl_into_4 = m.mk_eq(m_bv_util.mk_extract(2 * sbits + 1, 2 * sbits, sig_abs), - m_bv_util.mk_numeral(1, 2)); - dbg_decouple("fpa2bv_fma_ovfl_into_4", ovfl_into_4); - if (sbits > 5) { - expr_ref sticky_raw(m), sig_upper(m), sticky_redd(m), res_sig_norm(m); - sticky_raw = m_bv_util.mk_extract(sbits - 4, 0, sig_abs); - sig_upper = m_bv_util.mk_extract(2 * sbits, sbits - 3, sig_abs); - SASSERT(m_bv_util.get_bv_size(sig_upper) == sbits + 4); - sticky_redd = m.mk_app(bvfid, OP_BREDOR, sticky_raw.get()); - sticky = m_bv_util.mk_zero_extend(sbits + 3, sticky_redd); - expr * res_or_args[2] = { sig_upper, sticky }; - res_sig_norm = m_bv_util.mk_bv_or(2, res_or_args); + expr_ref extra(m), extra_is_zero(m); + extra = m_bv_util.mk_extract(2 * sbits + 1, 2 * sbits, sig_abs); + extra_is_zero = m.mk_eq(extra, m_bv_util.mk_numeral(0, 2)); + dbg_decouple("fpa2bv_fma_extra", extra); - expr_ref sticky_raw_ovfl(m), sig_upper_ovfl(m), sticky_redd_ovfl(m), sticky_ovfl(m), res_sig_ovfl(m); - sticky_raw_ovfl = m_bv_util.mk_extract(sbits - 4, 0, sig_abs); - sig_upper_ovfl = m_bv_util.mk_extract(2 * sbits, sbits - 3, sig_abs); - SASSERT(m_bv_util.get_bv_size(sig_upper_ovfl) == sbits + 4); - sticky_redd_ovfl = m.mk_app(bvfid, OP_BREDOR, sticky_raw_ovfl.get()); - sticky_ovfl = m_bv_util.mk_zero_extend(sbits + 3, sticky_redd_ovfl); - expr * res_or_args_ovfl[2] = { sig_upper_ovfl, sticky_ovfl }; - res_sig_ovfl = m_bv_util.mk_bv_or(2, res_or_args_ovfl); - - res_sig = m.mk_ite(ovfl_into_4, res_sig_ovfl, res_sig_norm); - res_exp = m.mk_ite(ovfl_into_4, m_bv_util.mk_bv_add(res_exp, m_bv_util.mk_numeral(1, ebits+2)), - res_exp); - } - else { - unsigned too_short = 6 - sbits; + unsigned too_short = 0; + if (sbits < 5) { + too_short = 6 - sbits + 1; sig_abs = m_bv_util.mk_concat(sig_abs, m_bv_util.mk_numeral(0, too_short)); - res_sig = m_bv_util.mk_extract(sbits + 3, 0, sig_abs); } - dbg_decouple("fpa2bv_fma_add_sum_sticky", sticky); - dbg_decouple("fpa2bv_fma_sig_abs", sig_abs); + + expr_ref sig_abs_h1(m), sticky_h1(m), sticky_h1_red(m), sig_abs_h1_f(m), res_sig_1(m); + sticky_h1 = m_bv_util.mk_extract(sbits - 5 + too_short, 0, sig_abs); + sig_abs_h1 = m_bv_util.mk_extract(2 * sbits + 1 + too_short, sbits - 4 + too_short, sig_abs); + sticky_h1_red = m_bv_util.mk_zero_extend(sbits + 5, m.mk_app(m_bv_util.get_fid(), OP_BREDOR, sticky_h1.get())); + expr * sticky_h1_red_args[2] = { sig_abs_h1, sticky_h1_red }; + sig_abs_h1_f = m_bv_util.mk_bv_or(2, sticky_h1_red_args); + res_sig_1 = m_bv_util.mk_extract(sbits + 3, 0, sig_abs_h1_f); + SASSERT(m_bv_util.get_bv_size(sticky_h1) == sbits - 4 + too_short); + SASSERT(m_bv_util.get_bv_size(sig_abs_h1) == sbits + 6); + SASSERT(m_bv_util.get_bv_size(sticky_h1_red) == sbits + 6); + SASSERT(m_bv_util.get_bv_size(sig_abs_h1_f) == sbits + 6); + SASSERT(m_bv_util.get_bv_size(res_sig_1) == sbits + 4); + + expr_ref sig_abs_h2(m), sticky_h2(m), sticky_h2_red(m), sig_abs_h2_f(m), res_sig_2(m); + sticky_h2 = m_bv_util.mk_extract(sbits - 4 + too_short, 0, sig_abs); + sig_abs_h2 = m_bv_util.mk_extract(2 * sbits + 1 + too_short, sbits - 3 + too_short, sig_abs); + sticky_h2_red = m_bv_util.mk_zero_extend(sbits + 4, m.mk_app(m_bv_util.get_fid(), OP_BREDOR, sticky_h1.get())); + expr * sticky_h2_red_args[2] = { sig_abs_h2, sticky_h2_red }; + sig_abs_h2_f = m_bv_util.mk_zero_extend(1, m_bv_util.mk_bv_or(2, sticky_h2_red_args)); + res_sig_2 = m_bv_util.mk_extract(sbits + 3, 0, sig_abs_h2_f); + SASSERT(m_bv_util.get_bv_size(sticky_h2) == sbits - 3 + too_short); + SASSERT(m_bv_util.get_bv_size(sig_abs_h2) == sbits + 5); + SASSERT(m_bv_util.get_bv_size(sticky_h2_red) == sbits + 5); + SASSERT(m_bv_util.get_bv_size(sig_abs_h2_f) == sbits + 6); + SASSERT(m_bv_util.get_bv_size(res_sig_2) == sbits + 4); + + res_sig = m.mk_ite(extra_is_zero, res_sig_1, res_sig_2); + res_exp = m.mk_ite(extra_is_zero, e_exp, m_bv_util.mk_bv_add(e_exp, m_bv_util.mk_numeral(1, ebits + 2))); + dbg_decouple("fpa2bv_fma_res_sig", res_sig); + dbg_decouple("fpa2bv_fma_res_exp", res_exp); SASSERT(m_bv_util.get_bv_size(res_sig) == sbits + 4); expr_ref is_zero_sig(m), nil_sbits4(m); @@ -3889,8 +3911,8 @@ void fpa2bv_converter::round(sort * s, expr_ref & rm, expr_ref & sgn, expr_ref & mk_min_exp(ebits, e_min); mk_max_exp(ebits, e_max); - TRACE("fpa2bv_dbg", tout << "e_min = " << mk_ismt2_pp(e_min, m) << std::endl << - "e_max = " << mk_ismt2_pp(e_max, m) << std::endl;); + TRACE("fpa2bv_dbg", tout << "e_min = " << mk_ismt2_pp(e_min, m) << std::endl; + tout << "e_max = " << mk_ismt2_pp(e_max, m) << std::endl;); expr_ref OVF1(m), e_top_three(m), sigm1(m), e_eq_emax_and_sigm1(m), e_eq_emax(m); expr_ref e3(m), ne3(m), e2(m), e1(m), e21(m), one_1(m), h_exp(m), sh_exp(m), th_exp(m); @@ -4106,7 +4128,6 @@ void fpa2bv_converter::round(sort * s, expr_ref & rm, expr_ref & sgn, expr_ref & dbg_decouple("fpa2bv_rnd_rm_is_to_neg", rm_is_to_neg); dbg_decouple("fpa2bv_rnd_rm_is_to_pos", rm_is_to_pos); - expr_ref sgn_is_zero(m), zero1(m); zero1 = m_bv_util.mk_numeral(0, 1); m_simp.mk_eq(sgn, zero1, sgn_is_zero); diff --git a/src/util/mpf.cpp b/src/util/mpf.cpp index e9d108cec..7b5c8d62d 100644 --- a/src/util/mpf.cpp +++ b/src/util/mpf.cpp @@ -886,7 +886,7 @@ void mpf_manager::fma(mpf_rounding_mode rm, mpf const & x, mpf const & y, mpf co 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))) + if (m_mpz_manager.ge(o.significand, m_powers2(2 * x.sbits + 1))) { extra++; o.exponent++; From 75b533f050c4029ad2d1b85532614a8749f1d9aa Mon Sep 17 00:00:00 2001 From: "Christoph M. Wintersteiger" Date: Tue, 25 Jul 2017 20:26:17 +0100 Subject: [PATCH 2/5] Fixed normalization shift in MPF rounder. Relates to #872. --- src/util/mpf.cpp | 46 ++++++++++++++++++---------------------------- 1 file changed, 18 insertions(+), 28 deletions(-) diff --git a/src/util/mpf.cpp b/src/util/mpf.cpp index 7b5c8d62d..fac314837 100644 --- a/src/util/mpf.cpp +++ b/src/util/mpf.cpp @@ -1947,35 +1947,25 @@ void mpf_manager::round(mpf_rounding_mode rm, mpf & o) { TRACE("mpf_dbg", tout << "Shift distance: " << m_mpz_manager.to_string(sigma) << " " << ((m_mpz_manager.is_nonneg(sigma))?"(LEFT)":"(RIGHT)") << std::endl;); - bool sticky = !m_mpz_manager.is_even(o.significand); - m_mpz_manager.machine_div2k(o.significand, 1); // Let f' = f_r/2 - - if (!m_mpz_manager.is_zero(sigma)) { - if (m_mpz_manager.is_neg(sigma)) { // Right shift - unsigned sigma_uint = (unsigned) -m_mpz_manager.get_int64(sigma); // sigma is capped, this is safe. - if (sticky) - m_mpz_manager.machine_div2k(o.significand, sigma_uint); - else - { - scoped_mpz sticky_rem(m_mpz_manager); - m_mpz_manager.machine_div_rem(o.significand, m_powers2(sigma_uint), o.significand, sticky_rem); - sticky = !m_mpz_manager.is_zero(sticky_rem); - } - } - else { // Left shift - unsigned sh_m = static_cast(m_mpz_manager.get_int64(sigma)); - m_mpz_manager.mul2k(o.significand, sh_m, o.significand); - m_mpz_manager.set(sigma, 0); - } + if (m_mpz_manager.le(sigma, 0)) { // Right shift + scoped_mpz sticky_rem(m_mpz_manager); + unsigned sigma_uint = (unsigned)-m_mpz_manager.get_int64(sigma); // sigma is capped, this is safe. + sigma_uint += 2; + m_mpz_manager.machine_div_rem(o.significand, m_powers2(sigma_uint), o.significand, sticky_rem); + bool sticky = !m_mpz_manager.is_zero(sticky_rem); + if (sticky && m_mpz_manager.is_even(o.significand)) + m_mpz_manager.inc(o.significand); + } + else { // Left shift + unsigned sigma_uint = static_cast(m_mpz_manager.get_int64(sigma)); + m_mpz_manager.mul2k(o.significand, sigma_uint - 1, o.significand); + bool sticky = !m_mpz_manager.is_even(o.significand); + m_mpz_manager.machine_div2k(o.significand, 1); + if (sticky && m_mpz_manager.is_even(o.significand)) + m_mpz_manager.inc(o.significand); } - TRACE("mpf_dbg", tout << "Before sticky: " << to_string(o) << std::endl;); - - // Make sure o.significand is a [sbits+2] bit number (i.e. f1[0:sbits+1] == f1[0:sbits-1][round][sticky]) - sticky = sticky || !m_mpz_manager.is_even(o.significand); - m_mpz_manager.machine_div2k(o.significand, 1); - if (sticky && m_mpz_manager.is_even(o.significand)) - m_mpz_manager.inc(o.significand); + TRACE("mpf_dbg", tout << "After sticky: " << to_string(o) << std::endl;); if (OVF1 && OVFen) { o.exponent = beta; @@ -1997,7 +1987,7 @@ void mpf_manager::round(mpf_rounding_mode rm, mpf & o) { // Significand rounding (sigrd) - sticky = !m_mpz_manager.is_even(o.significand); // new sticky bit! + bool sticky = !m_mpz_manager.is_even(o.significand); // new sticky bit! m_mpz_manager.machine_div2k(o.significand, 1); bool round = !m_mpz_manager.is_even(o.significand); m_mpz_manager.machine_div2k(o.significand, 1); From 33ebdc8adc0ac55ef2e071a15bdb3236ff0d01d7 Mon Sep 17 00:00:00 2001 From: "Christoph M. Wintersteiger" Date: Thu, 27 Jul 2017 23:08:35 +0100 Subject: [PATCH 3/5] Cleaned up mpf rounder. Rewrote mpf fma. Relates to #872. --- src/ast/fpa/fpa2bv_converter.cpp | 4 - src/util/mpf.cpp | 348 ++++++++++++++++--------------- src/util/mpf.h | 3 +- 3 files changed, 184 insertions(+), 171 deletions(-) diff --git a/src/ast/fpa/fpa2bv_converter.cpp b/src/ast/fpa/fpa2bv_converter.cpp index 5b71b2c2d..0e83d7fea 100644 --- a/src/ast/fpa/fpa2bv_converter.cpp +++ b/src/ast/fpa/fpa2bv_converter.cpp @@ -3904,9 +3904,6 @@ void fpa2bv_converter::round(sort * s, expr_ref & rm, expr_ref & sgn, expr_ref & SASSERT(m_bv_util.get_bv_size(sig) == sbits+4); SASSERT(m_bv_util.get_bv_size(exp) == ebits+2); - // bool UNFen = false; - // bool OVFen = false; - expr_ref e_min(m), e_max(m); mk_min_exp(ebits, e_min); mk_max_exp(ebits, e_max); @@ -4025,7 +4022,6 @@ void fpa2bv_converter::round(sort * s, expr_ref & rm, expr_ref & sgn, expr_ref & SASSERT(is_well_sorted(m, sig)); SASSERT(m_bv_util.get_bv_size(sig) == sbits+2); - // CMW: The (OVF1 && OVFen) and (TINY && UNFen) cases are never taken. expr_ref ext_emin(m); ext_emin = m_bv_util.mk_zero_extend(2, e_min); m_simp.mk_ite(TINY, ext_emin, beta, exp); diff --git a/src/util/mpf.cpp b/src/util/mpf.cpp index fac314837..36d9e49f1 100644 --- a/src/util/mpf.cpp +++ b/src/util/mpf.cpp @@ -805,7 +805,7 @@ void mpf_manager::fma(mpf_rounding_mode rm, mpf const & x, mpf const & y, mpf co o.ebits = x.ebits; o.sbits = x.sbits; - scoped_mpf mul_res(*this, x.ebits+2, 2*x.sbits); + scoped_mpf mul_res(*this, x.ebits+2, 2*x.sbits-1); scoped_mpf a(*this, x.ebits, x.sbits), b(*this, x.ebits, x.sbits), c(*this, x.ebits, x.sbits); set(a, x); set(b, y); @@ -814,9 +814,12 @@ void mpf_manager::fma(mpf_rounding_mode rm, mpf const & x, mpf const & y, mpf co unpack(b, true); unpack(c, true); - TRACE("mpf_dbg", tout << "A = " << to_string(a) << std::endl;); - TRACE("mpf_dbg", tout << "B = " << to_string(b) << std::endl;); - TRACE("mpf_dbg", tout << "C = " << to_string(c) << std::endl;); + TRACE("mpf_dbg", tout << "A = " << to_string(a) << std::endl; + tout << "B = " << to_string(b) << std::endl; + tout << "C = " << to_string(c) << std::endl; + tout << "A = " << to_string_binary(a, 1, 0) << std::endl; + tout << "B = " << to_string_binary(b, 1, 0) << std::endl; + tout << "C = " << to_string_binary(c, 1, 0) << std::endl;); SASSERT(m_mpz_manager.lt(a.significand(), m_powers2(x.sbits))); SASSERT(m_mpz_manager.lt(b.significand(), m_powers2(x.sbits))); @@ -824,98 +827,114 @@ void mpf_manager::fma(mpf_rounding_mode rm, mpf const & x, mpf const & y, mpf co SASSERT(m_mpz_manager.ge(a.significand(), m_powers2(x.sbits-1))); SASSERT(m_mpz_manager.ge(b.significand(), m_powers2(x.sbits-1))); + // A/B in _1.[sbits-1]. mul_res.get().sign = (a.sign() != b.sign()); mul_res.get().exponent = a.exponent() + b.exponent(); - m_mpz_manager.mul(a.significand(), b.significand(), mul_res.get().significand); + m_mpz_manager.mul(a.significand(), b.significand(), mul_res.significand()); - TRACE("mpf_dbg", tout << "PRODUCT = " << to_string(mul_res) << std::endl;); + TRACE("mpf_dbg", tout << "M = " << to_string(mul_res) << std::endl; + tout << "M = " << to_string_binary(mul_res, 1, 0) << std::endl;); - // mul_res is [-1][0].[2*sbits - 2], i.e., between 2*sbits-1 and 2*sbits. + // mul_res is [-1][0].[2*sbits - 2], i.e., >= 2^(2*sbits-2) and < 2^(2*sbits). SASSERT(m_mpz_manager.lt(mul_res.significand(), m_powers2(2*x.sbits))); SASSERT(m_mpz_manager.ge(mul_res.significand(), m_powers2(2*x.sbits - 2))); - // Introduce extra bits into c. - m_mpz_manager.mul2k(c.significand(), x.sbits-1, c.significand()); + // Introduce extra bits into c in _[0].[sbits-1] s.t. c in _[0].[2*sbits-2] + c.get().ebits = x.ebits + 2; + c.get().sbits = 2 * x.sbits - 1; + m_mpz_manager.mul2k(c.significand(), x.sbits - 1); - SASSERT(m_mpz_manager.lt(c.significand(), m_powers2(2 * x.sbits - 1))); + TRACE("mpf_dbg", tout << "C_= " << to_string(c) << std::endl; + tout << "C_= " << to_string_binary(c, 1, 0) << std::endl;); + + SASSERT(m_mpz_manager.lt(c.significand(), m_powers2(2*x.sbits))); SASSERT(m_mpz_manager.is_zero(c.significand()) || - m_mpz_manager.ge(c.significand(), m_powers2(2 * x.sbits - 2))); + m_mpz_manager.ge(c.significand(), m_powers2(2*x.sbits - 2))); - TRACE("mpf_dbg", tout << "C = " << to_string(c) << std::endl;); - - if (exp(c) > exp(mul_res)) + if (exp(c) > exp(mul_res)) { + TRACE("mpf_dbg", tout << "Swap!" << std::endl;); mul_res.swap(c); + } - mpf_exp_t exp_delta = exp(mul_res) - exp(c); + mpf_exp_t exp_delta_w = exp(mul_res) - exp(c); - SASSERT(exp(mul_res) >= exp(c) && exp_delta >= 0); + SASSERT(exp(mul_res) >= exp(c) && exp_delta_w >= 0); - if (exp_delta > 2 * x.sbits) - exp_delta = 2 * x.sbits; + if (exp_delta_w > 2 * x.sbits) + exp_delta_w = 2 * x.sbits; + unsigned exp_delta = (unsigned)exp_delta_w; TRACE("mpf_dbg", tout << "exp_delta = " << exp_delta << std::endl;); // Alignment shift with sticky bit computation. scoped_mpz sticky_rem(m_mpz_manager); - m_mpz_manager.machine_div_rem(c.significand(), m_powers2((int)exp_delta), c.significand(), sticky_rem); + m_mpz_manager.machine_div_rem(c.significand(), m_powers2(exp_delta), c.significand(), sticky_rem); TRACE("mpf_dbg", tout << "alignment shift -> sig = " << m_mpz_manager.to_string(c.significand()) << " sticky_rem = " << m_mpz_manager.to_string(sticky_rem) << std::endl;); if (!m_mpz_manager.is_zero(sticky_rem) && m_mpz_manager.is_even(c.significand())) m_mpz_manager.inc(c.significand()); - TRACE("mpf_dbg", tout << "M' = " << m_mpz_manager.to_string(mul_res.significand()) << std::endl;); - TRACE("mpf_dbg", tout << "C' = " << m_mpz_manager.to_string(c.significand()) << std::endl;); + TRACE("mpf_dbg", tout << "M'= " << m_mpz_manager.to_string(mul_res.significand()) << std::endl; + tout << "M'= " << to_string_binary(mul_res, 1, 0) << std::endl; ); + TRACE("mpf_dbg", tout << "C'= " << m_mpz_manager.to_string(c.significand()) << std::endl; + tout << "C'= " << to_string_binary(c, 1, 0) << std::endl; ); + + scoped_mpf res(c); + + res.get().sign = mul_res.sign(); + res.get().exponent = mul_res.exponent(); // Significand addition + if (sgn(mul_res) != sgn(c)) { TRACE("mpf_dbg", tout << "SUBTRACTING" << std::endl;); - m_mpz_manager.sub(mul_res.significand(), c.significand(), o.significand); + m_mpz_manager.sub(mul_res.significand(), c.significand(), res.significand()); + + if (m_mpz_manager.is_neg(res.significand())) { + m_mpz_manager.abs(res.significand()); + res.get().sign |= !mul_res.sign() && c.sign(); + } } else { TRACE("mpf_dbg", tout << "ADDING" << std::endl;); - m_mpz_manager.add(mul_res.significand(), c.significand(), o.significand); + m_mpz_manager.add(mul_res.significand(), c.significand(), res.significand()); } - TRACE("mpf_dbg", tout << "sum[-1:] = " << m_mpz_manager.to_string(o.significand) << std::endl;); + TRACE("mpf_dbg", tout << "S'= " << m_mpz_manager.to_string(res.significand()) << std::endl; + tout << "S'= " << to_string_binary(res, 1, 0) << std::endl; ); - bool neg = m_mpz_manager.is_neg(o.significand); - TRACE("mpf_dbg", tout << "NEG=" << neg << std::endl;); - if (neg) m_mpz_manager.abs(o.significand); - - o.exponent = mul_res.exponent(); - - 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))) + SASSERT(m_mpz_manager.lt(res.significand(), m_powers2(2 * x.sbits + 1))); + if (m_mpz_manager.ge(res.significand(), m_powers2(2 * x.sbits))) { - extra++; - o.exponent++; - TRACE("mpf_dbg", tout << "Addition overflew!" << std::endl;); + SASSERT(exp(res) < mk_max_exp(x.ebits)); // NYI. + + res.get().exponent++; + bool sticky = !m_mpz_manager.is_even(res.significand()); + m_mpz_manager.machine_div2k(res.significand(), 1); + if (sticky && m_mpz_manager.is_even(res.significand())) + m_mpz_manager.inc(res.significand()); + TRACE("mpf_dbg", tout << "Add/Sub overflew into 4.xxx!" << std::endl; + tout << "S*= " << to_string_binary(res, 2, 0) << std::endl;); } - // Remove the extra bits, keeping a sticky bit. - m_mpz_manager.set(sticky_rem, 0); - unsigned minbits = (4 + extra); - if (o.sbits >= minbits) - m_mpz_manager.machine_div_rem(o.significand, m_powers2(o.sbits - minbits), o.significand, sticky_rem); - else - m_mpz_manager.mul2k(o.significand, minbits - o.sbits, o.significand); + set(o, x.ebits, x.sbits, res.sign(), res.exponent(), mpz(0)); - if (!m_mpz_manager.is_zero(sticky_rem) && m_mpz_manager.is_even(o.significand)) - m_mpz_manager.inc(o.significand); + if (x.sbits >= 4) { + m_mpz_manager.set(sticky_rem, 0); + m_mpz_manager.machine_div_rem(res.significand(), m_powers2(x.sbits - 4), o.significand, sticky_rem); + if (!m_mpz_manager.is_zero(sticky_rem) && m_mpz_manager.is_even(o.significand)) + m_mpz_manager.inc(o.significand); + } + else { + m_mpz_manager.mul2k(res.significand(), 4 - x.sbits, o.significand); + } - TRACE("mpf_dbg", tout << "sum[-1:sbits+2] = " << m_mpz_manager.to_string(o.significand) << std::endl;); + TRACE("mpf_dbg", tout << "sum[-1:sbits+2] = " << m_mpz_manager.to_string(o.significand) << std::endl; + tout << "R = " << to_string_binary(o, 2, 3) << std::endl;); if (m_mpz_manager.is_zero(o.significand)) mk_zero(x.ebits, x.sbits, rm == MPF_ROUND_TOWARD_NEGATIVE, o); - else { - o.sign = ((!mul_res.sign() && c.sign() && neg) || - ( mul_res.sign() && !c.sign() && !neg) || - ( mul_res.sign() && c.sign())); - TRACE("mpf_dbg", tout << "before round = " << to_string(o) << std::endl << - "fs[-1:sbits+2] = " << m_mpz_manager.to_string(o.significand) << std::endl;); + else round(rm, o); - } } } @@ -1596,6 +1615,64 @@ std::string mpf_manager::to_string_hexfloat(mpf const & x) { return ss.str(); } +std::string mpf_manager::to_string_binary(mpf const & x, unsigned upper_extra, unsigned lower_extra) { + std::string res; + + if (is_nan(x)) + res = std::string("") + "#b0 " + + "#b" + std::string(x.ebits, '1') + " " + + "#b" + std::string(x.sbits-2, '0') + "1 " + + "(NaN)"; + else if (is_inf(x)) + res = std::string("") + "#b" + (sgn(x)?"1":"0") + " " + + "#b" + std::string(x.ebits, '1') + " " + + "#b" + std::string(x.sbits - 1, '0') + "1 " + + "(" + (sgn(x)?"-":"+") + "oo)"; + else if (is_zero(x)) + res = std::string("") + "#b" + (sgn(x) ? "1" : "0") + " " + + "#b" + std::string(x.ebits, '0') + " " + + "#b" + std::string(x.sbits - 1, '0') + " " + + "(" + (sgn(x) ? "-" : "+") + "zero)"; + else { + res = std::string("") + "#b" + (sgn(x) ? "1" : "0") + " "; + + scoped_mpz tmp(m_mpz_manager); + + if (is_denormal(x)) + m_mpz_manager.set(tmp, bias_exp(x.ebits, mk_min_exp(x.ebits))); + else { + m_mpz_manager.set(tmp, bias_exp(x.ebits, exp(x))); + } + + std::string tmp_str = ""; + for (unsigned i = 0; i < x.ebits; i++) { + tmp_str += m_mpz_manager.is_odd(tmp) ? "1" : "0"; + tmp /= 2; + } + std::reverse(tmp_str.begin(), tmp_str.end()); + res += "#b" + tmp_str + " "; + + tmp_str = ""; + m_mpz_manager.set(tmp, sig(x)); + unsigned num_bits = upper_extra + x.sbits + lower_extra; + for (unsigned i = 0; i < num_bits || !tmp.is_zero(); i++) { + tmp_str += m_mpz_manager.is_odd(tmp) ? "1" : "0"; + tmp /= 2; + if (i == lower_extra - 1) + tmp_str += ","; + if (i == x.sbits + lower_extra - 2) { + tmp_str += "."; + if (i == num_bits - 1) + tmp_str += " "; + } + } + std::reverse(tmp_str.begin(), tmp_str.end()); + res += "#b" + tmp_str; + } + + return res; +} + void mpf_manager::to_rational(mpf const & x, unsynch_mpq_manager & qm, mpq & o) { scoped_mpf a(*this); scoped_mpz n(m_mpq_manager), d(m_mpq_manager); @@ -1865,129 +1942,76 @@ void mpf_manager::mk_round_inf(mpf_rounding_mode rm, mpf & o) { } void mpf_manager::round(mpf_rounding_mode rm, mpf & o) { - // Assumptions: o.significand is of the form f[-1:0] . f[1:sbits-1] [guard,round,sticky], + // Assumptions: o.significand is of the form f[-1:0] . f[1:sbits-1] [round,extra,sticky], // i.e., it has 2 + (sbits-1) + 3 = sbits + 4 bits. - TRACE("mpf_dbg", tout << "RND: " << to_string(o) << std::endl;); + TRACE("mpf_dbg", tout << "RND: " << to_string(o) << std::endl; + tout << to_string_binary(o, 1, 3) << std::endl;); DEBUG_CODE({ - const mpz & p_m3 = m_powers2(o.sbits+5); + const mpz & p_m3 = m_powers2(o.sbits+4); SASSERT(m_mpz_manager.lt(o.significand, p_m3)); }); - // Structure of the rounder: - // (s, e_out, f_out) == (s, exprd(s, post(e, sigrd(s, f)))). + mpf_exp_t e_max = mk_max_exp(o.ebits); + mpf_exp_t e_min = mk_min_exp(o.ebits); - bool UNFen = false; // Are these supposed to be persistent flags accross calls? - bool OVFen = false; + unsigned sig_width = m_mpz_manager.prev_power_of_two(o.significand) + 1; + mpf_exp_t lz = o.sbits + 4 - sig_width; + mpf_exp_t beta = o.exponent - lz + 1; - mpf_exp_t e_max_norm = mk_max_exp(o.ebits); - mpf_exp_t e_min_norm = mk_min_exp(o.ebits); - scoped_mpz temporary(m_mpq_manager); + scoped_mpz sigma(m_mpz_manager); - TRACE("mpf_dbg", tout << "e_min_norm = " << e_min_norm << std::endl << - "e_max_norm = " << e_max_norm << std::endl;); - - const mpz & p_m1 = m_powers2(o.sbits+2); - const mpz & p_m2 = m_powers2(o.sbits+3); - (void)p_m1; - - TRACE("mpf_dbg", tout << "p_m1 = " << m_mpz_manager.to_string(p_m1) << std::endl << - "p_m2 = " << m_mpz_manager.to_string(p_m2) << std::endl;); - - bool OVF1 = o.exponent > e_max_norm || // Exponent OVF - (o.exponent == e_max_norm && m_mpz_manager.ge(o.significand, p_m2)); - - TRACE("mpf_dbg", tout << "OVF1 = " << OVF1 << std::endl;); - - int lz = 0; - scoped_mpz t(m_mpq_manager); - m_mpz_manager.set(t, p_m2); - while (m_mpz_manager.gt(t, o.significand)) { - m_mpz_manager.machine_div2k(t, 1); - lz++; + if (beta < e_min) { + // denormal significand/TINY + m_mpz_manager.set(sigma, o.exponent - e_min); + o.exponent = e_min; + } + else { + m_mpz_manager.set(sigma, lz - 1); + o.exponent = beta; } - TRACE("mpf_dbg", tout << "LZ = " << lz << std::endl;); + scoped_mpz sigma_cap(m_mpz_manager); + sigma_cap = o.sbits + 2; + m_mpz_manager.neg(sigma_cap); + if (m_mpz_manager.lt(sigma, sigma_cap)) + m_mpz_manager.set(sigma, sigma_cap); - m_mpz_manager.set(t, o.exponent); - m_mpz_manager.inc(t); - m_mpz_manager.sub(t, lz, t); - m_mpz_manager.set(temporary, e_min_norm); - m_mpz_manager.sub(t, temporary, t); - bool TINY = m_mpz_manager.is_neg(t); - - TRACE("mpf_dbg", tout << "TINY = " << TINY << std::endl;); - - mpf_exp_t alpha = 3 << (o.ebits-2); - mpf_exp_t beta = o.exponent - lz + 1; - - TRACE("mpf_dbg", tout << "alpha = " << alpha << std::endl << - "beta = " << beta << std::endl; ); - - scoped_mpz sigma(m_mpq_manager); - sigma = 0; - - if (TINY && !UNFen) { - m_mpz_manager.set(sigma, o.exponent); - m_mpz_manager.sub(sigma, temporary, sigma); - m_mpz_manager.inc(sigma); - } - else - m_mpz_manager.set(sigma, lz); - - scoped_mpz limit(m_mpq_manager); - limit = o.sbits + 2; - m_mpz_manager.neg(limit); - if (m_mpz_manager.lt(sigma, limit)) { - m_mpz_manager.set(sigma, limit); - } + TRACE("mpf_dbg", tout << "e_min_norm = " << e_min << std::endl; + tout << "e_max_norm = " << e_max << std::endl; + tout << "beta = " << beta << ", (beta < e_min) = " << (beta < e_min) << std::endl; + tout << "LZ = " << lz << std::endl; + tout << "sigma = " << m_mpz_manager.to_string(sigma) << std::endl; + tout << "sigma_cap = " << m_mpz_manager.to_string(sigma_cap) << std::endl;); // Normalization shift TRACE("mpf_dbg", tout << "Shift distance: " << m_mpz_manager.to_string(sigma) << " " << ((m_mpz_manager.is_nonneg(sigma))?"(LEFT)":"(RIGHT)") << std::endl;); - if (m_mpz_manager.le(sigma, 0)) { // Right shift + if (m_mpz_manager.le(sigma, -1)) { + // Right shift scoped_mpz sticky_rem(m_mpz_manager); - unsigned sigma_uint = (unsigned)-m_mpz_manager.get_int64(sigma); // sigma is capped, this is safe. - sigma_uint += 2; - m_mpz_manager.machine_div_rem(o.significand, m_powers2(sigma_uint), o.significand, sticky_rem); + unsigned nsigma_uint = (unsigned)-m_mpz_manager.get_int64(sigma); // sigma is capped, this is safe. + m_mpz_manager.machine_div_rem(o.significand, m_powers2(nsigma_uint), o.significand, sticky_rem); bool sticky = !m_mpz_manager.is_zero(sticky_rem); if (sticky && m_mpz_manager.is_even(o.significand)) m_mpz_manager.inc(o.significand); } - else { // Left shift + else { + // Left shift unsigned sigma_uint = static_cast(m_mpz_manager.get_int64(sigma)); - m_mpz_manager.mul2k(o.significand, sigma_uint - 1, o.significand); - bool sticky = !m_mpz_manager.is_even(o.significand); - m_mpz_manager.machine_div2k(o.significand, 1); - if (sticky && m_mpz_manager.is_even(o.significand)) - m_mpz_manager.inc(o.significand); + m_mpz_manager.mul2k(o.significand, sigma_uint, o.significand); } - TRACE("mpf_dbg", tout << "After sticky: " << to_string(o) << std::endl;); - - if (OVF1 && OVFen) { - o.exponent = beta; - o.exponent -= alpha; - } - else if (TINY && UNFen) { - o.exponent = beta; - o.exponent += alpha; - } - else if (TINY && !UNFen) - o.exponent = e_min_norm; - else - o.exponent = beta; - - TRACE("mpf_dbg", tout << "Shifted: " << to_string(o) << std::endl;); - - const mpz & p_sig = m_powers2(o.sbits); - SASSERT(TINY || (m_mpz_manager.ge(o.significand, p_sig))); + TRACE("mpf_dbg", tout << "Shifted: " << to_string(o) << std::endl; + tout << to_string_binary(o, 1, 3) << std::endl;); // Significand rounding (sigrd) - bool sticky = !m_mpz_manager.is_even(o.significand); // new sticky bit! + bool sticky = !m_mpz_manager.is_even(o.significand); + m_mpz_manager.machine_div2k(o.significand, 1); + sticky = sticky || !m_mpz_manager.is_even(o.significand); m_mpz_manager.machine_div2k(o.significand, 1); bool round = !m_mpz_manager.is_even(o.significand); m_mpz_manager.machine_div2k(o.significand, 1); @@ -2002,56 +2026,47 @@ void mpf_manager::round(mpf_rounding_mode rm, mpf & o) { bool inc = false; switch (rm) { case MPF_ROUND_NEAREST_TEVEN: inc = round && (last || sticky); break; - // case MPF_ROUND_NEAREST_TAWAY: inc = round; break; // CMW: Check - case MPF_ROUND_NEAREST_TAWAY: inc = round && (!last || sticky); break; // CMW: Fix ok? + case MPF_ROUND_NEAREST_TAWAY: inc = round && (!last || sticky); break; case MPF_ROUND_TOWARD_POSITIVE: inc = (!o.sign && (round || sticky)); break; case MPF_ROUND_TOWARD_NEGATIVE: inc = (o.sign && (round || sticky)); break; case MPF_ROUND_TOWARD_ZERO: inc = false; break; default: UNREACHABLE(); } - if (inc) { - TRACE("mpf_dbg", tout << "Rounding increment -> significand +1" << std::endl;); + TRACE("mpf_dbg", tout << "Rounding increment -> significand +" << (int)inc << std::endl;); + if (inc) m_mpz_manager.inc(o.significand); - } - else - TRACE("mpf_dbg", tout << "Rounding increment -> significand +0" << std::endl;); TRACE("mpf_dbg", tout << "Rounded significand: " << to_string(o) << std::endl;); - bool SIGovf = false; - // Post normalization (post) + const mpz & p_sig = m_powers2(o.sbits); if (m_mpz_manager.ge(o.significand, p_sig)) { m_mpz_manager.machine_div2k(o.significand, 1); o.exponent++; } - if (o.exponent > e_max_norm) - SIGovf = true; - + bool SIGovf = o.exponent > e_max; TRACE("mpf_dbg", tout << "Post-normalized: " << to_string(o) << std::endl;); - TRACE("mpf_dbg", tout << "SIGovf = " << SIGovf << std::endl;); // Exponent rounding (exprd) - bool o_has_max_exp = (o.exponent > e_max_norm); + bool o_has_max_exp = (o.exponent > e_max); bool OVF2 = SIGovf && o_has_max_exp; TRACE("mpf_dbg", tout << "OVF2 = " << OVF2 << std::endl;); TRACE("mpf_dbg", tout << "o_has_max_exp = " << o_has_max_exp << std::endl;); - if (!OVFen && OVF2) + if (OVF2) mk_round_inf(rm, o); else { const mpz & p = m_powers2(o.sbits-1); - TRACE("mpf_dbg", tout << "P: " << m_mpz_manager.to_string(p_m1) << std::endl;); if (m_mpz_manager.ge(o.significand, p)) { TRACE("mpf_dbg", tout << "NORMAL: " << m_mpz_manager.to_string(o.significand) << std::endl;); - m_mpz_manager.sub(o.significand, p, o.significand); + m_mpz_manager.sub(o.significand, p, o.significand); // Strips the hidden bit. } else { TRACE("mpf_dbg", tout << "DENORMAL: " << m_mpz_manager.to_string(o.significand) << std::endl;); @@ -2059,7 +2074,8 @@ void mpf_manager::round(mpf_rounding_mode rm, mpf & o) { } } - TRACE("mpf_dbg", tout << "ROUNDED = " << to_string(o) << std::endl;); + TRACE("mpf_dbg", tout << "ROUNDED = " << to_string(o) << std::endl; + tout << to_string_binary(o, -1, 0) << " (hidden bit, unbiased exp)." << std::endl;); } void mpf_manager::round_sqrt(mpf_rounding_mode rm, mpf & o) { diff --git a/src/util/mpf.h b/src/util/mpf.h index 31523c3ed..72c03f447 100644 --- a/src/util/mpf.h +++ b/src/util/mpf.h @@ -285,8 +285,9 @@ protected: }; std::string to_string_raw(mpf const & a); - std::string to_string_hexfloat(mpf const & a); + std::string to_string_hexfloat(mpf const & a); std::string to_string_hexfloat(bool sgn, mpf_exp_t exp, scoped_mpz const & sig, unsigned ebits, unsigned sbits, unsigned rbits); + std::string to_string_binary(mpf const & x, unsigned upper_extra, unsigned lower_extra); public: powers2 m_powers2; }; From 0610392a054f3e3ebed3a9a073c26b1ef961e406 Mon Sep 17 00:00:00 2001 From: "Christoph M. Wintersteiger" Date: Fri, 28 Jul 2017 20:16:13 +0100 Subject: [PATCH 4/5] Bugfix for fp.fma. Fixes #872. --- src/ast/fpa/fpa2bv_converter.cpp | 38 ++++++++------ src/util/mpf.cpp | 85 +++++++++++++++++--------------- src/util/mpf.h | 13 +++++ 3 files changed, 79 insertions(+), 57 deletions(-) diff --git a/src/ast/fpa/fpa2bv_converter.cpp b/src/ast/fpa/fpa2bv_converter.cpp index 0e83d7fea..979fed78e 100644 --- a/src/ast/fpa/fpa2bv_converter.cpp +++ b/src/ast/fpa/fpa2bv_converter.cpp @@ -1476,6 +1476,10 @@ void fpa2bv_converter::mk_fma(func_decl * f, unsigned num, expr * const * args, mk_ite(c71, zero_cond, z, v7); // else comes the fused multiplication. + expr_ref one_1(m), zero_1(m); + one_1 = m_bv_util.mk_numeral(1, 1); + zero_1 = m_bv_util.mk_numeral(0, 1); + unsigned ebits = m_util.get_ebits(f->get_range()); unsigned sbits = m_util.get_sbits(f->get_range()); @@ -1585,27 +1589,20 @@ void fpa2bv_converter::mk_fma(func_decl * f, unsigned num, expr * const * args, dbg_decouple("fpa2bv_fma_add_exp_delta_capped", exp_delta); // Alignment shift with sticky bit computation. - expr_ref shifted_big(m), shifted_f_sig(m), sticky_raw(m); + expr_ref shifted_big(m), shifted_f_sig(m); + expr_ref alignment_sticky_raw(m), alignment_sticky(m); shifted_big = m_bv_util.mk_bv_lshr( m_bv_util.mk_concat(f_sig, m_bv_util.mk_numeral(0, sbits)), m_bv_util.mk_zero_extend((3*sbits)-(ebits+2), exp_delta)); shifted_f_sig = m_bv_util.mk_extract(3*sbits-1, sbits, shifted_big); - sticky_raw = m_bv_util.mk_extract(sbits-1, 0, shifted_big); + alignment_sticky_raw = m_bv_util.mk_extract(sbits-1, 0, shifted_big); + alignment_sticky = m.mk_app(m_bv_util.get_fid(), OP_BREDOR, alignment_sticky_raw.get()); dbg_decouple("fpa2bv_fma_shifted_f_sig", shifted_f_sig); + dbg_decouple("fpa2bv_fma_f_sig_alignment_stickyw", alignment_sticky); + SASSERT(is_well_sorted(m, alignment_sticky)); SASSERT(m_bv_util.get_bv_size(shifted_f_sig) == 2 * sbits); SASSERT(is_well_sorted(m, shifted_f_sig)); - expr_ref sticky(m); - sticky = m_bv_util.mk_zero_extend(2*sbits-1, m.mk_app(m_bv_util.get_fid(), OP_BREDOR, sticky_raw.get())); - SASSERT(is_well_sorted(m, sticky)); - dbg_decouple("fpa2bv_fma_f_sig_sticky_raw", sticky_raw); - dbg_decouple("fpa2bv_fma_f_sig_sticky", sticky); - - expr * or_args[2] = { shifted_f_sig, sticky }; - shifted_f_sig = m_bv_util.mk_bv_or(2, or_args); - SASSERT(is_well_sorted(m, shifted_f_sig)); - dbg_decouple("fpa2bv_fma_f_sig_or_sticky", shifted_f_sig); - // Significant addition. // Two extra bits for the sign and for catching overflows. e_sig = m_bv_util.mk_zero_extend(2, e_sig); @@ -1621,9 +1618,19 @@ void fpa2bv_converter::mk_fma(func_decl * f, unsigned num, expr * const * args, dbg_decouple("fpa2bv_fma_add_e_sig", e_sig); dbg_decouple("fpa2bv_fma_add_shifted_f_sig", shifted_f_sig); - expr_ref sum(m), e_plus_f(m), e_minus_f(m); + expr_ref sum(m), e_plus_f(m), e_minus_f(m), sticky_wide(m); + sticky_wide = m_bv_util.mk_zero_extend(2*sbits+1, alignment_sticky); e_plus_f = m_bv_util.mk_bv_add(e_sig, shifted_f_sig); + e_plus_f = m.mk_ite(m.mk_eq(m_bv_util.mk_extract(0, 0, e_plus_f), zero_1), + m_bv_util.mk_bv_add(e_plus_f, sticky_wide), + e_plus_f); e_minus_f = m_bv_util.mk_bv_sub(e_sig, shifted_f_sig); + e_minus_f = m.mk_ite(m.mk_eq(m_bv_util.mk_extract(0, 0, e_minus_f), zero_1), + m_bv_util.mk_bv_sub(e_minus_f, sticky_wide), + e_minus_f); + SASSERT(is_well_sorted(m, shifted_f_sig)); + dbg_decouple("fpa2bv_fma_f_sig_or_sticky", shifted_f_sig); + m_simp.mk_ite(eq_sgn, e_plus_f, e_minus_f, sum); SASSERT(m_bv_util.get_bv_size(sum) == 2*sbits + 2); SASSERT(is_well_sorted(m, sum)); @@ -1635,8 +1642,7 @@ void fpa2bv_converter::mk_fma(func_decl * f, unsigned num, expr * const * args, dbg_decouple("fpa2bv_fma_add_sign_bv", sign_bv); dbg_decouple("fpa2bv_fma_add_n_sum", n_sum); - expr_ref res_sig_eq(m), sig_abs(m), one_1(m); - one_1 = m_bv_util.mk_numeral(1, 1); + expr_ref res_sig_eq(m), sig_abs(m); m_simp.mk_eq(sign_bv, one_1, res_sig_eq); m_simp.mk_ite(res_sig_eq, n_sum, sum, sig_abs); dbg_decouple("fpa2bv_fma_add_sig_abs", sig_abs); diff --git a/src/util/mpf.cpp b/src/util/mpf.cpp index 36d9e49f1..245186f52 100644 --- a/src/util/mpf.cpp +++ b/src/util/mpf.cpp @@ -556,8 +556,6 @@ void mpf_manager::add_sub(mpf_rounding_mode rm, mpf const & x, mpf const & y, mp SASSERT(exp_delta <= INT_MAX); scoped_mpz sticky_rem(m_mpz_manager); m_mpz_manager.machine_div_rem(b.significand(), m_powers2((int)exp_delta), b.significand(), sticky_rem); - if (!m_mpz_manager.is_zero(sticky_rem) && m_mpz_manager.is_even(b.significand())) - m_mpz_manager.inc(b.significand()); TRACE("mpf_dbg", tout << "A' = " << m_mpz_manager.to_string(a.significand()) << std::endl;); TRACE("mpf_dbg", tout << "B' = " << m_mpz_manager.to_string(b.significand()) << std::endl;); @@ -566,10 +564,14 @@ void mpf_manager::add_sub(mpf_rounding_mode rm, mpf const & x, mpf const & y, mp if (sgn(a) != sgn(b)) { TRACE("mpf_dbg", tout << "SUBTRACTING" << std::endl;); m_mpz_manager.sub(a.significand(), b.significand(), o.significand); + if (!sticky_rem.is_zero() && m_mpz_manager.is_even(o.significand)) + m_mpz_manager.dec(o.significand); } else { TRACE("mpf_dbg", tout << "ADDING" << std::endl;); m_mpz_manager.add(a.significand(), b.significand(), o.significand); + if (!sticky_rem.is_zero() && m_mpz_manager.is_even(o.significand)) + m_mpz_manager.inc(o.significand); } TRACE("mpf_dbg", tout << "sum[-2:sbits+2] = " << m_mpz_manager.to_string(o.significand) << std::endl;); @@ -805,7 +807,7 @@ void mpf_manager::fma(mpf_rounding_mode rm, mpf const & x, mpf const & y, mpf co o.ebits = x.ebits; o.sbits = x.sbits; - scoped_mpf mul_res(*this, x.ebits+2, 2*x.sbits-1); + scoped_mpf mr(*this); scoped_mpf a(*this, x.ebits, x.sbits), b(*this, x.ebits, x.sbits), c(*this, x.ebits, x.sbits); set(a, x); set(b, y); @@ -828,20 +830,18 @@ void mpf_manager::fma(mpf_rounding_mode rm, mpf const & x, mpf const & y, mpf co SASSERT(m_mpz_manager.ge(b.significand(), m_powers2(x.sbits-1))); // A/B in _1.[sbits-1]. - mul_res.get().sign = (a.sign() != b.sign()); - mul_res.get().exponent = a.exponent() + b.exponent(); - m_mpz_manager.mul(a.significand(), b.significand(), mul_res.significand()); + mr.set(x.ebits+2, 2*x.sbits-1, a.sign() != b.sign(), a.exponent() + b.exponent()); + m_mpz_manager.mul(a.significand(), b.significand(), mr.significand()); - TRACE("mpf_dbg", tout << "M = " << to_string(mul_res) << std::endl; - tout << "M = " << to_string_binary(mul_res, 1, 0) << std::endl;); + TRACE("mpf_dbg", tout << "M = " << to_string(mr) << std::endl; + tout << "M = " << to_string_binary(mr, 1, 0) << std::endl;); // mul_res is [-1][0].[2*sbits - 2], i.e., >= 2^(2*sbits-2) and < 2^(2*sbits). - SASSERT(m_mpz_manager.lt(mul_res.significand(), m_powers2(2*x.sbits))); - SASSERT(m_mpz_manager.ge(mul_res.significand(), m_powers2(2*x.sbits - 2))); + SASSERT(m_mpz_manager.lt(mr.significand(), m_powers2(2*x.sbits))); + SASSERT(m_mpz_manager.ge(mr.significand(), m_powers2(2*x.sbits - 2))); // Introduce extra bits into c in _[0].[sbits-1] s.t. c in _[0].[2*sbits-2] - c.get().ebits = x.ebits + 2; - c.get().sbits = 2 * x.sbits - 1; + c.set(x.ebits + 2, 2 * x.sbits - 1, c.sign(), c.exponent(), c.significand()); m_mpz_manager.mul2k(c.significand(), x.sbits - 1); TRACE("mpf_dbg", tout << "C_= " << to_string(c) << std::endl; @@ -851,67 +851,67 @@ void mpf_manager::fma(mpf_rounding_mode rm, mpf const & x, mpf const & y, mpf co SASSERT(m_mpz_manager.is_zero(c.significand()) || m_mpz_manager.ge(c.significand(), m_powers2(2*x.sbits - 2))); - if (exp(c) > exp(mul_res)) { + if (exp(c) > exp(mr)) { TRACE("mpf_dbg", tout << "Swap!" << std::endl;); - mul_res.swap(c); + mr.swap(c); } - mpf_exp_t exp_delta_w = exp(mul_res) - exp(c); - - SASSERT(exp(mul_res) >= exp(c) && exp_delta_w >= 0); + // Alignment shift with sticky bit computation. + mpf_exp_t exp_delta_w = exp(mr) - exp(c); + SASSERT(exp(mr) >= exp(c) && exp_delta_w >= 0); if (exp_delta_w > 2 * x.sbits) exp_delta_w = 2 * x.sbits; - unsigned exp_delta = (unsigned)exp_delta_w; - TRACE("mpf_dbg", tout << "exp_delta = " << exp_delta << std::endl;); - // Alignment shift with sticky bit computation. + unsigned exp_delta = (unsigned)exp_delta_w; + scoped_mpz sticky_rem(m_mpz_manager); m_mpz_manager.machine_div_rem(c.significand(), m_powers2(exp_delta), c.significand(), sticky_rem); - TRACE("mpf_dbg", tout << "alignment shift -> sig = " << m_mpz_manager.to_string(c.significand()) << + TRACE("mpf_dbg", tout << "alignment shift by delta=" << exp_delta << " -> sig = " << m_mpz_manager.to_string(c.significand()) << " sticky_rem = " << m_mpz_manager.to_string(sticky_rem) << std::endl;); - if (!m_mpz_manager.is_zero(sticky_rem) && m_mpz_manager.is_even(c.significand())) - m_mpz_manager.inc(c.significand()); + bool alignment_sticky = !m_mpz_manager.is_zero(sticky_rem); - TRACE("mpf_dbg", tout << "M'= " << m_mpz_manager.to_string(mul_res.significand()) << std::endl; - tout << "M'= " << to_string_binary(mul_res, 1, 0) << std::endl; ); + TRACE("mpf_dbg", tout << "M'= " << m_mpz_manager.to_string(mr.significand()) << std::endl; + tout << "M'= " << to_string_binary(mr, 1, 0) << std::endl; ); TRACE("mpf_dbg", tout << "C'= " << m_mpz_manager.to_string(c.significand()) << std::endl; tout << "C'= " << to_string_binary(c, 1, 0) << std::endl; ); - scoped_mpf res(c); - - res.get().sign = mul_res.sign(); - res.get().exponent = mul_res.exponent(); - // Significand addition + scoped_mpf res(mr); - if (sgn(mul_res) != sgn(c)) { + if (sgn(mr) != sgn(c)) { TRACE("mpf_dbg", tout << "SUBTRACTING" << std::endl;); - m_mpz_manager.sub(mul_res.significand(), c.significand(), res.significand()); + m_mpz_manager.sub(mr.significand(), c.significand(), res.significand()); + + if (alignment_sticky && m_mpz_manager.is_even(res.significand())) + m_mpz_manager.dec(res.get().significand); if (m_mpz_manager.is_neg(res.significand())) { m_mpz_manager.abs(res.significand()); - res.get().sign |= !mul_res.sign() && c.sign(); + res.get().sign |= c.sign(); } } else { TRACE("mpf_dbg", tout << "ADDING" << std::endl;); - m_mpz_manager.add(mul_res.significand(), c.significand(), res.significand()); + m_mpz_manager.add(mr.significand(), c.significand(), res.significand()); + + if (alignment_sticky && m_mpz_manager.is_even(res.significand())) + m_mpz_manager.inc(res.get().significand); } TRACE("mpf_dbg", tout << "S'= " << m_mpz_manager.to_string(res.significand()) << std::endl; tout << "S'= " << to_string_binary(res, 1, 0) << std::endl; ); + bool renorm_sticky = false; + SASSERT(m_mpz_manager.lt(res.significand(), m_powers2(2 * x.sbits + 1))); if (m_mpz_manager.ge(res.significand(), m_powers2(2 * x.sbits))) { SASSERT(exp(res) < mk_max_exp(x.ebits)); // NYI. res.get().exponent++; - bool sticky = !m_mpz_manager.is_even(res.significand()); + renorm_sticky = !m_mpz_manager.is_even(res.significand()); m_mpz_manager.machine_div2k(res.significand(), 1); - if (sticky && m_mpz_manager.is_even(res.significand())) - m_mpz_manager.inc(res.significand()); TRACE("mpf_dbg", tout << "Add/Sub overflew into 4.xxx!" << std::endl; tout << "S*= " << to_string_binary(res, 2, 0) << std::endl;); } @@ -919,17 +919,18 @@ void mpf_manager::fma(mpf_rounding_mode rm, mpf const & x, mpf const & y, mpf co set(o, x.ebits, x.sbits, res.sign(), res.exponent(), mpz(0)); if (x.sbits >= 4) { - m_mpz_manager.set(sticky_rem, 0); m_mpz_manager.machine_div_rem(res.significand(), m_powers2(x.sbits - 4), o.significand, sticky_rem); - if (!m_mpz_manager.is_zero(sticky_rem) && m_mpz_manager.is_even(o.significand)) - m_mpz_manager.inc(o.significand); + renorm_sticky |= !m_mpz_manager.is_zero(sticky_rem); } else { m_mpz_manager.mul2k(res.significand(), 4 - x.sbits, o.significand); } + if (renorm_sticky && m_mpz_manager.is_even(o.significand)) + m_mpz_manager.inc(o.significand); + TRACE("mpf_dbg", tout << "sum[-1:sbits+2] = " << m_mpz_manager.to_string(o.significand) << std::endl; - tout << "R = " << to_string_binary(o, 2, 3) << std::endl;); + tout << "R = " << to_string_binary(o, 1, 3) << std::endl;); if (m_mpz_manager.is_zero(o.significand)) mk_zero(x.ebits, x.sbits, rm == MPF_ROUND_TOWARD_NEGATIVE, o); @@ -1634,6 +1635,8 @@ std::string mpf_manager::to_string_binary(mpf const & x, unsigned upper_extra, u "#b" + std::string(x.sbits - 1, '0') + " " + "(" + (sgn(x) ? "-" : "+") + "zero)"; else { + SASSERT(m_mpz_manager.is_nonneg(sig(x))); + res = std::string("") + "#b" + (sgn(x) ? "1" : "0") + " "; scoped_mpz tmp(m_mpz_manager); diff --git a/src/util/mpf.h b/src/util/mpf.h index 72c03f447..07755e314 100644 --- a/src/util/mpf.h +++ b/src/util/mpf.h @@ -300,6 +300,19 @@ class scoped_mpf : public _scoped_numeral { mpf_exp_t exponent() const { return get().exponent; } unsigned sbits() const { return get().sbits; } void set(unsigned ebits, unsigned sbits) { get().set(ebits, sbits); } + void set(unsigned ebits, unsigned sbits, bool sign, mpf_exp_t exp, mpz & significand) { + get().set(ebits, sbits); + get().exponent = exp; + get().sign = sign; + if (&get().significand != &significand) + m().mpz_manager().set(get().significand, significand); + } + void set(unsigned ebits, unsigned sbits, bool sign, mpf_exp_t exp) { + get().set(ebits, sbits); + get().exponent = exp; + get().sign = sign; + m().mpz_manager().set(get().significand, 0); + } public: scoped_mpf(mpf_manager & m):_scoped_numeral(m) {} scoped_mpf(scoped_mpf const & n):_scoped_numeral(n) {} From e677030b7469933a14c02376a7935a307a919482 Mon Sep 17 00:00:00 2001 From: "Christoph M. Wintersteiger" Date: Fri, 28 Jul 2017 21:39:44 +0100 Subject: [PATCH 5/5] Fixed sign bug in mpf fp.fma. Relates to #872. --- src/util/mpf.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/util/mpf.cpp b/src/util/mpf.cpp index 245186f52..be19cbb9b 100644 --- a/src/util/mpf.cpp +++ b/src/util/mpf.cpp @@ -841,7 +841,7 @@ void mpf_manager::fma(mpf_rounding_mode rm, mpf const & x, mpf const & y, mpf co SASSERT(m_mpz_manager.ge(mr.significand(), m_powers2(2*x.sbits - 2))); // Introduce extra bits into c in _[0].[sbits-1] s.t. c in _[0].[2*sbits-2] - c.set(x.ebits + 2, 2 * x.sbits - 1, c.sign(), c.exponent(), c.significand()); + c.set(x.ebits+2, 2*x.sbits-1, c.sign(), c.exponent(), c.significand()); m_mpz_manager.mul2k(c.significand(), x.sbits - 1); TRACE("mpf_dbg", tout << "C_= " << to_string(c) << std::endl; @@ -888,7 +888,7 @@ void mpf_manager::fma(mpf_rounding_mode rm, mpf const & x, mpf const & y, mpf co if (m_mpz_manager.is_neg(res.significand())) { m_mpz_manager.abs(res.significand()); - res.get().sign |= c.sign(); + res.get().sign = !res.sign(); } } else {