Home | History | Annotate | Line # | Download | only in math 
      1 /* Compute remainder and a congruent to the quotient.
      2    Copyright (C) 1997-2018 Free Software Foundation, Inc.
      3    This file is part of the GNU C Library.
      4    Contributed by Ulrich Drepper <drepper (at) cygnus.com>, 1997 and
      5 		  Jakub Jelinek <jj (at) ultra.linux.cz>, 1999.
      6 
      7    The GNU C Library is free software; you can redistribute it and/or
      8    modify it under the terms of the GNU Lesser General Public
      9    License as published by the Free Software Foundation; either
     10    version 2.1 of the License, or (at your option) any later version.
     11 
     12    The GNU C Library is distributed in the hope that it will be useful,
     13    but WITHOUT ANY WARRANTY; without even the implied warranty of
     14    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
     15    Lesser General Public License for more details.
     16 
     17    You should have received a copy of the GNU Lesser General Public
     18    License along with the GNU C Library; if not, see
     19    <http://www.gnu.org/licenses/>.  */
     20 
     21 #include "quadmath-imp.h"
     22 
     23 static const __float128 zero = 0.0;
     24 
     25 
     26 __float128
     27 remquoq (__float128 x, __float128 y, int *quo)
     28 {
     29   int64_t hx,hy;
     30   uint64_t sx,lx,ly,qs;
     31   int cquo;
     32 
     33   GET_FLT128_WORDS64 (hx, lx, x);
     34   GET_FLT128_WORDS64 (hy, ly, y);
     35   sx = hx & 0x8000000000000000ULL;
     36   qs = sx ^ (hy & 0x8000000000000000ULL);
     37   hy &= 0x7fffffffffffffffLL;
     38   hx &= 0x7fffffffffffffffLL;
     39 
     40   /* Purge off exception values.  */
     41   if ((hy | ly) == 0)
     42     return (x * y) / (x * y); 			/* y = 0 */
     43   if ((hx >= 0x7fff000000000000LL)		/* x not finite */
     44       || ((hy >= 0x7fff000000000000LL)		/* y is NaN */
     45 	  && (((hy - 0x7fff000000000000LL) | ly) != 0)))
     46     return (x * y) / (x * y);
     47 
     48   if (hy <= 0x7ffbffffffffffffLL)
     49     x = fmodq (x, 8 * y);              /* now x < 8y */
     50 
     51   if (((hx - hy) | (lx - ly)) == 0)
     52     {
     53       *quo = qs ? -1 : 1;
     54       return zero * x;
     55     }
     56 
     57   x  = fabsq (x);
     58   y  = fabsq (y);
     59   cquo = 0;
     60 
     61   if (hy <= 0x7ffcffffffffffffLL && x >= 4 * y)
     62     {
     63       x -= 4 * y;
     64       cquo += 4;
     65     }
     66   if (hy <= 0x7ffdffffffffffffLL && x >= 2 * y)
     67     {
     68       x -= 2 * y;
     69       cquo += 2;
     70     }
     71 
     72   if (hy < 0x0002000000000000LL)
     73     {
     74       if (x + x > y)
     75 	{
     76 	  x -= y;
     77 	  ++cquo;
     78 	  if (x + x >= y)
     79 	    {
     80 	      x -= y;
     81 	      ++cquo;
     82 	    }
     83 	}
     84     }
     85   else
     86     {
     87       __float128 y_half = 0.5Q * y;
     88       if (x > y_half)
     89 	{
     90 	  x -= y;
     91 	  ++cquo;
     92 	  if (x >= y_half)
     93 	    {
     94 	      x -= y;
     95 	      ++cquo;
     96 	    }
     97 	}
     98     }
     99 
    100   *quo = qs ? -cquo : cquo;
    101 
    102   /* Ensure correct sign of zero result in round-downward mode.  */
    103   if (x == 0)
    104     x = 0;
    105   if (sx)
    106     x = -x;
    107   return x;
    108 }
    109