gonzui

    1: /* @(#)e_fmod.c 1.3 95/01/18 */
    2: /*-
    3:  * ====================================================
    4:  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
    5:  *
    6:  * Developed at SunSoft, a Sun Microsystems, Inc. business.
    7:  * Permission to use, copy, modify, and distribute this
    8:  * software is freely granted, provided that this notice 
    9:  * is preserved.
   10:  * ====================================================
   11:  */
   12: 
   13: #include <sys/cdefs.h>
   14: 
   15: #include "math.h"
   16: #include "math_private.h"
   17: 
   18: static const float Zero[] = {0.0, -0.0,};
   19: 
   20: /*
   21:  * Return the IEEE remainder and set *quo to the last n bits of the
   22:  * quotient, rounded to the nearest integer.  We choose n=31 because
   23:  * we wind up computing all the integer bits of the quotient anyway as
   24:  * a side-effect of computing the remainder by the shift and subtract
   25:  * method.  In practice, this is far more bits than are needed to use
   26:  * remquo in reduction algorithms.
   27:  */
   28: float
   29: remquof(float x, float y, int *quo)
   30: {
   31:         int32_t n,hx,hy,hz,ix,iy,sx,i;
   32:         u_int32_t q,sxy;
   33: 
   34:         GET_FLOAT_WORD(hx,x);
   35:         GET_FLOAT_WORD(hy,y);
   36:         sxy = (hx ^ hy) & 0x80000000;
   37:         sx = hx&0x80000000;            /* sign of x */
   38:         hx ^=sx;               /* |x| */
   39:         hy &= 0x7fffffff;      /* |y| */
   40: 
   41:     /* purge off exception values */
   42:         if(hy==0||hx>=0x7f800000||hy>0x7f800000) /* y=0,NaN;or x not finite */
   43:             return (x*y)/(x*y);
   44:         if(hx<hy) {
   45:             q = 0;
   46:             goto fixup;        /* |x|<|y| return x or x-y */
   47:         } else if(hx==hy) {
   48:             *quo = 1;
   49:             return Zero[(u_int32_t)sx>>31];    /* |x|=|y| return x*0*/
   50:         }
   51: 
   52:     /* determine ix = ilogb(x) */
   53:         if(hx<0x00800000) {    /* subnormal x */
   54:             for (ix = -126,i=(hx<<8); i>0; i<<=1) ix -=1;
   55:         } else ix = (hx>>23)-127;
   56: 
   57:     /* determine iy = ilogb(y) */
   58:         if(hy<0x00800000) {    /* subnormal y */
   59:             for (iy = -126,i=(hy<<8); i>0; i<<=1) iy -=1;
   60:         } else iy = (hy>>23)-127;
   61: 
   62:     /* set up {hx,lx}, {hy,ly} and align y to x */
   63:         if(ix >= -126)
   64:             hx = 0x00800000|(0x007fffff&hx);
   65:         else {         /* subnormal x, shift x to normal */
   66:             n = -126-ix;
   67:             hx <<= n;
   68:         }
   69:         if(iy >= -126)
   70:             hy = 0x00800000|(0x007fffff&hy);
   71:         else {         /* subnormal y, shift y to normal */
   72:             n = -126-iy;
   73:             hy <<= n;
   74:         }
   75: 
   76:     /* fix point fmod */
   77:         n = ix - iy;
   78:         q = 0;
   79:         while(n--) {
   80:             hz=hx-hy;
   81:             if(hz<0) hx = hx << 1;
   82:             else {hx = hz << 1; q++;}
   83:             q <<= 1;
   84:         }
   85:         hz=hx-hy;
   86:         if(hz>=0) {hx=hz;q++;}
   87: 
   88:     /* convert back to floating value and restore the sign */
   89:         if(hx==0) {                            /* return sign(x)*0 */
   90:             *quo = (sxy ? -q : q);
   91:             return Zero[(u_int32_t)sx>>31];
   92:         }
   93:         while(hx<0x00800000) {         /* normalize x */
   94:             hx <<= 1;
   95:             iy -= 1;
   96:         }
   97:         if(iy>= -126) {                /* normalize output */
   98:             hx = ((hx-0x00800000)|((iy+127)<<23));
   99:         } else {               /* subnormal output */
  100:             n = -126 - iy;
  101:             hx >>= n;
  102:         }
  103: fixup:
  104:         SET_FLOAT_WORD(x,hx);
  105:         y = fabsf(y);
  106:         if (y < 0x1p-125f) {
  107:             if (x+x>y || (x+x==y && (q & 1))) {
  108:                 q++;
  109:                 x-=y;
  110:             }
  111:         } else if (x>0.5f*y || (x==0.5f*y && (q & 1))) {
  112:             q++;
  113:             x-=y;
  114:         }
  115:         GET_FLOAT_WORD(hx,x);
  116:         SET_FLOAT_WORD(x,hx^sx);
  117:         q &= 0x7fffffff;
  118:         *quo = (sxy ? -q : q);
  119:         return x;
  120: }