libquadmath/math/remainderq.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 /* remainderq(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 fmodq() return x-[x/p]chopped*p exactlp.
  22  */
  23
  24 #include "quadmath-imp.h"
  25
  26 static const __float128 zero = 0.0Q;
  27
  28 __float128
  29 remainderq (__float128 x, __float128 p)
  30 {
  31   int64_t hx,hp;
  32   uint64_t sx,lx,lp;
  33   __float128 p_half;
  34
  35   GET_FLT128_WORDS64(hx,lx,x);
  36   GET_FLT128_WORDS64(hp,lp,p);
  37   sx = hx&0x8000000000000000ULL;
  38   hp &= 0x7fffffffffffffffLL;
  39   hx &= 0x7fffffffffffffffLL;
  40
  41   /* purge off exception values */
  42   if((hp|lp)==0) return (x*p)/(x*p);    /* p = 0 */
  43   if((hx>=0x7fff000000000000LL)||                       /* x not finite */
  44     ((hp>=0x7fff000000000000LL)&&                       /* p is NaN */
  45     (((hp-0x7fff000000000000LL)|lp)!=0)))
  46       return (x*p)/(x*p);
  47
  48   if (hp<=0x7ffdffffffffffffLL) x = fmodq (x,p+p);      /* now x < 2p */
  49   if (((hx-hp)|(lx-lp))==0) return zero*x;
  50   x  = fabsq(x);
  51   p  = fabsq(p);
  52   if (hp<0x0002000000000000LL) {
  53       if(x+x>p) {
  54           x-=p;
  55           if(x+x>=p) x -= p;
  56       }
  57   } else {
  58       p_half = 0.5Q*p;
  59       if(x>p_half) {
  60           x-=p;
  61           if(x>=p_half) x -= p;
  62       }
  63   }
  64   GET_FLT128_MSW64(hx,x);
  65   SET_FLT128_MSW64(x,hx^sx);
  66   return x;
  67 }