URL
https://opencores.org/ocsvn/openrisc/openrisc/trunk
Subversion Repositories openrisc
[/] [openrisc/] [trunk/] [gnu-old/] [newlib-1.17.0/] [newlib/] [libm/] [math/] [ef_fmod.c] - Rev 816
Compare with Previous | Blame | View Log
/* ef_fmod.c -- float version of e_fmod.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. * ==================================================== */ /* * __ieee754_fmodf(x,y) * Return x mod y in exact arithmetic * Method: shift and subtract */ #include "fdlibm.h" #ifdef __STDC__ static const float one = 1.0, Zero[] = {0.0, -0.0,}; #else static float one = 1.0, Zero[] = {0.0, -0.0,}; #endif #ifdef __STDC__ float __ieee754_fmodf(float x, float y) #else float __ieee754_fmodf(x,y) float x,y ; #endif { __int32_t n,hx,hy,hz,ix,iy,sx,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(FLT_UWORD_IS_ZERO(hy)|| !FLT_UWORD_IS_FINITE(hx)|| FLT_UWORD_IS_NAN(hy)) return (x*y)/(x*y); if(hx<hy) return x; /* |x|<|y| return x */ if(hx==hy) return Zero[(__uint32_t)sx>>31]; /* |x|=|y| return x*0*/ /* Note: y cannot be zero if we reach here. */ /* determine ix = ilogb(x) */ if(FLT_UWORD_IS_SUBNORMAL(hx)) { /* subnormal x */ for (ix = -126,i=(hx<<8); i>0; i<<=1) ix -=1; } else ix = (hx>>23)-127; /* determine iy = ilogb(y) */ if(FLT_UWORD_IS_SUBNORMAL(hy)) { /* 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; } if(iy >= -126) hy = 0x00800000|(0x007fffff&hy); else { /* subnormal y, shift y to normal */ n = -126-iy; hy = hy<<n; } /* 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; } } 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; } if(iy>= -126) { /* normalize output */ hx = ((hx-0x00800000)|((iy+127)<<23)); SET_FLOAT_WORD(x,hx|sx); } else { /* subnormal output */ /* If denormals are not supported, this code will generate a zero representation. */ n = -126 - iy; hx >>= n; SET_FLOAT_WORD(x,hx|sx); x *= one; /* create necessary signal */ } return x; /* exact output */ }