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 <libm-alias-double.h>
      23  #include <stdint.h>
      24  
      25  static const double zero = 0.0;
      26  
      27  
      28  double
      29  __remquo (double x, double y, int *quo)
      30  {
      31    int64_t hx, hy;
      32    uint64_t sx, qs;
      33    int cquo;
      34  
      35    EXTRACT_WORDS64 (hx, x);
      36    EXTRACT_WORDS64 (hy, y);
      37    sx = hx & UINT64_C(0x8000000000000000);
      38    qs = sx ^ (hy & UINT64_C(0x8000000000000000));
      39    hy &= UINT64_C(0x7fffffffffffffff);
      40    hx &= UINT64_C(0x7fffffffffffffff);
      41  
      42    /* Purge off exception values.  */
      43    if (__glibc_unlikely (hy == 0))
      44      return (x * y) / (x * y);			/* y = 0 */
      45    if (__builtin_expect (hx >= UINT64_C(0x7ff0000000000000) /* x not finite */
      46  			|| hy > UINT64_C(0x7ff0000000000000), 0))/* y is NaN */
      47      return (x * y) / (x * y);
      48  
      49    if (hy <= UINT64_C(0x7fbfffffffffffff))
      50      x = __ieee754_fmod (x, 8 * y);		/* now x < 8y */
      51  
      52    if (__glibc_unlikely (hx == hy))
      53      {
      54        *quo = qs ? -1 : 1;
      55        return zero * x;
      56      }
      57  
      58    x = fabs (x);
      59    INSERT_WORDS64 (y, hy);
      60    cquo = 0;
      61  
      62    if (hy <= UINT64_C(0x7fcfffffffffffff) && x >= 4 * y)
      63      {
      64        x -= 4 * y;
      65        cquo += 4;
      66      }
      67    if (hy <= UINT64_C(0x7fdfffffffffffff) && x >= 2 * y)
      68      {
      69        x -= 2 * y;
      70        cquo += 2;
      71      }
      72  
      73    if (hy < UINT64_C(0x0020000000000000))
      74      {
      75        if (x + x > y)
      76  	{
      77  	  x -= y;
      78  	  ++cquo;
      79  	  if (x + x >= y)
      80  	    {
      81  	      x -= y;
      82  	      ++cquo;
      83  	    }
      84  	}
      85      }
      86    else
      87      {
      88        double y_half = 0.5 * y;
      89        if (x > y_half)
      90  	{
      91  	  x -= y;
      92  	  ++cquo;
      93  	  if (x >= y_half)
      94  	    {
      95  	      x -= y;
      96  	      ++cquo;
      97  	    }
      98  	}
      99      }
     100  
     101    *quo = qs ? -cquo : cquo;
     102  
     103    /* Ensure correct sign of zero result in round-downward mode.  */
     104    if (x == 0.0)
     105      x = 0.0;
     106    if (sx)
     107      x = -x;
     108    return x;
     109  }
     110  libm_alias_double (__remquo, remquo)