diff options
Diffstat (limited to 'subprojects/exess/src/digits.c')
-rw-r--r-- | subprojects/exess/src/digits.c | 243 |
1 files changed, 243 insertions, 0 deletions
diff --git a/subprojects/exess/src/digits.c b/subprojects/exess/src/digits.c new file mode 100644 index 00000000..dd303269 --- /dev/null +++ b/subprojects/exess/src/digits.c @@ -0,0 +1,243 @@ +/* + Copyright 2019-2021 David Robillard <d@drobilla.net> + + Permission to use, copy, modify, and/or distribute this software for any + purpose with or without fee is hereby granted, provided that the above + copyright notice and this permission notice appear in all copies. + + THIS SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. +*/ + +#include "digits.h" + +#include "bigint.h" +#include "ieee_float.h" +#include "soft_float.h" +#include "warnings.h" + +#include <assert.h> +#include <math.h> +#include <stdbool.h> +#include <stddef.h> +#include <stdint.h> + +/* + This is more or less just an implementation of the classic rational number + based floating point print routine ("Dragon4"). See "How to Print + Floating-Point Numbers Accurately" by Guy L. Steele Jr. and Jon L White for + the canonical source. The basic idea is to find a big rational between 1 and + 10 where value = (numer / denom) * 10^e, then continuously divide it to + generate decimal digits. + + Unfortunately, this algorithm requires pretty massive bigints to work + correctly for all doubles, and isn't particularly fast. Something like + Grisu3 could be added to improve performance, but that has the annoying + property of needing a more precise fallback in some cases, meaning it would + only add more code, not replace any. Since this is already a pretty + ridiculous amount of code, I'll hold off on this until it becomes a problem, + or somebody comes up with a better algorithm. +*/ + +/// Return true if the number is within the lower boundary +static bool +within_lower(const ExessBigint* const numer, + const ExessBigint* const d_lower, + const bool is_even) +{ + return is_even ? exess_bigint_compare(numer, d_lower) <= 0 + : exess_bigint_compare(numer, d_lower) < 0; +} + +/// Return true if the number is within the upper boundary +static bool +within_upper(const ExessBigint* const numer, + const ExessBigint* const denom, + const ExessBigint* const d_upper, + const bool is_even) +{ + return is_even ? exess_bigint_plus_compare(numer, d_upper, denom) >= 0 + : exess_bigint_plus_compare(numer, d_upper, denom) > 0; +} + +/** + Find values so that 0.1 <= numer/denom < 1 or 1 <= numer/denom < 10. + + @param significand Double significand. + @param exponent Double exponent (base 2). + @param decimal_power Decimal exponent (log10 of the double). + @param[out] numer Numerator of rational number. + @param[out] denom Denominator of rational number. + @param[out] delta Distance to the lower and upper boundaries. +*/ +static void +calculate_initial_values(const uint64_t significand, + const int exponent, + const int decimal_power, + const bool lower_is_closer, + ExessBigint* const numer, + ExessBigint* const denom, + ExessBigint* const delta) +{ + /* Use a common denominator of 2^1 so that boundary distance is an integer. + If the lower boundary is closer, we need to scale everything but the + lower boundary to compensate, so add another factor of two here (this is + faster than shifting them again later as in the paper). */ + const unsigned lg_denom = 1u + lower_is_closer; + + if (exponent >= 0) { + // delta = 2^e + exess_bigint_set_u32(delta, 1); + exess_bigint_shift_left(delta, (unsigned)exponent); + + // numer = f * 2^e + exess_bigint_set_u64(numer, significand); + exess_bigint_shift_left(numer, (unsigned)exponent + lg_denom); + + // denom = 10^d + exess_bigint_set_pow10(denom, (unsigned)decimal_power); + exess_bigint_shift_left(denom, lg_denom); + } else if (decimal_power >= 0) { + // delta = 2^e, which is just 1 here since 2^-e is in the denominator + exess_bigint_set_u32(delta, 1); + + // numer = f + exess_bigint_set_u64(numer, significand); + exess_bigint_shift_left(numer, lg_denom); + + // denom = 10^d * 2^-e + exess_bigint_set_pow10(denom, (unsigned)decimal_power); + exess_bigint_shift_left(denom, (unsigned)-exponent + lg_denom); + } else { + // delta = 10^d + exess_bigint_set_pow10(delta, (unsigned)-decimal_power); + + // numer = f * 10^-d + exess_bigint_set(numer, delta); + exess_bigint_multiply_u64(numer, significand); + exess_bigint_shift_left(numer, lg_denom); + + // denom = 2^-exponent + exess_bigint_set_u32(denom, 1); + exess_bigint_shift_left(denom, (unsigned)-exponent + lg_denom); + } +} + +#ifndef NDEBUG +static bool +check_initial_values(const ExessBigint* const numer, + const ExessBigint* const denom, + const ExessBigint* const d_upper) +{ + ExessBigint upper = *numer; + exess_bigint_add(&upper, d_upper); + assert(exess_bigint_compare(&upper, denom) >= 0); + + const uint32_t div = exess_bigint_divmod(&upper, denom); + assert(div >= 1 && div < 10); + return true; +} +#endif + +static unsigned +emit_digits(ExessBigint* const numer, + const ExessBigint* const denom, + ExessBigint* const d_lower, + ExessBigint* const d_upper, + const bool is_even, + char* const buffer, + const size_t max_digits) +{ + unsigned length = 0; + for (size_t i = 0; i < max_digits; ++i) { + // Emit the next digit + const uint32_t digit = exess_bigint_divmod(numer, denom); + assert(digit <= 9); + buffer[length++] = (char)('0' + digit); + + // Check for termination + const bool within_low = within_lower(numer, d_lower, is_even); + const bool within_high = within_upper(numer, denom, d_upper, is_even); + if (!within_low && !within_high) { + exess_bigint_multiply_u32(numer, 10); + exess_bigint_multiply_u32(d_lower, 10); + if (d_lower != d_upper) { + exess_bigint_multiply_u32(d_upper, 10); + } + } else { + if (!within_low || (within_high && exess_bigint_plus_compare( + numer, numer, denom) >= 0)) { + // In high only, or halfway and the next digit is > 5, round up + assert(buffer[length - 1] != '9'); + buffer[length - 1]++; + } + + break; + } + } + + return length; +} + +ExessDigitCount +exess_digits(const double d, char* const buf, const unsigned max_digits) +{ + EXESS_DISABLE_CONVERSION_WARNINGS + assert(isfinite(d) && fpclassify(d) != FP_ZERO); + EXESS_RESTORE_WARNINGS + + const ExessSoftFloat value = soft_float_from_double(d); + const int power = (int)lrint(log10(d)); + const bool is_even = !(value.f & 1u); + const bool lower_is_closer = double_lower_boundary_is_closer(d); + + // Calculate initial values so that v = (numer / denom) * 10^power + ExessBigint numer; + ExessBigint denom; + ExessBigint d_lower; + calculate_initial_values( + value.f, value.e, power, lower_is_closer, &numer, &denom, &d_lower); + + ExessBigint d_upper_storage; + ExessBigint* d_upper = NULL; + if (lower_is_closer) { + // Scale upper boundary to account for the closer lower boundary + // (the numerator and denominator were already scaled above) + d_upper_storage = d_lower; + d_upper = &d_upper_storage; + exess_bigint_shift_left(d_upper, 1); + } else { + d_upper = &d_lower; // Boundaries are the same, reuse the lower + } + + // Scale if necessary to make 1 <= (numer + delta) / denom < 10 + ExessDigitCount count = {0, 0}; + if (within_upper(&numer, &denom, d_upper, is_even)) { + count.expt = power; + } else { + count.expt = power - 1; + exess_bigint_multiply_u32(&numer, 10); + exess_bigint_multiply_u32(&d_lower, 10); + if (d_upper != &d_lower) { + exess_bigint_multiply_u32(d_upper, 10); + } + } + + // Write digits to output + assert(check_initial_values(&numer, &denom, d_upper)); + count.count = + emit_digits(&numer, &denom, &d_lower, d_upper, is_even, buf, max_digits); + + // Trim trailing zeros + while (count.count > 1 && buf[count.count - 1] == '0') { + buf[--count.count] = 0; + } + + buf[count.count] = '\0'; + return count; +} |