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  
      23  
      24  static const float zero = 0.0;
      25  #include <libm-alias-float.h>
      26  
      27  
      28  float
      29  __remquof (float x, float y, int *quo)
      30  {
      31    int32_t hx,hy;
      32    uint32_t sx;
      33    int cquo, qs;
      34  
      35    GET_FLOAT_WORD (hx, x);
      36    GET_FLOAT_WORD (hy, y);
      37    sx = hx & 0x80000000;
      38    qs = sx ^ (hy & 0x80000000);
      39    hy &= 0x7fffffff;
      40    hx &= 0x7fffffff;
      41  
      42    /* Purge off exception values.  */
      43    if (hy == 0)
      44      return (x * y) / (x * y); 			/* y = 0 */
      45    if ((hx >= 0x7f800000)			/* x not finite */
      46        || (hy > 0x7f800000))			/* y is NaN */
      47      return (x * y) / (x * y);
      48  
      49    if (hy <= 0x7dffffff)
      50      x = __ieee754_fmodf (x, 8 * y);		/* now x < 8y */
      51  
      52    if ((hx - hy) == 0)
      53      {
      54        *quo = qs ? -1 : 1;
      55        return zero * x;
      56      }
      57  
      58    x  = fabsf (x);
      59    y  = fabsf (y);
      60    cquo = 0;
      61  
      62    if (hy <= 0x7e7fffff && x >= 4 * y)
      63      {
      64        x -= 4 * y;
      65        cquo += 4;
      66      }
      67    if (hy <= 0x7effffff && x >= 2 * y)
      68      {
      69        x -= 2 * y;
      70        cquo += 2;
      71      }
      72  
      73    if (hy < 0x01000000)
      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        float 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.0f)
     105      x = 0.0f;
     106    if (sx)
     107      x = -x;
     108    return x;
     109  }
     110  libm_alias_float (__remquo, remquo)