From 497b205caab21c2c42ebbf9851c7174133913e6c Mon Sep 17 00:00:00 2001 From: Eyal Rozenberg Date: Sat, 5 Feb 2022 00:00:49 +0200 Subject: [PATCH] * Factored out the code computing a power-of-10 (for the floored base-10 log of the number to print). * Increased the digits in the literals for the power-of-10 computation. * Beefed up the comments about how the computation is performed, including specific comments with the closed form of each magin number. --- src/printf/printf.c | 55 ++++++++++++++++++++++++++++++++------------- 1 file changed, 39 insertions(+), 16 deletions(-) diff --git a/src/printf/printf.c b/src/printf/printf.c index b8db9c82..a800234e 100644 --- a/src/printf/printf.c +++ b/src/printf/printf.c @@ -776,17 +776,47 @@ static int bastardized_floor(double x) return ( ((double) n) == x ) ? n : n-1; } -// Compute the base-10 logarithm of a positive number -static double log10_of_positive(double abs_number) +// Computes the base-10 logarithm of the input number - which must be an actual +// positive number (not infinity or NaN, nor a sub-normal) +static double log10_of_positive(double positive_number) { - double_with_bit_access dwba = get_bit_access(abs_number); + // The implementation follows David Gay (https://www.ampl.com/netlib/fp/dtoa.c). + // + // Since log_10 ( M * 2^x ) = log_10(M) + x , we can separate the components of + // our input number, and need only solve log_10(M) for M between 1 and 2 (as + // the base-2 mantissa is always 1-point-something). In that limited range, a + // Taylor series expansion of log10(x) should serve us well enough; and we'll + // take the mid-point, 1.5, as the point of expansion. + + double_with_bit_access dwba = get_bit_access(positive_number); // based on the algorithm by David Gay (https://www.ampl.com/netlib/fp/dtoa.c) int exp2 = get_exp2(dwba); // drop the exponent, so dwba.F comes into the range [1,2) dwba.U = (dwba.U & (((double_uint_t) (1) << DOUBLE_STORED_MANTISSA_BITS) - 1U)) | ((double_uint_t) DOUBLE_BASE_EXPONENT << DOUBLE_STORED_MANTISSA_BITS); - // now approximate log10 from the log2 integer part and an expansion of ln around 1.5 - return (0.1760912590558 + exp2 * 0.301029995663981 + (dwba.F - 1.5) * 0.289529654602168); + return ( + // Taylor expansion around 1.5: + 0.176091259055681242 // Expansion term 0: log_10(1.5) = ln(1.5)/ln(10) + + (dwba.F - 1.5) * 0.2895296546021678851 // Expansion term 1: (M - 1.5) * 1/(1.5 ln(10)) = (M/1.5 - 1) / ln(10) + // exact log_2 of the exponent x, with logarithm base change + + exp2 * 0.30102999566398119521 // = exp2 * log_10(2) = exp2 * ln(2)/ln(10) + ); + +} + + +static double pow10_of_int(int floored_exp10) +{ + // Compute 10^(floored_exp10) but (try to) make sure that doesn't overflow + double_with_bit_access dwba; + int exp2 = bastardized_floor(floored_exp10 * 3.321928094887362 + 0.5); + const double z = floored_exp10 * 2.302585092994046 - exp2 * 0.6931471805599453; + const double z2 = z * z; + dwba.U = ((double_uint_t)(exp2) + DOUBLE_BASE_EXPONENT) << DOUBLE_STORED_MANTISSA_BITS; + // compute exp(z) using continued fractions, + // see https://en.wikipedia.org/wiki/Exponential_function#Continued_fractions_for_ex + dwba.F *= 1 + 2 * z / (2 - z + (z2 / (6 + (z2 / (10 + z2 / 14))))); + return dwba.F; } // internal ftoa variant for exponential floating-point type, contributed by Martijn Jasperse @@ -809,21 +839,14 @@ static void print_exponential_number(output_gadget_t* output, double number, pri else { double exp10 = log10_of_positive(abs_number); floored_exp10 = bastardized_floor(exp10); - // now we want to compute 10^(floored_exp10) but we want to be sure it won't overflow - int exp2 = bastardized_floor(floored_exp10 * 3.321928094887362 + 0.5); - const double z = floored_exp10 * 2.302585092994046 - exp2 * 0.6931471805599453; - const double z2 = z * z; - double_with_bit_access dwba; - dwba.U = ((double_uint_t)(exp2) + DOUBLE_BASE_EXPONENT) << DOUBLE_STORED_MANTISSA_BITS; - // compute exp(z) using continued fractions, see https://en.wikipedia.org/wiki/Exponential_function#Continued_fractions_for_ex - dwba.F *= 1 + 2 * z / (2 - z + (z2 / (6 + (z2 / (10 + z2 / 14))))); + double p10 = pow10_of_int(floored_exp10); // correct for rounding errors - if (abs_number < dwba.F) { + if (abs_number < p10) { floored_exp10--; - dwba.F /= 10; + p10 /= 10; } abs_exp10_covered_by_powers_table = PRINTF_ABS(floored_exp10) < PRINTF_MAX_PRECOMPUTED_POWER_OF_10; - normalization.raw_factor = abs_exp10_covered_by_powers_table ? powers_of_10[PRINTF_ABS(floored_exp10)] : dwba.F; + normalization.raw_factor = abs_exp10_covered_by_powers_table ? powers_of_10[PRINTF_ABS(floored_exp10)] : p10; } // We now begin accounting for the widths of the two parts of our printed field: