0,0 → 1,113 |
/* 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 */ |
} |