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