lib/libm/src/e_remainder.c

   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 #include <sys/cdefs.h>
  14 #if defined(LIBM_SCCS) && !defined(lint)
  15 __RCSID("$NetBSD: e_remainder.c,v 1.11 2002/05/26 22:01:52 wiz Exp $");
  16 #endif
  17
  18 /* __ieee754_remainder(x,p)
  19  * Return :
  20  *      returns  x REM p  =  x - [x/p]*p as if in infinite
  21  *      precise arithmetic, where [x/p] is the (infinite bit)
  22  *      integer nearest x/p (in half way case choose the even one).
  23  * Method :
  24  *      Based on fmod() return x-[x/p]chopped*p exactlp.
  25  */
  26
  27 #include "math.h"
  28 #include "math_private.h"
  29
  30 static const double zero = 0.0;
  31
  32
  33 double
  34 __ieee754_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 = __ieee754_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 }