sysdeps/ieee754/ldbl-128/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 _Float128 zero = 0;
  28
  29
  30 _Float128
  31 __ieee754_remainderl(_Float128 x, _Float128 p)
  32 {
  33         int64_t hx,hp;
  34         u_int64_t sx,lx,lp;
  35         _Float128 p_half;
  36
  37         GET_LDOUBLE_WORDS64(hx,lx,x);
  38         GET_LDOUBLE_WORDS64(hp,lp,p);
  39         sx = hx&0x8000000000000000ULL;
  40         hp &= 0x7fffffffffffffffLL;
  41         hx &= 0x7fffffffffffffffLL;
  42
  43     /* purge off exception values */
  44         if((hp|lp)==0) return (x*p)/(x*p);      /* p = 0 */
  45         if((hx>=0x7fff000000000000LL)||                 /* x not finite */
  46           ((hp>=0x7fff000000000000LL)&&                 /* p is NaN */
  47           (((hp-0x7fff000000000000LL)|lp)!=0)))
  48             return (x*p)/(x*p);
  49
  50
  51         if (hp<=0x7ffdffffffffffffLL) x = __ieee754_fmodl(x,p+p);       /* now x < 2p */
  52         if (((hx-hp)|(lx-lp))==0) return zero*x;
  53         x  = fabsl(x);
  54         p  = fabsl(p);
  55         if (hp<0x0002000000000000LL) {
  56             if(x+x>p) {
  57                 x-=p;
  58                 if(x+x>=p) x -= p;
  59             }
  60         } else {
  61             p_half = L(0.5)*p;
  62             if(x>p_half) {
  63                 x-=p;
  64                 if(x>=p_half) x -= p;
  65             }
  66         }
  67         GET_LDOUBLE_MSW64(hx,x);
  68         SET_LDOUBLE_MSW64(x,hx^sx);
  69         return x;
  70 }
  71 strong_alias (__ieee754_remainderl, __remainderl_finite)