gonzui

    1: /* @(#)e_remainder.c 5.1 93/09/24 */
    2: /*
    3:  * ====================================================
    4:  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
    5:  *
    6:  * Developed at SunPro, 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: /* remainder(x,p)
   16:  * Return :                  
   17:  *      returns  x REM p  =  x - [x/p]*p as if in infinite 
   18:  *      precise arithmetic, where [x/p] is the (infinite bit) 
   19:  *      integer nearest x/p (in half way case choose the even one).
   20:  * Method : 
   21:  *      Based on fmod() return x-[x/p]chopped*p exactlp.
   22:  */
   23: 
   24: #include <sys/cdefs.h>
   25: #include <float.h>
   26: #include <math.h>
   27: 
   28: #include "math_private.h"
   29: 
   30: static const double zero = 0.0;
   31: 
   32: 
   33: double
   34: remainder(double x, double p)
   35: {
   36:         int32_t hx,hp;
   37:         u_int32_t sx,lx,lp;
   38:         double p_half;
   39: 
   40:         EXTRACT_WORDS(hx,lx,x);
   41:         EXTRACT_WORDS(hp,lp,p);
   42:         sx = hx&0x80000000;
   43:         hp &= 0x7fffffff;
   44:         hx &= 0x7fffffff;
   45: 
   46:     /* purge off exception values */
   47:         if((hp|lp)==0) return (x*p)/(x*p);     /* p = 0 */
   48:         if((hx>=0x7ff00000)||                  /* x not finite */
   49:           ((hp>=0x7ff00000)&&                  /* p is NaN */
   50:           (((hp-0x7ff00000)|lp)!=0)))
   51:             return (x*p)/(x*p);
   52: 
   53: 
   54:         if (hp<=0x7fdfffff) x = fmod(x,p+p);   /* now x < 2p */
   55:         if (((hx-hp)|(lx-lp))==0) return zero*x;
   56:         x  = fabs(x);
   57:         p  = fabs(p);
   58:         if (hp<0x00200000) {
   59:             if(x+x>p) {
   60:                 x-=p;
   61:                 if(x+x>=p) x -= p;
   62:             }
   63:         } else {
   64:             p_half = 0.5*p;
   65:             if(x>p_half) {
   66:                 x-=p;
   67:                 if(x>=p_half) x -= p;
   68:             }
   69:         }
   70:         GET_HIGH_WORD(hx,x);
   71:         SET_HIGH_WORD(x,hx^sx);
   72:         return x;
   73: }
   74: 
   75: #if     LDBL_MANT_DIG == 53
   76: #ifdef  lint
   77: /* PROTOLIB1 */
   78: long double remainderl(long double, long double);
   79: #else   /* lint */
   80: __weak_alias(remainderl, remainder);
   81: #endif  /* lint */
   82: #endif  /* LDBL_MANT_DIG == 53 */