sqrt is not correctly rounded, in the default mode or in the directed modes

[ibmibmibm]

TL;DR — IEEE 754-2019 5.4.1 lists squareRoot among the operations which give the
exact result rounded in the current mode. The 32 and 64 bit kernels of the library
truncate the root. An odd exponent multiplies the root by a rounded sqrt10 in the three
kernels. No kernel reads the mode. The decimal module of Python gives the nearest value.

fesetround(rounding_mode::fe_dec_to_nearest);
const auto a {sqrt(decimal32_t {2, 0})};    // a is 1.414213           it must be 1.414214
const auto b {sqrt(decimal64_t {26, 0})};   // b is 5.099019513592781  it must be 5.099019513592785
const auto c {sqrt(decimal32_t {1, -1})};   // c is 0.3162277          it must be 0.3162278
fesetround(rounding_mode::fe_dec_upward);
const auto d {sqrt(decimal64_t {2, 0})};    // d is 1.414213562373095  it must be 1.414213562373096

• Over 20000 random values, decimal32_t, decimal64_t and decimal128_t miss the
  nearest value 57%, 74% and 37% of the time, by up to five steps. The fast types give
  the same values.
• In the three directed modes the value of the mode is missed 33% to 99% of the time.
• No test compares sqrt against a decimal reference value in any mode.

The parts that follow give the values, the cause, the correction, and the measurements.

The values

The default mode is set. The nearest value comes from the decimal module of Python
3.14.7 at 7, 16 and 34 digits with ROUND_HALF_EVEN. "steps" is the distance in units of
the last digit. A negative value is too low.

input	type	sqrt of the library	the nearest value	steps
2	decimal32_t	1.414213	1.414214	-1
3	decimal32_t	1.732050	1.732051	-1
5	decimal32_t	2.236067	2.236068	-1
5	decimal64_t	2.236067977499789	2.236067977499790	-1
7	decimal64_t	2.645751311064590	2.645751311064591	-1
11	decimal32_t	3.316622	3.316625	-3
11	decimal64_t	3.316624790355398	3.316624790355400	-2
11	decimal128_t	3.316624790355399849114932736670688	3.316624790355399849114932736670687	+1
26	decimal32_t	5.099018	5.099020	-2
26	decimal64_t	5.099019513592781	5.099019513592785	-4
48	decimal64_t	6.928203230275507	6.928203230275509	-2
48	decimal128_t	6.928203230275509174109785366023492	6.928203230275509174109785366023489	+3

The inputs 2, 3, 5 and 7 have an even exponent. The inputs 11, 26 and 48 have an odd
exponent. A power of ten takes a separate path. An odd power over 1 is right, because the
path multiplies by the nearest sqrt10_v<T>. An odd power under 1 divides by it, and the
rounded quotient of the rounded constant is one step off for the three types:

input	type	sqrt of the library	the nearest value	steps
0.1	decimal32_t	0.3162277	0.3162278	-1
0.1	decimal64_t	0.3162277660168380	0.3162277660168379	+1
0.1	decimal128_t	0.3162277660168379331998893544432718	0.3162277660168379331998893544432719	-1

The misses over the integers 2 to 1000, split by the parity of the exponent. 908 inputs
have an even exponent, and 91 have an odd one:

type	even exponent	odd exponent	all
decimal32_t	470, each -1	61, from -3 to +1	531
decimal64_t	441, each -1	88, from -4 to -1	529
decimal128_t	0	70, from -1 to +3	70

decimal_fast32_t, decimal_fast64_t and decimal_fast128_t give the same three rows.

A sweep of 20000 random values per type gives the same picture. Each value has a random
16 digit significand and a random exponent across the range of the type. The steps are the
distance from the nearest value:

type	misses of 20000	steps
decimal32_t	11344	-3 to +2
decimal64_t	14825	-5 to -1
decimal128_t	7332	-1 to +3

The same sweep in the five modes, against the value of each mode. The floor, the ceiling
and the nearest value of each root come from an integer square root. The sqrt of Python
is not the reference here, because it rounds half even in every mode. The two nearest
modes give the same value, because a root is never a tie:

type	nearest, nearest from zero	toward zero, downward	upward
decimal32_t	59%, -3 to +2	33%, -3 to +1	85%, -3 to +2
decimal64_t	74%, -5 to -1	49%, -5 to -1	99%, -5 to -1
decimal128_t	37%, -1 to +3	61%, -2 to +3	62%, -2 to +3

The truncated root of the 32 and 64 bit kernels is the toward zero value by accident, for
an even exponent only. This is why the second column is the lowest.

Show the program of the random sweep

// Print x and sqrt(x) for random values with random exponents, per type and mode.
#include <boost/decimal.hpp>
#include <iostream>
#include <iomanip>
#include <random>

using namespace boost::decimal;

template <typename T>
static void dump(const char* tag, int digits, int emin, int emax, int count, const char* mname, rounding_mode mode)
{
    fesetround(mode);
    std::mt19937_64 gen(42);
    std::uniform_int_distribution<int> de(emin, emax);
    std::uniform_int_distribution<std::uint64_t> ds(1, 9999999999999999ULL);
    for (int i = 0; i < count; ++i)
    {
        const T x {ds(gen), de(gen)};
        std::cout << tag << ' ' << mname << ' ' << std::scientific << std::setprecision(digits - 1) << x << ' ' << sqrt(x) << '\n';
    }
}

template <typename T>
static void all(const char* tag, int digits, int emin, int emax, int count)
{
    dump<T>(tag, digits, emin, emax, count, "nearest", rounding_mode::fe_dec_to_nearest);
    dump<T>(tag, digits, emin, emax, count, "nearest_from_zero", rounding_mode::fe_dec_to_nearest_from_zero);
    dump<T>(tag, digits, emin, emax, count, "toward_zero", rounding_mode::fe_dec_toward_zero);
    dump<T>(tag, digits, emin, emax, count, "downward", rounding_mode::fe_dec_downward);
    dump<T>(tag, digits, emin, emax, count, "upward", rounding_mode::fe_dec_upward);
}

int main()
{
    all<decimal32_t>("d32", 7, -95, 90, 20000);
    all<decimal64_t>("d64", 16, -380, 370, 20000);
    all<decimal128_t>("d128", 34, -6100, 6100, 20000);
    return 0;
}

# The exact floor, ceiling and nearest value of sqrt(x) at p digits, from an integer root.
import sys
from math import isqrt
from decimal import Decimal, Context, ROUND_HALF_EVEN, getcontext
from collections import Counter, defaultdict

getcontext().prec = 200
getcontext().Emax = 999999
getcontext().Emin = -999999
DIGITS = {"d32": 7, "d64": 16, "d128": 34}
WANT = {"nearest": 2, "nearest_from_zero": 2, "toward_zero": 0, "downward": 0, "upward": 1}

def refs(x, p):
    t = x.as_tuple()
    s = int("".join(map(str, t.digits)))
    e = t.exponent
    f = Context(prec=p + 10, rounding=ROUND_HALF_EVEN).sqrt(x).adjusted() - (p - 1)
    while True:
        n = s * 10 ** (e - 2 * f)
        r = isqrt(n)
        if r >= 10 ** p:
            f += 1
        elif r < 10 ** (p - 1):
            f -= 1
        else:
            break
    exact = r * r == n
    lo = Decimal(r).scaleb(f)
    hi = lo if exact else Decimal(r + 1).scaleb(f)
    near = lo if (exact or 4 * n <= (2 * r + 1) ** 2) else hi
    return lo, hi, near, f

stats, total = defaultdict(Counter), Counter()
for line in open(sys.argv[1]):
    tag, mode, xs, rs = line.split()
    x, lib = Decimal(xs), Decimal(rs)
    if not x.is_finite() or x == 0:
        continue
    lo, hi, near, f = refs(x, DIGITS[tag])
    ref = (lo, hi, near)[WANT[mode]]
    total[(tag, mode)] += 1
    if lib != ref:
        stats[(tag, mode)][int((lib - ref).scaleb(-f))] += 1

for tag in DIGITS:
    for mode in WANT:
        c = stats[(tag, mode)]
        miss = sum(c.values())
        print(f"{tag:5} {mode:18} miss={miss:6}/{total[(tag, mode)]}  steps {min(c) if c else '-'} to {max(c) if c else '-'}")

Output on develop at 628546d2, GCC 16.2.1, -std=c++17 -O2. The counts under 20000 are
the finite inputs, because a random exponent at the top of the range overflows:

d32   nearest            miss= 11344/19062  steps -3 to 2
d32   nearest_from_zero  miss= 11344/19062  steps -3 to 2
d32   toward_zero        miss=  6684/20000  steps -3 to 1
d32   downward           miss=  6684/20000  steps -3 to 1
d32   upward             miss= 16222/19062  steps -3 to 2
d64   nearest            miss= 14825/19975  steps -5 to -1
d64   nearest_from_zero  miss= 14825/19975  steps -5 to -1
d64   toward_zero        miss=  9854/20000  steps -5 to -1
d64   downward           miss=  9854/20000  steps -5 to -1
d64   upward             miss= 19807/19975  steps -5 to -1
d128  nearest            miss=  7332/20000  steps -1 to 3
d128  nearest_from_zero  miss=  7332/20000  steps -1 to 3
d128  toward_zero        miss= 12270/20000  steps -2 to 3
d128  downward           miss= 12270/20000  steps -2 to 3
d128  upward             miss= 12310/20000  steps -2 to 3

Show the program of the library column

// Print sqrt(n) at full precision for n = 2..1000, per type.
#include <boost/decimal.hpp>
#include <iostream>
#include <iomanip>

using namespace boost::decimal;

template <typename T>
static void dump(const char* tag, int digits)
{
    fesetround(rounding_mode::fe_dec_to_nearest);
    for (int n = 2; n <= 1000; ++n)
    {
        const T r {sqrt(T{n, 0})};
        std::cout << tag << ' ' << n << ' ' << std::setprecision(digits) << r << '\n';
    }
}

int main()
{
    dump<decimal32_t>("d32", 7);
    dump<decimal64_t>("d64", 16);
    dump<decimal128_t>("d128", 34);
    dump<decimal_fast32_t>("f32", 7);
    dump<decimal_fast64_t>("f64", 16);
    dump<decimal_fast128_t>("f128", 34);
    return 0;
}

Show the program of the nearest column and the counts

# Compare the dump of the library against the nearest value, per type and parity.
import sys
from decimal import Decimal, Context, ROUND_HALF_EVEN
from collections import Counter

DIGITS = {"d32": 7, "d64": 16, "d128": 34, "f32": 7, "f64": 16, "f128": 34}
stats = {t: {"even": Counter(), "odd": Counter()} for t in DIGITS}

for line in open(sys.argv[1]):
    tag, ns, val = line.split()
    n = int(ns)
    prec = DIGITS[tag]
    ref = Context(prec=prec, rounding=ROUND_HALF_EVEN).sqrt(Decimal(n))
    lib = Decimal(val)
    step = Decimal(1).scaleb(ref.adjusted() - (prec - 1))
    ulps = int(((lib - ref) / step).to_integral_value())
    parity = "even" if Decimal(n).adjusted() % 2 == 0 else "odd"
    stats[tag][parity][ulps] += 1

for t in DIGITS:
    for parity in ("even", "odd"):
        cnt = stats[t][parity]
        miss = sum(cnt.values()) - cnt.get(0, 0)
        dist = " ".join(f"{k:+d}:{cnt[k]}" for k in sorted(cnt) if k != 0)
        print(f"{t:5} {parity:4} miss={miss:3}  steps: {dist}")

Output on develop at 628546d2, GCC 16.2.1, -std=c++17 -O2:

d32   even miss=470  steps: -1:470
d32   odd  miss= 61  steps: -3:2 -2:31 -1:19 +1:9
d64   even miss=441  steps: -1:441
d64   odd  miss= 88  steps: -4:8 -3:34 -2:24 -1:22
d128  even miss=  0  steps:
d128  odd  miss= 70  steps: -1:5 +1:33 +2:25 +3:7
f32   even miss=470  steps: -1:470
f32   odd  miss= 61  steps: -3:2 -2:31 -1:19 +1:9
f64   even miss=441  steps: -1:441
f64   odd  miss= 88  steps: -4:8 -3:34 -2:24 -1:22
f128  even miss=  0  steps:
f128  odd  miss= 70  steps: -1:5 +1:33 +2:25 +3:7

Environment

• Boost.Decimal at commit 628546d2 of develop
• GCC 16.2.1 and clang 22.1.8, with -std=c++17 -O2, on x86-64 Linux
• The nearest values come from the decimal module of Python 3.14.7, whose sqrt is
  correctly rounded with ROUND_HALF_EVEN
• The values of the directed modes come from an integer square root in Python, because
  the decimal module rounds sqrt half even in every mode

The cause

IEEE 754-2019 4.3: "Except where stated otherwise, every operation shall be performed as
if it first produced an intermediate result correct to infinite precision and with
unbounded range, and then rounded that result according to one of the attributes in this
clause." 4.3.1: in roundTiesToEven "the floating-point number nearest to the infinitely
precise result shall be delivered". 4.3.2: roundTowardPositive delivers the number
"closest to and no less than the infinitely precise result", and roundTowardNegative
the one "closest to and no greater than". 4.3.3: "The rounding-direction attribute affects
all computational operations that might be inexact." 5.4.1 lists squareRoot among the
arithmetic operations, next to the four basic operations and fusedMultiplyAdd.

Each kernel gets the significand gx in [1, 10), computes an integer root sig_z of
gx at the precision of the type, and rescales. There are four defects.

The 32 and 64 bit kernels only make sure that sig_z is a lower bound. They never move up
to the nearest integer. sqrt32_impl.hpp:86, and the same lines at sqrt64_impl.hpp:99:

    // ---------- Final rounding check ----------
    // Ensure z² ≤ gx (z is a lower bound)
    if (rem < 0)
    {
        --sig_z;
    }

For sqrt(decimal32_t{2}) the target is 2 * 10^12 and sig_z is 1414213. Its
square is 1999998409369, and the remainder is 1590631. The remainder is more than
sig_z, thus 1414214 is nearer. The kernel keeps 1414213. This is the source of each
-1 of the even exponents. The 128 bit kernel has the missing step at
sqrt128_impl.hpp:178:

        // Step 2: Round-to-nearest check
        // If (sig_z + 0.5)² < target, then sig_z+1 is closer
        // Equivalent: sig_z² + sig_z + 0.25 < target
        // Since we work with integers: if target - sig_z² > sig_z, round up
        u256 rem;
        i256_sub(target, sig_z_sq, rem);  // rem = target - sig_z², guaranteed non-negative

        if (rem > sig_z)
        {
            u256 one{1};
            sig_z = sig_z + one;
        }

The three kernels rescale an odd exponent with a multiply by the rounded constant
sqrt10_v<T>. sqrt32_impl.hpp:102, sqrt64_impl.hpp:119 and sqrt128_impl.hpp:212:

    if ((exp10val & 1) != 0)
    {
        z *= numbers::sqrt10_v<T>;
    }

The root of gx is rounded once, the constant is rounded, and the product is rounded a
third time. For sqrt(decimal32_t{11}) the root of 1.1 is 1.048808 after the
truncation. sqrt10_v<decimal32_t> is 3.162278. The product 3.316622 is three steps
under the nearest value 3.316625. This is the source of the misses of the odd exponents,
and of every miss of decimal128_t.

The fast path of a power of ten in sqrt.hpp:98 has the same double rounding for an odd
power under 1. sqrt(0.1) is T{1, 0} / sqrt10_v<T>, a rounded quotient of a rounded
constant, and it is one step off for the three types:

        if (p10_mod2 == 1)
        {
            result *= numbers::sqrt10_v<T>;
        }
        else if (p10_mod2 == -1)
        {
            result /= numbers::sqrt10_v<T>;
        }

No kernel reads the rounding mode. Whatever the mode is, the final step above picks the
lower bound in the 32 and 64 bit kernels. It picks the nearest integer in the 128 bit
kernel. The only step which reads the mode is the multiply by sqrt10_v<T>, which rounds
a value that is wrong already. In the upward mode the result is under the ceiling 85% to
99% of the time for the 32 and 64 bit types.

A fourth defect hides behind the first. The estimate of the 32 bit kernel comes from
approx_recip_sqrt32. That function works at the scale 10^5 and multiplies by 100,
thus the estimate has five digits. The Newton step truncates its correction toward zero.
If the root ends near .999, one step can then land two above the floor. The single step
down of the old final step stops one above the floor. This is the source of the +1 of
926 random values of decimal32_t.

The correction

For an odd exponent, the kernels put the factor 10 into the integer before the root. The
multiply by sqrt10_v<T> goes away. A helper reads the mode, and the last step of each
kernel keeps the floor or steps up as the mode asks. The change is +92/-58 lines in five
headers. After it, 20000 random values per type give the value of the mode every time, in
the five modes. An even exponent costs 5 to 25 ns more, and an odd exponent costs 10 to
220 ns less.

A helper in detail/fenv_rounding.hpp, next to overflow_is_finite of #1460, reads the
mode and decides the last step from two facts about the exact remainder:

// IEEE 754-2019 4.3: the floor of an inexact root steps up in the upward mode, past the half
// point in the nearest modes, and never in the modes toward zero. A root is never a tie.
BOOST_DECIMAL_CUDA_CONSTEXPR auto sqrt_steps_up(const bool inexact, const bool past_half) noexcept -> bool
{
    auto round {_boost_decimal_global_rounding_mode};
    #ifndef BOOST_DECIMAL_NO_CONSTEVAL_DETECTION
    if (!BOOST_DECIMAL_IS_CONSTANT_EVALUATED(inexact))
    {
        round = _boost_decimal_global_runtime_rounding_mode;
    }
    #endif
    if (round == rounding_mode::fe_dec_upward)
    {
        return inexact;
    }
    if (round == rounding_mode::fe_dec_toward_zero || round == rounding_mode::fe_dec_downward)
    {
        return false;
    }
    return past_half;
}

The final step of each kernel goes down to the floor in a loop. One Newton step lands at
or above the floor, and at most a few steps above it. Then the helper decides the step up.
Here sqrt32_impl.hpp:

     // ---------- Final rounding check ----------
-    // Ensure z² ≤ gx (z is a lower bound)
-    if (rem < 0)
+    // Ensure z² ≤ gx (z is a lower bound); Newton can land a few steps above it
+    while (rem < 0)
     {
         --sig_z;
+        rem += 2 * static_cast<std::int64_t>(sig_z) + 1;
+    }
+    // Round up if the mode asks for it: target - sig_z² > sig_z means (sig_z + 0.5)² < target
+    if (sqrt_steps_up(rem != 0, rem > static_cast<std::int64_t>(sig_z)))
+    {
+        ++sig_z;
     }

For an odd exponent the three kernels multiply the target by 10 before the integer root,
in place of the multiply by sqrt10_v<T> after it. The root is then sqrt(10 * gx),
which is rounded once. The estimate of the root is scaled by an integer sqrt(10), and
the Newton steps correct it. The widths do not change: the target is less than 10^14,
10^32 and 10^68 for the three kernels. No step uses more digits than the type has.

-    // Precompute target = sig_gx * 10^6 (avoids recomputing in Newton and rounding)
-    const std::uint64_t target = static_cast<std::uint64_t>(sig_gx) * scale6;
+    // If exp is odd, take √(10 × gx) instead so the result is rounded once (no ×√10 later)
+    // sig_z ≈ √gx × 10^6 × √10, then Newton corrects it
+    const bool odd = (exp10val & 1) != 0;
+    if (odd)
+    {
+        sig_z = static_cast<std::uint32_t>(static_cast<std::uint64_t>(sig_z) * 31622777ULL / scale7);
+    }
+
+    // Precompute target = sig_gx * 10^6, ×10 if exp was odd (avoids recomputing in Newton and rounding)
+    const std::uint64_t target = static_cast<std::uint64_t>(sig_gx) * (odd ? scale7 : scale6);

     if (half_exp != 0)
     {
         z *= T{1, half_exp};
     }
-    if ((exp10val & 1) != 0)
-    {
-        z *= numbers::sqrt10_v<T>;
-    }

The 64 and 128 bit kernels get the same three changes with their own widths and scales.
The 128 bit kernel already has the loop down to the floor.

The fast path of sqrt.hpp keeps an even power of ten, which has an exact root. An odd
power now takes the kernel with gx at 1, thus sqrt(10) is rounded once, in the mode:

-        if (p10 == 0)
+        if ((p10 % 2) == 0)
         {
-            return T{1};
+            return T{1, p10 / 2};
         }
-
-        const int p10_mod2 = (p10 % 2);
-        T result = T{1, p10 / 2};
-
-        if (p10_mod2 == 1)
-        {
-            result *= numbers::sqrt10_v<T>;
-        }
-        else if (p10_mod2 == -1)
-        {
-            result /= numbers::sqrt10_v<T>;
-        }
-        return result;

A new test compares sqrt of 29 integers and three random values against the floor and
the ceiling of the root. It covers the six types in the five modes, with two
static_assert for the compile time path. A sibling test does the same under the
compile-time upward mode without constant evaluation detection. Before the correction the
tests fail 482 of 900 and 112 of 180 checks, and the two static_assert fail.
test_sqrt.cpp expected the wrong sqrt(0.1). Its control is now sqrt10_v<T> times an
exact power of ten for both signs of the power. After the correction, the three tests and
178 other tests of the suite pass with GCC 16.2.1. The build uses C++17, the warning set
of the Jamfile and -Werror. The new tests also pass with GCC in C++14 and C++20, and
with clang 22.1.8 in C++14, C++17 and C++20. The sweep of 20000 random values and 83
powers of ten per type gives no miss in any mode.

What the correction costs

ns per sqrt on 20000 random full width operands, best of 10, one core, GCC 16.2.1
-O2, in the default mode. The even column gets the mode read and one more remainder
compare. The odd column loses a decimal multiply, and the integer scale of the estimate is
cheaper.

type	even exponent, before / after	odd exponent, before / after
decimal32_t	80 / 88	104 / 102
decimal_fast32_t	81 / 84	101 / 82
decimal64_t	167 / 181	192 / 197
decimal_fast64_t	194 / 218	225 / 235
decimal128_t	776 / 798	893 / 825
decimal_fast128_t	760 / 775	877 / 804

The same measurement for the three IEEE types in the five runtime modes. Before, the
multiply by sqrt10_v<T> was the one step which read the mode, and it made the upward
mode the slowest. After, the mode read is the same in every mode:

type	mode	even, before / after	odd, before / after
decimal32_t	nearest	81 / 86	111 / 91
decimal32_t	nearest from zero	85 / 91	115 / 96
decimal32_t	toward zero	84 / 89	104 / 93
decimal32_t	downward	85 / 88	107 / 92
decimal32_t	upward	88 / 92	111 / 97
decimal64_t	nearest	158 / 184	177 / 197
decimal64_t	nearest from zero	171 / 189	204 / 200
decimal64_t	toward zero	162 / 177	188 / 189
decimal64_t	downward	163 / 174	189 / 185
decimal64_t	upward	173 / 178	200 / 189
decimal128_t	nearest	796 / 796	912 / 822
decimal128_t	nearest from zero	818 / 823	998 / 848
decimal128_t	toward zero	766 / 753	891 / 775
decimal128_t	downward	767 / 753	893 / 776
decimal128_t	upward	863 / 842	1081 / 865

Why the tests do not find it

test_sqrt.cpp compares the decimal result against the binary std::sqrt of a float,
a double or a long double. The random tests use a tolerance of 16 epsilon, and 64 for
the 128 bit types. The dense scan of [1.01, 9.99] uses 32 epsilon, and 1024 for
decimal128_t. A miss of one to four decimal steps is inside these tolerances. The test
of sqrt(2) and sqrt(5) only makes sure that the result is between 1 and 2, and between
2 and 3. The test of the powers of ten expects T{1, np / 2} / sqrt10_v<T> for an odd
negative power, which is the wrong value that the fast path computed. No test compares
sqrt against a decimal reference value, and no test runs sqrt in a directed mode.