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