a2ddfa5ea7
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().
66 lines
1.4 KiB
C
66 lines
1.4 KiB
C
/* e_remainderf.c -- float version of e_remainder.c.
|
|
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
|
*/
|
|
|
|
/*
|
|
* ====================================================
|
|
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
|
*
|
|
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
|
* Permission to use, copy, modify, and distribute this
|
|
* software is freely granted, provided that this notice
|
|
* is preserved.
|
|
* ====================================================
|
|
*/
|
|
|
|
#include <sys/cdefs.h>
|
|
__FBSDID("$FreeBSD$");
|
|
|
|
#include "math.h"
|
|
#include "math_private.h"
|
|
|
|
static const float zero = 0.0;
|
|
|
|
|
|
float
|
|
__ieee754_remainderf(float x, float p)
|
|
{
|
|
int32_t hx,hp;
|
|
u_int32_t sx;
|
|
float p_half;
|
|
|
|
GET_FLOAT_WORD(hx,x);
|
|
GET_FLOAT_WORD(hp,p);
|
|
sx = hx&0x80000000;
|
|
hp &= 0x7fffffff;
|
|
hx &= 0x7fffffff;
|
|
|
|
/* purge off exception values */
|
|
if(hp==0) return (x*p)/(x*p); /* p = 0 */
|
|
if((hx>=0x7f800000)|| /* x not finite */
|
|
((hp>0x7f800000))) /* p is NaN */
|
|
return ((long double)x*p)/((long double)x*p);
|
|
|
|
|
|
if (hp<=0x7effffff) x = __ieee754_fmodf(x,p+p); /* now x < 2p */
|
|
if ((hx-hp)==0) return zero*x;
|
|
x = fabsf(x);
|
|
p = fabsf(p);
|
|
if (hp<0x01000000) {
|
|
if(x+x>p) {
|
|
x-=p;
|
|
if(x+x>=p) x -= p;
|
|
}
|
|
} else {
|
|
p_half = (float)0.5*p;
|
|
if(x>p_half) {
|
|
x-=p;
|
|
if(x>=p_half) x -= p;
|
|
}
|
|
}
|
|
GET_FLOAT_WORD(hx,x);
|
|
if ((hx&0x7fffffff)==0) hx = 0;
|
|
SET_FLOAT_WORD(x,hx^sx);
|
|
return x;
|
|
}
|