sysdeps/ieee754/ldbl-128ibm/e_remainderl.c

   1 /* e_fmodl.c -- long double version of e_fmod.c.
   2  * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
   3  */
   4 /*
   5  * ====================================================
   6  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   7  *
   8  * Developed at SunPro, a Sun Microsystems, Inc. business.
   9  * Permission to use, copy, modify, and distribute this
  10  * software is freely granted, provided that this notice
  11  * is preserved.
  12  * ====================================================
  13  */
  14
  15 /* __ieee754_remainderl(x,p)
  16  * Return :
  17  *      returns  x REM p  =  x - [x/p]*p as if in infinite
  18  *      precise arithmetic, where [x/p] is the (infinite bit)
  19  *      integer nearest x/p (in half way case choose the even one).
  20  * Method :
  21  *      Based on fmodl() return x-[x/p]chopped*p exactlp.
  22  */
  23
  24 #include <math.h>
  25 #include <math_private.h>
  26
  27 static const long double zero = 0.0L;
  28
  29
  30 long double
  31 __ieee754_remainderl(long double x, long double p)
  32 {
  33         int64_t hx,hp;
  34         u_int64_t sx,lx,lp;
  35         long double p_half;
  36         double xhi, xlo, phi, plo;
  37
  38         ldbl_unpack (x, &xhi, &xlo);
  39         EXTRACT_WORDS64 (hx, xhi);
  40         EXTRACT_WORDS64 (lx, xlo);
  41         ldbl_unpack (p, &phi, &plo);
  42         EXTRACT_WORDS64 (hp, phi);
  43         EXTRACT_WORDS64 (lp, plo);
  44         sx = hx&0x8000000000000000ULL;
  45         hp &= 0x7fffffffffffffffLL;
  46         hx &= 0x7fffffffffffffffLL;
  47
  48     /* purge off exception values */
  49         if(hp==0) return (x*p)/(x*p);   /* p = 0 */
  50         if((hx>=0x7ff0000000000000LL)||                 /* x not finite */
  51            (hp>0x7ff0000000000000LL))                   /* p is NaN */
  52             return (x*p)/(x*p);
  53
  54
  55         if (hp<=0x7fdfffffffffffffLL) x = __ieee754_fmodl(x,p+p);       /* now x < 2p */
  56         if (((hx-hp)|(lx-lp))==0) return zero*x;
  57         x  = fabsl(x);
  58         p  = fabsl(p);
  59         if (hp<0x0020000000000000LL) {
  60             if(x+x>p) {
  61                 x-=p;
  62                 if(x+x>=p) x -= p;
  63             }
  64         } else {
  65             p_half = 0.5L*p;
  66             if(x>p_half) {
  67                 x-=p;
  68                 if(x>=p_half) x -= p;
  69             }
  70         }
  71         if (sx)
  72           x = -x;
  73         return x;
  74 }
  75 strong_alias (__ieee754_remainderl, __remainderl_finite)