A divide under a non-default rounding mode drops the remainder of its quotient

**TL;DR** — Under a mode which is not `fe_dec_to_nearest`, the divide truncates its wide
quotient and gives it to the constructor without the remainder. A mode which must move
away from the truncated quotient cannot.

```cpp
fesetround(rounding_mode::fe_dec_upward);
const auto x {1_DL / 3_DL}; // x is 0.3333333333333333333333333333333333 it must end in 334
fesetround(rounding_mode::fe_dec_downward);
const auto y {-1_DL / 3_DL}; // y is -0.3333333333333333333333333333333333 it must end in 334
fesetround(rounding_mode::fe_dec_upward);
const auto z {1_DF / 1.000001_DF}; // z is 0.999999 it must be 0.9999991
```

The wide quotient of the 128 bit types holds zero or one digit past the precision, thus
the constructor has almost nothing to round. The wide quotient of the 32 and 64 bit types
holds five or six digits past the precision, thus the defect shows only when all of those
digits are zero and the remainder is not.

- Of 10000 random quotients of `decimal128_t`, 2793 are wrong under the upward mode, 2797
under the downward mode and 2563 under `fe_dec_to_nearest_from_zero`. `decimal_fast128_t`
gives 2751, 2687 and 2410. The toward zero mode asks for the truncated quotient, thus it
is correct.
- Of 200000 random quotients of `decimal32_t` and of `decimal_fast32_t`, 2 and 1 are wrong
under the downward mode. `1 / 1.000001` under the upward mode shows the shape directly,
and so does `1 / 1.000000000000001` for the 64 bit types.
- No test of the suite divides under a directed mode.

**The correction** — the remainder goes in as one sticky digit after the wide quotient,
thus the constructor rounds with the exact sticky bit in every mode. The 128 bit divide
also takes one digit more, thus the constructor has a quotient digit to round before the
sticky digit. The correction is +21/-19 lines in `div_impl.hpp`, and 0 of 240000 random
quotients are wrong after it.

The parts that follow give the measurements, the exact values, and the programs.

The quotients of the six types

The fast types give the same values as their plain types, thus the table shows the three
widths. The correct values come from the `decimal` module of Python at the precision of the
type.

| type | mode | quotient | library | correct | agreement |
| --- | --- | --- | --- | --- | --- |
| `decimal128_t` | upward | `1 / 3` | 3.333333333333333333333333333333333e-01 | 3.333333333333333333333333333333334e-01 | **no** |
| `decimal128_t` | downward | `-1 / 3` | -3.333333333333333333333333333333333e-01 | -3.333333333333333333333333333333334e-01 | **no** |
| `decimal128_t` | to nearest from zero | `2 / 3` | 6.666666666666666666666666666666666e-01 | 6.666666666666666666666666666666667e-01 | **no** |
| `decimal128_t` | upward | `1 / 1.000000000000000000000000000000001` | 9.99999999999999999999999999999999e-01 | 9.999999999999999999999999999999991e-01 | **no** |
| `decimal128_t` | toward zero | `-1 / 3` | -3.333333333333333333333333333333333e-01 | -3.333333333333333333333333333333333e-01 | yes |
| `decimal64_t` | upward | `1 / 3` | 3.333333333333334e-01 | 3.333333333333334e-01 | yes |
| `decimal64_t` | downward | `-1 / 3` | -3.333333333333334e-01 | -3.333333333333334e-01 | yes |
| `decimal64_t` | to nearest from zero | `2 / 3` | 6.666666666666667e-01 | 6.666666666666667e-01 | yes |
| `decimal64_t` | upward | `1 / 1.000000000000001` | 9.99999999999999e-01 | 9.999999999999991e-01 | **no** |
| `decimal32_t` | upward | `1 / 3` | 3.333334e-01 | 3.333334e-01 | yes |
| `decimal32_t` | downward | `-1 / 3` | -3.333334e-01 | -3.333334e-01 | yes |
| `decimal32_t` | to nearest from zero | `2 / 3` | 6.666667e-01 | 6.666667e-01 | yes |
| `decimal32_t` | upward | `1 / 1.000001` | 9.99999e-01 | 9.999991e-01 | **no** |
| `decimal32_t` | downward | `-2.554295e+01 / 3.328344e+01` | -7.674372e-01 | -7.674373e-01 | **no** |

`1 / 1.000001` is `0.999999000000999999...`, thus the digits 8 to 12 are zero and the
remainder is not. The wide quotient of `decimal32_t` holds 12 digits, thus the constructor
sees `999999000000` and it has no reason to go up. The last row is a random pair of the
same shape: `-0.76743720000096...` has zeros in the digits 8 to 12.

Show the program of the table

```cpp
// Prints a/b for the six types under the four non-default modes.
#include
#include
#include
using namespace boost::decimal;

template
void show(const char* name, const char* a, const char* b)
{
static const rounding_mode modes[] {rounding_mode::fe_dec_upward, rounding_mode::fe_dec_downward,
rounding_mode::fe_dec_toward_zero, rounding_mode::fe_dec_to_nearest_from_zero};
static const char* mode_names[] {"upward", "downward", "toward_zero", "nearest_from_zero"};
T x {}, y {};
from_chars(a, a + std::strlen(a), x);
from_chars(b, b + std::strlen(b), y);
for (int i = 0; i < 4; ++i)
{
fesetround(modes[i]);
const T q {x / y};
char buf[64] {};
auto r = to_chars(buf, buf + sizeof(buf), q, chars_format::scientific);
*r.ptr = 0;
std::printf("%-6s %-18s %s / %s = %s\n", name, mode_names[i], a, b, buf);
}
fesetround(rounding_mode::fe_dec_to_nearest);
}

int main(int argc, char** argv)
{
const char* a = argc > 1 ? argv[1] : "1";
const char* b = argc > 2 ? argv[2] : "3";
show("d32", a, b);
show("f32", a, b);
show("d64", a, b);
show("f64", a, b);
show("d128", a, b);
show("f128", a, b);
return 0;
}
```

```
$ ./probe 1 3 | grep upward
d32 upward 1 / 3 = 3.333334e-01
f32 upward 1 / 3 = 3.333334e-01
d64 upward 1 / 3 = 3.333333333333334e-01
f64 upward 1 / 3 = 3.333333333333334e-01
d128 upward 1 / 3 = 3.333333333333333333333333333333333e-01
f128 upward 1 / 3 = 3.333333333333333333333333333333333e-01
$ ./probe 1 1.000001 | grep '^d32'
d32 upward 1 / 1.000001 = 9.99999e-01
d32 downward 1 / 1.000001 = 9.99999e-01
d32 toward_zero 1 / 1.000001 = 9.99999e-01
d32 nearest_from_zero 1 / 1.000001 = 9.99999e-01
```

Environment

- Boost.Decimal at commit `bad17ca9` of `develop`
- GCC 16.2.1, with `-std=c++17 -O2`, on x86-64 Linux
- The correct values come from the `decimal` module of Python 3.14.7, which follows the
same arithmetic standard as this library, and from the `_Decimal64` and `_Decimal128`
of GCC, which are the BID library of Intel

The cause

Each of the three divide drivers in `include/boost/decimal/detail/div_impl.hpp` has a
branch for the modes which are not `fe_dec_to_nearest`. The branch of the 128 bit types:

```cpp
if (BOOST_DECIMAL_UNLIKELY(!impl::div_default_rounding(lhs_c.sig)))
{
const auto wide_tens {pow10(int128::uint128_t{static_cast(precision_v)})};
const auto wide_sig {detail::umul256(lhs_c.sig, wide_tens)};
auto wide_q {wide_sig / rhs_c.sig};
auto wide_exp {lhs_c.exp - rhs_c.exp - static_cast(precision_v)};

if (wide_q[3] != 0U || wide_q[2] != 0U)
{
...
}

BOOST_DECIMAL_ASSERT((wide_q[3] | wide_q[2]) == 0U);
return DecimalType{int128::uint128_t{wide_q[1], wide_q[0]}, wide_exp, sign};
}
```

`wide_sig / rhs_c.sig` drops its remainder, and nothing reads it. The constructor gets
`wide_q` and it rounds the digits past the precision with `fenv_round`, which reads the
mode. The remainder is not one of those digits, thus the constructor sees an exact
quotient when the digits are zero.

Both operands have 34 digits after `expand_significand`, thus `wide_q` has 34 or 35
digits. With 34 digits the constructor rounds nothing, and `1 / 3` is `0.333...333` in
every mode. With 35 digits the constructor rounds one digit, and the modes which move away
from zero are correct only when that digit is not zero.

The branches of the 32 and 64 bit types are the same with `uint64_t` and `uint128_t`.
Their `wide_offset` is `digits10 - precision`, thus `wide_q` has 12 or 13 digits for
`decimal32_t` and 22 or 23 for `decimal64_t`. The constructor rounds five or six digits,
and the remainder matters only when all of them are zero.

The block for a `wide_q` of more than 128 bits is dead. `lhs_c.sig` is below `10^34` and
`rhs_c.sig` is at least `10^33`, thus `wide_q` is below `10^35`.

The correction

The remainder goes in as one sticky digit after the wide quotient, and the exponent goes
down by one. The constructor then drops the sticky digit with the other digits, thus
`fenv_round` gets the exact sticky bit. The branch of the 32 bit types:

```cpp
const auto wide_sig {static_cast(lhs_c.sig) * wide_tens};
const auto wide_div {static_cast(rhs_c.sig)};
const auto wide_q {wide_sig / wide_div};
// The constructor rounds by the digits it drops, thus the remainder goes in as one digit
const auto sticky {static_cast(wide_sig - wide_q * wide_div != 0U)};
const auto wide_exp {(lhs_c.exp - static_cast(wide_offset) - 1) - rhs_c.exp};
return DecimalType{wide_q * 10U + sticky, wide_exp, sign};
```

The branch of the 64 bit types is the same with `uint128_t`. The branch of the 128 bit
types scales by `10^35` in place of `10^34`, thus the quotient has one or two digits past
the precision before the sticky digit. A sticky digit alone cannot decide a tie of
`fe_dec_to_nearest_from_zero`, thus that quotient digit is necessary. `impl::div_mod`
gives the quotient and the remainder at the cost of the one divide, and the dead block
goes away:

```cpp
// One digit past the precision for the constructor to round, then the remainder goes in
// as one more digit. Both operands have 34 digits, thus the quotient has at most 37.
const auto wide_tens {pow10(int128::uint128_t{static_cast(precision_v + 1)})};
const auto wide_sig {detail::umul256(lhs_c.sig, wide_tens)};
const auto wide_dr {impl::div_mod(wide_sig, rhs_c.sig)};
const auto wide_exp {lhs_c.exp - rhs_c.exp - static_cast(precision_v) - 2};

BOOST_DECIMAL_ASSERT((wide_dr.quotient[3] | wide_dr.quotient[2]) == 0U);
const int128::uint128_t wide_q {wide_dr.quotient[1], wide_dr.quotient[0]};
const auto sticky {static_cast(wide_dr.remainder != u256{})};
return DecimalType{wide_q * 10U + sticky, wide_exp, sign};
```

The `decimal` module of Python and the BID library of Intel give the values of the
correction:

| mode | quotient | Boost.Decimal before | after | Python `decimal` | Intel BID |
| --- | --- | --- | --- | --- | --- |
| upward | `1 / 3` at 34 digits | 0.3333333333333333333333333333333333 | 0.3333333333333333333333333333333334 | 0.3333333333333333333333333333333334 | 0.3333333333333333333333333333333334 |
| downward | `-1 / 3` at 34 digits | -0.3333333333333333333333333333333333 | -0.3333333333333333333333333333333334 | -0.3333333333333333333333333333333334 | -0.3333333333333333333333333333333334 |
| to nearest from zero | `2 / 3` at 34 digits | 0.6666666666666666666666666666666666 | 0.6666666666666666666666666666666667 | 0.6666666666666666666666666666666667 | 0.6666666666666666666666666666666667 |
| upward | `1 / 1.000000000000000000000000000000001` | 0.999999999999999999999999999999999 | 0.9999999999999999999999999999999991 | 0.9999999999999999999999999999999991 | 0.9999999999999999999999999999999991 |
| upward | `1 / 1.000000000000001` at 16 digits | 0.999999999999999 | 0.9999999999999991 | 0.9999999999999991 | 0.9999999999999991 |

Show the program of the BID column

```c
/* GCC's _Decimal64 and _Decimal128 are Intel BID from libgcc. Its rounding mode is the global below. */
#include
extern __thread int __bid_IDEC_glbround;
enum { BID_NEAREST = 0, BID_DOWN = 1, BID_UP = 2, BID_TO_ZERO = 3, BID_NEAREST_AWAY = 4 };
static const char* names[] = {"nearest", "downward", "upward", "toward_zero", "nearest_from_zero"};

/* Prints the 34 digits of q in (-1, 1). The operands are volatile, thus GCC cannot fold the divide. */
static void show128(const char* label, _Decimal128 q)
{
char sign = q < 0 ? '-' : '+';
_Decimal128 m = q < 0 ? -q : q;
unsigned long long hi = (unsigned long long)(m * 1e17DL);
unsigned long long lo = (unsigned long long)((m * 1e17DL - hi) * 1e17DL);
printf("%-28s %c0.%017llu%017llu\n", label, sign, hi, lo);
}

static void show64(const char* label, _Decimal64 q)
{
char sign = q < 0 ? '-' : '+';
_Decimal64 m = q < 0 ? -q : q;
unsigned long long d = (unsigned long long)(m * 1e16DD);
printf("%-28s %c0.%016llu\n", label, sign, d);
}

int main(void)
{
volatile _Decimal128 one = 1.DL, two = 2.DL, three = 3.DL, near_one = 1.000000000000000000000000000000001DL;
volatile _Decimal64 one64 = 1.DD, near_one64 = 1.000000000000001DD;
const int modes[] = {BID_UP, BID_DOWN, BID_TO_ZERO, BID_NEAREST_AWAY};
for (int i = 0; i < 4; ++i)
{
__bid_IDEC_glbround = modes[i];
printf("== %s\n", names[modes[i]]);
show128("1/3", one / three);
show128("-1/3", -one / three);
show128("2/3", two / three);
show128("1/1.000...001 (34 digits)", one / near_one);
show64("1/1.000000000000001", one64 / near_one64);
}
return 0;
}
```

```
== upward
1/3 +0.3333333333333333333333333333333334
-1/3 -0.3333333333333333333333333333333333
2/3 +0.6666666666666666666666666666666667
1/1.000...001 (34 digits) +0.9999999999999999999999999999999991
1/1.000000000000001 +0.9999999999999991
== downward
1/3 +0.3333333333333333333333333333333333
-1/3 -0.3333333333333333333333333333333334
2/3 +0.6666666666666666666666666666666666
1/1.000...001 (34 digits) +0.9999999999999999999999999999999990
1/1.000000000000001 +0.9999999999999990
== toward_zero
1/3 +0.3333333333333333333333333333333333
-1/3 -0.3333333333333333333333333333333333
2/3 +0.6666666666666666666666666666666666
1/1.000...001 (34 digits) +0.9999999999999999999999999999999990
1/1.000000000000001 +0.9999999999999990
== nearest_from_zero
1/3 +0.3333333333333333333333333333333333
-1/3 -0.3333333333333333333333333333333333
2/3 +0.6666666666666666666666666666666667
1/1.000...001 (34 digits) +0.9999999999999999999999999999999990
1/1.000000000000001 +0.9999999999999990
```

What the correction changes: 10000 pairs of each type

Random operands with seven, sixteen or 34 digits and with exponents from `10^-3` to
`10^3`, 10000 pairs of each type, and each pair under the four modes which are not
`fe_dec_to_nearest`. The counts are the quotients which are not equal to the value of the
`decimal` module of Python at the precision of the type:

| type | upward | downward | toward zero | to nearest from zero |
| --- | --- | --- | --- | --- |
| `decimal32_t` | 0 → 0 | 0 → 0 | 0 → 0 | 0 → 0 |
| `decimal64_t` | 0 → 0 | 0 → 0 | 0 → 0 | 0 → 0 |
| `decimal128_t` | 2793 → 0 | 2797 → 0 | 0 → 0 | 2563 → 0 |
| `decimal_fast32_t` | 0 → 0 | 0 → 0 | 0 → 0 | 0 → 0 |
| `decimal_fast64_t` | 0 → 0 | 0 → 0 | 0 → 0 | 0 → 0 |
| `decimal_fast128_t` | 2751 → 0 | 2687 → 0 | 0 → 0 | 2410 → 0 |

The run holds 240000 quotients, of which 16001 are wrong before the correction and 0
after it.

The 32 and 64 bit types need five or six zero digits past the precision, thus a random
pair shows the defect about one time in 100000. A second run of 200000 pairs of each 32
bit type gives 2 wrong quotients of `decimal32_t` and 1 of `decimal_fast32_t`, all under
the downward mode, and 0 after the correction. One of them is
`-3.28626e-03 / 6.292506e-01`, which is `-5.222498e-03` before and `-5.222499e-03`
after.

A C++ filter prints the quotient of each pair under the four modes, and a Python driver
compares against the `decimal` module. I can add both programs in a comment.

What the correction costs

`perf stat` of 400000 quotients of random operands with seven digits, under the default
mode as a control and under the upward mode. The mean of five runs:

| type | mode | instructions before | instructions after | cycles before | cycles after |
| --- | --- | --- | --- | --- | --- |
| `decimal32_t` | to nearest | 128.8 M | 128.8 M | 60.1 M | 60.3 M |
| `decimal32_t` | upward | 187.4 M | 189.0 M | 82.3 M | 87.4 M |
| `decimal64_t` | to nearest | 250.1 M | 250.1 M | 113.2 M | 119.3 M |
| `decimal64_t` | upward | 366.0 M | 375.3 M | 137.8 M | 148.8 M |
| `decimal128_t` | to nearest | 827.3 M | 819.7 M | 383.7 M | 382.3 M |
| `decimal128_t` | upward | 797.5 M | 966.1 M | 376.1 M | 420.5 M |
| `decimal_fast32_t` | to nearest | 109.0 M | 109.0 M | 66.9 M | 67.9 M |
| `decimal_fast32_t` | upward | 125.4 M | 126.8 M | 76.5 M | 78.3 M |
| `decimal_fast64_t` | to nearest | 217.3 M | 217.9 M | 133.6 M | 134.7 M |
| `decimal_fast64_t` | upward | 286.1 M | 296.6 M | 144.3 M | 142.0 M |
| `decimal_fast128_t` | to nearest | 725.3 M | 725.3 M | 325.6 M | 342.3 M |
| `decimal_fast128_t` | upward | 676.7 M | 856.8 M | 332.8 M | 376.0 M |

The default mode does not run the branch, and its instruction counts move by 1% or less.
The 32 and 64 bit types run 1% to 4% more instructions under the upward mode, which is the
multiply for the remainder and the sticky digit. The 128 bit types run 21% to 27% more
instructions and 12% to 13% more cycles under the upward mode. That is the rounding which
the constructor now does on every quotient: before the correction it did nothing for a
quotient with 34 digits.

Why the tests do not find it

`test_upward_rounding.cpp` and `test_downward_rounding.cpp` add and subtract. The random
tests of the arithmetic run under the default mode, which takes the other path with
`div_finalize_*` and the remainder. No test divides under a directed mode.

The 32 and 64 bit types are correct at `1 / 3` and at most other quotients, thus a test of
those types needs the shape with zero digits past the precision.

Start in include/boost/decimal/detail/div_impl.hpp, focusing on the non-default rounding branches of the three divide drivers and the existing division tests. Reproduce the directed-mode examples from the issue, then verify that quotient remainders affect rounding and that the reported random cases produce the expected values.

Written by the indexing model from the issue text.