Fix etc__parse_number_f64 goto-crosses-init
diff --git a/internal/cgen/base/strconv-impl.c b/internal/cgen/base/strconv-impl.c
index 1d1c3af..06710a1 100644
--- a/internal/cgen/base/strconv-impl.c
+++ b/internal/cgen/base/strconv-impl.c
@@ -1107,178 +1107,184 @@
wuffs_base__private_implementation__medium_prec_bin* m,
const wuffs_base__private_implementation__high_prec_dec* h,
bool skip_fast_path_for_tests) {
- // m->mantissa is a uint64_t, which is an integer approximation to a rational
- // value - h's underlying digits after m's normalization. This error is an
- // upper bound on the difference between the approximate and actual value.
- //
- // The DiyFpStrtod function in https://github.com/google/double-conversion
- // uses a finer grain (1/8th of the ULP, Unit in the Last Place) when
- // tracking error. This implementation is coarser (1 ULP) but simpler.
- //
- // It is an error in the "numerical approximation" sense, not in the typical
- // programming sense (as in "bad input" or "a result type").
- uint64_t error = 0;
-
- // Convert up to 19 decimal digits (in h->digits) to 64 binary digits (in
- // m->mantissa): (1e19 < (1<<64)) and ((1<<64) < 1e20). If we have more than
- // 19 digits, we're truncating (with error).
- uint32_t i;
- uint32_t i_end = h->num_digits;
- if (i_end > 19) {
- i_end = 19;
- error = 1;
- }
- uint64_t mantissa = 0;
- for (i = 0; i < i_end; i++) {
- mantissa = (10 * mantissa) + h->digits[i];
- }
- m->mantissa = mantissa;
- m->exp2 = 0;
-
- // Check that exp10 lies in the (big_powers_of_10 + small_powers_of_10)
- // range, -348 ..= +347, stepping big_powers_of_10 by 8 (which is 87 triples)
- // and small_powers_of_10 by 1 (which is 8 triples).
- int32_t exp10 = h->decimal_point - ((int32_t)(i_end));
- if (exp10 < -348) {
- goto fail;
- }
- uint32_t bpo10 = ((uint32_t)(exp10 + 348)) / 8;
- uint32_t spo10 = ((uint32_t)(exp10 + 348)) % 8;
- if (bpo10 >= 87) {
- goto fail;
- }
-
- // Try a fast path, if float64 math would be exact.
- //
- // 15 is such that 1e15 can be losslessly represented in a float64 mantissa:
- // (1e15 < (1<<53)) and ((1<<53) < 1e16).
- //
- // 22 is the maximum valid index for the
- // wuffs_base__private_implementation__f64_powers_of_10 array.
do {
- if (skip_fast_path_for_tests || ((mantissa >> 52) != 0)) {
- break;
+ // m->mantissa is a uint64_t, which is an integer approximation to a
+ // rational value - h's underlying digits after m's normalization. This
+ // error is an upper bound on the difference between the approximate and
+ // actual value.
+ //
+ // The DiyFpStrtod function in https://github.com/google/double-conversion
+ // uses a finer grain (1/8th of the ULP, Unit in the Last Place) when
+ // tracking error. This implementation is coarser (1 ULP) but simpler.
+ //
+ // It is an error in the "numerical approximation" sense, not in the
+ // typical programming sense (as in "bad input" or "a result type").
+ uint64_t error = 0;
+
+ // Convert up to 19 decimal digits (in h->digits) to 64 binary digits (in
+ // m->mantissa): (1e19 < (1<<64)) and ((1<<64) < 1e20). If we have more
+ // than 19 digits, we're truncating (with error).
+ uint32_t i;
+ uint32_t i_end = h->num_digits;
+ if (i_end > 19) {
+ i_end = 19;
+ error = 1;
}
- double d = (double)mantissa;
+ uint64_t mantissa = 0;
+ for (i = 0; i < i_end; i++) {
+ mantissa = (10 * mantissa) + h->digits[i];
+ }
+ m->mantissa = mantissa;
+ m->exp2 = 0;
- if (exp10 == 0) {
- wuffs_base__result_f64 ret;
- ret.status.repr = NULL;
- ret.value = h->negative ? -d : +d;
- return ret;
+ // Check that exp10 lies in the (big_powers_of_10 + small_powers_of_10)
+ // range, -348 ..= +347, stepping big_powers_of_10 by 8 (which is 87
+ // triples) and small_powers_of_10 by 1 (which is 8 triples).
+ int32_t exp10 = h->decimal_point - ((int32_t)(i_end));
+ if (exp10 < -348) {
+ goto fail;
+ }
+ uint32_t bpo10 = ((uint32_t)(exp10 + 348)) / 8;
+ uint32_t spo10 = ((uint32_t)(exp10 + 348)) % 8;
+ if (bpo10 >= 87) {
+ goto fail;
+ }
- } else if (exp10 > 0) {
- if (exp10 > 22) {
- if (exp10 > (15 + 22)) {
- break;
- }
- // If exp10 is in the range 23 ..= 37, try moving a few of the zeroes
- // from the exponent to the mantissa. If we're still under 1e15, we
- // haven't truncated any mantissa bits.
- if (exp10 > 22) {
- d *= wuffs_base__private_implementation__f64_powers_of_10[exp10 - 22];
- exp10 = 22;
- if (d >= 1e15) {
- break;
- }
- }
- }
- d *= wuffs_base__private_implementation__f64_powers_of_10[exp10];
- wuffs_base__result_f64 ret;
- ret.status.repr = NULL;
- ret.value = h->negative ? -d : +d;
- return ret;
-
- } else { // "if (exp10 < 0)" is effectively "if (true)" here.
- if (exp10 < -22) {
+ // Try a fast path, if float64 math would be exact.
+ //
+ // 15 is such that 1e15 can be losslessly represented in a float64
+ // mantissa: (1e15 < (1<<53)) and ((1<<53) < 1e16).
+ //
+ // 22 is the maximum valid index for the
+ // wuffs_base__private_implementation__f64_powers_of_10 array.
+ do {
+ if (skip_fast_path_for_tests || ((mantissa >> 52) != 0)) {
break;
}
- d /= wuffs_base__private_implementation__f64_powers_of_10[-exp10];
- wuffs_base__result_f64 ret;
- ret.status.repr = NULL;
- ret.value = h->negative ? -d : +d;
- return ret;
+ double d = (double)mantissa;
+
+ if (exp10 == 0) {
+ wuffs_base__result_f64 ret;
+ ret.status.repr = NULL;
+ ret.value = h->negative ? -d : +d;
+ return ret;
+
+ } else if (exp10 > 0) {
+ if (exp10 > 22) {
+ if (exp10 > (15 + 22)) {
+ break;
+ }
+ // If exp10 is in the range 23 ..= 37, try moving a few of the zeroes
+ // from the exponent to the mantissa. If we're still under 1e15, we
+ // haven't truncated any mantissa bits.
+ if (exp10 > 22) {
+ d *= wuffs_base__private_implementation__f64_powers_of_10[exp10 -
+ 22];
+ exp10 = 22;
+ if (d >= 1e15) {
+ break;
+ }
+ }
+ }
+ d *= wuffs_base__private_implementation__f64_powers_of_10[exp10];
+ wuffs_base__result_f64 ret;
+ ret.status.repr = NULL;
+ ret.value = h->negative ? -d : +d;
+ return ret;
+
+ } else { // "if (exp10 < 0)" is effectively "if (true)" here.
+ if (exp10 < -22) {
+ break;
+ }
+ d /= wuffs_base__private_implementation__f64_powers_of_10[-exp10];
+ wuffs_base__result_f64 ret;
+ ret.status.repr = NULL;
+ ret.value = h->negative ? -d : +d;
+ return ret;
+ }
+ } while (0);
+
+ // Normalize (and scale the error).
+ error <<= wuffs_base__private_implementation__medium_prec_bin__normalize(m);
+
+ // Multiplying two MPB values nominally multiplies two mantissas, call them
+ // A and B, which are integer approximations to the precise values (A+a)
+ // and (B+b) for some error terms a and b.
+ //
+ // MPB multiplication calculates (((A+a) * (B+b)) >> 64) to be ((A*B) >>
+ // 64). Shifting (truncating) and rounding introduces further error. The
+ // difference between the calculated result:
+ // ((A*B ) >> 64)
+ // and the true result:
+ // ((A*B + A*b + a*B + a*b) >> 64) + rounding_error
+ // is:
+ // (( A*b + a*B + a*b) >> 64) + rounding_error
+ // which can be re-grouped as:
+ // ((A*b) >> 64) + ((a*(B+b)) >> 64) + rounding_error
+ //
+ // Now, let A and a be "m->mantissa" and "error", and B and b be the
+ // pre-calculated power of 10. A and B are both less than (1 << 64), a is
+ // the "error" local variable and b is less than 1.
+ //
+ // An upper bound (in absolute value) on ((A*b) >> 64) is therefore 1.
+ //
+ // An upper bound on ((a*(B+b)) >> 64) is a, also known as error.
+ //
+ // Finally, the rounding_error is at most 1.
+ //
+ // In total, calling mpb__mul_pow_10 will raise the worst-case error by 2.
+ // The subsequent re-normalization can multiply that by a further factor.
+
+ // Multiply by small_powers_of_10[etc].
+ wuffs_base__private_implementation__medium_prec_bin__mul_pow_10(
+ m, &wuffs_base__private_implementation__small_powers_of_10[3 * spo10]);
+ error += 2;
+ error <<= wuffs_base__private_implementation__medium_prec_bin__normalize(m);
+
+ // Multiply by big_powers_of_10[etc].
+ wuffs_base__private_implementation__medium_prec_bin__mul_pow_10(
+ m, &wuffs_base__private_implementation__big_powers_of_10[3 * bpo10]);
+ error += 2;
+ error <<= wuffs_base__private_implementation__medium_prec_bin__normalize(m);
+
+ // We have a good approximation of h, but we still have to check whether
+ // the error is small enough. Equivalently, whether the number of surplus
+ // mantissa bits (the bits dropped when going from m's 64 mantissa bits to
+ // the smaller number of double-precision mantissa bits) would always round
+ // up or down, even when perturbed by ±error. We start at 11 surplus bits
+ // (m has 64, double-precision has 1+52), but it can be higher for
+ // subnormals.
+ //
+ // In many cases, the error is small enough and we return true.
+ const int32_t f64_bias = -1023;
+ int32_t subnormal_exp2 = f64_bias - 63;
+ uint32_t surplus_bits = 11;
+ if (subnormal_exp2 >= m->exp2) {
+ surplus_bits += 1 + ((uint32_t)(subnormal_exp2 - m->exp2));
}
+
+ uint64_t surplus_mask =
+ (((uint64_t)1) << surplus_bits) - 1; // e.g. 0x07FF.
+ uint64_t surplus = m->mantissa & surplus_mask;
+ uint64_t halfway = ((uint64_t)1) << (surplus_bits - 1); // e.g. 0x0400.
+
+ // Do the final calculation in *signed* arithmetic.
+ int64_t i_surplus = (int64_t)surplus;
+ int64_t i_halfway = (int64_t)halfway;
+ int64_t i_error = (int64_t)error;
+
+ if ((i_surplus > (i_halfway - i_error)) &&
+ (i_surplus < (i_halfway + i_error))) {
+ goto fail;
+ }
+
+ wuffs_base__result_f64 ret;
+ ret.status.repr = NULL;
+ ret.value = wuffs_base__private_implementation__medium_prec_bin__as_f64(
+ m, h->negative);
+ return ret;
} while (0);
- // Normalize (and scale the error).
- error <<= wuffs_base__private_implementation__medium_prec_bin__normalize(m);
-
- // Multiplying two MPB values nominally multiplies two mantissas, call them A
- // and B, which are integer approximations to the precise values (A+a) and
- // (B+b) for some error terms a and b.
- //
- // MPB multiplication calculates (((A+a) * (B+b)) >> 64) to be ((A*B) >> 64).
- // Shifting (truncating) and rounding introduces further error. The
- // difference between the calculated result:
- // ((A*B ) >> 64)
- // and the true result:
- // ((A*B + A*b + a*B + a*b) >> 64) + rounding_error
- // is:
- // (( A*b + a*B + a*b) >> 64) + rounding_error
- // which can be re-grouped as:
- // ((A*b) >> 64) + ((a*(B+b)) >> 64) + rounding_error
- //
- // Now, let A and a be "m->mantissa" and "error", and B and b be the
- // pre-calculated power of 10. A and B are both less than (1 << 64), a is the
- // "error" local variable and b is less than 1.
- //
- // An upper bound (in absolute value) on ((A*b) >> 64) is therefore 1.
- //
- // Similarly, an upper bound on ((a*(B+b)) >> 64) is a, also known as error.
- //
- // Finally, the rounding_error is at most 1.
- //
- // In total, calling mpb__mul_pow_10 will increase the worst-case error by 2.
- // The subsequent re-normalization can multiply that by a further factor.
-
- // Multiply by small_powers_of_10[etc].
- wuffs_base__private_implementation__medium_prec_bin__mul_pow_10(
- m, &wuffs_base__private_implementation__small_powers_of_10[3 * spo10]);
- error += 2;
- error <<= wuffs_base__private_implementation__medium_prec_bin__normalize(m);
-
- // Multiply by big_powers_of_10[etc].
- wuffs_base__private_implementation__medium_prec_bin__mul_pow_10(
- m, &wuffs_base__private_implementation__big_powers_of_10[3 * bpo10]);
- error += 2;
- error <<= wuffs_base__private_implementation__medium_prec_bin__normalize(m);
-
- // We have a good approximation of h, but we still have to check whether the
- // error is small enough. Equivalently, whether the number of surplus
- // mantissa bits (the bits dropped when going from m's 64 mantissa bits to
- // the smaller number of double-precision mantissa bits) would always round
- // up or down, even when perturbed by ±error. We start at 11 surplus bits (m
- // has 64, double-precision has 1+52), but it can be higher for subnormals.
- //
- // In many cases, the error is small enough and we return true.
- const int32_t f64_bias = -1023;
- int32_t subnormal_exp2 = f64_bias - 63;
- uint32_t surplus_bits = 11;
- if (subnormal_exp2 >= m->exp2) {
- surplus_bits += 1 + ((uint32_t)(subnormal_exp2 - m->exp2));
- }
-
- uint64_t surplus_mask = (((uint64_t)1) << surplus_bits) - 1; // e.g. 0x07FF.
- uint64_t surplus = m->mantissa & surplus_mask;
- uint64_t halfway = ((uint64_t)1) << (surplus_bits - 1); // e.g. 0x0400.
-
- // Do the final calculation in *signed* arithmetic.
- int64_t i_surplus = (int64_t)surplus;
- int64_t i_halfway = (int64_t)halfway;
- int64_t i_error = (int64_t)error;
-
- if ((i_surplus > (i_halfway - i_error)) &&
- (i_surplus < (i_halfway + i_error))) {
- goto fail;
- }
-
- wuffs_base__result_f64 ret;
- ret.status.repr = NULL;
- ret.value = wuffs_base__private_implementation__medium_prec_bin__as_f64(
- m, h->negative);
- return ret;
-
fail:
do {
wuffs_base__result_f64 ret;
diff --git a/internal/cgen/data.go b/internal/cgen/data.go
index da600e3..437a056 100644
--- a/internal/cgen/data.go
+++ b/internal/cgen/data.go
@@ -118,14 +118,14 @@
"normalized).\n//\n// The etc_powers_of_10 triple is already normalized.\nstatic void //\nwuffs_base__private_implementation__medium_prec_bin__mul_pow_10(\n wuffs_base__private_implementation__medium_prec_bin* m,\n const uint32_t* p) {\n uint64_t p_mantissa = ((uint64_t)p[0]) | (((uint64_t)p[1]) << 32);\n int32_t p_exp2 = (int32_t)p[2];\n\n wuffs_base__multiply_u64__output o =\n wuffs_base__multiply_u64(m->mantissa, p_mantissa);\n // Round the mantissa up. It cannot overflow because the maximum possible\n // value of o.hi is 0xFFFFFFFFFFFFFFFE.\n m->mantissa = o.hi + (o.lo >> 63);\n m->exp2 = m->exp2 + p_exp2 + 64;\n}\n\n// wuffs_base__private_implementation__medium_prec_bin__as_f64 converts m to a\n// double (what C calls a double-precision float64).\n//\n// Preconditions:\n// - m is non-NULL.\n// - m->mantissa is non-zero.\n// - m->mantissa's high bit is set (i.e. m is normalized).\nstatic double //\nwuffs_base__private_implementation__medium_prec_bin__as_f64(\n const wuffs_base__private_implementation__mediu" +
"m_prec_bin* m,\n bool negative) {\n uint64_t mantissa64 = m->mantissa;\n // An mpb's mantissa has the implicit (binary) decimal point at the right\n // hand end of the mantissa's explicit digits. A double-precision's mantissa\n // has that decimal point near the left hand end. There's also an explicit\n // versus implicit leading 1 bit (binary digit). Together, the difference in\n // semantics corresponds to adding 63.\n int32_t exp2 = m->exp2 + 63;\n\n // Ensure that exp2 is at least -1022, the minimum double-precision exponent\n // for normal (as opposed to subnormal) numbers.\n if (-1022 > exp2) {\n uint32_t n = (uint32_t)(-1022 - exp2);\n mantissa64 >>= n;\n exp2 += (int32_t)n;\n }\n\n // Extract the (1 + 52) bits from the 64-bit mantissa64. 52 is the number of\n // explicit mantissa bits in a double-precision f64.\n //\n // Before, we have 64 bits and due to normalization, the high bit 'H' is 1.\n // 63 55 47 etc 15 7\n // H210_9876_5432_1098_7654_etc_etc_etc_5432_109" +
"8_7654_3210\n // ++++_++++_++++_++++_++++_etc_etc_etc_++++_+..._...._.... Kept bits.\n // ...._...._...H_2109_8765_etc_etc_etc_6543_2109_8765_4321 After shifting.\n // After, we have 53 bits (and bit #52 is this 'H' bit).\n uint64_t mantissa53 = mantissa64 >> 11;\n\n // Round up if the old bit #10 (the highest bit dropped by shifting) was set.\n // We also fix any overflow from rounding up.\n if (mantissa64 & 1024) {\n mantissa53++;\n if ((mantissa53 >> 53) != 0) {\n mantissa53 >>= 1;\n exp2++;\n }\n }\n\n // Handle double-precision infinity (a nominal exponent of 1024) and\n // subnormals (an exponent of -1023 and no implicit mantissa bit, bit #52).\n if (exp2 >= 1024) {\n mantissa53 = 0;\n exp2 = 1024;\n } else if ((mantissa53 >> 52) == 0) {\n exp2 = -1023;\n }\n\n // Pack the bits and return.\n const int32_t f64_bias = -1023;\n uint64_t exp2_bits =\n (uint64_t)((exp2 - f64_bias) & 0x07FF); // (1 << 11) - 1.\n uint64_t bits = (mantissa53 & 0x000FFFFFFFFFFFFF) | // (1 << 52" +
- ") - 1.\n (exp2_bits << 52) | //\n (negative ? 0x8000000000000000 : 0); // (1 << 63).\n return wuffs_base__ieee_754_bit_representation__to_f64(bits);\n}\n\n// wuffs_base__private_implementation__medium_prec_bin__parse_number_f64\n// converts from an HPD to a double, using an MPB as scratch space. It returns\n// a NULL status.repr if there is no ambiguity in the truncation or rounding to\n// a float64 (an IEEE 754 double-precision floating point value).\n//\n// It may modify m even if it returns a non-NULL status.repr.\nstatic wuffs_base__result_f64 //\nwuffs_base__private_implementation__medium_prec_bin__parse_number_f64(\n wuffs_base__private_implementation__medium_prec_bin* m,\n const wuffs_base__private_implementation__high_prec_dec* h,\n bool skip_fast_path_for_tests) {\n // m->mantissa is a uint64_t, which is an integer approximation to a rational\n // value - h's underlying digits after m's normalization. This error is an\n // upper bound on the difference " +
- "between the approximate and actual value.\n //\n // The DiyFpStrtod function in https://github.com/google/double-conversion\n // uses a finer grain (1/8th of the ULP, Unit in the Last Place) when\n // tracking error. This implementation is coarser (1 ULP) but simpler.\n //\n // It is an error in the \"numerical approximation\" sense, not in the typical\n // programming sense (as in \"bad input\" or \"a result type\").\n uint64_t error = 0;\n\n // Convert up to 19 decimal digits (in h->digits) to 64 binary digits (in\n // m->mantissa): (1e19 < (1<<64)) and ((1<<64) < 1e20). If we have more than\n // 19 digits, we're truncating (with error).\n uint32_t i;\n uint32_t i_end = h->num_digits;\n if (i_end > 19) {\n i_end = 19;\n error = 1;\n }\n uint64_t mantissa = 0;\n for (i = 0; i < i_end; i++) {\n mantissa = (10 * mantissa) + h->digits[i];\n }\n m->mantissa = mantissa;\n m->exp2 = 0;\n\n // Check that exp10 lies in the (big_powers_of_10 + small_powers_of_10)\n // range, -348 ..= +347, stepping big_powers_of_10 by " +
- "8 (which is 87 triples)\n // and small_powers_of_10 by 1 (which is 8 triples).\n int32_t exp10 = h->decimal_point - ((int32_t)(i_end));\n if (exp10 < -348) {\n goto fail;\n }\n uint32_t bpo10 = ((uint32_t)(exp10 + 348)) / 8;\n uint32_t spo10 = ((uint32_t)(exp10 + 348)) % 8;\n if (bpo10 >= 87) {\n goto fail;\n }\n\n // Try a fast path, if float64 math would be exact.\n //\n // 15 is such that 1e15 can be losslessly represented in a float64 mantissa:\n // (1e15 < (1<<53)) and ((1<<53) < 1e16).\n //\n // 22 is the maximum valid index for the\n // wuffs_base__private_implementation__f64_powers_of_10 array.\n do {\n if (skip_fast_path_for_tests || ((mantissa >> 52) != 0)) {\n break;\n }\n double d = (double)mantissa;\n\n if (exp10 == 0) {\n wuffs_base__result_f64 ret;\n ret.status.repr = NULL;\n ret.value = h->negative ? -d : +d;\n return ret;\n\n } else if (exp10 > 0) {\n if (exp10 > 22) {\n if (exp10 > (15 + 22)) {\n break;\n }\n // If exp10 is in the " +
- "range 23 ..= 37, try moving a few of the zeroes\n // from the exponent to the mantissa. If we're still under 1e15, we\n // haven't truncated any mantissa bits.\n if (exp10 > 22) {\n d *= wuffs_base__private_implementation__f64_powers_of_10[exp10 - 22];\n exp10 = 22;\n if (d >= 1e15) {\n break;\n }\n }\n }\n d *= wuffs_base__private_implementation__f64_powers_of_10[exp10];\n wuffs_base__result_f64 ret;\n ret.status.repr = NULL;\n ret.value = h->negative ? -d : +d;\n return ret;\n\n } else { // \"if (exp10 < 0)\" is effectively \"if (true)\" here.\n if (exp10 < -22) {\n break;\n }\n d /= wuffs_base__private_implementation__f64_powers_of_10[-exp10];\n wuffs_base__result_f64 ret;\n ret.status.repr = NULL;\n ret.value = h->negative ? -d : +d;\n return ret;\n }\n } while (0);\n\n // Normalize (and scale the error).\n error <<= wuffs_base__private_implementation__medium_prec_bin__normalize(m);\n\n" +
- " // Multiplying two MPB values nominally multiplies two mantissas, call them A\n // and B, which are integer approximations to the precise values (A+a) and\n // (B+b) for some error terms a and b.\n //\n // MPB multiplication calculates (((A+a) * (B+b)) >> 64) to be ((A*B) >> 64).\n // Shifting (truncating) and rounding introduces further error. The\n // difference between the calculated result:\n // ((A*B ) >> 64)\n // and the true result:\n // ((A*B + A*b + a*B + a*b) >> 64) + rounding_error\n // is:\n // (( A*b + a*B + a*b) >> 64) + rounding_error\n // which can be re-grouped as:\n // ((A*b) >> 64) + ((a*(B+b)) >> 64) + rounding_error\n //\n // Now, let A and a be \"m->mantissa\" and \"error\", and B and b be the\n // pre-calculated power of 10. A and B are both less than (1 << 64), a is the\n // \"error\" local variable and b is less than 1.\n //\n // An upper bound (in absolute value) on ((A*b) >> 64) is therefore 1.\n //\n // Similarly, an upper bound on ((a*(B+b)) >> 64) is a, " +
- "also known as error.\n //\n // Finally, the rounding_error is at most 1.\n //\n // In total, calling mpb__mul_pow_10 will increase the worst-case error by 2.\n // The subsequent re-normalization can multiply that by a further factor.\n\n // Multiply by small_powers_of_10[etc].\n wuffs_base__private_implementation__medium_prec_bin__mul_pow_10(\n m, &wuffs_base__private_implementation__small_powers_of_10[3 * spo10]);\n error += 2;\n error <<= wuffs_base__private_implementation__medium_prec_bin__normalize(m);\n\n // Multiply by big_powers_of_10[etc].\n wuffs_base__private_implementation__medium_prec_bin__mul_pow_10(\n m, &wuffs_base__private_implementation__big_powers_of_10[3 * bpo10]);\n error += 2;\n error <<= wuffs_base__private_implementation__medium_prec_bin__normalize(m);\n\n // We have a good approximation of h, but we still have to check whether the\n // error is small enough. Equivalently, whether the number of surplus\n // mantissa bits (the bits dropped when going from m's 64 mantissa bits to\n /" +
- "/ the smaller number of double-precision mantissa bits) would always round\n // up or down, even when perturbed by ±error. We start at 11 surplus bits (m\n // has 64, double-precision has 1+52), but it can be higher for subnormals.\n //\n // In many cases, the error is small enough and we return true.\n const int32_t f64_bias = -1023;\n int32_t subnormal_exp2 = f64_bias - 63;\n uint32_t surplus_bits = 11;\n if (subnormal_exp2 >= m->exp2) {\n surplus_bits += 1 + ((uint32_t)(subnormal_exp2 - m->exp2));\n }\n\n uint64_t surplus_mask = (((uint64_t)1) << surplus_bits) - 1; // e.g. 0x07FF.\n uint64_t surplus = m->mantissa & surplus_mask;\n uint64_t halfway = ((uint64_t)1) << (surplus_bits - 1); // e.g. 0x0400.\n\n // Do the final calculation in *signed* arithmetic.\n int64_t i_surplus = (int64_t)surplus;\n int64_t i_halfway = (int64_t)halfway;\n int64_t i_error = (int64_t)error;\n\n if ((i_surplus > (i_halfway - i_error)) &&\n (i_surplus < (i_halfway + i_error))) {\n goto fail;\n }\n\n wuffs_base__result_f64" +
- " ret;\n ret.status.repr = NULL;\n ret.value = wuffs_base__private_implementation__medium_prec_bin__as_f64(\n m, h->negative);\n return ret;\n\nfail:\n do {\n wuffs_base__result_f64 ret;\n ret.status.repr = \"#base: mpb__parse_number_f64 failed\";\n ret.value = 0;\n return ret;\n } while (0);\n}\n\n" +
+ ") - 1.\n (exp2_bits << 52) | //\n (negative ? 0x8000000000000000 : 0); // (1 << 63).\n return wuffs_base__ieee_754_bit_representation__to_f64(bits);\n}\n\n// wuffs_base__private_implementation__medium_prec_bin__parse_number_f64\n// converts from an HPD to a double, using an MPB as scratch space. It returns\n// a NULL status.repr if there is no ambiguity in the truncation or rounding to\n// a float64 (an IEEE 754 double-precision floating point value).\n//\n// It may modify m even if it returns a non-NULL status.repr.\nstatic wuffs_base__result_f64 //\nwuffs_base__private_implementation__medium_prec_bin__parse_number_f64(\n wuffs_base__private_implementation__medium_prec_bin* m,\n const wuffs_base__private_implementation__high_prec_dec* h,\n bool skip_fast_path_for_tests) {\n do {\n // m->mantissa is a uint64_t, which is an integer approximation to a\n // rational value - h's underlying digits after m's normalization. This\n // error is an upper bound on th" +
+ "e difference between the approximate and\n // actual value.\n //\n // The DiyFpStrtod function in https://github.com/google/double-conversion\n // uses a finer grain (1/8th of the ULP, Unit in the Last Place) when\n // tracking error. This implementation is coarser (1 ULP) but simpler.\n //\n // It is an error in the \"numerical approximation\" sense, not in the\n // typical programming sense (as in \"bad input\" or \"a result type\").\n uint64_t error = 0;\n\n // Convert up to 19 decimal digits (in h->digits) to 64 binary digits (in\n // m->mantissa): (1e19 < (1<<64)) and ((1<<64) < 1e20). If we have more\n // than 19 digits, we're truncating (with error).\n uint32_t i;\n uint32_t i_end = h->num_digits;\n if (i_end > 19) {\n i_end = 19;\n error = 1;\n }\n uint64_t mantissa = 0;\n for (i = 0; i < i_end; i++) {\n mantissa = (10 * mantissa) + h->digits[i];\n }\n m->mantissa = mantissa;\n m->exp2 = 0;\n\n // Check that exp10 lies in the (big_powers_of_10 + small_po" +
+ "wers_of_10)\n // range, -348 ..= +347, stepping big_powers_of_10 by 8 (which is 87\n // triples) and small_powers_of_10 by 1 (which is 8 triples).\n int32_t exp10 = h->decimal_point - ((int32_t)(i_end));\n if (exp10 < -348) {\n goto fail;\n }\n uint32_t bpo10 = ((uint32_t)(exp10 + 348)) / 8;\n uint32_t spo10 = ((uint32_t)(exp10 + 348)) % 8;\n if (bpo10 >= 87) {\n goto fail;\n }\n\n // Try a fast path, if float64 math would be exact.\n //\n // 15 is such that 1e15 can be losslessly represented in a float64\n // mantissa: (1e15 < (1<<53)) and ((1<<53) < 1e16).\n //\n // 22 is the maximum valid index for the\n // wuffs_base__private_implementation__f64_powers_of_10 array.\n do {\n if (skip_fast_path_for_tests || ((mantissa >> 52) != 0)) {\n break;\n }\n double d = (double)mantissa;\n\n if (exp10 == 0) {\n wuffs_base__result_f64 ret;\n ret.status.repr = NULL;\n ret.value = h->negative ? -d : +d;\n return ret;\n\n } else if (e" +
+ "xp10 > 0) {\n if (exp10 > 22) {\n if (exp10 > (15 + 22)) {\n break;\n }\n // If exp10 is in the range 23 ..= 37, try moving a few of the zeroes\n // from the exponent to the mantissa. If we're still under 1e15, we\n // haven't truncated any mantissa bits.\n if (exp10 > 22) {\n d *= wuffs_base__private_implementation__f64_powers_of_10[exp10 -\n 22];\n exp10 = 22;\n if (d >= 1e15) {\n break;\n }\n }\n }\n d *= wuffs_base__private_implementation__f64_powers_of_10[exp10];\n wuffs_base__result_f64 ret;\n ret.status.repr = NULL;\n ret.value = h->negative ? -d : +d;\n return ret;\n\n } else { // \"if (exp10 < 0)\" is effectively \"if (true)\" here.\n if (exp10 < -22) {\n break;\n }\n d /= wuffs_base__private_implementation__f64_powers_of_10[-exp10];\n wuffs_bas" +
+ "e__result_f64 ret;\n ret.status.repr = NULL;\n ret.value = h->negative ? -d : +d;\n return ret;\n }\n } while (0);\n\n // Normalize (and scale the error).\n error <<= wuffs_base__private_implementation__medium_prec_bin__normalize(m);\n\n // Multiplying two MPB values nominally multiplies two mantissas, call them\n // A and B, which are integer approximations to the precise values (A+a)\n // and (B+b) for some error terms a and b.\n //\n // MPB multiplication calculates (((A+a) * (B+b)) >> 64) to be ((A*B) >>\n // 64). Shifting (truncating) and rounding introduces further error. The\n // difference between the calculated result:\n // ((A*B ) >> 64)\n // and the true result:\n // ((A*B + A*b + a*B + a*b) >> 64) + rounding_error\n // is:\n // (( A*b + a*B + a*b) >> 64) + rounding_error\n // which can be re-grouped as:\n // ((A*b) >> 64) + ((a*(B+b)) >> 64) + rounding_error\n //\n // Now, let A and a be \"m->mantissa\" and \"erro" +
+ "r\", and B and b be the\n // pre-calculated power of 10. A and B are both less than (1 << 64), a is\n // the \"error\" local variable and b is less than 1.\n //\n // An upper bound (in absolute value) on ((A*b) >> 64) is therefore 1.\n //\n // An upper bound on ((a*(B+b)) >> 64) is a, also known as error.\n //\n // Finally, the rounding_error is at most 1.\n //\n // In total, calling mpb__mul_pow_10 will raise the worst-case error by 2.\n // The subsequent re-normalization can multiply that by a further factor.\n\n // Multiply by small_powers_of_10[etc].\n wuffs_base__private_implementation__medium_prec_bin__mul_pow_10(\n m, &wuffs_base__private_implementation__small_powers_of_10[3 * spo10]);\n error += 2;\n error <<= wuffs_base__private_implementation__medium_prec_bin__normalize(m);\n\n // Multiply by big_powers_of_10[etc].\n wuffs_base__private_implementation__medium_prec_bin__mul_pow_10(\n m, &wuffs_base__private_implementation__big_powers_of_10[3 * bpo10]);\n err" +
+ "or += 2;\n error <<= wuffs_base__private_implementation__medium_prec_bin__normalize(m);\n\n // We have a good approximation of h, but we still have to check whether\n // the error is small enough. Equivalently, whether the number of surplus\n // mantissa bits (the bits dropped when going from m's 64 mantissa bits to\n // the smaller number of double-precision mantissa bits) would always round\n // up or down, even when perturbed by ±error. We start at 11 surplus bits\n // (m has 64, double-precision has 1+52), but it can be higher for\n // subnormals.\n //\n // In many cases, the error is small enough and we return true.\n const int32_t f64_bias = -1023;\n int32_t subnormal_exp2 = f64_bias - 63;\n uint32_t surplus_bits = 11;\n if (subnormal_exp2 >= m->exp2) {\n surplus_bits += 1 + ((uint32_t)(subnormal_exp2 - m->exp2));\n }\n\n uint64_t surplus_mask =\n (((uint64_t)1) << surplus_bits) - 1; // e.g. 0x07FF.\n uint64_t surplus = m->mantissa & surplus_mask;\n uint64_t" +
+ " halfway = ((uint64_t)1) << (surplus_bits - 1); // e.g. 0x0400.\n\n // Do the final calculation in *signed* arithmetic.\n int64_t i_surplus = (int64_t)surplus;\n int64_t i_halfway = (int64_t)halfway;\n int64_t i_error = (int64_t)error;\n\n if ((i_surplus > (i_halfway - i_error)) &&\n (i_surplus < (i_halfway + i_error))) {\n goto fail;\n }\n\n wuffs_base__result_f64 ret;\n ret.status.repr = NULL;\n ret.value = wuffs_base__private_implementation__medium_prec_bin__as_f64(\n m, h->negative);\n return ret;\n } while (0);\n\nfail:\n do {\n wuffs_base__result_f64 ret;\n ret.status.repr = \"#base: mpb__parse_number_f64 failed\";\n ret.value = 0;\n return ret;\n } while (0);\n}\n\n" +
"" +
"// --------\n\nwuffs_base__result_f64 //\nwuffs_base__parse_number_f64_special(wuffs_base__slice_u8 s,\n const char* fallback_status_repr) {\n do {\n uint8_t* p = s.ptr;\n uint8_t* q = s.ptr + s.len;\n\n for (; (p < q) && (*p == '_'); p++) {\n }\n if (p >= q) {\n goto fallback;\n }\n\n // Parse sign.\n bool negative = false;\n do {\n if (*p == '+') {\n p++;\n } else if (*p == '-') {\n negative = true;\n p++;\n } else {\n break;\n }\n for (; (p < q) && (*p == '_'); p++) {\n }\n } while (0);\n if (p >= q) {\n goto fallback;\n }\n\n bool nan = false;\n switch (p[0]) {\n case 'I':\n case 'i':\n if (((q - p) < 3) || //\n ((p[1] != 'N') && (p[1] != 'n')) || //\n ((p[2] != 'F') && (p[2] != 'f'))) {\n goto fallback;\n }\n p += 3;\n\n if ((p >= q) || (*p == '_')) {\n break;\n } else if (((q - p) < 5) || " +
" //\n ((p[0] != 'I') && (p[0] != 'i')) || //\n ((p[1] != 'N') && (p[1] != 'n')) || //\n ((p[2] != 'I') && (p[2] != 'i')) || //\n ((p[3] != 'T') && (p[3] != 't')) || //\n ((p[4] != 'Y') && (p[4] != 'y'))) {\n goto fallback;\n }\n p += 5;\n\n if ((p >= q) || (*p == '_')) {\n break;\n }\n goto fallback;\n\n case 'N':\n case 'n':\n if (((q - p) < 3) || //\n ((p[1] != 'A') && (p[1] != 'a')) || //\n ((p[2] != 'N') && (p[2] != 'n'))) {\n goto fallback;\n }\n p += 3;\n\n if ((p >= q) || (*p == '_')) {\n nan = true;\n break;\n }\n goto fallback;\n\n default:\n goto fallback;\n }\n\n // Finish.\n for (; (p < q) && (*p == '_'); p++) {\n }\n if (p != q) {\n goto fallback;\n }\n wuffs_base__result_f64 ret;\n ret.status.repr = NULL;\n ret.value = w" +
diff --git a/release/c/wuffs-unsupported-snapshot.c b/release/c/wuffs-unsupported-snapshot.c
index 042f052..8262813 100644
--- a/release/c/wuffs-unsupported-snapshot.c
+++ b/release/c/wuffs-unsupported-snapshot.c
@@ -9814,178 +9814,184 @@
wuffs_base__private_implementation__medium_prec_bin* m,
const wuffs_base__private_implementation__high_prec_dec* h,
bool skip_fast_path_for_tests) {
- // m->mantissa is a uint64_t, which is an integer approximation to a rational
- // value - h's underlying digits after m's normalization. This error is an
- // upper bound on the difference between the approximate and actual value.
- //
- // The DiyFpStrtod function in https://github.com/google/double-conversion
- // uses a finer grain (1/8th of the ULP, Unit in the Last Place) when
- // tracking error. This implementation is coarser (1 ULP) but simpler.
- //
- // It is an error in the "numerical approximation" sense, not in the typical
- // programming sense (as in "bad input" or "a result type").
- uint64_t error = 0;
-
- // Convert up to 19 decimal digits (in h->digits) to 64 binary digits (in
- // m->mantissa): (1e19 < (1<<64)) and ((1<<64) < 1e20). If we have more than
- // 19 digits, we're truncating (with error).
- uint32_t i;
- uint32_t i_end = h->num_digits;
- if (i_end > 19) {
- i_end = 19;
- error = 1;
- }
- uint64_t mantissa = 0;
- for (i = 0; i < i_end; i++) {
- mantissa = (10 * mantissa) + h->digits[i];
- }
- m->mantissa = mantissa;
- m->exp2 = 0;
-
- // Check that exp10 lies in the (big_powers_of_10 + small_powers_of_10)
- // range, -348 ..= +347, stepping big_powers_of_10 by 8 (which is 87 triples)
- // and small_powers_of_10 by 1 (which is 8 triples).
- int32_t exp10 = h->decimal_point - ((int32_t)(i_end));
- if (exp10 < -348) {
- goto fail;
- }
- uint32_t bpo10 = ((uint32_t)(exp10 + 348)) / 8;
- uint32_t spo10 = ((uint32_t)(exp10 + 348)) % 8;
- if (bpo10 >= 87) {
- goto fail;
- }
-
- // Try a fast path, if float64 math would be exact.
- //
- // 15 is such that 1e15 can be losslessly represented in a float64 mantissa:
- // (1e15 < (1<<53)) and ((1<<53) < 1e16).
- //
- // 22 is the maximum valid index for the
- // wuffs_base__private_implementation__f64_powers_of_10 array.
do {
- if (skip_fast_path_for_tests || ((mantissa >> 52) != 0)) {
- break;
+ // m->mantissa is a uint64_t, which is an integer approximation to a
+ // rational value - h's underlying digits after m's normalization. This
+ // error is an upper bound on the difference between the approximate and
+ // actual value.
+ //
+ // The DiyFpStrtod function in https://github.com/google/double-conversion
+ // uses a finer grain (1/8th of the ULP, Unit in the Last Place) when
+ // tracking error. This implementation is coarser (1 ULP) but simpler.
+ //
+ // It is an error in the "numerical approximation" sense, not in the
+ // typical programming sense (as in "bad input" or "a result type").
+ uint64_t error = 0;
+
+ // Convert up to 19 decimal digits (in h->digits) to 64 binary digits (in
+ // m->mantissa): (1e19 < (1<<64)) and ((1<<64) < 1e20). If we have more
+ // than 19 digits, we're truncating (with error).
+ uint32_t i;
+ uint32_t i_end = h->num_digits;
+ if (i_end > 19) {
+ i_end = 19;
+ error = 1;
}
- double d = (double)mantissa;
+ uint64_t mantissa = 0;
+ for (i = 0; i < i_end; i++) {
+ mantissa = (10 * mantissa) + h->digits[i];
+ }
+ m->mantissa = mantissa;
+ m->exp2 = 0;
- if (exp10 == 0) {
- wuffs_base__result_f64 ret;
- ret.status.repr = NULL;
- ret.value = h->negative ? -d : +d;
- return ret;
+ // Check that exp10 lies in the (big_powers_of_10 + small_powers_of_10)
+ // range, -348 ..= +347, stepping big_powers_of_10 by 8 (which is 87
+ // triples) and small_powers_of_10 by 1 (which is 8 triples).
+ int32_t exp10 = h->decimal_point - ((int32_t)(i_end));
+ if (exp10 < -348) {
+ goto fail;
+ }
+ uint32_t bpo10 = ((uint32_t)(exp10 + 348)) / 8;
+ uint32_t spo10 = ((uint32_t)(exp10 + 348)) % 8;
+ if (bpo10 >= 87) {
+ goto fail;
+ }
- } else if (exp10 > 0) {
- if (exp10 > 22) {
- if (exp10 > (15 + 22)) {
- break;
- }
- // If exp10 is in the range 23 ..= 37, try moving a few of the zeroes
- // from the exponent to the mantissa. If we're still under 1e15, we
- // haven't truncated any mantissa bits.
- if (exp10 > 22) {
- d *= wuffs_base__private_implementation__f64_powers_of_10[exp10 - 22];
- exp10 = 22;
- if (d >= 1e15) {
- break;
- }
- }
- }
- d *= wuffs_base__private_implementation__f64_powers_of_10[exp10];
- wuffs_base__result_f64 ret;
- ret.status.repr = NULL;
- ret.value = h->negative ? -d : +d;
- return ret;
-
- } else { // "if (exp10 < 0)" is effectively "if (true)" here.
- if (exp10 < -22) {
+ // Try a fast path, if float64 math would be exact.
+ //
+ // 15 is such that 1e15 can be losslessly represented in a float64
+ // mantissa: (1e15 < (1<<53)) and ((1<<53) < 1e16).
+ //
+ // 22 is the maximum valid index for the
+ // wuffs_base__private_implementation__f64_powers_of_10 array.
+ do {
+ if (skip_fast_path_for_tests || ((mantissa >> 52) != 0)) {
break;
}
- d /= wuffs_base__private_implementation__f64_powers_of_10[-exp10];
- wuffs_base__result_f64 ret;
- ret.status.repr = NULL;
- ret.value = h->negative ? -d : +d;
- return ret;
+ double d = (double)mantissa;
+
+ if (exp10 == 0) {
+ wuffs_base__result_f64 ret;
+ ret.status.repr = NULL;
+ ret.value = h->negative ? -d : +d;
+ return ret;
+
+ } else if (exp10 > 0) {
+ if (exp10 > 22) {
+ if (exp10 > (15 + 22)) {
+ break;
+ }
+ // If exp10 is in the range 23 ..= 37, try moving a few of the zeroes
+ // from the exponent to the mantissa. If we're still under 1e15, we
+ // haven't truncated any mantissa bits.
+ if (exp10 > 22) {
+ d *= wuffs_base__private_implementation__f64_powers_of_10[exp10 -
+ 22];
+ exp10 = 22;
+ if (d >= 1e15) {
+ break;
+ }
+ }
+ }
+ d *= wuffs_base__private_implementation__f64_powers_of_10[exp10];
+ wuffs_base__result_f64 ret;
+ ret.status.repr = NULL;
+ ret.value = h->negative ? -d : +d;
+ return ret;
+
+ } else { // "if (exp10 < 0)" is effectively "if (true)" here.
+ if (exp10 < -22) {
+ break;
+ }
+ d /= wuffs_base__private_implementation__f64_powers_of_10[-exp10];
+ wuffs_base__result_f64 ret;
+ ret.status.repr = NULL;
+ ret.value = h->negative ? -d : +d;
+ return ret;
+ }
+ } while (0);
+
+ // Normalize (and scale the error).
+ error <<= wuffs_base__private_implementation__medium_prec_bin__normalize(m);
+
+ // Multiplying two MPB values nominally multiplies two mantissas, call them
+ // A and B, which are integer approximations to the precise values (A+a)
+ // and (B+b) for some error terms a and b.
+ //
+ // MPB multiplication calculates (((A+a) * (B+b)) >> 64) to be ((A*B) >>
+ // 64). Shifting (truncating) and rounding introduces further error. The
+ // difference between the calculated result:
+ // ((A*B ) >> 64)
+ // and the true result:
+ // ((A*B + A*b + a*B + a*b) >> 64) + rounding_error
+ // is:
+ // (( A*b + a*B + a*b) >> 64) + rounding_error
+ // which can be re-grouped as:
+ // ((A*b) >> 64) + ((a*(B+b)) >> 64) + rounding_error
+ //
+ // Now, let A and a be "m->mantissa" and "error", and B and b be the
+ // pre-calculated power of 10. A and B are both less than (1 << 64), a is
+ // the "error" local variable and b is less than 1.
+ //
+ // An upper bound (in absolute value) on ((A*b) >> 64) is therefore 1.
+ //
+ // An upper bound on ((a*(B+b)) >> 64) is a, also known as error.
+ //
+ // Finally, the rounding_error is at most 1.
+ //
+ // In total, calling mpb__mul_pow_10 will raise the worst-case error by 2.
+ // The subsequent re-normalization can multiply that by a further factor.
+
+ // Multiply by small_powers_of_10[etc].
+ wuffs_base__private_implementation__medium_prec_bin__mul_pow_10(
+ m, &wuffs_base__private_implementation__small_powers_of_10[3 * spo10]);
+ error += 2;
+ error <<= wuffs_base__private_implementation__medium_prec_bin__normalize(m);
+
+ // Multiply by big_powers_of_10[etc].
+ wuffs_base__private_implementation__medium_prec_bin__mul_pow_10(
+ m, &wuffs_base__private_implementation__big_powers_of_10[3 * bpo10]);
+ error += 2;
+ error <<= wuffs_base__private_implementation__medium_prec_bin__normalize(m);
+
+ // We have a good approximation of h, but we still have to check whether
+ // the error is small enough. Equivalently, whether the number of surplus
+ // mantissa bits (the bits dropped when going from m's 64 mantissa bits to
+ // the smaller number of double-precision mantissa bits) would always round
+ // up or down, even when perturbed by ±error. We start at 11 surplus bits
+ // (m has 64, double-precision has 1+52), but it can be higher for
+ // subnormals.
+ //
+ // In many cases, the error is small enough and we return true.
+ const int32_t f64_bias = -1023;
+ int32_t subnormal_exp2 = f64_bias - 63;
+ uint32_t surplus_bits = 11;
+ if (subnormal_exp2 >= m->exp2) {
+ surplus_bits += 1 + ((uint32_t)(subnormal_exp2 - m->exp2));
}
+
+ uint64_t surplus_mask =
+ (((uint64_t)1) << surplus_bits) - 1; // e.g. 0x07FF.
+ uint64_t surplus = m->mantissa & surplus_mask;
+ uint64_t halfway = ((uint64_t)1) << (surplus_bits - 1); // e.g. 0x0400.
+
+ // Do the final calculation in *signed* arithmetic.
+ int64_t i_surplus = (int64_t)surplus;
+ int64_t i_halfway = (int64_t)halfway;
+ int64_t i_error = (int64_t)error;
+
+ if ((i_surplus > (i_halfway - i_error)) &&
+ (i_surplus < (i_halfway + i_error))) {
+ goto fail;
+ }
+
+ wuffs_base__result_f64 ret;
+ ret.status.repr = NULL;
+ ret.value = wuffs_base__private_implementation__medium_prec_bin__as_f64(
+ m, h->negative);
+ return ret;
} while (0);
- // Normalize (and scale the error).
- error <<= wuffs_base__private_implementation__medium_prec_bin__normalize(m);
-
- // Multiplying two MPB values nominally multiplies two mantissas, call them A
- // and B, which are integer approximations to the precise values (A+a) and
- // (B+b) for some error terms a and b.
- //
- // MPB multiplication calculates (((A+a) * (B+b)) >> 64) to be ((A*B) >> 64).
- // Shifting (truncating) and rounding introduces further error. The
- // difference between the calculated result:
- // ((A*B ) >> 64)
- // and the true result:
- // ((A*B + A*b + a*B + a*b) >> 64) + rounding_error
- // is:
- // (( A*b + a*B + a*b) >> 64) + rounding_error
- // which can be re-grouped as:
- // ((A*b) >> 64) + ((a*(B+b)) >> 64) + rounding_error
- //
- // Now, let A and a be "m->mantissa" and "error", and B and b be the
- // pre-calculated power of 10. A and B are both less than (1 << 64), a is the
- // "error" local variable and b is less than 1.
- //
- // An upper bound (in absolute value) on ((A*b) >> 64) is therefore 1.
- //
- // Similarly, an upper bound on ((a*(B+b)) >> 64) is a, also known as error.
- //
- // Finally, the rounding_error is at most 1.
- //
- // In total, calling mpb__mul_pow_10 will increase the worst-case error by 2.
- // The subsequent re-normalization can multiply that by a further factor.
-
- // Multiply by small_powers_of_10[etc].
- wuffs_base__private_implementation__medium_prec_bin__mul_pow_10(
- m, &wuffs_base__private_implementation__small_powers_of_10[3 * spo10]);
- error += 2;
- error <<= wuffs_base__private_implementation__medium_prec_bin__normalize(m);
-
- // Multiply by big_powers_of_10[etc].
- wuffs_base__private_implementation__medium_prec_bin__mul_pow_10(
- m, &wuffs_base__private_implementation__big_powers_of_10[3 * bpo10]);
- error += 2;
- error <<= wuffs_base__private_implementation__medium_prec_bin__normalize(m);
-
- // We have a good approximation of h, but we still have to check whether the
- // error is small enough. Equivalently, whether the number of surplus
- // mantissa bits (the bits dropped when going from m's 64 mantissa bits to
- // the smaller number of double-precision mantissa bits) would always round
- // up or down, even when perturbed by ±error. We start at 11 surplus bits (m
- // has 64, double-precision has 1+52), but it can be higher for subnormals.
- //
- // In many cases, the error is small enough and we return true.
- const int32_t f64_bias = -1023;
- int32_t subnormal_exp2 = f64_bias - 63;
- uint32_t surplus_bits = 11;
- if (subnormal_exp2 >= m->exp2) {
- surplus_bits += 1 + ((uint32_t)(subnormal_exp2 - m->exp2));
- }
-
- uint64_t surplus_mask = (((uint64_t)1) << surplus_bits) - 1; // e.g. 0x07FF.
- uint64_t surplus = m->mantissa & surplus_mask;
- uint64_t halfway = ((uint64_t)1) << (surplus_bits - 1); // e.g. 0x0400.
-
- // Do the final calculation in *signed* arithmetic.
- int64_t i_surplus = (int64_t)surplus;
- int64_t i_halfway = (int64_t)halfway;
- int64_t i_error = (int64_t)error;
-
- if ((i_surplus > (i_halfway - i_error)) &&
- (i_surplus < (i_halfway + i_error))) {
- goto fail;
- }
-
- wuffs_base__result_f64 ret;
- ret.status.repr = NULL;
- ret.value = wuffs_base__private_implementation__medium_prec_bin__as_f64(
- m, h->negative);
- return ret;
-
fail:
do {
wuffs_base__result_f64 ret;
