gonzui

    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: /* 
   17:  * fmodf(x,y)
   18:  * Return x mod y in exact arithmetic
   19:  * Method: shift and subtract
   20:  */
   21: 
   22: #include "math.h"
   23: #include "math_private.h"
   24: 
   25: static const float one = 1.0, Zero[] = {0.0, -0.0,};
   26: 
   27: float
   28: fmodf(float x, float y)
   29: {
   30:         int32_t n,hx,hy,hz,ix,iy,sx,i;
   31: 
   32:         GET_FLOAT_WORD(hx,x);
   33:         GET_FLOAT_WORD(hy,y);
   34:         sx = hx&0x80000000;            /* sign of x */
   35:         hx ^=sx;               /* |x| */
   36:         hy &= 0x7fffffff;      /* |y| */
   37: 
   38:     /* purge off exception values */
   39:         if(hy==0||(hx>=0x7f800000)||           /* y=0,or x not finite */
   40:            (hy>0x7f800000))                    /* or y is NaN */
   41:             return (x*y)/(x*y);
   42:         if(hx<hy) return x;                    /* |x|<|y| return x */
   43:         if(hx==hy)
   44:             return Zero[(u_int32_t)sx>>31];    /* |x|=|y| return x*0*/
   45: 
   46:     /* determine ix = ilogb(x) */
   47:         if(hx<0x00800000) {    /* subnormal x */
   48:             for (ix = -126,i=(hx<<8); i>0; i<<=1) ix -=1;
   49:         } else ix = (hx>>23)-127;
   50: 
   51:     /* determine iy = ilogb(y) */
   52:         if(hy<0x00800000) {    /* subnormal y */
   53:             for (iy = -126,i=(hy<<8); i>=0; i<<=1) iy -=1;
   54:         } else iy = (hy>>23)-127;
   55: 
   56:     /* set up {hx,lx}, {hy,ly} and align y to x */
   57:         if(ix >= -126) 
   58:             hx = 0x00800000|(0x007fffff&hx);
   59:         else {         /* subnormal x, shift x to normal */
   60:             n = -126-ix;
   61:             hx = hx<<n;
   62:         }
   63:         if(iy >= -126) 
   64:             hy = 0x00800000|(0x007fffff&hy);
   65:         else {         /* subnormal y, shift y to normal */
   66:             n = -126-iy;
   67:             hy = hy<<n;
   68:         }
   69: 
   70:     /* fix point fmod */
   71:         n = ix - iy;
   72:         while(n--) {
   73:             hz=hx-hy;
   74:             if(hz<0){hx = hx+hx;}
   75:             else {
   76:                if(hz==0)                 /* return sign(x)*0 */
   77:                     return Zero[(u_int32_t)sx>>31];
   78:                hx = hz+hz;
   79:             }
   80:         }
   81:         hz=hx-hy;
   82:         if(hz>=0) {hx=hz;}
   83: 
   84:     /* convert back to floating value and restore the sign */
   85:         if(hx==0)                      /* return sign(x)*0 */
   86:             return Zero[(u_int32_t)sx>>31];    
   87:         while(hx<0x00800000) {         /* normalize x */
   88:             hx = hx+hx;
   89:             iy -= 1;
   90:         }
   91:         if(iy>= -126) {                /* normalize output */
   92:             hx = ((hx-0x00800000)|((iy+127)<<23));
   93:             SET_FLOAT_WORD(x,hx|sx);
   94:         } else {               /* subnormal output */
   95:             n = -126 - iy;
   96:             hx >>= n;
   97:             SET_FLOAT_WORD(x,hx|sx);
   98:             x *= one;          /* create necessary signal */
   99:         }
  100:         return x;              /* exact output */
  101: }