compiler/stdc/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 #ifndef lint
  15 static char rcsid[] = "$FreeBSD: src/lib/msun/src/e_remainder.c,v 1.12 2008/03/30 20:47:42 das 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 <float.h>
  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         uint32_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 AROS_MAKE_ASM_SYM(typeof(remainderl), remainderl, AROS_CSYM_FROM_ASM_NAME(remainderl), AROS_CSYM_FROM_ASM_NAME(remainder));
  79 AROS_EXPORT_ASM_SYM(AROS_CSYM_FROM_ASM_NAME(remainderl));
  80 #endif