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