diff options
Diffstat (limited to 'newlib/libm/common/sf_remquo.c')
| -rw-r--r-- | newlib/libm/common/sf_remquo.c | 150 |
1 files changed, 115 insertions, 35 deletions
diff --git a/newlib/libm/common/sf_remquo.c b/newlib/libm/common/sf_remquo.c index 41c2038d4..7a756af8d 100644 --- a/newlib/libm/common/sf_remquo.c +++ b/newlib/libm/common/sf_remquo.c @@ -1,50 +1,130 @@ -/* Copyright (C) 2002 by Red Hat, Incorporated. All rights reserved. +/* Adapted for Newlib, 2009. (Allow for int < 32 bits; return *quo=0 during + * errors to make test scripts easier.) */ +/* @(#)e_fmod.c 1.3 95/01/18 */ +/*- + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * - * Permission to use, copy, modify, and distribute this software - * is freely granted, provided that this notice is preserved. + * 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 <math.h> #include "fdlibm.h" -#ifdef __STDC__ - float remquof(float x, float y, int *quo) /* wrapper remquof */ -#else - float remquof(x,y,quo) /* wrapper remquof */ - float x,y; - int *quo; +/* For quotient, return either all 31 bits that can from calculation (using + * int32_t), or as many as can fit into an int that is smaller than 32 bits. */ +#if INT_MAX > 0x7FFFFFFFL + #define QUO_MASK 0x7FFFFFFF +# else + #define QUO_MASK INT_MAX #endif -{ - int signx, signy, signres; - int wx; - int wy; - float x_over_y; - GET_FLOAT_WORD(wx, x); - GET_FLOAT_WORD(wy, y); +static const float Zero[] = {0.0, -0.0,}; - signx = (wx & 0x80000000) >> 31; - signy = (wy & 0x80000000) >> 31; +/* + * Return the IEEE remainder and set *quo to the last n bits of the + * quotient, rounded to the nearest integer. We choose n=31--if that many fit-- + * 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; + __uint32_t q,sxy; - signres = (signx ^ signy) ? -1 : 1; + 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| */ - x_over_y = fabsf(x / y); + /* purge off exception values */ + if(hy==0||hx>=0x7f800000||hy>0x7f800000) { /* y=0,NaN;or x not finite */ + *quo = 0; /* Not necessary, but return consistent value */ + return (x*y)/(x*y); + } + if(hx<hy) { + q = 0; + goto fixup; /* |x|<|y| return x or x-y */ + } else if(hx==hy) { + *quo = 1; + return Zero[(__uint32_t)sx>>31]; /* |x|=|y| return x*0*/ + } - *quo = signres * (lrintf(x_over_y) & 0x7f); + /* 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; - return remainderf(x,y); -} + /* 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; -#ifdef _DOUBLE_IS_32BITS + /* 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; + } -#ifdef __STDC__ - double remquo(double x, double y, int *quo) /* wrapper remquof */ -#else - double remquo(x,y,quo) /* wrapper remquof */ - double x,y; - int *quo; -#endif -{ - return (double) remquof((float) x, (float) y, quo); -} + /* 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++;} -#endif /* defined(_DOUBLE_IS_32BITS) */ + /* convert back to floating value and restore the sign */ + if(hx==0) { /* return sign(x)*0 */ + *quo = (sxy ? -q : q); + return Zero[(__uint32_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; +} |