diff options
Diffstat (limited to 'src/math/fmodf.c')
| -rw-r--r-- | src/math/fmodf.c | 139 |
1 files changed, 50 insertions, 89 deletions
diff --git a/src/math/fmodf.c b/src/math/fmodf.c index 837e9371..ff58f933 100644 --- a/src/math/fmodf.c +++ b/src/math/fmodf.c @@ -1,104 +1,65 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/e_fmodf.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. - * ==================================================== - */ -/* - * fmodf(x,y) - * Return x mod y in exact arithmetic - * Method: shift and subtract - */ - -#include "libm.h" - -static const float Zero[] = {0.0, -0.0,}; +#include <math.h> +#include <stdint.h> float fmodf(float x, float y) { - int32_t n,hx,hy,hz,ix,iy,sx,i; + union {float f; uint32_t i;} ux = {x}, uy = {y}; + int ex = ux.i>>23 & 0xff; + int ey = uy.i>>23 & 0xff; + uint32_t sx = ux.i & 0x80000000; + uint32_t i; + uint32_t uxi = ux.i; - GET_FLOAT_WORD(hx, x); - GET_FLOAT_WORD(hy, y); - sx = hx & 0x80000000; /* sign of x */ - hx ^= sx; /* |x| */ - hy &= 0x7fffffff; /* |y| */ - - /* purge off exception values */ - if (hy == 0 || hx >= 0x7f800000 || /* y=0,or x not finite */ - hy > 0x7f800000) /* or y is NaN */ + if (uy.i<<1 == 0 || isnan(y) || ex == 0xff) return (x*y)/(x*y); - if (hx < hy) /* |x| < |y| */ + if (uxi<<1 <= uy.i<<1) { + if (uxi<<1 == uy.i<<1) + return 0*x; return x; - if (hx == hy) /* |x| = |y|, return x*0 */ - return Zero[(uint32_t)sx>>31]; - - /* 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 = 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 = 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; - while (n--) { - hz = hx-hy; - if (hz<0) - hx = hx+hx; - else { - if(hz == 0) /* return sign(x)*0 */ - return Zero[(uint32_t)sx>>31]; - hx = hz+hz; + /* x mod y */ + for (; ex > ey; ex--) { + i = uxi - uy.i; + if (i >> 31 == 0) { + if (i == 0) + return 0*x; + uxi = i; } + uxi <<= 1; } - hz = hx-hy; - if (hz >= 0) - hx = hz; - - /* convert back to floating value and restore the sign */ - if (hx == 0) /* return sign(x)*0 */ - return Zero[(uint32_t)sx>>31]; - while (hx < 0x00800000) { /* normalize x */ - hx = hx+hx; - iy -= 1; + i = uxi - uy.i; + if (i >> 31 == 0) { + if (i == 0) + return 0*x; + uxi = i; } - if (iy >= -126) { /* normalize output */ - hx = ((hx-0x00800000)|((iy+127)<<23)); - SET_FLOAT_WORD(x, hx|sx); - } else { /* subnormal output */ - n = -126 - iy; - hx >>= n; - SET_FLOAT_WORD(x, hx|sx); + for (; uxi>>23 == 0; uxi <<= 1, ex--); + + /* scale result up */ + if (ex > 0) { + uxi -= 1U << 23; + uxi |= (uint32_t)ex << 23; + } else { + uxi >>= -ex + 1; } - return x; /* exact output */ + uxi |= sx; + ux.i = uxi; + return ux.f; } |