diff options
Diffstat (limited to 'src/math/remquof.c')
| -rw-r--r-- | src/math/remquof.c | 166 |
1 files changed, 61 insertions, 105 deletions
diff --git a/src/math/remquof.c b/src/math/remquof.c index 536a050a..2f41ff70 100644 --- a/src/math/remquof.c +++ b/src/math/remquof.c @@ -1,126 +1,82 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_remquof.c */ -/*- - * ==================================================== - * 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. - * ==================================================== - */ -/* - * 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. - */ - -#include "libm.h" - -static const float Zero[] = {0.0, -0.0,}; +#include <math.h> +#include <stdint.h> float remquof(float x, float y, int *quo) { - int32_t n,hx,hy,hz,ix,iy,sx,i; - uint32_t q,sxy; + union {float f; uint32_t i;} ux = {x}, uy = {y}; + int ex = ux.i>>23 & 0xff; + int ey = uy.i>>23 & 0xff; + int sx = ux.i>>31; + int sy = uy.i>>31; + uint32_t q; + uint32_t i; + uint32_t uxi = ux.i; - 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 */ + *quo = 0; + if (uy.i<<1 == 0 || isnan(y) || ex == 0xff) return (x*y)/(x*y); - if (hx < hy) { /* |x| < |y| return x or x-y */ - q = 0; - goto fixup; - } else if(hx==hy) { /* |x| = |y| return x*0*/ - *quo = sxy ? -1 : 1; - return Zero[(uint32_t)sx>>31]; - } + if (ux.i<<1 == 0) + return x; - /* determine ix = ilogb(x) */ - if (hx < 0x00800000) { /* subnormal x */ - for (ix = -126, i=hx<<8; i>0; i<<=1) ix--; - } 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--; - } 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; + /* normalize x and y */ + if (!ex) { + for (i = uxi<<9; i>>31 == 0; ex--, i <<= 1); + uxi <<= -ex + 1; + } else { + uxi &= -1U >> 9; + uxi |= 1U << 23; } - if (iy >= -126) - hy = 0x00800000|(0x007fffff&hy); - else { /* subnormal y, shift y to normal */ - n = -126 - iy; - hy <<= n; + if (!ey) { + for (i = uy.i<<9; i>>31 == 0; ey--, i <<= 1); + uy.i <<= -ey + 1; + } else { + uy.i &= -1U >> 9; + uy.i |= 1U << 23; } - /* fix point fmod */ - n = ix - iy; q = 0; - while (n--) { - hz = hx - hy; - if (hz < 0) - hx = hx << 1; - else { - hx = hz << 1; + if (ex < ey) { + if (ex+1 == ey) + goto end; + return x; + } + + /* x mod y */ + for (; ex > ey; ex--) { + i = uxi - uy.i; + if (i >> 31 == 0) { + uxi = i; q++; } + uxi <<= 1; q <<= 1; } - hz = hx - hy; - if (hz >= 0) { - hx = hz; + i = uxi - uy.i; + if (i >> 31 == 0) { + uxi = i; 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[(uint32_t)sx>>31]; - } - while (hx < 0x00800000) { /* normalize x */ - hx <<= 1; - iy--; + if (uxi == 0) + ex = -30; + else + for (; uxi>>23 == 0; uxi <<= 1, ex--); +end: + /* scale result and decide between |x| and |x|-|y| */ + if (ex > 0) { + uxi -= 1U << 23; + uxi |= (uint32_t)ex << 23; + } else { + uxi >>= -ex + 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++; + ux.i = uxi; + x = ux.f; + if (sy) + y = -y; + if (ex == ey || (ex+1 == ey && (2*x > y || (2*x == y && q%2)))) { x -= y; + q++; } - GET_FLOAT_WORD(hx, x); - SET_FLOAT_WORD(x, hx ^ sx); q &= 0x7fffffff; - *quo = sxy ? -q : q; - return x; + *quo = sx^sy ? -(int)q : (int)q; + return sx ? -x : x; } |