1 /* e_fmodf.c -- float version of e_fmod.c. 2 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. 3 */ 4 5 /* 6 * ==================================================== 7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 8 * 9 * Developed at SunPro, a Sun Microsystems, Inc. business. 10 * Permission to use, copy, modify, and distribute this 11 * software is freely granted, provided that this notice 12 * is preserved. 13 * ==================================================== 14 */ 15 16 #ifndef lint 17 static char rcsid[] = "$Id: e_fmodf.c,v 1.1.1.1 1994/08/19 09:39:54 jkh Exp $"; 18 #endif 19 20 /* 21 * __ieee754_fmodf(x,y) 22 * Return x mod y in exact arithmetic 23 * Method: shift and subtract 24 */ 25 26 #include "math.h" 27 #include "math_private.h" 28 29 #ifdef __STDC__ 30 static const float one = 1.0, Zero[] = {0.0, -0.0,}; 31 #else 32 static float one = 1.0, Zero[] = {0.0, -0.0,}; 33 #endif 34 35 #ifdef __STDC__ 36 float __ieee754_fmodf(float x, float y) 37 #else 38 float __ieee754_fmodf(x,y) 39 float x,y ; 40 #endif 41 { 42 int32_t n,hx,hy,hz,ix,iy,sx,i; 43 44 GET_FLOAT_WORD(hx,x); 45 GET_FLOAT_WORD(hy,y); 46 sx = hx&0x80000000; /* sign of x */ 47 hx ^=sx; /* |x| */ 48 hy &= 0x7fffffff; /* |y| */ 49 50 /* purge off exception values */ 51 if(hy==0||(hx>=0x7f800000)|| /* y=0,or x not finite */ 52 (hy>0x7f800000)) /* or y is NaN */ 53 return (x*y)/(x*y); 54 if(hx<hy) return x; /* |x|<|y| return x */ 55 if(hx==hy) 56 return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/ 57 58 /* determine ix = ilogb(x) */ 59 if(hx<0x00800000) { /* subnormal x */ 60 for (ix = -126,i=(hx<<8); i>0; i<<=1) ix -=1; 61 } else ix = (hx>>23)-127; 62 63 /* determine iy = ilogb(y) */ 64 if(hy<0x00800000) { /* subnormal y */ 65 for (iy = -126,i=(hy<<8); i>=0; i<<=1) iy -=1; 66 } else iy = (hy>>23)-127; 67 68 /* set up {hx,lx}, {hy,ly} and align y to x */ 69 if(ix >= -126) 70 hx = 0x00800000|(0x007fffff&hx); 71 else { /* subnormal x, shift x to normal */ 72 n = -126-ix; 73 hx = hx<<n; 74 } 75 if(iy >= -126) 76 hy = 0x00800000|(0x007fffff&hy); 77 else { /* subnormal y, shift y to normal */ 78 n = -126-iy; 79 hy = hy<<n; 80 } 81 82 /* fix point fmod */ 83 n = ix - iy; 84 while(n--) { 85 hz=hx-hy; 86 if(hz<0){hx = hx+hx;} 87 else { 88 if(hz==0) /* return sign(x)*0 */ 89 return Zero[(u_int32_t)sx>>31]; 90 hx = hz+hz; 91 } 92 } 93 hz=hx-hy; 94 if(hz>=0) {hx=hz;} 95 96 /* convert back to floating value and restore the sign */ 97 if(hx==0) /* return sign(x)*0 */ 98 return Zero[(u_int32_t)sx>>31]; 99 while(hx<0x00800000) { /* normalize x */ 100 hx = hx+hx; 101 iy -= 1; 102 } 103 if(iy>= -126) { /* normalize output */ 104 hx = ((hx-0x00800000)|((iy+127)<<23)); 105 SET_FLOAT_WORD(x,hx|sx); 106 } else { /* subnormal output */ 107 n = -126 - iy; 108 hx >>= n; 109 SET_FLOAT_WORD(x,hx|sx); 110 x *= one; /* create necessary signal */ 111 } 112 return x; /* exact output */ 113 } 114