The branch main has been updated by fuz: URL: https://cgit.FreeBSD.org/src/commit/?id=c3f6dcea199289329c1d3b91b69e5a4821fc3dff
commit c3f6dcea199289329c1d3b91b69e5a4821fc3dff Author: Steve Kargl <[email protected]> AuthorDate: 2026-06-07 19:12:16 +0000 Commit: Robert Clausecker <[email protected]> CommitDate: 2026-06-07 20:59:19 +0000 msun: Fix up for recent rsqrt[fl] functions Paul Zimmermann (of Core-Math and MPFR fame) graciously tested the recently committed rsqrt[fl]() functions. He identified 127 incorrectly rounded values for rsqrtf() in round-to-nearest mode. This patch fixes the rounding in RN. Exhaustive testing now shows that rsqrtf() is corrected rounded for RN. He also tested rsqrt() and rsqrtl() in the interval [1,4). Both appear to be correctly rounded. Finally, the patch includes small changes to comments. A concise list of changes is * lib/msun/src/s_rsqrt.c: . Fix comments. * lib/msun/src/s_rsqrtf.c . Fix comments. . Exhaustive testing by Paul Zimmermann found 127 incorrectly rounded values in round-to-nearests. These gave have the form 0x1.13e07pN with N an odd integer. With this patch, all values are now correctly rounded in round-to-nearest. * lib/msun/src/s_rsqrtl.c . Fix comments. . Move all variable declarations to top of function and sort. PR: 295706 MFC after: 1 week Fixes: 3085fc9d97bd83785ba3ba43e0378d7d67987d1f --- lib/msun/src/s_rsqrt.c | 4 ++-- lib/msun/src/s_rsqrtf.c | 37 ++++++++++++++++++++++--------------- lib/msun/src/s_rsqrtl.c | 8 +++----- 3 files changed, 27 insertions(+), 22 deletions(-) diff --git a/lib/msun/src/s_rsqrt.c b/lib/msun/src/s_rsqrt.c index 0a513afe8ed2..cac16f837a05 100644 --- a/lib/msun/src/s_rsqrt.c +++ b/lib/msun/src/s_rsqrt.c @@ -108,14 +108,14 @@ rsqrt(double x) hx = (hx & 0x000fffff) | 0x3fe00000; SET_HIGH_WORD(x, hx); - /* m is odd. Put x into [0.25,5) and increase m. */ + /* m is odd. Put x into [0.25,0.5) and increase m. */ if (m & 1) { x /= 2; m += 1; } m = -(m >> 1); /* Prepare for 2^(-m/2). */ - y = 1 / sqrt(x); /* ~52-bit estimate. */ + y = 1 / sqrt(x); /* ~52-bit estimate. */ h = y / 2; diff --git a/lib/msun/src/s_rsqrtf.c b/lib/msun/src/s_rsqrtf.c index b71f7baf5657..faeb1c300e18 100644 --- a/lib/msun/src/s_rsqrtf.c +++ b/lib/msun/src/s_rsqrtf.c @@ -108,7 +108,7 @@ rsqrtf(float x) ix = (ix & 0x007fffff) | 0x3f000000; SET_FLOAT_WORD(x, ix); /* x is in [0.5,1). */ - /* m is odd. Put x into [0.25,5) and increase m. */ + /* m is odd. Put x into [0.25,0.5) and increase m. */ if (m & 1) { x /= 2; m += 1; @@ -122,33 +122,40 @@ rsqrtf(float x) * a max ulp of ~1.49. */ y = 1 / sqrtf(x); - h = y / 2; /* - * For values of x with a representation of 0x1.fffffcpN with - * N an odd integer, the computed rsqrtf() is not correctly - * rounded in round-to-nearest without toggling the rounding - * mode to FE_TOWARDZERO. Note, FE_DOWNWARD also works. - * However, messing with the rounding mode is expensive, so - * only do it when necessary. Example, x = 0x1.fffffcp3 gives - * y --> 0x3f800001. + * For values of x with representations of the forms 0x1.fffffcpN + * or 0x1.13e07pN with N an odd integer, the computed rsqrtf() is + * not correctly rounded in round-to-nearest without fiddling with + * the rounding mode and/or bits of x. */ GET_FLOAT_WORD(ix, y); - if ((ix & 0x000fffff) == 1) { + if ((ix & 0x000fffff) == 1) { /* 0x1.fffffcpN */ + rnd = fegetround(); + fesetround(FE_DOWNWARD); + _MUL(x, y, zh, zl); + _XMUL(zh, zl, h, 0, ph, pl); + fesetround(rnd); + _XADD(-ph, -pl, half, 0, rh, rl); + y = h + h * rh; + } else if((ix & 0x00ffffff) == 0x00ae6055) { /* 0x1.13e07pN */ + GET_FLOAT_WORD(ix, x); + SET_FLOAT_WORD(x, ix - 1); rnd = fegetround(); - fesetround(FE_TOWARDZERO); + fesetround(FE_DOWNWARD); _MUL(x, y, zh, zl); _XMUL(zh, zl, h, 0, ph, pl); - fesetround(rnd); + _XADD(-ph, -pl, half, 0, rh, rl); + y = h + h * rh; + fesetround(rnd); } else { _MUL(x, y, zh, zl); _XMUL(zh, zl, h, 0, ph, pl); + _XADD(-ph, -pl, half, 0, rh, rl); + y = h * rh + h; } - _XADD(-ph, -pl, half, 0, rh, rl); - y = h * rh + h; - ix = (uint32_t)(m + 128) << 23; SET_FLOAT_WORD(x, ix); return (y *= x); diff --git a/lib/msun/src/s_rsqrtl.c b/lib/msun/src/s_rsqrtl.c index 34d5cdc48173..bf3d32481e3f 100644 --- a/lib/msun/src/s_rsqrtl.c +++ b/lib/msun/src/s_rsqrtl.c @@ -108,7 +108,7 @@ rsqrtl(long double x) } u.bits.exp = 0x3ffe; - /* m is odd. Put x into [0.25,5) and increase m. */ + /* m is odd. Put x into [0.25,0.5) and increase m. */ if (m & 1) { u.e /= 2; m += 1; @@ -141,8 +141,9 @@ long double rsqrtl(long double x) { volatile static const double vzero = 0; + static const double half = 0.5; int hx, m, rnd; - long double y; + long double h, ph, pl, rh, rl, y, zh, zl; /* x = +-0. Raise exception. */ if (x == 0) @@ -183,9 +184,6 @@ rsqrtl(long double x) y = 1 / sqrt((double)x); /* ~52-bit estimate. */ y -= y * (x * y * y - 1) / 2; /* ~104-bit estimate. */ - static const double half = 0.5; - long double h, ph, pl, rh, rl, zh, zl; - h = y / 2; rnd = fegetround();
