1  /* Compute remainder and a congruent to the quotient.
       2     Copyright (C) 1997-2023 Free Software Foundation, Inc.
       3     This file is part of the GNU C Library.
       4  
       5     The GNU C Library is free software; you can redistribute it and/or
       6     modify it under the terms of the GNU Lesser General Public
       7     License as published by the Free Software Foundation; either
       8     version 2.1 of the License, or (at your option) any later version.
       9  
      10     The GNU C Library is distributed in the hope that it will be useful,
      11     but WITHOUT ANY WARRANTY; without even the implied warranty of
      12     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      13     Lesser General Public License for more details.
      14  
      15     You should have received a copy of the GNU Lesser General Public
      16     License along with the GNU C Library; if not, see
      17     <https://www.gnu.org/licenses/>.  */
      18  
      19  #include <math.h>
      20  
      21  #include <math_private.h>
      22  #include <math_ldbl_opt.h>
      23  
      24  
      25  static const long double zero = 0.0;
      26  
      27  
      28  long double
      29  __remquol (long double x, long double y, int *quo)
      30  {
      31    int64_t hx,hy;
      32    uint64_t sx,lx,ly,qs;
      33    int cquo;
      34    double xhi, xlo, yhi, ylo;
      35  
      36    ldbl_unpack (x, &xhi, &xlo);
      37    EXTRACT_WORDS64 (hx, xhi);
      38    EXTRACT_WORDS64 (lx, xlo);
      39    ldbl_unpack (y, &yhi, &ylo);
      40    EXTRACT_WORDS64 (hy, yhi);
      41    EXTRACT_WORDS64 (ly, ylo);
      42    sx = hx & 0x8000000000000000ULL;
      43    qs = sx ^ (hy & 0x8000000000000000ULL);
      44    ly ^= hy & 0x8000000000000000ULL;
      45    hy &= 0x7fffffffffffffffLL;
      46    lx ^= sx;
      47    hx &= 0x7fffffffffffffffLL;
      48  
      49    /* Purge off exception values.  */
      50    if (hy == 0)
      51      return (x * y) / (x * y); 			/* y = 0 */
      52    if ((hx >= 0x7ff0000000000000LL)		/* x not finite */
      53        || (hy > 0x7ff0000000000000LL))		/* y is NaN */
      54      return (x * y) / (x * y);
      55  
      56    if (hy <= 0x7fbfffffffffffffLL)
      57      x = __ieee754_fmodl (x, 8 * y);              /* now x < 8y */
      58  
      59    if (((hx - hy) | (lx - ly)) == 0)
      60      {
      61        *quo = qs ? -1 : 1;
      62        return zero * x;
      63      }
      64  
      65    x  = fabsl (x);
      66    y  = fabsl (y);
      67    cquo = 0;
      68  
      69    if (hy <= 0x7fcfffffffffffffLL && x >= 4 * y)
      70      {
      71        x -= 4 * y;
      72        cquo += 4;
      73      }
      74    if (hy <= 0x7fdfffffffffffffLL && x >= 2 * y)
      75      {
      76        x -= 2 * y;
      77        cquo += 2;
      78      }
      79  
      80    if (hy < 0x0020000000000000LL)
      81      {
      82        if (x + x > y)
      83  	{
      84  	  x -= y;
      85  	  ++cquo;
      86  	  if (x + x >= y)
      87  	    {
      88  	      x -= y;
      89  	      ++cquo;
      90  	    }
      91  	}
      92      }
      93    else
      94      {
      95        long double y_half = 0.5L * y;
      96        if (x > y_half)
      97  	{
      98  	  x -= y;
      99  	  ++cquo;
     100  	  if (x >= y_half)
     101  	    {
     102  	      x -= y;
     103  	      ++cquo;
     104  	    }
     105  	}
     106      }
     107  
     108    *quo = qs ? -cquo : cquo;
     109  
     110    /* Ensure correct sign of zero result in round-downward mode.  */
     111    if (x == 0.0L)
     112      x = 0.0L;
     113    if (sx)
     114      x = -x;
     115    return x;
     116  }
     117  long_double_symbol (libm, __remquol, remquol);