nsz Git - musl/blob - src/math/remainder.c

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