A decimal128_t product of short operands near the smallest exponent is wrong

[ibmibmibm]

TL;DR

Some decimal128_t products have wrong digits and a wrong exponent. The error can be a
factor of 10 or more. Two conditions cause it. The exponents of the operands add to near the
smallest exponent. The product of the two significands has 35 to 38 digits. The defect is in
each rounding mode.

const auto a {34e-6168_DL * 6304777571276319608815211782449451e-6_DL};
// a is 6.670563082001761558497347434477494e-6141, it must be 2.143624374233948666997172006032813e-6139
const auto b {9999e-6160_DL * 9999999999999999999999999999999999e-10_DL};
// b is 7.131692053359185016762684537821425e-6137, it must be 9.998999999999999999999999999999999e-6133

• The multiply has a separate path for an exponent sum below -6142. It does not expand the
  significands. It rounds only a product with more than 38 digits. The limit must be 34,
  the precision of the type. A product of 35 to 38 digits goes to the encoder with too many
  digits. The encoder then writes digits into bits outside the significand field.
• In a sweep of 20000 such products, 962 are wrong in each of the five rounding modes.
• The correction changes one line. After it, no product of the sweep is wrong.
• The correction needs the correction of The constructor rounds some values below the smallest subnormal wrongly, and some of them stop the program #1476 (Round a value below the smallest subnormal in each rounding mode #1477), which is in develop. Without it,
  a few products of 35 to 38 digits reach a defect of coefficient_rounding.

The values

The decimal module of Python at 34 digits gives the correct values.

operation	library	correct
34e-6168_DL * 6304777571276319608815211782449451e-6_DL	6.670563082001761558497347434477494e-6141	2.143624374233948666997172006032813e-6139
99e-6150_DL * 1234567890123456789012345678901234e-20_DL	7.991690234456014284551302968380054e-6137	1.222222211222222221122222222112222e-6135
9999e-6160_DL * 9999999999999999999999999999999999e-10_DL	7.131692053359185016762684537821425e-6137	9.998999999999999999999999999999999e-6133

The exact products have 36, 36 and 38 digits.

Show the program of the table

#include <boost/decimal.hpp>
#include <iomanip>
#include <iostream>

using namespace boost::decimal;
using namespace boost::decimal::literals;

int main()
{
    std::cout << std::scientific << std::setprecision(33);
    std::cout << 34e-6168_DL * 6304777571276319608815211782449451e-6_DL << '\n'
              << "  it must be 2.143624374233948666997172006032813e-6139\n";
    std::cout << 99e-6150_DL * 1234567890123456789012345678901234e-20_DL << '\n'
              << "  it must be 1.222222211222222221122222222112222e-6135\n";
    std::cout << 9999e-6160_DL * 9999999999999999999999999999999999e-10_DL << '\n'
              << "  it must be 9.998999999999999999999999999999999e-6133\n";
}

Its output with develop:

6.670563082001761558497347434477494e-6141
  it must be 2.143624374233948666997172006032813e-6139
7.991690234456014284551302968380054e-6137
  it must be 1.222222211222222221122222222112222e-6135
7.131692053359185016762684537821425e-6137
  it must be 9.998999999999999999999999999999999e-6133

A sweep takes operands with 1 to 34 random digits and random signs. The exponents add to a
value from -6222 to -6143. The decimal module of Python gives the correct values.

mode	wrong with develop	wrong with the correction
fe_dec_to_nearest	962 of 20000	0
fe_dec_to_nearest_from_zero	962 of 20000	0
fe_dec_toward_zero	962 of 20000	0
fe_dec_upward	962 of 20000	0
fe_dec_downward	962 of 20000	0

Show the program of the sweep

// decimal128_t products whose exponent sum is below -bias + 34 (the subnormal bail-out of the
// multiply). Operands have 1 to 34 random digits. Prints "mode a b result".
#include <boost/decimal.hpp>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <random>
#include <string>

using namespace boost::decimal;
using S = decimal128_t::significand_type;

static std::string digits_of(S v)
{
    std::string s;
    do
    {
        s.insert(s.begin(), static_cast<char>('0' + static_cast<int>(v % 10U)));
        v /= 10U;
    } while (v != 0U);
    return s;
}

static std::string str(decimal128_t x)
{
    int e {};
    const auto sig = frexp10(x, &e);
    return std::string(signbit(x) ? "-" : "") + digits_of(sig) + "e" + std::to_string(e);
}

static S rand_sig(std::mt19937_64& g, int d)
{
    std::uniform_int_distribution<int> digit(0, 9);
    S v {static_cast<unsigned>(std::uniform_int_distribution<int>(1, 9)(g))};
    for (int i = 1; i < d; ++i) v = v * 10U + static_cast<unsigned>(digit(g));
    return v;
}

int main(int argc, char** argv)
{
    const char* m = argc > 1 ? argv[1] : "near";
    if (std::strcmp(m, "from0") == 0) fesetround(rounding_mode::fe_dec_to_nearest_from_zero);
    else if (std::strcmp(m, "up") == 0) fesetround(rounding_mode::fe_dec_upward);
    else if (std::strcmp(m, "down") == 0) fesetround(rounding_mode::fe_dec_downward);
    else if (std::strcmp(m, "zero") == 0) fesetround(rounding_mode::fe_dec_toward_zero);
    std::mt19937_64 g(99);
    std::uniform_int_distribution<int> dd(1, 34), ed(-6176, -20), sgn(0, 1);
    const int n {argc > 2 ? std::atoi(argv[2]) : 20000};
    for (int i = 0; i < n; ++i)
    {
        const int da {dd(g)}, db {dd(g)};
        const int ea {ed(g)};
        // The exponent sum is below -6176 + 34, and the product is near the subnormal range
        const int eb {std::uniform_int_distribution<int>(-6176 + 34 - ea - 80, -6176 + 34 - ea - 1)(g)};
        const decimal128_t a {rand_sig(g, da), ea, sgn(g) == 1};
        const decimal128_t b {rand_sig(g, db), eb, sgn(g) == 1};
        std::printf("%s %s %s %s\n", m, str(a).c_str(), str(b).c_str(), str(a * b).c_str());
    }
}

#!/usr/bin/env python3
# mul128_check.py < mul128_sweep output (any modes): wrong products per mode, with examples
import sys
from collections import Counter
from decimal import (Context, Decimal, ROUND_CEILING, ROUND_DOWN, ROUND_FLOOR,
                     ROUND_HALF_EVEN, ROUND_HALF_UP)

rnd = {"near": ROUND_HALF_EVEN, "from0": ROUND_HALF_UP, "up": ROUND_CEILING,
       "down": ROUND_FLOOR, "zero": ROUND_DOWN}
bad = Counter()
total = Counter()
ex = {}
for line in sys.stdin:
    mode, a, b, r = line.split()
    c = Context(prec=34, Emin=-6143, Emax=6144, rounding=rnd[mode], clamp=1, traps=[])
    exact = c.multiply(Decimal(a), Decimal(b))
    got = Decimal(r)
    total[mode] += 1
    if exact != got or exact.is_signed() != got.is_signed():
        bad[mode] += 1
        ex.setdefault(mode, []).append(f"{a} * {b} = {r}, correct {exact}")
for m in ["near", "from0", "zero", "up", "down"]:
    if m in total:
        print(m, "wrong:", bad[m], "of", total[m])
        for e in ex.get(m, [])[:2]:
            print("   ", e)

Environment

• Boost.Decimal at commit 7c79789a of develop
• GCC 16.2.1 and clang 22.1.8, on x86-64 Linux
• The decimal module of Python 3.14.7

The cause

mul_impl.hpp:498 sends a product with an exponent sum below -bias + 34 to a separate path.
The comment says that it saves time for products that become zero. This path does not
expand the significands to 34 digits. The product can have 1 to 68 digits.
mul_impl.hpp:503:

        const auto sig_dig {detail::num_digits(res_sig)};
        const auto digit_delta {sig_dig - std::numeric_limits<sig_type>::digits10};
        if (BOOST_DECIMAL_LIKELY(digit_delta > 0))
        {
            auto biased_exp {res_exp_mut + detail::bias_v<ReturnType>};
            detail::coefficient_rounding<ReturnType>(res_sig, res_exp_mut, biased_exp, sign, sig_dig);
        }

sig_type is the 128-bit significand type, and its digits10 is 38. The type holds 34
digits. A product of 35 to 38 digits does not go to coefficient_rounding. If its exponent
is in the range of the type, pack_in_range then encodes it directly. The encoder writes the
digits above the 34th into bits outside the significand field.

mul_impl.hpp:543 has the same expression. On that path the two significands have 34 digits
each, and the product has 67 or 68 digits. The limit then does not change the result.

The correction

The limit is the precision of the type:

diff --git a/include/boost/decimal/detail/mul_impl.hpp b/include/boost/decimal/detail/mul_impl.hpp
index 6c085767..7322109f 100644
--- a/include/boost/decimal/detail/mul_impl.hpp
+++ b/include/boost/decimal/detail/mul_impl.hpp
@@ -500,7 +500,7 @@ BOOST_DECIMAL_CUDA_CONSTEXPR auto d128_mul_impl(const T1& lhs_sig_in, const U1 l
         auto res_sig {detail::umul256(lhs_sig_in, rhs_sig_in)};
         auto res_exp_mut {static_cast<typename ReturnType::biased_exponent_type>(lhs_exp_in + rhs_exp_in)};
         const auto sig_dig {detail::num_digits(res_sig)};
-        const auto digit_delta {sig_dig - std::numeric_limits<sig_type>::digits10};
+        const auto digit_delta {sig_dig - detail::precision_v<ReturnType>};
         if (BOOST_DECIMAL_LIKELY(digit_delta > 0))
         {
             auto biased_exp {res_exp_mut + detail::bias_v<ReturnType>};

After this change, products of 35 to 38 digits go to coefficient_rounding. That function has
a defect for a 128-bit value with few digits and a large shift (#1476). For example,
9999999999999999999e-3112_DL * 9999999999999999999e-3113_DL gives 1e-6176 with this change
alone, but 0 is correct. #1477 corrected that defect, and it is in develop.

A new test has the three values of the table, a tie, and the product above. It also
multiplies such products in each of the four other rounding modes. The expected values come
from the decimal module of Python. Before the correction, each check fails, except the
product above: the old path does not round it at all, and it gives 0. The test passes with
GCC 16 and clang 22.

What the correction changes

A second program compares develop at 1e058ab3 and this correction bit by bit.
mul_impl.hpp did not change after 1e058ab3. It uses 400000 pairs for
each of the six types, and the four operations in the five modes, thus 48 million results.
No result changes. The operands of that program do not give products of 35 to 38 digits near
the smallest exponent.

What the correction costs

The change is in the path for an exponent sum below -6142. The other products do not use
it. The fast types also use this path. Their significands always have 34 digits, thus the
product always has more than 38 digits, and the change does not affect them.

Google Benchmark 1.9.5, GCC 16 with -O3, one core of a Ryzen 9 3900X. Each value is the
median time of one operation over 7 runs, with develop at 1e058ab3. "p" is the U test of
compare.py. The sets are:

• stock: the operands of test/benchmarks.cpp, with random exponents over the whole range.
• full: random operands with all digits, and exponents 0 to 4 apart.
• sameexp: random operands with all digits and equal exponents.
• band: exponents 4 to 37 apart.
• short: operands with 1 to 6 digits and two decimals.

Other programs used many cores of the computer during the runs. Changes of a few percent are
noise.

op	type	set	mode	develop ns	new ns	change	p
mul	decimal128_t	stock	near	40.6	40.5	0%	0.053
mul	decimal_fast128_t	stock	near	54.6	56.8	+4%	0.001
mul	decimal128_t	full	near	36.2	35.5	-2%	0.001
mul	decimal_fast128_t	full	near	32.9	32.1	-2%	0.001
mul	decimal128_t	full	up	175.8	172.0	-2%	0.026
mul	decimal_fast128_t	full	up	166.1	174.7	+5%	0.001
mul	decimal128_t	sameexp	near	36.1	35.9	0%	0.018
mul	decimal_fast128_t	sameexp	near	33.0	32.8	-1%	0.097
mul	decimal128_t	sameexp	up	177.3	172.0	-3%	0.001
mul	decimal_fast128_t	sameexp	up	167.7	174.1	+4%	0.001
mul	decimal128_t	band	near	36.3	35.5	-2%	0.026
mul	decimal_fast128_t	band	near	33.0	32.3	-2%	0.001
mul	decimal128_t	band	up	177.7	173.4	-2%	0.001
mul	decimal_fast128_t	band	up	168.1	175.5	+4%	0.001
mul	decimal128_t	short	near	6.7	6.4	-3%	0.001
mul	decimal_fast128_t	short	near	29.6	29.5	0%	0.620
mul	decimal128_t	short	up	6.5	6.5	0%	0.710
mul	decimal_fast128_t	short	up	115.8	124.2	+7%	0.001

Why the tests do not find it

The tests of the multiply use operands with all 34 digits, or products far from the smallest
exponent. No test multiplies short decimal128_t operands whose exponents add to near the
smallest exponent.