crt: arm: Use bit twiddling implementations of remquo* This avoids issues when the quotient exceeds the range for integer variables. Signed-off-by: Martin Storsjö <martin@martin.st>
diff --git a/mingw-w64-crt/Makefile.am b/mingw-w64-crt/Makefile.am index e07feba..d499e66 100644 --- a/mingw-w64-crt/Makefile.am +++ b/mingw-w64-crt/Makefile.am
@@ -303,9 +303,9 @@ math/arm-common/remainder.c \ math/arm-common/remainderf.c \ math/arm-common/remainderl.c \ - math/arm-common/remquo.c \ - math/arm-common/remquof.c \ math/arm-common/remquol.c \ + math/arm-common/s_remquo.c \ + math/arm-common/s_remquof.c \ math/arm-common/scalbn.c endif @@ -340,9 +340,9 @@ math/arm-common/remainder.c \ math/arm-common/remainderf.c \ math/arm-common/remainderl.c \ - math/arm-common/remquo.c \ - math/arm-common/remquof.c \ math/arm-common/remquol.c \ + math/arm-common/s_remquo.c \ + math/arm-common/s_remquof.c \ math/arm-common/scalbn.c \ math/arm64/exp2.S \ math/arm64/exp2f.S \
diff --git a/mingw-w64-crt/math/arm-common/remquo.c b/mingw-w64-crt/math/arm-common/remquo.c deleted file mode 100644 index d8ec534..0000000 --- a/mingw-w64-crt/math/arm-common/remquo.c +++ /dev/null
@@ -1,35 +0,0 @@ -/** - * This file has no copyright assigned and is placed in the Public Domain. - * This file is part of the mingw-w64 runtime package. - * No warranty is given; refer to the file DISCLAIMER.PD within this package. - */ - -#include <math.h> -#include <errno.h> - -double remquo(double x, double y, int *quo) -{ - if (isnan(x)) - return x; - if (isnan(y)) - return y; - if (isinf(x) || y == 0) - return NAN; - double div = x/y; - double integral; - double frac = modf(div, &integral); - int iintegral = (int)integral; - if (frac == 0.5) { - if (iintegral & 1) - *quo = iintegral + 1; - else - *quo = iintegral; - } else if (frac == -0.5) { - if (iintegral & 1) - *quo = iintegral - 1; - else - *quo = iintegral; - } else - *quo = round(div); - return x - *quo * y; -}
diff --git a/mingw-w64-crt/math/arm-common/remquof.c b/mingw-w64-crt/math/arm-common/remquof.c deleted file mode 100644 index d928002..0000000 --- a/mingw-w64-crt/math/arm-common/remquof.c +++ /dev/null
@@ -1,35 +0,0 @@ -/** - * This file has no copyright assigned and is placed in the Public Domain. - * This file is part of the mingw-w64 runtime package. - * No warranty is given; refer to the file DISCLAIMER.PD within this package. - */ - -#include <math.h> -#include <errno.h> - -float remquof(float x, float y, int *quo) -{ - if (isnan(x)) - return x; - if (isnan(y)) - return y; - if (isinf(x) || y == 0) - return NAN; - float div = x/y; - float integral; - float frac = modff(div, &integral); - int iintegral = (int)integral; - if (frac == 0.5) { - if (iintegral & 1) - *quo = iintegral + 1; - else - *quo = iintegral; - } else if (frac == -0.5) { - if (iintegral & 1) - *quo = iintegral - 1; - else - *quo = iintegral; - } else - *quo = roundf(div); - return x - *quo * y; -}
diff --git a/mingw-w64-crt/math/arm-common/s_remquo.c b/mingw-w64-crt/math/arm-common/s_remquo.c new file mode 100644 index 0000000..48b7317 --- /dev/null +++ b/mingw-w64-crt/math/arm-common/s_remquo.c
@@ -0,0 +1,154 @@ +/* @(#)e_fmod.c 1.3 95/01/18 */ +/*- + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, 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> + +#include <float.h> + +#include <math.h> +#include "../bsd_private_base.h" + +static const double Zero[] = {0.0, -0.0,}; + +/* + * Return the IEEE remainder and set *quo to the last n bits of the + * quotient, rounded to the nearest integer. We choose n=31 because + * we wind up computing all the integer bits of the quotient anyway as + * a side-effect of computing the remainder by the shift and subtract + * method. In practice, this is far more bits than are needed to use + * remquo in reduction algorithms. + */ +double +remquo(double x, double y, int *quo) +{ + int32_t n,hx,hy,hz,ix,iy,sx,i; + u_int32_t lx,ly,lz,q,sxy; + + EXTRACT_WORDS(hx,lx,x); + EXTRACT_WORDS(hy,ly,y); + sxy = (hx ^ hy) & 0x80000000; + sx = hx&0x80000000; /* sign of x */ + hx ^=sx; /* |x| */ + hy &= 0x7fffffff; /* |y| */ + + /* purge off exception values */ + if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */ + ((hy|((ly|-ly)>>31))>0x7ff00000)) /* or y is NaN */ + return nan_mix_op(x, y, *)/nan_mix_op(x, y, *); + if(hx<=hy) { + if((hx<hy)||(lx<ly)) { + q = 0; + goto fixup; /* |x|<|y| return x or x-y */ + } + if(lx==ly) { + *quo = (sxy ? -1 : 1); + return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/ + } + } + + /* determine ix = ilogb(x) */ + if(hx<0x00100000) { /* subnormal x */ + if(hx==0) { + for (ix = -1043, i=lx; i>0; i<<=1) ix -=1; + } else { + for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1; + } + } else ix = (hx>>20)-1023; + + /* determine iy = ilogb(y) */ + if(hy<0x00100000) { /* subnormal y */ + if(hy==0) { + for (iy = -1043, i=ly; i>0; i<<=1) iy -=1; + } else { + for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1; + } + } else iy = (hy>>20)-1023; + + /* set up {hx,lx}, {hy,ly} and align y to x */ + if(ix >= -1022) + hx = 0x00100000|(0x000fffff&hx); + else { /* subnormal x, shift x to normal */ + n = -1022-ix; + if(n<=31) { + hx = (hx<<n)|(lx>>(32-n)); + lx <<= n; + } else { + hx = lx<<(n-32); + lx = 0; + } + } + if(iy >= -1022) + hy = 0x00100000|(0x000fffff&hy); + else { /* subnormal y, shift y to normal */ + n = -1022-iy; + if(n<=31) { + hy = (hy<<n)|(ly>>(32-n)); + ly <<= n; + } else { + hy = ly<<(n-32); + ly = 0; + } + } + + /* fix point fmod */ + n = ix - iy; + q = 0; + while(n--) { + hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; + if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;} + else {hx = hz+hz+(lz>>31); lx = lz+lz; q++;} + q <<= 1; + } + hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; + if(hz>=0) {hx=hz;lx=lz;q++;} + + /* convert back to floating value and restore the sign */ + if((hx|lx)==0) { /* return sign(x)*0 */ + q &= 0x7fffffff; + *quo = (sxy ? -q : q); + return Zero[(u_int32_t)sx>>31]; + } + while(hx<0x00100000) { /* normalize x */ + hx = hx+hx+(lx>>31); lx = lx+lx; + iy -= 1; + } + if(iy>= -1022) { /* normalize output */ + hx = ((hx-0x00100000)|((iy+1023)<<20)); + } else { /* subnormal output */ + n = -1022 - iy; + if(n<=20) { + lx = (lx>>n)|((u_int32_t)hx<<(32-n)); + hx >>= n; + } else if (n<=31) { + lx = (hx<<(32-n))|(lx>>n); hx = 0; + } else { + lx = hx>>(n-32); hx = 0; + } + } +fixup: + INSERT_WORDS(x,hx,lx); + y = fabs(y); + if (y < 0x1p-1021) { + if (x+x>y || (x+x==y && (q & 1))) { + q++; + x-=y; + } + } else if (x>0.5*y || (x==0.5*y && (q & 1))) { + q++; + x-=y; + } + GET_HIGH_WORD(hx,x); + SET_HIGH_WORD(x,hx^sx); + q &= 0x7fffffff; + *quo = (sxy ? -q : q); + return x; +}
diff --git a/mingw-w64-crt/math/arm-common/s_remquof.c b/mingw-w64-crt/math/arm-common/s_remquof.c new file mode 100644 index 0000000..f2a4e0d --- /dev/null +++ b/mingw-w64-crt/math/arm-common/s_remquof.c
@@ -0,0 +1,121 @@ +/* @(#)e_fmod.c 1.3 95/01/18 */ +/*- + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, 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> + +#include <math.h> +#include "../bsd_private_base.h" + +static const float Zero[] = {0.0, -0.0,}; + +/* + * Return the IEEE remainder and set *quo to the last n bits of the + * quotient, rounded to the nearest integer. We choose n=31 because + * we wind up computing all the integer bits of the quotient anyway as + * a side-effect of computing the remainder by the shift and subtract + * method. In practice, this is far more bits than are needed to use + * remquo in reduction algorithms. + */ +float +remquof(float x, float y, int *quo) +{ + int32_t n,hx,hy,hz,ix,iy,sx,i; + u_int32_t q,sxy; + + GET_FLOAT_WORD(hx,x); + GET_FLOAT_WORD(hy,y); + sxy = (hx ^ hy) & 0x80000000; + sx = hx&0x80000000; /* sign of x */ + hx ^=sx; /* |x| */ + hy &= 0x7fffffff; /* |y| */ + + /* purge off exception values */ + if(hy==0||hx>=0x7f800000||hy>0x7f800000) /* y=0,NaN;or x not finite */ + return nan_mix_op(x, y, *)/nan_mix_op(x, y, *); + if(hx<hy) { + q = 0; + goto fixup; /* |x|<|y| return x or x-y */ + } else if(hx==hy) { + *quo = (sxy ? -1 : 1); + return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/ + } + + /* determine ix = ilogb(x) */ + if(hx<0x00800000) { /* subnormal x */ + for (ix = -126,i=(hx<<8); i>0; i<<=1) ix -=1; + } else ix = (hx>>23)-127; + + /* determine iy = ilogb(y) */ + if(hy<0x00800000) { /* subnormal y */ + for (iy = -126,i=(hy<<8); i>0; i<<=1) iy -=1; + } else iy = (hy>>23)-127; + + /* set up {hx,lx}, {hy,ly} and align y to x */ + if(ix >= -126) + hx = 0x00800000|(0x007fffff&hx); + else { /* subnormal x, shift x to normal */ + n = -126-ix; + hx <<= n; + } + if(iy >= -126) + hy = 0x00800000|(0x007fffff&hy); + else { /* subnormal y, shift y to normal */ + n = -126-iy; + hy <<= n; + } + + /* fix point fmod */ + n = ix - iy; + q = 0; + while(n--) { + hz=hx-hy; + if(hz<0) hx = hx << 1; + else {hx = hz << 1; q++;} + q <<= 1; + } + hz=hx-hy; + if(hz>=0) {hx=hz;q++;} + + /* convert back to floating value and restore the sign */ + if(hx==0) { /* return sign(x)*0 */ + q &= 0x7fffffff; + *quo = (sxy ? -q : q); + return Zero[(u_int32_t)sx>>31]; + } + while(hx<0x00800000) { /* normalize x */ + hx <<= 1; + iy -= 1; + } + if(iy>= -126) { /* normalize output */ + hx = ((hx-0x00800000)|((iy+127)<<23)); + } else { /* subnormal output */ + n = -126 - iy; + hx >>= n; + } +fixup: + SET_FLOAT_WORD(x,hx); + y = fabsf(y); + if (y < 0x1p-125f) { + if (x+x>y || (x+x==y && (q & 1))) { + q++; + x-=y; + } + } else if (x>0.5f*y || (x==0.5f*y && (q & 1))) { + q++; + x-=y; + } + GET_FLOAT_WORD(hx,x); + SET_FLOAT_WORD(x,hx^sx); + q &= 0x7fffffff; + *quo = (sxy ? -q : q); + return x; +}
diff --git a/mingw-w64-crt/math/bsd_private_base.h b/mingw-w64-crt/math/bsd_private_base.h index de5b935..598b815 100644 --- a/mingw-w64-crt/math/bsd_private_base.h +++ b/mingw-w64-crt/math/bsd_private_base.h
@@ -121,3 +121,28 @@ } while (0) #endif #endif /* FLT_EVAL_METHOD */ + +/* + * Mix 0, 1 or 2 NaNs. First add 0 to each arg. This normally just turns + * signaling NaNs into quiet NaNs by setting a quiet bit. We do this + * because we want to never return a signaling NaN, and also because we + * don't want the quiet bit to affect the result. Then mix the converted + * args using the specified operation. + * + * When one arg is NaN, the result is typically that arg quieted. When both + * args are NaNs, the result is typically the quietening of the arg whose + * mantissa is largest after quietening. When neither arg is NaN, the + * result may be NaN because it is indeterminate, or finite for subsequent + * construction of a NaN as the indeterminate 0.0L/0.0L. + * + * Technical complications: the result in bits after rounding to the final + * precision might depend on the runtime precision and/or on compiler + * optimizations, especially when different register sets are used for + * different precisions. Try to make the result not depend on at least the + * runtime precision by always doing the main mixing step in long double + * precision. Try to reduce dependencies on optimizations by adding the + * the 0's in different precisions (unless everything is in long double + * precision). + */ +#define nan_mix(x, y) (nan_mix_op((x), (y), +)) +#define nan_mix_op(x, y, op) (((x) + 0.0L) op ((y) + 0))