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@cygnus.com>, 1997 and
       5  		  Jakub Jelinek <jj@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  }