From 19ec6ba298580bb1c650778388a18c85da7c410a Mon Sep 17 00:00:00 2001 From: Dhruv Maroo Date: Tue, 25 Jun 2024 22:06:24 +0530 Subject: [PATCH] Use SoftFloat 3e for implementing arithmetic operations in RzFloat (#4535) * Add a new test for checking 80-bit floating point operations * New test `f80_ieee_div_test` tests the division of two 80-bit floats * Add SoftFloat 2c as a meson subproject * Add softfloat code to make the failing test case pass * Update the hash for the latest softfloat revision * Implement `rz_float_sqrt` using SoftFloat * Run the `f80_ieee_div_test` only for x86 * Replace SoftFloat version 2c with 3e * 3e has less bugs and more features * Modify the implementation in accordance * Update SoftFloat revision and add a guard around the 80-bit div test * Use SoftFloat for add, sub, mul operations as well * Make rem and mod also use SoftFloat functions * Also add test for mod and rem, and fix behavior of rem * Add comment about behavior of mod and rem * Add comments for tests which have different results for mod and rem * Simplify usage of loop variable as suggested in review * Remove unused macro from float.c * Implement `FMA` and `ROUND` using SoftFloat API * Add more tests for 80-bit floats * Change remote to a repository under rizinorg * Add info about the rounding mode in the Doxygen for rem and mod * Add comments in tests for rem and mod in `test_float.c` * Use bitvectors to initialize 80-bit soft floats * This makes the tests more portable and hence they can be run on any platform * Remove guards for f80 tests since they are portable now --- .gitignore | 1 + doc/PACKAGERS.md | 1 + librz/util/float/float.c | 1563 ++++++----------------------- librz/util/float/float_internal.c | 20 - librz/util/meson.build | 2 +- meson.build | 13 + meson_options.txt | 1 + subprojects/softfloat.wrap | 5 + test/unit/test_float.c | 185 +++- 9 files changed, 489 insertions(+), 1302 deletions(-) create mode 100644 subprojects/softfloat.wrap diff --git a/.gitignore b/.gitignore index bf258e2681..2447d8df65 100644 --- a/.gitignore +++ b/.gitignore @@ -130,6 +130,7 @@ subprojects/libmspack/ subprojects/blake3/ subprojects/xz-*/ subprojects/zstd-*/ +subprojects/softfloat/ dist/windows/Output # Core files generated by OpenBSD *.core diff --git a/doc/PACKAGERS.md b/doc/PACKAGERS.md index 5330f5bfd0..fa3dca038f 100644 --- a/doc/PACKAGERS.md +++ b/doc/PACKAGERS.md @@ -81,6 +81,7 @@ At the time of writing, these are: * `use_sys_libmspack` * `use_sys_pcre2` * `use_sys_tree_sitter` +* `use_sys_softfloat` See [meson_options.txt][] for a complete list of compile-time options. diff --git a/librz/util/float/float.c b/librz/util/float/float.c index 96ec3ef65e..66e904c7c3 100644 --- a/librz/util/float/float.c +++ b/librz/util/float/float.c @@ -20,6 +20,7 @@ #include #include #include +#include /** * \defgroup Generate Nan and infinite for float/double/long double @@ -54,6 +55,134 @@ define_types_gen_inf(f64, double); define_types_gen_inf(f128, long double); /**@}*/ +/** \defgroup Helper utilities to inter-operate with SoftFloat. + * @ { + */ +static inline float32_t to_float32(RzFloat *f32) { + rz_warn_if_fail(f32->r == RZ_FLOAT_IEEE754_BIN_32); + float32_t ret = { + .v = rz_bv_to_ut32(f32->s) + }; + return ret; +} + +static inline float64_t to_float64(RzFloat *f64) { + rz_warn_if_fail(f64->r == RZ_FLOAT_IEEE754_BIN_64); + float64_t ret = { + .v = rz_bv_to_ut64(f64->s) + }; + return ret; +} + +static inline extFloat80_t to_float80(RzFloat *f80) { + rz_warn_if_fail(f80->r == RZ_FLOAT_IEEE754_BIN_80); + + extFloat80_t ret; + ret.signif = rz_bv_to_ut64(f80->s); + + ut16 upper = 0; + for (ut8 i = 79; i >= 64; i--) { + upper <<= 1; + upper |= rz_bv_get(f80->s, i); + } + ret.signExp = upper; + + return ret; +} + +static inline float128_t to_float128(RzFloat *f128) { + rz_warn_if_fail(f128->r != RZ_FLOAT_IEEE754_BIN_128); + + float128_t ret; + ret.v[0] = rz_bv_to_ut64(f128->s); + + ut64 upper = 0; + for (ut8 i = 127; i >= 64; i--) { + upper <<= 1; + upper |= rz_bv_get(f128->s, i); + } + ret.v[1] = upper; + + return ret; +} + +static inline RzFloat *set_exception_flags(RzFloat *f) { + if (softfloat_exceptionFlags & softfloat_flag_inexact) { + f->exception |= RZ_FLOAT_E_INEXACT; + } + if (softfloat_exceptionFlags & softfloat_flag_underflow) { + f->exception |= RZ_FLOAT_E_UNDERFLOW; + } + if (softfloat_exceptionFlags & softfloat_flag_overflow) { + f->exception |= RZ_FLOAT_E_OVERFLOW; + } + if (softfloat_exceptionFlags & softfloat_flag_infinite) { + f->exception |= RZ_FLOAT_E_DIV_ZERO; + } + if (softfloat_exceptionFlags & softfloat_flag_invalid) { + f->exception |= RZ_FLOAT_E_INVALID_OP; + } + + softfloat_exceptionFlags = 0; + return f; +} + +static inline RzFloat *of_float32(float32_t f32) { + RzFloat *ret = rz_float_new(RZ_FLOAT_IEEE754_BIN_32); + + rz_bv_set_from_ut64(ret->s, f32.v); + return set_exception_flags(ret); +} + +static inline RzFloat *of_float64(float64_t f64) { + RzFloat *ret = rz_float_new(RZ_FLOAT_IEEE754_BIN_64); + + rz_bv_set_from_ut64(ret->s, f64.v); + return set_exception_flags(ret); +} + +static inline RzFloat *of_float80(extFloat80_t f80) { + RzFloat *ret = rz_float_new(RZ_FLOAT_IEEE754_BIN_80); + + rz_bv_set_from_ut64(ret->s, f80.signif); + + ut16 upper = f80.signExp; + for (ut8 i = 0; i < 16; i++) { + rz_bv_set(ret->s, 64 + i, upper & 1); + upper >>= 1; + } + + return set_exception_flags(ret); +} + +static inline RzFloat *of_float128(float128_t f128) { + RzFloat *ret = rz_float_new(RZ_FLOAT_IEEE754_BIN_128); + + rz_bv_set_from_ut64(ret->s, f128.v[0]); + + ut64 upper = f128.v[1]; + for (ut8 i = 0; i < 64; i++) { + rz_bv_set(ret->s, 64 + i, upper & 1); + upper >>= 1; + } + + return set_exception_flags(ret); +} + +static int8_t rounding_mode_mapping[] = { + [RZ_FLOAT_RMODE_RNE] = softfloat_round_near_even, + [RZ_FLOAT_RMODE_RNA] = softfloat_round_near_maxMag, + [RZ_FLOAT_RMODE_RTP] = softfloat_round_max, + [RZ_FLOAT_RMODE_RTN] = softfloat_round_min, + [RZ_FLOAT_RMODE_RTZ] = softfloat_round_minMag, + [RZ_FLOAT_RMODE_UNK] = 6, +}; + +static inline void set_float_rounding_mode(RzFloatRMode mode) { + softfloat_roundingMode = rounding_mode_mapping[mode]; +} +/**@}*/ + /** * \brief return the bitvector string of a float * \param f float @@ -231,23 +360,6 @@ RZ_API RZ_OWN char *rz_float_as_dec_string(RZ_NULLABLE RzFloat *f) { return rz_str_newf("%" LDBLFMTg, result); } -/* - * Common NaN and Inf detection - * */ -#define PROC_SPECIAL_FLOAT_START(left, right) \ - { \ - RzFloatSpec l_type, r_type; \ - l_type = rz_float_detect_spec((left)); \ - r_type = rz_float_detect_spec((right)); \ - bool l_is_inf = (l_type == RZ_FLOAT_SPEC_PINF || l_type == RZ_FLOAT_SPEC_NINF); \ - bool r_is_inf = (r_type == RZ_FLOAT_SPEC_PINF || r_type == RZ_FLOAT_SPEC_NINF); \ - bool l_is_nan = (l_type == RZ_FLOAT_SPEC_SNAN || l_type == RZ_FLOAT_SPEC_QNAN); \ - bool r_is_nan = (r_type == RZ_FLOAT_SPEC_SNAN || r_type == RZ_FLOAT_SPEC_QNAN); \ - bool l_is_zero = l_type == RZ_FLOAT_SPEC_ZERO; \ - bool r_is_zero = r_type == RZ_FLOAT_SPEC_ZERO; - -#define PROC_SPECIAL_FLOAT_END } - /** * \brief Get const attributes from float * \param format RzFloatFormat, format of a float @@ -1041,355 +1153,6 @@ RZ_API RZ_OWN RzFloat *rz_float_new_snan(RzFloatFormat format) { return ret; } -/** - * \brief propagate NaN and trigger signal (set exception for a NaN), - * used in float arithmetic to deal with NaN operand - */ -static RZ_OWN RzFloat *propagate_float_nan(RZ_NONNULL RzFloat *left, RzFloatSpec ltype, RZ_NONNULL RzFloat *right, RzFloatSpec rtype) { - bool l_is_sig_nan = ltype == RZ_FLOAT_SPEC_SNAN; - bool r_is_sig_nan = rtype == RZ_FLOAT_SPEC_SNAN; - - // gen a quiet NaN for return - RzFloatFormat format = left->r; - RzFloat *ret = rz_float_new(left->r); - RzBitVector *bv = ret->s; - ut32 exp_start = rz_float_get_format_info(format, RZ_FLOAT_INFO_MAN_LEN); - ut32 exp_end = exp_start + rz_float_get_format_info(format, RZ_FLOAT_INFO_EXP_LEN); - - // set exponent part to all 1 - rz_bv_set_range(bv, exp_start, exp_end - 1, true); - - // set is_quiet to 1 - rz_bv_set(bv, exp_start - 1, true); - - // signal an exception - if (l_is_sig_nan || r_is_sig_nan) { - ret->exception |= RZ_FLOAT_E_INVALID_OP; - } - - return ret; -} - -/** - * \brief add magnitude (absolute value) - */ -static RZ_OWN RzFloat *fadd_mag(RZ_NONNULL RzFloat *left, RZ_NONNULL RzFloat *right, bool sign, RzFloatRMode mode) { - RzFloat *result = NULL; - - /* Process NaN and Inf cases */ - PROC_SPECIAL_FLOAT_START(left, right) - // propagate NaN - if (l_is_nan || r_is_nan) { - return propagate_float_nan(left, l_type, right, r_type); - } - - if (l_is_inf || r_is_inf) { - // inf + inf = inf - return rz_float_new_inf(left->r, sign); - } - - if (l_is_zero || r_is_zero) { - return rz_float_dup(l_is_zero ? right : left); - } - PROC_SPECIAL_FLOAT_END - - /* Process normal cases */ - // Extract attribute from format - RzFloatFormat format = left->r; - ut32 exp_len = rz_float_get_format_info(format, RZ_FLOAT_INFO_EXP_LEN); - ut32 total_len = rz_float_get_format_info(format, RZ_FLOAT_INFO_TOTAL_LEN); - - // Extract fields from num - RzBitVector *l_exp_squashed = get_exp_squashed(left->s, left->r); - RzBitVector *r_exp_squashed = get_exp_squashed(right->s, right->r); - RzBitVector *l_mantissa = get_man(left->s, left->r); - RzBitVector *r_mantissa = get_man(right->s, right->r); - - if (!l_exp_squashed || !r_exp_squashed || !l_mantissa || !r_mantissa) { - RZ_LOG_ERROR("float: fadd: Error when parsing RzFloat\n"); - return NULL; - } - - RzBitVector *l_borrowed_sig = l_mantissa; - RzBitVector *r_borrowed_sig = r_mantissa; - RzBitVector *result_sig = NULL; - RzBitVector *exp_one = rz_bv_new_one(exp_len); - bool unused; - - // Handle normal float add - ut32 l_exp_val = rz_bv_to_ut32(l_exp_squashed); - ut32 r_exp_val = rz_bv_to_ut32(r_exp_squashed); - st32 exp_diff = (st32)(l_exp_val - r_exp_val); - ut32 abs_exp_diff = exp_diff; - ut32 l_borrow_exp_val = l_exp_val; - ut32 r_borrow_exp_val = r_exp_val; - - // left shift to prevent some tail bits being discard during calculating - // should reserve 3 bits before mantissa : ABCM MMMM MMMM MMMM ... - // C : for the hidden significant bit - // B : carry bit - // A : a space for possible overflow during rounding - // M : represent for mantissa bits - ut32 shift_dist = (exp_len + 1) - 3; // mantissa have (exp_len + sign_len) free bits, and then reserve 3 bits - ut32 hidden_bit_pos = total_len - 3; // the 3rd bit counted from MSB - ut32 carry_bit_pos = total_len - 2; // the 2nd bit counted from MSB - - if (exp_diff == 0) { - // normalized float, hidden bit is 1, recover it in significant - // 1.MMMM MMMM ... - if (l_borrow_exp_val != 0) { - rz_bv_lshift(l_mantissa, shift_dist); - rz_bv_lshift(r_mantissa, shift_dist); - rz_bv_set(l_borrowed_sig, hidden_bit_pos, true); - rz_bv_set(r_borrowed_sig, hidden_bit_pos, true); - } else { - // sub-normal + sub-normal - // sub-normal float, hidden bit is 0, so we do nothing to sigs - // 0.MMMM MMMM ... - // calculate and then pack to return - result = RZ_NEW0(RzFloat); - result->r = format; - result->s = rz_bv_add(left->s, r_mantissa, &unused); - goto clean; - } - } else { // exp_diff != 0 - rz_bv_lshift(l_mantissa, shift_dist); - rz_bv_lshift(r_mantissa, shift_dist); - // should align exponent, chose the max(l_exp, r_exp) as final exp - if (exp_diff < 0) { - // swap to keep l_exp > r_exp - l_borrowed_sig = r_mantissa; - r_borrowed_sig = l_mantissa; - l_borrow_exp_val = r_exp_val; - r_borrow_exp_val = l_exp_val; - abs_exp_diff = -exp_diff; - } - - // check if the small one (right) is normalized ? - if (r_borrow_exp_val != 0) { - // normalized, and then we recover the leading bit 1 - // 1.MMMM MMMM ... - rz_bv_set(r_borrowed_sig, hidden_bit_pos, true); - } else { - // sub-normal (or denormalized float) case - // in IEEE, the value of exp is (1 - bias) for sub-normal, instead of (0 - bias) - // but we considered it as (0 - bias) when calculate the exp_diff = l_exp_field - r_exp_field - // we should r-shift (l_exp_field - bias) - (1 - bias) = l_exp_field - 1, - // but we r-shift (l_exp_field - bias) - (0 - bias) = l_exp_filed - // thus we need to l-shift 1 bit to fix this incompatible - rz_bv_lshift(r_borrowed_sig, 1); - } - - // revealed the hidden bit of the bigger one : 1.MMMM - rz_bv_set(l_borrowed_sig, hidden_bit_pos, true); - // aligned exponent, and generate sticky bit - rz_bv_shift_right_jammed(r_borrowed_sig, abs_exp_diff); - } - - // set result exponent - ut32 result_exp_val = l_borrow_exp_val; - - // now l_exp == r_exp - // calculate significant - result_sig = rz_bv_add(l_borrowed_sig, r_borrowed_sig, &unused); - - if (rz_bv_get(result_sig, carry_bit_pos)) { - result_exp_val += 1; - } - - // round - result = round_float_bv_new( - sign, - result_exp_val, - result_sig, - format, - format, - mode); - -// clean -clean: - rz_bv_free(l_exp_squashed); - rz_bv_free(l_mantissa); - rz_bv_free(r_exp_squashed); - rz_bv_free(r_mantissa); - rz_bv_free(result_sig); - rz_bv_free(exp_one); - return result; -} - -/** - * \brief sub magnitude (absolute value) - */ -static RZ_OWN RzFloat *fsub_mag(RZ_NONNULL RzFloat *left, RZ_NONNULL RzFloat *right, bool sign, RzFloatRMode mode) { - RzFloat *result = NULL; - - /* Process NaN and Inf cases */ - PROC_SPECIAL_FLOAT_START(left, right) - // propagate NaN - if (l_is_nan || r_is_nan) { - return propagate_float_nan(left, l_type, right, r_type); - } - - bool l_sign = rz_float_is_negative(left); - bool r_sign = rz_float_is_negative(right); - if (l_is_inf || r_is_inf) { - if (l_is_inf && r_is_inf) { - // +inf - inf = NaN - return rz_float_new_qnan(left->r); - } - return l_is_inf ? rz_float_new_inf(left->r, l_sign) : rz_float_new_inf(left->r, r_sign); - } - - if (l_is_zero || r_is_zero) { - RzFloat *ret_spec = rz_float_dup(l_is_zero ? right : left); - if (l_is_zero) { - rz_bv_set(ret_spec->s, ret_spec->s->len - 1, !r_sign); - } - return ret_spec; - } - PROC_SPECIAL_FLOAT_END - - // Extract attribute from format - RzFloatFormat format = left->r; - ut32 exp_len = rz_float_get_format_info(format, RZ_FLOAT_INFO_EXP_LEN); - ut32 total_len = rz_float_get_format_info(format, RZ_FLOAT_INFO_TOTAL_LEN); - - // Extract fields from num - RzBitVector *l_exp_squashed = get_exp_squashed(left->s, left->r); - RzBitVector *r_exp_squashed = get_exp_squashed(right->s, right->r); - RzBitVector *l_mantissa = get_man(left->s, left->r); - RzBitVector *r_mantissa = get_man(right->s, right->r); - - if (!l_exp_squashed || !r_exp_squashed || !l_mantissa || !r_mantissa) { - RZ_LOG_ERROR("float: fsub: Error when parsing RzFloat\n"); - rz_bv_free(l_exp_squashed); - rz_bv_free(r_exp_squashed); - rz_bv_free(l_mantissa); - rz_bv_free(r_mantissa); - return NULL; - } - - RzBitVector *l_borrowed_sig = l_mantissa; - RzBitVector *r_borrowed_sig = r_mantissa; - RzBitVector *result_sig = NULL, *result_exp_squashed = NULL; - bool unused; - - // Handle normal float add - ut32 l_exp_val = rz_bv_to_ut32(l_exp_squashed); - ut32 r_exp_val = rz_bv_to_ut32(r_exp_squashed); - st32 exp_diff = (st32)(l_exp_val - r_exp_val); - ut32 abs_exp_diff = exp_diff; - ut32 l_borrow_exp_val = l_exp_val; - ut32 r_borrow_exp_val = r_exp_val; - st32 res_exp_val; - - // similar to `add`, but remember that sub would never produce a carry bit - // we create ABMM MMMM MMMM MMMM ... - // B : for the leading significant bit - // A : space - ut32 shift_dist = (exp_len + 1) - 2; // mantissa have (exp_len + sign_len) free bits, and then reserve 2 bits - ut32 hidden_bit_pos = total_len - 2; // the 2nd bit counted from MSB - - // if l_exp = r_exp - if (exp_diff == 0) { - // compare result - ut8 sdiff_neg = rz_bv_ule(l_mantissa, r_mantissa); - ut8 sdiff_pos = rz_bv_ule(r_mantissa, l_mantissa); - ut8 sig_diff_is_zero = sdiff_neg && sdiff_pos; - RzBitVector *sig_diff = NULL; - if (sig_diff_is_zero) { - // pack to return, exp = 0, sig = 0 - result = RZ_NEW0(RzFloat); - result->r = format; - result->s = rz_bv_new_zero(total_len); - rz_bv_set(result->s, total_len - 1, mode == RZ_FLOAT_RMODE_RTN); - goto clean; - } - - // calculate the correct sig diff - if (sdiff_neg) { - sign = !sign; - sig_diff = rz_bv_sub(r_mantissa, l_mantissa, &unused); - } else { - sig_diff = rz_bv_sub(l_mantissa, r_mantissa, &unused); - } - - // normalize sig - // clz - exp_len - sign_len + 1 (reserve the leading bit) = clz - exp_len - shift_dist = rz_bv_clz(sig_diff) - exp_len; - res_exp_val = (st32)(l_exp_val - shift_dist); - if (res_exp_val < 0) { - // too tiny after shifting, limit to exp_A - shift_dist = l_exp_val; - res_exp_val = 0; - } - // normalize sig diff, reveal the hidden bit pos - rz_bv_lshift(sig_diff, shift_dist); - - result_exp_squashed = rz_bv_new_from_ut64(l_exp_squashed->len, res_exp_val); - result = RZ_NEW0(RzFloat); - result->r = format; - result->s = pack_float_bv(sign, result_exp_squashed, sig_diff, format); - - rz_bv_free(sig_diff); - goto clean; - } else { - rz_bv_lshift(l_mantissa, shift_dist); - rz_bv_lshift(r_mantissa, shift_dist); - // l_exp != r_exp - if (exp_diff < 0) { - // swap to keep l_exp > r_exp - l_borrow_exp_val = r_exp_val; - r_borrow_exp_val = l_exp_val; - l_borrowed_sig = r_mantissa; - r_borrowed_sig = l_mantissa; - abs_exp_diff = -exp_diff; - sign = !sign; - } - - // check if the small one (right) is sub-normal ? - if (r_borrow_exp_val != 0) { - // normalized, and then we recover the leading bit 1 - // 1.MMMM MMMM ... - rz_bv_set(r_borrowed_sig, hidden_bit_pos, true); - } - - // revealed the hidden bit of the bigger one : 1.MMMM - rz_bv_set(l_borrowed_sig, hidden_bit_pos, true); - // aligned exponent, and generate sticky bit - rz_bv_shift_right_jammed(r_borrowed_sig, abs_exp_diff); - } - - // result_exp = bigger_exp - res_exp_val = l_borrow_exp_val; - // result_sig = bigger_sig - small_sig - result_sig = rz_bv_sub(l_borrowed_sig, r_borrowed_sig, &unused); - - ut32 borrow_pos = hidden_bit_pos; - if (!rz_bv_get(result_sig, borrow_pos)) { - // borrow happens - res_exp_val -= 1; - } - - result = round_float_bv_new( - sign, - res_exp_val, - result_sig, - format, - format, - mode); - -clean: - rz_bv_free(l_exp_squashed); - rz_bv_free(l_mantissa); - rz_bv_free(r_exp_squashed); - rz_bv_free(r_mantissa); - rz_bv_free(result_exp_squashed); - rz_bv_free(result_sig); - - return result; -} - /** * \defgroup rz_float_arithmetic_group Arithmetic Operations * implements add, sub, mul, div, fma, rem, sqrt for binary32/binary64/binary128 @@ -1402,12 +1165,24 @@ clean: * \return result of arithmetic operation */ RZ_API RZ_OWN RzFloat *rz_float_add_ieee_bin(RZ_NONNULL RzFloat *left, RZ_NONNULL RzFloat *right, RzFloatRMode mode) { - bool l_sign = rz_float_is_negative(left); - bool r_sign = rz_float_is_negative(right); - if (l_sign == r_sign) { - return fadd_mag(left, right, l_sign, mode); + rz_return_val_if_fail(left && right && left->r == right->r, NULL); + + RzFloatFormat format = left->r; + set_float_rounding_mode(mode); + + switch (format) { + case RZ_FLOAT_IEEE754_BIN_32: + return of_float32(f32_add(to_float32(left), to_float32(right))); + case RZ_FLOAT_IEEE754_BIN_64: + return of_float64(f64_add(to_float64(left), to_float64(right))); + case RZ_FLOAT_IEEE754_BIN_80: + return of_float80(extF80_add(to_float80(left), to_float80(right))); + case RZ_FLOAT_IEEE754_BIN_128: + return of_float128(f128_add(to_float128(left), to_float128(right))); + default: + RZ_LOG_ERROR("float: ADD operation unimplemented for format %d\n", format); + return NULL; } - return fsub_mag(left, right, l_sign, mode); } /** @@ -1416,12 +1191,24 @@ RZ_API RZ_OWN RzFloat *rz_float_add_ieee_bin(RZ_NONNULL RzFloat *left, RZ_NONNUL * \return result of arithmetic operation */ RZ_API RZ_OWN RzFloat *rz_float_sub_ieee_bin(RZ_NONNULL RzFloat *left, RZ_NONNULL RzFloat *right, RzFloatRMode mode) { - bool l_sign = rz_float_is_negative(left); - bool r_sign = rz_float_is_negative(right); - if (l_sign == r_sign) { - return fsub_mag(left, right, l_sign, mode); + rz_return_val_if_fail(left && right && left->r == right->r, NULL); + + RzFloatFormat format = left->r; + set_float_rounding_mode(mode); + + switch (format) { + case RZ_FLOAT_IEEE754_BIN_32: + return of_float32(f32_sub(to_float32(left), to_float32(right))); + case RZ_FLOAT_IEEE754_BIN_64: + return of_float64(f64_sub(to_float64(left), to_float64(right))); + case RZ_FLOAT_IEEE754_BIN_80: + return of_float80(extF80_sub(to_float80(left), to_float80(right))); + case RZ_FLOAT_IEEE754_BIN_128: + return of_float128(f128_sub(to_float128(left), to_float128(right))); + default: + RZ_LOG_ERROR("float: SUB operation unimplemented for format %d\n", format); + return NULL; } - return fadd_mag(left, right, l_sign, mode); } /** @@ -1430,129 +1217,24 @@ RZ_API RZ_OWN RzFloat *rz_float_sub_ieee_bin(RZ_NONNULL RzFloat *left, RZ_NONNUL * \return result of arithmetic operation */ RZ_API RZ_OWN RzFloat *rz_float_mul_ieee_bin(RZ_NONNULL RzFloat *left, RZ_NONNULL RzFloat *right, RzFloatRMode mode) { - RzFloat *result = NULL; + rz_return_val_if_fail(left && right && left->r == right->r, NULL); - /* Process NaN and Inf cases */ - PROC_SPECIAL_FLOAT_START(left, right) - // propagate NaN - if (l_is_nan || r_is_nan) { - return propagate_float_nan(left, l_type, right, r_type); - } - - bool l_sign = rz_float_is_negative(left); - bool r_sign = rz_float_is_negative(right); - bool spec_sign = l_sign ^ r_sign; - - if (l_is_inf) { - return r_is_zero ? rz_float_new_qnan(left->r) : rz_float_new_inf(left->r, spec_sign); - } - - if (r_is_inf) { - return l_is_zero ? rz_float_new_qnan(left->r) : rz_float_new_inf(left->r, spec_sign); - } - - if (l_is_zero || r_is_zero) { - // 0 * x = 0 - return rz_float_new(left->r); - } - PROC_SPECIAL_FLOAT_END - - // Extract attribute from format RzFloatFormat format = left->r; - ut32 exp_len = rz_float_get_format_info(format, RZ_FLOAT_INFO_EXP_LEN); - ut32 total_len = rz_float_get_format_info(format, RZ_FLOAT_INFO_TOTAL_LEN); - ut32 bias = rz_float_get_format_info(format, RZ_FLOAT_INFO_BIAS); - ut32 extra_len = total_len; + set_float_rounding_mode(mode); - // Extract fields from num - RzBitVector *l_exp_squashed = get_exp_squashed(left->s, left->r); - RzBitVector *r_exp_squashed = get_exp_squashed(right->s, right->r); - RzBitVector *l_mantissa = get_man_stretched(left->s, left->r); - RzBitVector *r_mantissa = get_man_stretched(right->s, right->r); - RzBitVector *result_sig = NULL, *result_exp_squashed = NULL; - bool l_sign = get_sign(left->s, left->r); - bool r_sign = get_sign(right->s, right->r); - bool result_sign = l_sign ^ r_sign; - - // Handle normal float multiply - ut32 l_exp_val = rz_bv_to_ut32(l_exp_squashed); - ut32 r_exp_val = rz_bv_to_ut32(r_exp_squashed); - ut32 shift_dist; - - // no biased one - rz_bv_lshift(l_mantissa, exp_len - 1); - rz_bv_lshift(r_mantissa, exp_len - 1); - - st32 lexp_nobias = rz_float_get_exponent_val_no_bias(left); - st32 rexp_nobias = rz_float_get_exponent_val_no_bias(right); - st32 result_exp_val = lexp_nobias + rexp_nobias; - - // remember we would like to make 01.MM MMMM ... (but leave higher extra bits empty) - ut32 hidden_bit_pos = total_len - 2; - - // set leading bit - if (l_exp_val != 0) { - rz_bv_set(l_mantissa, hidden_bit_pos, true); + switch (format) { + case RZ_FLOAT_IEEE754_BIN_32: + return of_float32(f32_mul(to_float32(left), to_float32(right))); + case RZ_FLOAT_IEEE754_BIN_64: + return of_float64(f64_mul(to_float64(left), to_float64(right))); + case RZ_FLOAT_IEEE754_BIN_80: + return of_float80(extF80_mul(to_float80(left), to_float80(right))); + case RZ_FLOAT_IEEE754_BIN_128: + return of_float128(f128_mul(to_float128(left), to_float128(right))); + default: + RZ_LOG_ERROR("float: MUL operation unimplemented for format %d\n", format); + return NULL; } - - if (r_exp_val != 0) { - rz_bv_set(r_mantissa, hidden_bit_pos, true); - } - - // multiplication - // since operands have 0H.MMMM... form, and 0H.MMMMM... - // result would be 00XX.MMMM... - result_sig = rz_bv_mul(l_mantissa, r_mantissa); - - // check if a carry happen, if not, l-shift to force a leading 1 - // check MSB and the bit after MSB - if (rz_bv_get(result_sig, total_len + extra_len - 3)) { - // carry case, think about 01.10 * 01.10 => 0001.0010 - // 001X.MMMM... -> 001.0MMMMM.. - result_exp_val += 1; - rz_bv_shift_right_jammed(result_sig, 1); - } - - // check result and normalize it if needed - ut32 clz = rz_bv_clz(result_sig); - if (clz > 3) { - // means there are sub normal as factor - // try shift - shift_dist = clz - 3; - if (result_exp_val - (st32)shift_dist < 1 - bias) { - // too small, represent as sub-normal - shift_dist = result_exp_val - (1 - bias); - } - rz_bv_lshift(result_sig, shift_dist); - - // biased one - result_exp_val = 0; - - // for those who may be sub-normal, use fake hidden bit for rounding - // note that result sig has 000H.MMMM... form - rz_bv_set(result_sig, rz_bv_len(result_sig) - 4, true); - } - // others has 0001.MMMM... - else { - result_exp_val += bias; - } - - result = round_float_bv_new( - result_sign, - result_exp_val, - result_sig, - format, - format, - mode); - - rz_bv_free(l_exp_squashed); - rz_bv_free(r_exp_squashed); - rz_bv_free(l_mantissa); - rz_bv_free(r_mantissa); - rz_bv_free(result_exp_squashed); - rz_bv_free(result_sig); - - return result; } /** @@ -1567,438 +1249,98 @@ RZ_API RZ_OWN RzFloat *rz_float_mul_ieee_bin(RZ_NONNULL RzFloat *left, RZ_NONNUL * \return result of arithmetic operation */ RZ_API RZ_OWN RzFloat *rz_float_div_ieee_bin(RZ_NONNULL RzFloat *left, RZ_NONNULL RzFloat *right, RzFloatRMode mode) { - RzFloat *result = NULL; + rz_return_val_if_fail(left && right && left->r == right->r, NULL); - PROC_SPECIAL_FLOAT_START(left, right) - bool l_sign = rz_float_is_negative(left); - bool r_sign = rz_float_is_negative(right); - bool sign = l_sign ^ r_sign; - RzFloat *spec_ret = NULL; - - if (l_is_nan || r_is_nan) { - return rz_float_new_qnan(left->r); - } - - if (l_is_inf) { - if (!r_is_inf) { - return rz_float_new_inf(left->r, sign); - } else { - spec_ret = rz_float_new_qnan(left->r); - spec_ret->exception |= RZ_FLOAT_E_INVALID_OP; - return spec_ret; - } - } else { - if (r_is_inf) { - return rz_float_new_zero(left->r); - } - } - - if (l_is_zero) { - if (r_is_zero) { - spec_ret = rz_float_new_qnan(left->r); - spec_ret->exception |= RZ_FLOAT_E_INVALID_OP; - return spec_ret; - } else { - return rz_float_new(left->r); - } - } else { - if (r_is_zero) { - return rz_float_new_inf(left->r, sign); - } - } - PROC_SPECIAL_FLOAT_END - - // Extract attribute from format RzFloatFormat format = left->r; - ut32 exp_len = rz_float_get_format_info(format, RZ_FLOAT_INFO_EXP_LEN); - ut32 total_len = rz_float_get_format_info(format, RZ_FLOAT_INFO_TOTAL_LEN); - ut32 bias = rz_float_get_format_info(format, RZ_FLOAT_INFO_BIAS); - ut32 extra_len = total_len; + set_float_rounding_mode(mode); - // Extract fields from num - RzBitVector *l_exp_squashed = get_exp_squashed(left->s, left->r); - RzBitVector *r_exp_squashed = get_exp_squashed(right->s, right->r); - RzBitVector *l_mantissa = get_man_stretched(left->s, left->r); - RzBitVector *r_mantissa = get_man_stretched(right->s, right->r); - RzBitVector *result_sig = NULL, *result_exp_squashed = NULL; - bool l_sign = get_sign(left->s, left->r); - bool r_sign = get_sign(right->s, right->r); - bool result_sign = l_sign ^ r_sign; - - // Handle normal float multiply - ut32 l_exp_val = rz_bv_to_ut32(l_exp_squashed); - ut32 r_exp_val = rz_bv_to_ut32(r_exp_squashed); - ut32 shift_dist; - - // normalize sub-normal num - // similar to multiplication - if (l_exp_val == 0) { - // is sub-normal - shift_dist = rz_bv_clz(l_mantissa) - (1 + exp_len) + 1 - extra_len; - l_exp_val = 1 - shift_dist; - rz_bv_lshift(l_mantissa, shift_dist); + switch (format) { + case RZ_FLOAT_IEEE754_BIN_32: + return of_float32(f32_div(to_float32(left), to_float32(right))); + case RZ_FLOAT_IEEE754_BIN_64: + return of_float64(f64_div(to_float64(left), to_float64(right))); + case RZ_FLOAT_IEEE754_BIN_80: + return of_float80(extF80_div(to_float80(left), to_float80(right))); + case RZ_FLOAT_IEEE754_BIN_128: + return of_float128(f128_div(to_float128(left), to_float128(right))); + default: + RZ_LOG_ERROR("float: DIV operation unimplemented for format %d\n", format); + return NULL; } - - if (r_exp_val == 0) { - // is sub-normal - shift_dist = rz_bv_clz(r_mantissa) - (1 + exp_len) + 1 - extra_len; - r_exp_val = 1 - shift_dist; - rz_bv_lshift(r_mantissa, shift_dist); - } - - ut32 result_exp_val = l_exp_val - r_exp_val + bias; - - // remember we would like to make the pattern 01.MM MMMM ... - shift_dist = (exp_len + 1) - 2; - ut32 hiddent_bit_pos = total_len - (1 + exp_len); - - // set leading bit - rz_bv_set(l_mantissa, hiddent_bit_pos, true); - rz_bv_set(r_mantissa, hiddent_bit_pos, true); - - // shift to make sure left is large enough to div - // Fx = Mx * 2^x, Fy = My * 2^y - // we have Mx as 01MM MMMM MMMM ... - // now expand left operand to have more bits - // dividend 01MM ..MM 0000 0000 0000 ... - // divisor 00...0000 01MM MMMM MMMM ... - rz_bv_lshift(l_mantissa, shift_dist + extra_len); - rz_bv_lshift(r_mantissa, shift_dist); - - // both dividend and divisor have the form 1.MM... - // and thus the first bit-1 must be set in - // a. LSB of extra bits (dividend sig >= divisor sig) - // b. MSB of original bits (dividend sig < divisor sig) - // the clz should be 31 or 32 respectively - result_sig = rz_bv_div(l_mantissa, r_mantissa); - ut32 clz = rz_bv_clz(result_sig); - - // check if normalization needed - shift_dist = clz == extra_len ? 1 : 0; - - // Convert to original length bitvector - // normalize it - // and make 01MM MMMM MMMM ... format - rz_bv_shift_right_jammed(result_sig, 2 - shift_dist); - - // dec exp according to normalization - // exp -= shift - result_exp_val -= shift_dist; - - if ((st32)result_exp_val < 0) { - // underflow ? - result_exp_val = 0; - } - - result = round_float_bv_new( - result_sign, - result_exp_val, - result_sig, - format, - format, - mode); - - rz_bv_free(l_exp_squashed); - rz_bv_free(r_exp_squashed); - rz_bv_free(l_mantissa); - rz_bv_free(r_mantissa); - rz_bv_free(result_exp_squashed); - rz_bv_free(result_sig); - - return result; } /** - * \brief calculate remainder of \p left % \p right and round the result after + * \brief Returns the value of \p left % \p right, with quotient rounded to an integer with rounding mode RNE * \details * Any % 0 => NaN * Inf % Any => NaN, invalid * Any % Inf -> Any * 0 % Any -> 0 - * \param quo_rnd quotient round mode, fmod use RTZ, frem use RNE - * \param mode rounding mode - * \return result of arithmetic operation - */ -static RZ_OWN RzFloat *rz_float_rem_internal(RZ_NONNULL RzFloat *left, RZ_NONNULL RzFloat *right, RzFloatRMode quo_rnd, RzFloatRMode mode) { - PROC_SPECIAL_FLOAT_START(left, right) - RzFloat *spec_ret = NULL; - - if (l_is_nan || r_is_nan) { - return rz_float_new_qnan(left->r); - } - - if (l_is_inf || r_is_zero) { - spec_ret = rz_float_new_qnan(left->r); - spec_ret->exception |= RZ_FLOAT_E_INVALID_OP; - return spec_ret; - } - - if (r_is_inf) { - return rz_float_dup(left); - } - - if (l_is_zero) { - return rz_float_new_zero(left->r); - } - PROC_SPECIAL_FLOAT_END - - // extract info from args - // left = mx * 2^(ex), right = my * 2^(ey) - RzBitVector *mx = rz_float_get_mantissa(left); - RzBitVector *my = rz_float_get_mantissa(right); - RzBitVector *exp_x = rz_float_get_exponent(left); - RzBitVector *exp_y = rz_float_get_exponent(right); - ut32 bias = rz_float_get_format_info(left->r, RZ_FLOAT_INFO_BIAS); - st32 ex = (st32)(rz_bv_to_ut32(exp_x) - bias); - st32 ey = (st32)(rz_bv_to_ut32(exp_y) - bias); - rz_bv_free(exp_x); - rz_bv_free(exp_y); - - bool sign_x = rz_float_is_negative(left); - - /* quo(-x,-y) = quo(x,y), rem(-x,-y) = -rem(x,y) - * quo(-x,y) = -quo(x,y), rem(-x,y) = -rem(x,y) - * thus quo = sign(x/y)*quo(|x|,|y|), rem = sign(x)*rem(|x|,|y|) */ - bool sign_z = sign_x; - - // reveal the hidden bit in IEEE, adjust exponent and mantissa - ut32 man_len = rz_float_get_format_info(left->r, RZ_FLOAT_INFO_MAN_LEN); - - rz_bv_set(mx, man_len, true); - ex -= man_len; - rz_bv_set(my, man_len, true); - ey -= man_len; - - // every mantissa would become an big integer with clz(num) = 0 - ex -= rz_bv_clz(mx); - ey -= rz_bv_clz(my); - rz_bv_lshift(mx, rz_bv_clz(mx)); - rz_bv_lshift(my, rz_bv_clz(my)); - - // help flag - bool tiny = 0; - st32 compare = false; - bool quo_is_odd = false; - - // result of rem(x, y) - RzBitVector *mz = NULL; - ut32 ez; - RzFloat *z; - - // make last bit of mantissa is 1 - // TODO : add a scan function to bitvector lib (like clz but cnted from LSB to MSB) - ut32 k; - for (k = 0; k < my->len; ++k) { - if (rz_bv_get(my, k)) { - break; - } - } - - ey += k; - rz_bv_rshift(my, k); - - // q = x/y = mx/(my*2^(ey-ex)) - if (ex <= ey) { - // detect magnitude - ut32 sx = mx->len - rz_bv_clz(mx); - ut32 sy = my->len - rz_bv_clz(my); - ut32 mag_level_mx = sx + ex; - ut32 mag_level_my = sy + ey; - - if (mag_level_mx < mag_level_my) { - // tiny, quotient = 0, remainder = mx - tiny = 1; - z = rz_float_dup(left); - goto clean; - } else { - // mx mod my*2^(ey-ex) - // construct real number real_my = 2^(ey - ex) * my - RzBitVector *real_my = rz_bv_prepend_zero(my, my->len); - rz_bv_lshift(real_my, ey - ex); - - // stretch mx to have the same length for calculation - RzBitVector *stretched_mx = rz_bv_prepend_zero(mx, mx->len); - RzBitVector *stretched_mz = rz_bv_mod(stretched_mx, real_my); - mz = rz_bv_cut_head(stretched_mz, my->len); - - rz_bv_free(real_my); - rz_bv_free(stretched_mx); - rz_bv_free(stretched_mz); - } - } else { - // ex > ey - // preprocess for rounding - if (quo_rnd == RZ_FLOAT_RMODE_RTN) { - // let my = my * 2 - rz_bv_lshift(my, 1); - } - - // r = mx * (2^(ex - ey) mod my) mod my - // 1. build 2^(ex - ey) bv - ut32 aligned_length = ex - ey + 1; - - RzBitVector *two_exponent_fact; - RzBitVector *stretched_my; - bool is_stretched = false; - if (aligned_length < my->len) { - two_exponent_fact = rz_bv_new(my->len); - stretched_my = rz_bv_dup(my); - } else { - is_stretched = true; - two_exponent_fact = rz_bv_new(aligned_length); - stretched_my = rz_bv_prepend_zero(my, aligned_length - my->len); - } - rz_bv_set(two_exponent_fact, aligned_length - 1, true); - - // 2. mod my for the 1st time - RzBitVector *fact_mod = rz_bv_mod(two_exponent_fact, stretched_my); - - RzBitVector *mx_fact; - mx_fact = is_stretched ? rz_bv_cut_head(fact_mod, aligned_length - my->len) : rz_bv_dup(fact_mod); - - // 3. mul with mx, and then mod my - // mul maybe overflow, so stretch both - RzBitVector *mx_ext = rz_bv_prepend_zero(mx, mx->len); - RzBitVector *mx_fact_ext = rz_bv_prepend_zero(mx_fact, mx_fact->len); - RzBitVector *my_ext = rz_bv_prepend_zero(my, my->len); - RzBitVector *mul_ext = rz_bv_mul(mx_ext, mx_fact_ext); - RzBitVector *mz_ext; - mz_ext = rz_bv_mod(mul_ext, my_ext); - mz = rz_bv_cut_head(mz_ext, my->len); - - // free temp bv - rz_bv_free(two_exponent_fact); - rz_bv_free(stretched_my); - rz_bv_free(fact_mod); - rz_bv_free(mx_fact); - rz_bv_free(mx_ext); - rz_bv_free(mul_ext); - rz_bv_free(my_ext); - rz_bv_free(mz_ext); - rz_bv_free(mx_fact_ext); - - // rounding - if (quo_rnd == RZ_FLOAT_RMODE_RTN) { - // let my = my / 2 - rz_bv_shift_right_jammed(my, 1); - quo_is_odd = rz_bv_ule(my, mz); - if (quo_is_odd) { - // mz = mz - my - RzBitVector *tmp = rz_bv_sub(mz, my, NULL); - rz_bv_free(mz); - mz = tmp; - tmp = NULL; - } - } - } - - // r == 0, return 0 - if (rz_bv_is_zero_vector(mz)) { - z = rz_float_new_zero(left->r); - rz_bv_set(z->s, z->s->len, sign_z); - goto clean; - } - - // 2r < y ? round(r) : round(r-my) - if (quo_rnd == RZ_FLOAT_RMODE_RTN) { - // r = 2 * r - rz_bv_lshift(mz, 1); - - if (tiny) { - // detect magnitude - ut32 sz = mx->len - rz_bv_clz(mx); - ut32 sy = my->len - rz_bv_clz(my); - ut32 mag_level_mz = sz + ex; - ut32 mag_level_my = sy + ey; - - if (mag_level_mz > mag_level_my) { - // equal - compare = 0; - } else { - // sz >= ey + sr - ex, shift is safe - // my * 2^(ey - ex) - rz_bv_lshift(my, ey - ex); - compare = rz_bv_cmp(mz, my); - } - } else { - // cmp mz with my - compare = rz_bv_cmp(mz, my); - } - - rz_bv_shift_right_jammed(mz, 1); - if ((compare > 0) || - ((mode == RZ_FLOAT_RMODE_RTN) && (compare == 0) && (quo_is_odd))) { - // r = mz - my - RzBitVector *tmp = rz_bv_sub(mz, my, NULL); - rz_bv_free(mz); - mz = tmp; - tmp = NULL; - } - } - - // result exponent - ez = ex > ey ? ey : ex; - - // normalize - // make total - clz = man_len + 1, a normalized mz with hidden bit set - ut32 exp_len = rz_float_get_format_info(left->r, RZ_FLOAT_INFO_EXP_LEN); - st32 shift_dist = (st32)(rz_bv_clz(mz) - exp_len); - ez -= shift_dist; - if (shift_dist < 0) { - rz_bv_shift_right_jammed(mz, -shift_dist); - } else { - rz_bv_lshift(mz, shift_dist); - } - - // recover IEEE mantissa and exponent - ez += man_len; - ez = ez == 1 - bias ? 0 : ez + bias; - - // apply to round_float_bv required format - // 01 MMMM MMMM ... - shift_dist = (st32)(exp_len - 1); - rz_bv_lshift(mz, shift_dist); - - z = round_float_bv_new( - sign_z, - ez, - mz, - left->r, - left->r, - mode); -clean: - rz_bv_free(mx); - rz_bv_free(my); - rz_bv_free(mz); - return z; -} - -/** - * \brief calculate \p left % \p right and round the result after, return the result - * \details - * Any % 0 => NaN - * Inf % Any => NaN, invalid - * Any % Inf -> Any - * 0 % Any -> 0 - * \param mode rounding mode + * \param mode rounding mode used for calculating the quotient * \return result of arithmetic operation + * + * Can be positive or negative. Range: [ -abs(right)/2, abs(right)/2 ] */ RZ_API RZ_OWN RzFloat *rz_float_rem_ieee_bin(RZ_NONNULL RzFloat *left, RZ_NONNULL RzFloat *right, RzFloatRMode mode) { - return rz_float_rem_internal(left, right, RZ_FLOAT_RMODE_RNE, mode); + rz_return_val_if_fail(left && right && left->r == right->r, NULL); + + RzFloatFormat format = left->r; + set_float_rounding_mode(mode); + + switch (format) { + case RZ_FLOAT_IEEE754_BIN_32: + return of_float32(f32_rem(to_float32(left), to_float32(right))); + case RZ_FLOAT_IEEE754_BIN_64: + return of_float64(f64_rem(to_float64(left), to_float64(right))); + case RZ_FLOAT_IEEE754_BIN_80: + return of_float80(extF80_rem(to_float80(left), to_float80(right))); + case RZ_FLOAT_IEEE754_BIN_128: + return of_float128(f128_rem(to_float128(left), to_float128(right))); + default: + RZ_LOG_ERROR("float: REM operation unimplemented for format %d\n", format); + return NULL; + } } /** - * \brief calculate \p left % \p right and round the result after, return the result + * \brief Returns the value of \p left % \p right, with quotient rounded to an integer with rounding mode RTZ * \details * Any % 0 => NaN * Inf % Any => NaN, invalid * Any % Inf -> Any * 0 % Any -> 0 - * \param mode rounding mode + * \param mode rounding mode used for calculating the quotient * \return result of arithmetic operation + * + * Mod is guaranteed to be of the same sign as \p left. + * Range: + * - [ 0, abs(right) ) if left >= 0 + * - ( -abs(right), 0 ] if left <= 0 */ RZ_API RZ_OWN RzFloat *rz_float_mod_ieee_bin(RZ_NONNULL RzFloat *left, RZ_NONNULL RzFloat *right, RzFloatRMode mode) { - return rz_float_rem_internal(left, right, RZ_FLOAT_RMODE_RTZ, mode); + rz_return_val_if_fail(left && right && left->r == right->r, NULL); + + RzFloat *ret = rz_float_rem_ieee_bin(left, right, mode); + if (rz_float_get_sign(ret) != rz_float_get_sign(left)) { + if (rz_float_is_zero(ret)) { + /* If a zero is returned, it should still have the same sign as the dividend. */ + rz_float_set_sign(ret, rz_float_get_sign(left)); + } else { + RzFloat *same_sign = NULL; + RzFloat *right_abs = rz_float_abs(right); + + if (rz_float_is_negative(ret)) { + same_sign = rz_float_add(ret, right_abs, mode); + } else { + same_sign = rz_float_sub(ret, right_abs, mode); + } + + rz_float_free(ret); + ret = same_sign; + } + } + + return ret; } /** @@ -2007,244 +1349,34 @@ RZ_API RZ_OWN RzFloat *rz_float_mod_ieee_bin(RZ_NONNULL RzFloat *left, RZ_NONNUL * \return result of arithmetic operation */ RZ_API RZ_OWN RzFloat *rz_float_fma_ieee_bin(RZ_NONNULL RzFloat *a, RZ_NONNULL RzFloat *b, RZ_NONNULL RzFloat *c, RzFloatRMode mode) { - // process NaN / Inf - { - RzFloatSpec a_type, b_type, c_type; - a_type = rz_float_detect_spec(a); - b_type = rz_float_detect_spec(b); - c_type = rz_float_detect_spec(c); - bool a_is_inf = (a_type == RZ_FLOAT_SPEC_PINF || a_type == RZ_FLOAT_SPEC_NINF); - bool b_is_inf = (b_type == RZ_FLOAT_SPEC_PINF || b_type == RZ_FLOAT_SPEC_NINF); - bool c_is_inf = (c_type == RZ_FLOAT_SPEC_PINF || c_type == RZ_FLOAT_SPEC_NINF); - bool a_is_nan = (a_type == RZ_FLOAT_SPEC_SNAN || a_type == RZ_FLOAT_SPEC_QNAN); - bool b_is_nan = (b_type == RZ_FLOAT_SPEC_SNAN || b_type == RZ_FLOAT_SPEC_QNAN); - bool c_is_nan = (c_type == RZ_FLOAT_SPEC_SNAN || c_type == RZ_FLOAT_SPEC_QNAN); + rz_return_val_if_fail(a && b && c && a->r == b->r && b->r == c->r, NULL); - bool a_sign = get_sign(a->s, a->r); - bool b_sign = get_sign(b->s, b->r); - bool c_sign = get_sign(c->s, c->r); - - // simplified, may not be exactly correct - if (a_is_nan || b_is_nan || c_is_nan) { - return rz_float_new_qnan(a->r); - } - - if (a_is_inf || b_is_inf || c_is_inf) { - return rz_float_new_inf(a->r, a_is_inf ? a_sign : b_is_inf ? b_sign - : c_sign); - } - } - - // Extract attribute from format RzFloatFormat format = a->r; - ut32 exp_len = rz_float_get_format_info(format, RZ_FLOAT_INFO_EXP_LEN); - ut32 total_len = rz_float_get_format_info(format, RZ_FLOAT_INFO_TOTAL_LEN); - ut32 bias = rz_float_get_format_info(format, RZ_FLOAT_INFO_BIAS); - ut32 extra_len = total_len; + set_float_rounding_mode(mode); - // extra fields from a and b for multiply - RzBitVector *a_exp_squashed = get_exp_squashed(a->s, a->r); - RzBitVector *b_exp_squashed = get_exp_squashed(b->s, b->r); - RzBitVector *a_mantissa = get_man_stretched(a->s, a->r); - RzBitVector *b_mantissa = get_man_stretched(b->s, b->r); - RzBitVector *mul_sig = NULL; - bool a_sign = get_sign(a->s, a->r); - bool b_sign = get_sign(b->s, b->r); - bool mul_sign = a_sign ^ b_sign; - bool res_sign; - ut32 res_exp_val; - RzBitVector *res_sig; - RzFloat *ret_f; + switch (format) { + case RZ_FLOAT_IEEE754_BIN_32: + return of_float32(f32_mulAdd(to_float32(a), to_float32(b), to_float32(c))); + case RZ_FLOAT_IEEE754_BIN_64: + return of_float64(f64_mulAdd(to_float64(a), to_float64(b), to_float64(c))); + case RZ_FLOAT_IEEE754_BIN_80: { + /* We don't have a 80-bit FMA available in SoftFloat, so we cast the + * float to 128-bit, perform FMA and cast it back. This should be fine + * since th 80-bit and the 128-bit format differ only in the size of + * their mantissa. */ + float128_t a_resized = extF80_to_f128(to_float80(a)); + float128_t b_resized = extF80_to_f128(to_float80(b)); + float128_t c_resized = extF80_to_f128(to_float80(c)); - // Handle normal float multiply - ut32 a_exp_val = rz_bv_to_ut32(a_exp_squashed); - ut32 b_exp_val = rz_bv_to_ut32(b_exp_squashed); - ut32 shift_dist; - - // remember we would like to make 01.MM MMMM ... (but leave higher extra bits empty) - shift_dist = (exp_len + 1) - 2; - ut32 hidden_bit_pos = total_len - 2; - - rz_bv_lshift(a_mantissa, shift_dist); - rz_bv_lshift(b_mantissa, shift_dist); - - st32 aexp_nobias = rz_float_get_exponent_val_no_bias(a); - st32 bexp_nobias = rz_float_get_exponent_val_no_bias(b); - st32 mul_exp_val = aexp_nobias + bexp_nobias; - - // set leading bit - if (a_exp_val != 0) { - rz_bv_set(a_mantissa, hidden_bit_pos, true); + float128_t fma_resized = f128_mulAdd(a_resized, b_resized, c_resized); + return of_float80(f128_to_extF80(fma_resized)); } - - if (b_exp_val != 0) { - rz_bv_set(b_mantissa, hidden_bit_pos, true); + case RZ_FLOAT_IEEE754_BIN_128: + return of_float128(f128_mulAdd(to_float128(a), to_float128(b), to_float128(c))); + default: + RZ_LOG_ERROR("float: FMA operation unimplemented for format %d\n", format); + return NULL; } - - // multiplication - mul_sig = rz_bv_mul(a_mantissa, b_mantissa); - - // check if a carry happen, if not, l-shift to force a leading 1 - // check MSB and the bit after MSB - if (rz_bv_get(mul_sig, total_len + extra_len - 3)) { - // carry case, think about 01.10 * 01.10 => 0001.0010 - // 001X.MMMM... -> 001.0MMMMM.. - mul_exp_val += 1; - rz_bv_shift_right_jammed(mul_sig, 1); - } - - // check result and normalize it if needed - ut32 clz = rz_bv_clz(mul_sig); - if (clz > 3) { - // means there are sub normal as factor - // try shift - shift_dist = clz - 3; - if (mul_exp_val - (st32)shift_dist < 1 - bias) { - // too small, represent as sub-normal - shift_dist = mul_exp_val - (1 - bias); - } - rz_bv_lshift(mul_sig, shift_dist); - - // biased one - mul_exp_val = 0; - - // for those who may be sub-normal, use fake hidden bit for rounding - // note that result sig has 000H.MMMM... form - rz_bv_set(mul_sig, rz_bv_len(mul_sig) - 4, true); - } - // others has 0001.MMMM... - else { - mul_exp_val += bias; - } - - // note that mul sig has 000H.MMMM form - // addition we have 00H.MMMM form - rz_bv_lshift(mul_sig, 1); - - // calculating addition - RzBitVector *c_exp_squashed = get_exp_squashed(c->s, c->r); - ut32 c_exp_val = rz_bv_to_ut32(c_exp_squashed); - bool c_sign = get_sign(c->s, c->r); - RzBitVector *c_mantissa = get_man_stretched(c->s, c->r); - - res_sign = mul_sign; - if (!c_exp_val) { - if (rz_bv_is_zero_vector(c_mantissa)) { - res_exp_val = mul_exp_val - 1; - res_sig = mul_sig; - mul_sig = NULL; - goto round; - } - - // normalize sub-normal c - // TODO : create a function - normalize_subnorm - shift_dist = rz_bv_clz(c_mantissa) - (1 + exp_len) + 1; - res_exp_val = 1 - shift_dist; - rz_bv_lshift(c_mantissa, shift_dist); - } - - // prepare c_sig for addition - // set hidden bit 1 and shift to construct (00H.M MMMM ...) - hidden_bit_pos = total_len - 3; - rz_bv_lshift(c_mantissa, exp_len - 2); - rz_bv_set(c_mantissa, hidden_bit_pos, true); - rz_bv_lshift(c_mantissa, extra_len); - - st32 exp_diff_val = (st32)(mul_exp_val - c_exp_val); - st32 abs_exp_diff_val = exp_diff_val > 0 ? exp_diff_val : -exp_diff_val; - if (mul_sign == c_sign) { - // addition - if (exp_diff_val <= 0) { - res_exp_val = c_exp_val; - rz_bv_shift_right_jammed(mul_sig, abs_exp_diff_val); - } else { - res_exp_val = mul_exp_val; - rz_bv_shift_right_jammed(c_mantissa, abs_exp_diff_val); - } - - // calc - res_sig = rz_bv_add(mul_sig, c_mantissa, NULL); - - // check if we should normalize when carry - ut32 new_total_len = rz_bv_len(res_sig); - if (rz_bv_get(res_sig, new_total_len - 2)) { - res_exp_val += 1; - rz_bv_shift_right_jammed(res_sig, 1); - } - } else { - // sub - if (exp_diff_val < 0) { - res_sign = c_sign; - res_exp_val = c_exp_val; - rz_bv_shift_right_jammed(mul_sig, abs_exp_diff_val); - res_sig = rz_bv_sub(c_mantissa, mul_sig, NULL); - } else if (exp_diff_val == 0) { - res_exp_val = mul_exp_val; - res_sig = rz_bv_sub(mul_sig, c_mantissa, NULL); - if (rz_bv_is_zero_vector(res_sig)) { - goto zero; - } - if (rz_bv_msb(res_sig)) { - // if negative, turn to (+/- absolute val) from 2's complement - res_sign = !res_sign; - RzBitVector *tmp = rz_bv_complement_2(res_sig); - rz_bv_free(res_sig); - res_sig = tmp; - tmp = NULL; - } - - } else { - // exp_diff > 0 - res_exp_val = mul_exp_val; - rz_bv_shift_right_jammed(c_mantissa, abs_exp_diff_val); - res_sig = rz_bv_sub(mul_sig, c_mantissa, NULL); - } - - // note that we have 00H.MMMMM... form - shift_dist = rz_bv_clz(res_sig) - 2; - res_exp_val -= shift_dist; - if (shift_dist < 0) { - rz_bv_shift_right_jammed(res_sig, -shift_dist); - } else { - rz_bv_lshift(res_sig, shift_dist); - } - } - - // drop extra length - // recovered to original length - rz_bv_shift_right_jammed(res_sig, extra_len); - RzBitVector *tmp = rz_bv_cut_head(res_sig, extra_len); - rz_bv_free(res_sig); - res_sig = tmp; - tmp = NULL; - - goto round; - -zero: - // complete zero - ret_f = rz_float_new(format); - ret_f->s = rz_bv_new(total_len); - rz_bv_set(ret_f->s, total_len - 1, mode == RZ_FLOAT_RMODE_RTN); - goto clean; -round: - ret_f = round_float_bv_new( - res_sign, - res_exp_val, - res_sig, - format, - format, - mode); -clean: - rz_bv_free(a_mantissa); - rz_bv_free(a_exp_squashed); - rz_bv_free(b_mantissa); - rz_bv_free(b_exp_squashed); - rz_bv_free(mul_sig); - rz_bv_free(c_exp_squashed); - rz_bv_free(c_mantissa); - rz_bv_free(res_sig); - - return ret_f; } /** @@ -2253,47 +1385,24 @@ clean: * \return result of arithmetic operation */ RZ_API RZ_OWN RzFloat *rz_float_sqrt_ieee_bin(RZ_NONNULL RzFloat *n, RzFloatRMode mode) { - // Use Newton method now, May Optimize - RzFloat *eps = rz_float_new_zero(n->r); - ut32 bias = rz_float_get_format_info(n->r, RZ_FLOAT_INFO_BIAS); - ut32 man_len = rz_float_get_format_info(n->r, RZ_FLOAT_INFO_MAN_LEN); - ut32 eps_magic = bias - man_len; + rz_return_val_if_fail(n, NULL); - RzBitVector *eps_bv = rz_bv_new_from_ut64(n->s->len, eps_magic); - rz_bv_lshift(eps_bv, man_len); - RzFloat *x = rz_float_new_from_bv(eps_bv); - rz_bv_free(eps_bv); + RzFloatFormat format = n->r; + set_float_rounding_mode(mode); - while (true) { - RzFloat *q = rz_float_div_ieee_bin(n, x, mode); - RzFloat *sum = rz_float_add_ieee_bin(x, q, mode); - RzFloat *sum_half = rz_half_float(sum); - RzFloat *abs = rz_float_sub_ieee_bin(x, sum_half, mode); - rz_make_fabs(abs); - - // abs <= eps, both are positive - if (rz_bv_ule(abs->s, eps->s)) { - rz_float_free(q); - rz_float_free(abs); - rz_float_free(sum); - rz_float_free(sum_half); - break; - } - - rz_float_free(x); - x = sum_half; - sum_half = NULL; - - rz_float_free(q); - rz_float_free(abs); - rz_float_free(sum); - sum = NULL; - q = NULL; - abs = NULL; + switch (format) { + case RZ_FLOAT_IEEE754_BIN_32: + return of_float32(f32_sqrt(to_float32(n))); + case RZ_FLOAT_IEEE754_BIN_64: + return of_float64(f64_sqrt(to_float64(n))); + case RZ_FLOAT_IEEE754_BIN_80: + return of_float80(extF80_sqrt(to_float80(n))); + case RZ_FLOAT_IEEE754_BIN_128: + return of_float128(f128_sqrt(to_float128(n))); + default: + RZ_LOG_ERROR("float: SQRT operation unimplemented for format %d\n", format); + return NULL; } - - rz_float_free(eps); - return x; } /** \} */ // end rz_float_arithmetic_group @@ -2353,105 +1462,23 @@ RZ_API RZ_OWN RzFloat *rz_float_trunc(RZ_NONNULL RzFloat *f) { */ RZ_API RZ_OWN RzFloat *rz_float_round_to_integral(RZ_NONNULL RzFloat *f, RzFloatRMode mode) { rz_return_val_if_fail(f, NULL); - RzFloat *ret; - RzBitVector *tmp, *rounded; - ut32 exp = float_exponent(f); + RzFloatFormat format = f->r; - bool sign = get_sign(f->s, format); - ut32 bias = rz_float_get_format_info(format, RZ_FLOAT_INFO_BIAS); - bool is_subnormal = exp == 0; - st32 exp_no_bias = is_subnormal ? (1 - bias) : (exp - bias); - ut32 total_len = rz_float_get_format_info(format, RZ_FLOAT_INFO_TOTAL_LEN); - ut32 man_len = rz_float_get_format_info(format, RZ_FLOAT_INFO_MAN_LEN); + set_float_rounding_mode(mode); - // rounding float to get an integer means - // we should try to reserve `exponent` bits of mantissa - // drop extra bits or append zeros - // 1.MM..M * 2^exp = 1MM..M * 2^0 (integer) - bool should_inc = false; - RzBitVector *sig = rz_float_get_mantissa(f); - - // sub normal one has no hidden bit, others should set to 1 - if (!is_subnormal) { - rz_bv_set(sig, man_len, true); + switch (format) { + case RZ_FLOAT_IEEE754_BIN_32: + return of_float32(f32_roundToInt(to_float32(f), softfloat_roundingMode, false)); + case RZ_FLOAT_IEEE754_BIN_64: + return of_float64(f64_roundToInt(to_float64(f), softfloat_roundingMode, false)); + case RZ_FLOAT_IEEE754_BIN_80: + return of_float80(extF80_roundToInt(to_float80(f), softfloat_roundingMode, false)); + case RZ_FLOAT_IEEE754_BIN_128: + return of_float128(f128_roundToInt(to_float128(f), softfloat_roundingMode, false)); + default: + RZ_LOG_ERROR("float: ROUND operation unimplemented for format %d\n", format); + return NULL; } - - if (exp_no_bias >= 0) { - // has `exp_no_bias` + 3 + 1 length - tmp = round_significant(sign, sig, exp_no_bias, mode, &should_inc); - } else { - // float 1.M..M * 2^exp, when exp < 0 - // flatten it and we have 0.0..1M..M (|exp|+1 zeros before 1MMM...) - // set a fake 1 before radix point, and we can use round_significant to round - ut32 remained_zeros = total_len - man_len - 1; - RzBitVector *fake_f; - if (-exp_no_bias > remained_zeros) { - // prepend - fake_f = rz_bv_prepend_zero(sig, -exp_no_bias - remained_zeros); - } else { - fake_f = rz_bv_dup(sig); - } - rz_bv_set(fake_f, rz_bv_len(fake_f) - 1, true); - tmp = round_significant(sign, fake_f, 0, mode, &should_inc); - - // unset the fake 1 in tmp - // tmp has 3 + 1 + precision = 4 - rz_bv_set(tmp, 0, false); - rz_bv_free(fake_f); - } - rz_bv_free(sig); - sig = NULL; - - // rounded result, rounded has (3 + 1 + precision) length - if (should_inc) { - // WARN: possible overflow => no enough length - RzBitVector *bv_one; - bv_one = rz_bv_new_one(rz_bv_len(tmp)); - rounded = rz_bv_add(tmp, bv_one, NULL); - rz_bv_free(bv_one); - } else { - rounded = rz_bv_dup(tmp); - } - rz_bv_free(tmp); - tmp = NULL; - - // now we have an integer bitv, convert it to significant - // 0001 MMMM = 1.MMMM * 2^4 - st32 integral_exp_val = rz_bv_len(rounded) - rz_bv_clz(rounded) - 1; - if (integral_exp_val < 0) { - // -1, means rounded is all zero - ret = rz_float_new_zero(format); - rz_float_set_sign(ret, sign); - - rz_bv_free(rounded); - return ret; - } - RzBitVector *integeral_exp = rz_bv_new_from_ut64(32, integral_exp_val + bias); - - if (man_len > integral_exp_val) { - sig = rz_bv_append_zero(rounded, man_len - integral_exp_val); - rz_bv_free(rounded); - rounded = NULL; - } else { - // right shift zero bits - rz_bv_rshift(rounded, integral_exp_val - man_len); - sig = rounded; - rounded = NULL; - } - - ret = RZ_NEW0(RzFloat); - if (!ret) { - rz_bv_free(integeral_exp); - rz_bv_free(sig); - return ret; - } - - ret->r = format; - ret->s = pack_float_bv(sign, integeral_exp, sig, format); - - rz_bv_free(integeral_exp); - rz_bv_free(sig); - return ret; } /** @@ -2525,7 +1552,7 @@ RZ_API RZ_OWN RzBitVector *rz_float_cast_int(RZ_NONNULL RzFloat *f, ut32 length, * cast_sint s rm x returns an integer closest to x. * The resulting bitvector should be interpreted as a signed two-complement integer. * \param f float - * \param length length of returnded bitvector + * \param length length of returned bitvector * \param mode rounding mode * \return signed bitvector in 2's complement */ diff --git a/librz/util/float/float_internal.c b/librz/util/float/float_internal.c index f3ef83901a..76d9d611cc 100644 --- a/librz/util/float/float_internal.c +++ b/librz/util/float/float_internal.c @@ -246,26 +246,6 @@ static bool rz_make_fabs(RzFloat *f) { return rz_bv_set(f->s, f->s->len - 1, false); } -/** - * get the half value of a float (by decreasing exponent value) - * \param f float - * \return half value of a float - */ -static RzFloat *rz_half_float(RzFloat *f) { - ut32 total = rz_float_get_format_info(f->r, RZ_FLOAT_INFO_TOTAL_LEN); - ut32 exp_start = rz_float_get_format_info(f->r, RZ_FLOAT_INFO_MAN_LEN); - - // for exp sub 1 - RzBitVector *sub = rz_bv_new(total); - rz_bv_set(sub, exp_start, true); - RzFloat *half = rz_float_new(f->r); - rz_bv_free(half->s); - half->s = rz_bv_sub(f->s, sub, NULL); - - rz_bv_free(sub); - return half; -} - /** * Pack sign, exponent, and significant together to float bv * \param sign sign of float diff --git a/librz/util/meson.build b/librz/util/meson.build index 8e9ac6ba51..f022ed8f27 100644 --- a/librz/util/meson.build +++ b/librz/util/meson.build @@ -93,7 +93,7 @@ rz_util_common_sources = [ ] rz_util_sources = rz_util_common_sources -rz_util_deps = [ldl, lrt, mth, th, utl, pcre2_dep] + platform_deps +rz_util_deps = [ldl, lrt, mth, th, utl, pcre2_dep, softfloat_dep] + platform_deps if zlib_dep.found() rz_util_deps += [zlib_dep] endif diff --git a/meson.build b/meson.build index b6b1f0730a..6bf5a32423 100644 --- a/meson.build +++ b/meson.build @@ -608,6 +608,18 @@ if get_option('use_lzma') endif endif +# handle softfloat dependency +r = run_command(py3_exe, check_meson_subproject_py, 'softfloat', check: false) +if r.returncode() == 1 and get_option('subprojects_check') + error(subproject_clean_error_msg) +endif + +softfloat_dep = dependency('softfloat', required: get_option('use_sys_softfloat'), static: is_static_build) +if not softfloat_dep.found() + softfloat_proj = subproject('softfloat', default_options: ['default_library=static', 'warning_level=0']) + softfloat_dep = softfloat_proj.get_variable('softfloat_dep') +endif + # handle zip dependency r = run_command(py3_exe, check_meson_subproject_py, 'libzip', check: false) if r.returncode() == 1 and get_option('subprojects_check') @@ -798,6 +810,7 @@ summary({ 'System tree-sitter library': tree_sitter_dep.found() and tree_sitter_dep.type_name() != 'internal', 'System lz4 library': lz4_dep.found() and lz4_dep.type_name() != 'internal', 'System lzma library': liblzma_dep.found() and liblzma_dep.type_name() != 'internal', + 'System softfloat library': softfloat_dep.found() and softfloat_dep.type_name() != 'internal', 'System zlib library': zlib_dep.found() and zlib_dep.type_name() != 'internal', 'System zstd library': libzstd_dep.found() and libzstd_dep.type_name() != 'internal', 'System zip library': libzip_dep.found() and libzip_dep.type_name() != 'internal', diff --git a/meson_options.txt b/meson_options.txt index 5cba4c4ad0..b774ad26a4 100644 --- a/meson_options.txt +++ b/meson_options.txt @@ -36,6 +36,7 @@ option('use_sys_openssl', type: 'feature', value: 'disabled') option('use_sys_libmspack', type: 'feature', value: 'disabled') option('use_sys_tree_sitter', type: 'feature', value: 'disabled') option('use_sys_pcre2', type: 'feature', value: 'disabled') +option('use_sys_softfloat', type: 'feature', value: 'disabled') option('use_swift_demangler', type: 'boolean', value: true, description: 'If false, disables the swift demangler') option('use_gpl', type: 'boolean', value: true, description: 'Set to false when you want to disable gpl code') option('install_sigdb', type: 'boolean', value: false, description: 'Downloads and installs rizin sigdb') diff --git a/subprojects/softfloat.wrap b/subprojects/softfloat.wrap new file mode 100644 index 0000000000..a41cb5507b --- /dev/null +++ b/subprojects/softfloat.wrap @@ -0,0 +1,5 @@ +[wrap-git] +url = https://github.com/rizinorg/softfloat +revision = e06f4fcbdc6b3545c0624ac70725daf95eda53e8 +directory = softfloat +depth = 1 diff --git a/test/unit/test_float.c b/test/unit/test_float.c index 3f3c129591..f0bfb7efbe 100644 --- a/test/unit/test_float.c +++ b/test/unit/test_float.c @@ -558,9 +558,15 @@ bool f32_ieee_special_num_test(void) { } bool f32_ieee_rem_test(void) { + /* mod(x, y) = x - round(x/y, RNE) * y */ + + /* This test should return a different result in mod. */ RzFloat *a1 = rz_float_new_from_f32(4.0f); RzFloat *b1 = rz_float_new_from_f32(1.5f); - RzFloat *expect1 = rz_float_new_from_f32(1.0f); + RzFloat *expect1 = rz_float_new_from_f32(-0.5f); + /* quot = x/y = 4.0/1.5 = 2.666... + * rounded_quot (RNE) = 3 + */ RzFloat *rem1 = rz_float_rem(a1, b1, RZ_FLOAT_RMODE_RNE); mu_assert_true(is_equal_float(rem1, expect1), "rem test 1"); rz_float_free(a1); @@ -578,9 +584,13 @@ bool f32_ieee_rem_test(void) { rz_float_free(expect2); rz_float_free(rem2); - RzFloat *a3 = rz_float_new_from_ut32_as_f32(0x3F7FFF3F); - RzFloat *b3 = rz_float_new_from_ut32_as_f32(0x957CE0B6); - RzFloat *expect3 = rz_float_new_from_ut32_as_f32(0x145F53B0); + /* This test should return a different result in mod. */ + RzFloat *a3 = rz_float_new_from_ut32_as_f32(0xCBF83FFF); + RzFloat *b3 = rz_float_new_from_ut32_as_f32(0x44800FF0); + RzFloat *expect3 = rz_float_new_from_ut32_as_f32(0x43E63BC0); + /* quot = x/y = -32538622.0/1024.498046875 = -31760.550543997346 + * rounded_quot (RNE) = -31761 + */ RzFloat *rem3 = rz_float_rem(a3, b3, RZ_FLOAT_RMODE_RNE); mu_assert_true(is_equal_float(rem3, expect3), "rem test 3"); rz_float_free(a3); @@ -588,6 +598,70 @@ bool f32_ieee_rem_test(void) { rz_float_free(expect3); rz_float_free(rem3); + RzFloat *a4 = rz_float_new_from_ut32_as_f32(0x3F7FFF3F); + RzFloat *b4 = rz_float_new_from_ut32_as_f32(0x957CE0B6); + RzFloat *expect4 = rz_float_new_from_ut32_as_f32(0x145F53B0); + RzFloat *rem4 = rz_float_rem(a4, b4, RZ_FLOAT_RMODE_RNE); + mu_assert_true(is_equal_float(rem4, expect4), "rem test 4"); + rz_float_free(a4); + rz_float_free(b4); + rz_float_free(expect4); + rz_float_free(rem4); + + mu_end; +} + +bool f32_ieee_mod_test(void) { + /* mod(x, y) = x - round(x/y, RTZ) * y */ + + /* This test should return a different result in rem. */ + RzFloat *a1 = rz_float_new_from_f32(4.0f); + RzFloat *b1 = rz_float_new_from_f32(1.5f); + RzFloat *expect1 = rz_float_new_from_f32(1.0f); + /* quot = x/y = 4.0/1.5 = 2.666... + * rounded_quot (RTZ) = 2 + */ + RzFloat *rem1 = rz_float_mod(a1, b1, RZ_FLOAT_RMODE_RNE); + mu_assert_true(is_equal_float(rem1, expect1), "rem test 1"); + rz_float_free(a1); + rz_float_free(b1); + rz_float_free(expect1); + rz_float_free(rem1); + + RzFloat *a2 = rz_float_new_from_ut32_as_f32(0xCBF83FFF); + RzFloat *b2 = rz_float_new_from_ut32_as_f32(0x44801003); + RzFloat *expect2 = rz_float_new_from_ut32_as_f32(0xC3F52F40); + RzFloat *rem2 = rz_float_mod(a2, b2, RZ_FLOAT_RMODE_RNE); + mu_assert_true(is_equal_float(rem2, expect2), "rem test 2"); + rz_float_free(a2); + rz_float_free(b2); + rz_float_free(expect2); + rz_float_free(rem2); + + /* This test should return a different result in rem. */ + RzFloat *a3 = rz_float_new_from_ut32_as_f32(0xCBF83FFF); + RzFloat *b3 = rz_float_new_from_ut32_as_f32(0x44801002); + RzFloat *expect3 = rz_float_new_from_ut32_as_f32(0xC3F71F80); + /* quot = x/y = -32538622.0/1024.498046875 = -31760.550543997346 + * rounded_quot (RTZ) = -31760 + */ + RzFloat *rem3 = rz_float_mod(a3, b3, RZ_FLOAT_RMODE_RNE); + mu_assert_true(is_equal_float(rem3, expect3), "rem test 3"); + rz_float_free(a3); + rz_float_free(b3); + rz_float_free(expect3); + rz_float_free(rem3); + + RzFloat *a4 = rz_float_new_from_ut32_as_f32(0x3F7FFF3F); + RzFloat *b4 = rz_float_new_from_ut32_as_f32(0x957CE0B6); + RzFloat *expect4 = rz_float_new_from_ut32_as_f32(0x145F53B0); + RzFloat *rem4 = rz_float_mod(a4, b4, RZ_FLOAT_RMODE_RNE); + mu_assert_true(is_equal_float(rem4, expect4), "rem test 4"); + rz_float_free(a4); + rz_float_free(b4); + rz_float_free(expect4); + rz_float_free(rem4); + mu_end; } @@ -1414,10 +1488,91 @@ bool f32_ieee_cast_test(void) { mu_end; } -bool f80_round_test(void) { +static RzFloat *new_f80_from_bytes(const char *bytes) { + RzBitVector *bv = rz_bv_new_from_bytes_be((const unsigned char *)bytes, 0, 80); + RzFloat *ret = rz_float_new_from_bv(bv); + rz_bv_free(bv); + + return ret; +} + +bool f80_ieee_add_test(void) { + RzFloat *x_f80 = new_f80_from_bytes("\x40\x00\xa6\x5f\x8f\x48\x12\x44\xca\x0b"); + RzFloat *y_f80 = new_f80_from_bytes("\x3f\xfd\xc4\xf4\x67\x6e\x7d\x93\x40\x00"); + RzFloat *sum_f80 = rz_float_add(x_f80, y_f80, RZ_FLOAT_RMODE_RNE); + RzFloat *expected_f80 = new_f80_from_bytes("\x40\x00\xbe\xfe\x1c\x35\xe1\xf7\x32\x0b"); + mu_assert_false(rz_float_cmp(sum_f80, expected_f80), "Add 80-bit floats"); + + rz_float_free(x_f80); + rz_float_free(y_f80); + rz_float_free(sum_f80); + rz_float_free(expected_f80); + + mu_end; +} + +bool f80_ieee_sub_test(void) { + RzFloat *x_f80 = new_f80_from_bytes("\x40\x00\xa6\x5f\x8f\x48\x12\x44\xca\x0b"); + RzFloat *y_f80 = new_f80_from_bytes("\x3f\xfd\xc4\xf4\x67\x6e\x7d\x93\x40\x00"); + RzFloat *diff_f80 = rz_float_sub(x_f80, y_f80, RZ_FLOAT_RMODE_RNE); + RzFloat *expected_f80 = new_f80_from_bytes("\x40\x00\x8d\xc1\x02\x5a\x42\x92\x62\x0b"); + mu_assert_false(rz_float_cmp(diff_f80, expected_f80), "Subtract 80-bit floats"); + + rz_float_free(x_f80); + rz_float_free(y_f80); + rz_float_free(diff_f80); + rz_float_free(expected_f80); + + mu_end; +} + +bool f80_ieee_mul_test(void) { + RzFloat *x_f80 = new_f80_from_bytes("\x40\x00\xa6\x5f\x8f\x48\x12\x44\xca\x0b"); + RzFloat *y_f80 = new_f80_from_bytes("\x3f\xfd\xc4\xf4\x67\x6e\x7d\x93\x40\x00"); + RzFloat *prod_f80 = rz_float_mul(x_f80, y_f80, RZ_FLOAT_RMODE_RNE); + RzFloat *expected_f80 = new_f80_from_bytes("\x3f\xff\x80\x00\x00\x00\x00\x00\x00\x00"); + mu_assert_false(rz_float_cmp(prod_f80, expected_f80), "Multiply 80-bit floats"); + + rz_float_free(x_f80); + rz_float_free(y_f80); + rz_float_free(prod_f80); + rz_float_free(expected_f80); + + mu_end; +} + +bool f80_ieee_div_test(void) { + RzFloat *x_f80 = new_f80_from_bytes("\x3f\xff\x80\x00\x00\x00\x00\x00\x00\x00"); + RzFloat *y_f80 = new_f80_from_bytes("\xbf\xfd\xc4\xf4\x67\x6e\x7d\x93\x40\x00"); + RzFloat *quot_f80 = rz_float_div(x_f80, y_f80, RZ_FLOAT_RMODE_RNE); + RzFloat *expected_f80 = new_f80_from_bytes("\xc0\x00\xa6\x5f\x8f\x48\x12\x44\xca\x0b"); + mu_assert_false(rz_float_cmp(quot_f80, expected_f80), "Divide 80-bit floats"); + + rz_float_free(x_f80); + rz_float_free(y_f80); + rz_float_free(quot_f80); + rz_float_free(expected_f80); + + mu_end; +} + +bool f80_ieee_sqrt_test(void) { + RzFloat *x_f80 = new_f80_from_bytes("\x3f\xff\xab\x27\x32\x90\xa7\x8b\x0c\x29"); + RzFloat *sqrt_f80 = rz_float_sqrt(x_f80, RZ_FLOAT_RMODE_RNE); + RzFloat *expected_f80 = new_f80_from_bytes("\x3f\xff\x94\x03\x1c\xc0\x8d\xdc\xfb\xb5"); + mu_assert_false(rz_float_cmp(sqrt_f80, expected_f80), "Square root 80-bit floats"); + + rz_float_free(x_f80); + rz_float_free(sqrt_f80); + rz_float_free(expected_f80); + + mu_end; +} + +bool f80_ieee_cast_test(void) { /* To 80-bit */ RzFloat *old_f = rz_float_new_from_f64(14.285714285714286); - RzFloat *expect_f = rz_float_new_from_f80(14.2857142857142864756l); + RzFloat *expect_f = new_f80_from_bytes("\x40\x02\xe4\x92\x49\x24\x92\x49\x28\x00"); RzFloat *new_cast = rz_float_convert(old_f, RZ_FLOAT_IEEE754_BIN_80, RZ_FLOAT_RMODE_RNE); mu_assert_false(rz_float_cmp(expect_f, new_cast), "test convert 14.285714285714286d to 14.2857142857142864756l"); rz_float_free(old_f); @@ -1425,23 +1580,21 @@ bool f80_round_test(void) { rz_float_free(new_cast); /* From 80-bit */ - old_f = rz_float_new_from_f80(13.37l); + old_f = new_f80_from_bytes("\x40\x02\xd5\xeb\x85\x1e\xb8\x51\xeb\x85"); expect_f = rz_float_new_from_f32(13.37f); new_cast = rz_float_convert(old_f, RZ_FLOAT_IEEE754_BIN_32, RZ_FLOAT_RMODE_RNE); mu_assert_false(rz_float_cmp(expect_f, new_cast), "test convert 13.37l to 13.37f"); rz_float_free(old_f); rz_float_free(expect_f); rz_float_free(new_cast); - mu_end; /* From 80-bit to 80-bit (should lead to the same value) */ - old_f = rz_float_new_from_f80(66668466788774.6870804l); - expect_f = rz_float_new_from_f80(66668466788774.6870804l); + old_f = new_f80_from_bytes("\x40\x2c\xf2\x89\xd9\x1f\x66\x9a\xbf\x92"); new_cast = rz_float_convert(old_f, RZ_FLOAT_IEEE754_BIN_80, RZ_FLOAT_RMODE_RNE); - mu_assert_false(rz_float_cmp(expect_f, new_cast), "test convert 66668466788774.6870804l to itself"); + mu_assert_false(rz_float_cmp(old_f, new_cast), "test convert 66668466788774.6870804l to itself"); rz_float_free(old_f); - rz_float_free(expect_f); rz_float_free(new_cast); + mu_end; } @@ -1459,6 +1612,7 @@ bool all_tests() { mu_run_test(f32_ieee_round_test); mu_run_test(f32_ieee_sqrt_test); mu_run_test(f32_ieee_rem_test); + mu_run_test(f32_ieee_mod_test); mu_run_test(f32_ieee_special_num_test); mu_run_test(float_load_from_bitvector); mu_run_test(float_print_num); @@ -1468,7 +1622,12 @@ bool all_tests() { mu_run_test(f32_new_round_test); mu_run_test(f32_ieee_fround_test); mu_run_test(f32_ieee_cast_test); - mu_run_test(f80_round_test); + mu_run_test(f80_ieee_add_test); + mu_run_test(f80_ieee_sub_test); + mu_run_test(f80_ieee_mul_test); + mu_run_test(f80_ieee_div_test); + mu_run_test(f80_ieee_sqrt_test); + mu_run_test(f80_ieee_cast_test); return tests_passed != tests_run; }