Compute frexp with integer arithmetic

[ibmibmibm]

frexp now gives the correctly rounded result for all finite inputs in the current rounding mode, and it is 16 to 39 times faster.

The old frexp_impl has four defects:

• The result is not correctly rounded. frexp_impl multiplies by 2^-t from pow_2_impl and then divides by powers of 2. Each step rounds, so the error is up to 852 ulp.
• The loop does not stop below 1e-91 (decimal32), 1e-369 (decimal64) and 1e-6111 (decimal128). The fast types return 0.
• A subnormal input gives 0.
• -0 gives +0.

The new frexp_impl

• Writes v * 10^n as s * 5^K * 2^K, and gets s * 5^K as a binary value of 64, 128 or 256 bits. The value comes from a table of 5^(+-d * 16^j), so that |K| needs at most four products.
• Gets the digits with one shift. fenv_round rounds them once in the current rounding mode, so that the constructor gets p digits.
• Uses n = p + 1 digits for decimal32 and decimal64. decimal128 uses n = p digits. fenv_round then gets the parity of the last digit and a digit for the bits below, so that it does not divide 128-bit digits by 10.
• First computes s * 5^K in a word of 64 or 128 bits for decimal64 and decimal128. This is half of the width that an exact result can need. After nm products, that M is below the exact value by less than 6 * nm + 1 units. If the bits below the p digits are not that close to 0, to one half or to a carry, they give the result. Else the wide pass runs.
• Knows that M is exact if s and the table entries fit in the word together.
• Gets the products of decimal128 from umul256_hi and umul512_hi. For K > 0, the first narrow product of decimal128 uses 64 bits of the entry, because 5^1 to 5^15 are exact in 64 bits. A static_assert stops a radix with a digit above 27.
• Removes the factors of 5 only when K >= -1.5 * p, because for a smaller K no result is exact. A product with the inverse of 5 modulo 2^w removes each of them without a division.
• Finds NaN, infinity and zero without fpclassify, which counts the digits. A zero is returned as it is, so that it keeps its sign.
• Makes the result with pack_in_range, because the result has p digits and an exponent in range. pack_in_range does not do the checks of the constructor.
• Ends all finite paths in one pack_in_range. With one return for each path, GCC merges the results of the fast types with a memory copy, and that copy stalls on the stores of the single members.
• Cuts the table entries to the word width at compile time, so that a program keeps only the widths that it uses. decimal32 keeps a table of 736 bytes and decimal64 a table of 1472 bytes instead of 2944 bytes.

Other changes

• countl_zero gives the correct count for types of fewer than 64 bits on compilers without __builtin_clz, for example MSVC. Before, the fallback counted the leading zeros of a 64-bit value, so frexp of decimal32 gave 0.
• frexp computes in the type of the value, also when BOOST_DECIMAL_DEC_EVAL_METHOD is 1 or 2. Before, frexp computed in the wider evaluation type, and the cast back to the type of the value can round the fraction to 1. For example, frexp(9223372e12_DF) gave 1 * 2^63 and not 0.5 * 2^64.
• A first commit adds umul256_hi, the high 128 bits of the product of two uint128_t, and computes umul512_hi with a schoolbook product of 64-bit limbs. The old umul512_hi used four umul256 products, with branches and carry compares. The value is the same. frexp uses both functions.

Results

• Accuracy: a GMP reference computes the correctly rounded result in the five rounding modes. The inputs are random values over the full range, powers of 2 and their neighbors, subnormals, short significands and significands with factors of 5. Seeds 7 and 11 give 1104710 results each. All results are correctly rounded, in the range [0.5, 1), and have the correct exponent.
• Proof for the wide pass: its M is below the exact value by less than 6 * nm units. A search over all K and all s finds the exact values nearest to a rounding boundary. The nearest one is 38.7 units away for decimal32 (bound 13), 703.6 units for decimal64 (bound 19) and 169830 units for decimal128 (bound 25). A brute-force scan gives the same cases for all s of decimal32 and for s < 10^9 of decimal64.
• Speed (GCC 15, -O3, |v| in [1e-20, 1e20], ns for each call, develop to this PR): decimal32 176.9 to 6.8, decimal64 272.0 to 9.8, decimal128 736.0 to 18.7, decimal_fast32 130.4 to 8.0, decimal_fast64 323.1 to 16.1, decimal_fast128 693.4 to 27.1. This is 16 to 39 times faster.
• Speed of umul512_hi (GCC 15, -O3, ns for each call): divmod_pow10_u256, which rounds decimal128 results, goes from 74 to 78 ns to 63 to 64 ns. Multiplication, division and fma of decimal128 and decimal_fast128 do not change more than the noise.
• Code size (GCC, -O2, .text of one frexp call): decimal32 10449 to 986 bytes, decimal64 15692 to 2243 bytes, decimal128 33508 to 3484 bytes.

Fixes #1487

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 98.8%. Comparing base (c157871) to head (d44c69e).

Additional details and impacted files

@@            Coverage Diff            @@
##           develop   #1488     +/-   ##
=========================================
+ Coverage     98.8%   98.8%   +0.1%     
=========================================
  Files          318     319      +1     
  Lines        26496   26596    +100     
  Branches      2252    2258      +6     
=========================================
+ Hits         26160   26260    +100     
  Misses         336     336             

Files with missing lines	Coverage Δ
include/boost/decimal/detail/cmath/frexp.hpp	100.0% <100.0%> (ø)
include/boost/decimal/detail/countl.hpp	100.0% <ø> (ø)
include/boost/decimal/detail/u256.hpp	95.3% <100.0%> (-0.1%)	⬇️
test/github_issue_1487_eval_method.cpp	100.0% <100.0%> (ø)
test/test_frexp_ldexp.cpp	100.0% <100.0%> (ø)

---

Continue to review full report in Codecov by Harness.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update c157871...d44c69e. Read the comment docs.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[ckormanyos]

Hi @ibmibmibm I believe you can find the 128 and 256 bit integer operations in the existing codebase already without writing them from scratch. I did not check the exact operations, but normal ones should be available. If you need some specialized optimized ones, however, they might not be already present in the repo.

[ibmibmibm]

Duplicated functions removed, thanks