lib/libm/src/e_remainder.c

   1
   2 /* @(#)e_remainder.c 1.3 95/01/18 */
   3 /*
   4  * ====================================================
   5  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   6  *
   7  * Developed at SunSoft, a Sun Microsystems, Inc. business.
   8  * Permission to use, copy, modify, and distribute this
   9  * software is freely granted, provided that this notice
  10  * is preserved.
  11  * ====================================================
  12  */
  13
  14
  15 /* __ieee754_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 <float.h>
  25
  26 #include "math.h"
  27 #include "math_private.h"
  28
  29 static const double zero = 0.0;
  30
  31
  32 double
  33 __ieee754_remainder(double x, double p)
  34 {
  35         int32_t hx,hp;
  36         u_int32_t sx,lx,lp;
  37         double p_half;
  38
  39         EXTRACT_WORDS(hx,lx,x);
  40         EXTRACT_WORDS(hp,lp,p);
  41         sx = hx&0x80000000;
  42         hp &= 0x7fffffff;
  43         hx &= 0x7fffffff;
  44
  45     /* purge off exception values */
  46         if((hp|lp)==0) return (x*p)/(x*p);      /* p = 0 */
  47         if((hx>=0x7ff00000)||                   /* x not finite */
  48           ((hp>=0x7ff00000)&&                   /* p is NaN */
  49           (((hp-0x7ff00000)|lp)!=0)))
  50             return ((long double)x*p)/((long double)x*p);
  51
  52
  53         if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p);  /* now x < 2p */
  54         if (((hx-hp)|(lx-lp))==0) return zero*x;
  55         x  = fabs(x);
  56         p  = fabs(p);
  57         if (hp<0x00200000) {
  58             if(x+x>p) {
  59                 x-=p;
  60                 if(x+x>=p) x -= p;
  61             }
  62         } else {
  63             p_half = 0.5*p;
  64             if(x>p_half) {
  65                 x-=p;
  66                 if(x>=p_half) x -= p;
  67             }
  68         }
  69         GET_HIGH_WORD(hx,x);
  70         if ((hx&0x7fffffff)==0) hx = 0;
  71         SET_HIGH_WORD(x,hx^sx);
  72         return x;
  73 }
  74
  75 #if LDBL_MANT_DIG == 53
  76 __weak_reference(remainder, remainderl);
  77 #endif