git: c3f6dcea1992 - main - msun: Fix up for recent rsqrt[fl] functions

The branch main has been updated by fuz:

URL: https://cgit.FreeBSD.org/src/commit/?id=c3f6dcea199289329c1d3b91b69e5a4821fc3dff

commit c3f6dcea199289329c1d3b91b69e5a4821fc3dff
Author:     Steve Kargl <kargl@FreeBSD.org>
AuthorDate: 2026-06-07 19:12:16 +0000
Commit:     Robert Clausecker <fuz@FreeBSD.org>
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();