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