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: /* LINTLIBRARY */
   14: 
   15: #include <sys/cdefs.h>
   16: #include <float.h>
   17: #include <math.h>
   18: 
   19: #include "math_private.h"
   20: 
   21: static const double Zero[] = {0.0, -0.0,};
   22: 
   23: /*
   24:  * Return the IEEE remainder and set *quo to the last n bits of the
   25:  * quotient, rounded to the nearest integer.  We choose n=31 because
   26:  * we wind up computing all the integer bits of the quotient anyway as
   27:  * a side-effect of computing the remainder by the shift and subtract
   28:  * method.  In practice, this is far more bits than are needed to use
   29:  * remquo in reduction algorithms.
   30:  */
   31: double
   32: remquo(double x, double y, int *quo)
   33: {
   34:         int32_t n,hx,hy,hz,ix,iy,sx,i;
   35:         u_int32_t lx,ly,lz,q,sxy;
   36: 
   37:         EXTRACT_WORDS(hx,lx,x);
   38:         EXTRACT_WORDS(hy,ly,y);
   39:         sxy = (hx ^ hy) & 0x80000000;
   40:         sx = hx&0x80000000;            /* sign of x */
   41:         hx ^=sx;               /* |x| */
   42:         hy &= 0x7fffffff;      /* |y| */
   43: 
   44:     /* purge off exception values */
   45:         if((hy|ly)==0||(hx>=0x7ff00000)||      /* y=0,or x not finite */
   46:           ((hy|((ly|-ly)>>31))>0x7ff00000))    /* or y is NaN */
   47:             return (x*y)/(x*y);
   48:         if(hx<=hy) {
   49:             if((hx<hy)||(lx<ly)) {
   50:                 q = 0;
   51:                 goto fixup;   /* |x|<|y| return x or x-y */
   52:             }
   53:             if(lx==ly) {
   54:                 *quo = 1;
   55:                 return Zero[(u_int32_t)sx>>31];       /* |x|=|y| return x*0*/
   56:             }
   57:         }
   58: 
   59:     /* determine ix = ilogb(x) */
   60:         if(hx<0x00100000) {    /* subnormal x */
   61:             if(hx==0) {
   62:                 for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
   63:             } else {
   64:                 for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
   65:             }
   66:         } else ix = (hx>>20)-1023;
   67: 
   68:     /* determine iy = ilogb(y) */
   69:         if(hy<0x00100000) {    /* subnormal y */
   70:             if(hy==0) {
   71:                 for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
   72:             } else {
   73:                 for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
   74:             }
   75:         } else iy = (hy>>20)-1023;
   76: 
   77:     /* set up {hx,lx}, {hy,ly} and align y to x */
   78:         if(ix >= -1022) 
   79:             hx = 0x00100000|(0x000fffff&hx);
   80:         else {         /* subnormal x, shift x to normal */
   81:             n = -1022-ix;
   82:             if(n<=31) {
   83:                 hx = (hx<<n)|(lx>>(32-n));
   84:                 lx <<= n;
   85:             } else {
   86:                 hx = lx<<(n-32);
   87:                 lx = 0;
   88:             }
   89:         }
   90:         if(iy >= -1022) 
   91:             hy = 0x00100000|(0x000fffff&hy);
   92:         else {         /* subnormal y, shift y to normal */
   93:             n = -1022-iy;
   94:             if(n<=31) {
   95:                 hy = (hy<<n)|(ly>>(32-n));
   96:                 ly <<= n;
   97:             } else {
   98:                 hy = ly<<(n-32);
   99:                 ly = 0;
  100:             }
  101:         }
  102: 
  103:     /* fix point fmod */
  104:         n = ix - iy;
  105:         q = 0;
  106:         while(n--) {
  107:             hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
  108:             if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;}
  109:             else {hx = hz+hz+(lz>>31); lx = lz+lz; q++;}
  110:             q <<= 1;
  111:         }
  112:         hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
  113:         if(hz>=0) {hx=hz;lx=lz;q++;}
  114: 
  115:     /* convert back to floating value and restore the sign */
  116:         if((hx|lx)==0) {                       /* return sign(x)*0 */
  117:             *quo = (sxy ? -q : q);
  118:             return Zero[(u_int32_t)sx>>31];
  119:         }
  120:         while(hx<0x00100000) {         /* normalize x */
  121:             hx = hx+hx+(lx>>31); lx = lx+lx;
  122:             iy -= 1;
  123:         }
  124:         if(iy>= -1022) {       /* normalize output */
  125:             hx = ((hx-0x00100000)|((iy+1023)<<20));
  126:         } else {               /* subnormal output */
  127:             n = -1022 - iy;
  128:             if(n<=20) {
  129:                 lx = (lx>>n)|((u_int32_t)hx<<(32-n));
  130:                 hx >>= n;
  131:             } else if (n<=31) {
  132:                 lx = (hx<<(32-n))|(lx>>n); hx = sx;
  133:             } else {
  134:                 lx = hx>>(n-32); hx = sx;
  135:             }
  136:         }
  137: fixup:
  138:         INSERT_WORDS(x,hx,lx);
  139:         y = fabs(y);
  140:         if (y < 0x1p-1021) {
  141:             if (x+x>y || (x+x==y && (q & 1))) {
  142:                 q++;
  143:                 x-=y;
  144:             }
  145:         } else if (x>0.5*y || (x==0.5*y && (q & 1))) {
  146:             q++;
  147:             x-=y;
  148:         }
  149:         GET_HIGH_WORD(hx,x);
  150:         SET_HIGH_WORD(x,hx^sx);
  151:         q &= 0x7fffffff;
  152:         *quo = (sxy ? -q : q);
  153:         return x;
  154: }
  155: 
  156: #if     LDBL_MANT_DIG == 53
  157: #ifdef  lint
  158: /* PROTOLIB1 */
  159: long double remquol(long double, long double, int *);
  160: #else   /* lint */
  161: __weak_alias(remquol, remquo);
  162: #endif  /* lint */
  163: #endif  /* LDBL_MANT_DIG == 53 */