diff options
Diffstat (limited to 'src/MSL_C.PPCEABI.bare.H/e_fmod.c')
| -rw-r--r-- | src/MSL_C.PPCEABI.bare.H/e_fmod.c | 167 |
1 files changed, 167 insertions, 0 deletions
diff --git a/src/MSL_C.PPCEABI.bare.H/e_fmod.c b/src/MSL_C.PPCEABI.bare.H/e_fmod.c new file mode 100644 index 0000000..b0f2524 --- /dev/null +++ b/src/MSL_C.PPCEABI.bare.H/e_fmod.c @@ -0,0 +1,167 @@ +/* @(#)e_fmod.c 1.2 95/01/04 */ +/* + * ==================================================== + * 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. + * ==================================================== + */ + +/* + * __ieee754_fmod(x,y) + * Return x mod y in exact arithmetic + * Method: shift and subtract + */ + +#include "fdlibm.h" + +#ifdef __STDC__ +static const double one = 1.0, Zero[] = { + 0.0, + -0.0, +}; +#else +static double one = 1.0, Zero[] = { + 0.0, + -0.0, +}; +#endif + +#ifdef __STDC__ +double __ieee754_fmod(double x, double y) +#else +double __ieee754_fmod(x, y) +double x, y; +#endif +{ + _INT32 n, hx, hy, hz, ix, iy, sx, i; + _UINT32 lx, ly, lz; + + hx = __HI(x); /* high word of x */ + lx = __LO(x); /* low word of x */ + hy = __HI(y); /* high word of y */ + ly = __LO(y); /* low word of y */ + 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 (x * y) / (x * y); + if (hx <= hy) { + if ((hx < hy) || (lx < ly)) + return x; /* |x|<|y| return x */ + if (lx == ly) + return Zero[(_UINT32)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; + 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 { + if ((hz | lz) == 0) /* return sign(x)*0 */ + return Zero[(_UINT32)sx >> 31]; + hx = hz + hz + (lz >> 31); + lx = lz + lz; + } + } + hz = hx - hy; + lz = lx - ly; + if (lx < ly) + hz -= 1; + if (hz >= 0) { + hx = hz; + lx = lz; + } + + /* convert back to floating value and restore the sign */ + if ((hx | lx) == 0) /* return sign(x)*0 */ + return Zero[(_UINT32)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)); + __HI(x) = hx | sx; + __LO(x) = lx; + } else { /* subnormal output */ + n = -1022 - iy; + if (n <= 20) { + lx = (lx >> n) | ((_UINT32)hx << (32 - n)); + hx >>= n; + } else if (n <= 31) { + lx = (hx << (32 - n)) | (lx >> n); + hx = sx; + } else { + lx = hx >> (n - 32); + hx = sx; + } + __HI(x) = hx | sx; + __LO(x) = lx; + x *= one; /* create necessary signal */ + } + return x; /* exact output */ +} |