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:		3085fc9d97

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;
 

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_TOWARDZERO);
+	    fesetround(FE_DOWNWARD);
 	    _MUL(x, y, zh, zl);
 	    _XMUL(zh, zl, h, 0, ph, pl);
-	   fesetround(rnd);
+	    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_DOWNWARD);
+	    _MUL(x, y, zh, zl);
+	    _XMUL(zh, zl, h, 0, ph, pl);
+            _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);

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();