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().
This commit is contained in:
parent
54f8dbc48f
commit
a2ddfa5ea7
Notes:
svn2git
2020-12-20 02:59:44 +00:00
svn path=/head/; revision=176207
@ -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;
|
||||
}
|
||||
|
@ -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;
|
||||
}
|
||||
|
Loading…
Reference in New Issue
Block a user