gonzui

    1: /* e_remainderf.c -- float version of e_remainder.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: #include "math.h"
   17: #include "math_private.h"
   18: 
   19: static const float zero = 0.0;
   20: 
   21: float
   22: remainderf(float x, float p)
   23: {
   24:         int32_t hx,hp;
   25:         u_int32_t sx;
   26:         float p_half;
   27: 
   28:         GET_FLOAT_WORD(hx,x);
   29:         GET_FLOAT_WORD(hp,p);
   30:         sx = hx&0x80000000;
   31:         hp &= 0x7fffffff;
   32:         hx &= 0x7fffffff;
   33: 
   34:     /* purge off exception values */
   35:         if(hp==0) return (x*p)/(x*p);           /* p = 0 */
   36:         if((hx>=0x7f800000)||                  /* x not finite */
   37:           ((hp>0x7f800000)))                   /* p is NaN */
   38:             return (x*p)/(x*p);
   39: 
   40: 
   41:         if (hp<=0x7effffff) x = fmodf(x,p+p);  /* now x < 2p */
   42:         if ((hx-hp)==0) return zero*x;
   43:         x  = fabsf(x);
   44:         p  = fabsf(p);
   45:         if (hp<0x01000000) {
   46:             if(x+x>p) {
   47:                 x-=p;
   48:                 if(x+x>=p) x -= p;
   49:             }
   50:         } else {
   51:             p_half = (float)0.5*p;
   52:             if(x>p_half) {
   53:                 x-=p;
   54:                 if(x>=p_half) x -= p;
   55:             }
   56:         }
   57:         GET_FLOAT_WORD(hx,x);
   58:         SET_FLOAT_WORD(x,hx^sx);
   59:         return x;
   60: }