Fix remainder() and remainderf() in round-towards-minus-infinity mode

when the result is +-0.  IEEE754 requires (in all rounding modes) that
if the result is +-0 then its sign is the same as that of the first
arg, but in round-towards-minus-infinity mode an uncorrected implementation
detail always reversed the sign.  (The detail is that x-x with x's
sign positive gives -0 in this mode only, but the algorithm assumed
that x-x always has positive sign for such x.)

remquo() and remquof() seem to need the same fix, but I cannot test them
yet.

Use long doubles when mixing NaN args.  This trick improves consistency
of results on at least amd64, so that more serious problems like the
above aren't hidden in simple regression tests by noise for the NaNs.
On amd64, hardware remainder should be used since it is about 10 times
faster than software remainder and is already used for remquo(), but
it involves using the i387 even for floats and doubles, and the i387
does NaN mixing which is better than but inconsistent with SSE NaN mixing.
Software remainder() would probably have been inconsistent with
software remainderl() for the same reason if the latter existed.

Signaling NaNs cause further inconsistencies on at least ia64 and i386.

Use __FBSDID().

2 files changed, 8 insertions, 8 deletions

diff --git a/lib/msun/src/e_remainder.c b/lib/msun/src/e_remainder.c
index 1cea6c1..2d285d6 100644
--- a/lib/msun/src/e_remainder.c
+++ b/lib/msun/src/e_remainder.c
@@ -11,9 +11,8 @@
  * ====================================================
  */
 
-#ifndef lint
-static char rcsid[] = "$FreeBSD$";
-#endif
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
 
 /* __ieee754_remainder(x,p)
  * Return :                  
@@ -48,7 +47,7 @@ __ieee754_remainder(double x, double p)
 	if((hx>=0x7ff00000)||			/* x not finite */
 	  ((hp>=0x7ff00000)&&			/* p is NaN */
 	  (((hp-0x7ff00000)|lp)!=0)))
-	    return (x*p)/(x*p);
+	    return ((long double)x*p)/((long double)x*p);
 
 
 	if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p);	/* now x < 2p */
@@ -68,6 +67,7 @@ __ieee754_remainder(double x, double p)
 	    }
 	}
 	GET_HIGH_WORD(hx,x);
+	if ((hx&0x7fffffff)==0) hx = 0;
 	SET_HIGH_WORD(x,hx^sx);
 	return x;
 }
diff --git a/lib/msun/src/e_remainderf.c b/lib/msun/src/e_remainderf.c
index 7841bdc..b0014ae 100644
--- a/lib/msun/src/e_remainderf.c
+++ b/lib/msun/src/e_remainderf.c
@@ -13,9 +13,8 @@
  * ====================================================
  */
 
-#ifndef lint
-static char rcsid[] = "$FreeBSD$";
-#endif
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
 
 #include "math.h"
 #include "math_private.h"
@@ -40,7 +39,7 @@ __ieee754_remainderf(float x, float p)
 	if(hp==0) return (x*p)/(x*p);	 	/* p = 0 */
 	if((hx>=0x7f800000)||			/* x not finite */
 	  ((hp>0x7f800000)))			/* p is NaN */
-	    return (x*p)/(x*p);
+	    return ((long double)x*p)/((long double)x*p);
 
 
 	if (hp<=0x7effffff) x = __ieee754_fmodf(x,p+p);	/* now x < 2p */
@@ -60,6 +59,7 @@ __ieee754_remainderf(float x, float p)
 	    }
 	}
 	GET_FLOAT_WORD(hx,x);
+	if ((hx&0x7fffffff)==0) hx = 0;
 	SET_FLOAT_WORD(x,hx^sx);
 	return x;
 }