sysdeps/ieee754/flt-32/e_fmodf.c

   1 /* Floating-point remainder function.
   2    Copyright (C) 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 <libm-alias-finite.h>
  20 #include <libm-alias-float.h>
  21 #include <math-svid-compat.h>
  22 #include <math.h>
  23 #include "math_config.h"
  24
  25 /* With x = mx * 2^ex and y = my * 2^ey (mx, my, ex, ey being integers), the
  26    simplest implementation is:
  27
  28    mx * 2^ex == 2 * mx * 2^(ex - 1)
  29
  30    or
  31
  32    while (ex > ey)
  33      {
  34        mx *= 2;
  35        --ex;
  36        mx %= my;
  37      }
  38
  39    With the mathematical equivalence of:
  40
  41    r == x % y == (x % (N * y)) % y
  42
  43    And with mx/my being mantissa of a single floating point number (which uses
  44    less bits than the storage type), on each step the argument reduction can
  45    be improved by 8 (which is the size of uint32_t minus MANTISSA_WIDTH plus
  46    the implicit one bit):
  47
  48    mx * 2^ex == 2^8 * mx * 2^(ex - 8)
  49
  50    or
  51
  52    while (ex > ey)
  53      {
  54        mx << 8;
  55        ex -= 8;
  56        mx %= my;
  57      }
  58
  59    Special cases:
  60      - If x or y is a NaN, a NaN is returned.
  61      - If x is an infinity, or y is zero, a NaN is returned and EDOM is set.
  62      - If x is +0/-0, and y is not zero, +0/-0 is returned.  */
  63
  64 float
  65 __fmodf (float x, float y)
  66 {
  67   uint32_t hx = asuint (x);
  68   uint32_t hy = asuint (y);
  69
  70   uint32_t sx = hx & SIGN_MASK;
  71   /* Get |x| and |y|.  */
  72   hx ^= sx;
  73   hy &= ~SIGN_MASK;
  74
  75   if (__glibc_likely (hx < hy))
  76     {
  77       /* If y is a NaN, return a NaN.  */
  78       if (__glibc_unlikely (hy > EXPONENT_MASK))
  79         return x * y;
  80       return x;
  81     }
  82
  83   int ex = hx >> MANTISSA_WIDTH;
  84   int ey = hy >> MANTISSA_WIDTH;
  85   int exp_diff = ex - ey;
  86
  87   /* Common case where exponents are close: |x/y| < 2^9, x not inf/NaN
  88      and |x%y| not denormal.  */
  89   if (__glibc_likely (ey < (EXPONENT_MASK >> MANTISSA_WIDTH) - EXPONENT_WIDTH
  90                      && ey > MANTISSA_WIDTH
  91                      && exp_diff <= EXPONENT_WIDTH))
  92     {
  93       uint32_t mx = (hx << EXPONENT_WIDTH) | SIGN_MASK;
  94       uint32_t my = (hy << EXPONENT_WIDTH) | SIGN_MASK;
  95
  96       mx %= (my >> exp_diff);
  97
  98       if (__glibc_unlikely (mx == 0))
  99         return asfloat (sx);
 100       int shift = __builtin_clz (mx);
 101       ex -= shift + 1;
 102       mx <<= shift;
 103       mx = sx | (mx >> EXPONENT_WIDTH);
 104       return asfloat (mx + ((uint32_t)ex << MANTISSA_WIDTH));
 105     }
 106
 107   if (__glibc_unlikely (hy == 0 || hx >= EXPONENT_MASK))
 108     {
 109       /* If x is a NaN, return a NaN.  */
 110       if (hx > EXPONENT_MASK)
 111         return x * y;
 112
 113       /* If x is an infinity or y is zero, return a NaN and set EDOM.  */
 114       return __math_edomf ((x * y) / (x * y));
 115     }
 116
 117   /* Special case, both x and y are denormal.  */
 118   if (__glibc_unlikely (ex == 0))
 119     return asfloat (sx | hx % hy);
 120
 121   /* Extract normalized mantissas - hx is not denormal and hy != 0.  */
 122   uint32_t mx = get_mantissa (hx) | (MANTISSA_MASK + 1);
 123   uint32_t my = get_mantissa (hy) | (MANTISSA_MASK + 1);
 124   int lead_zeros_my = EXPONENT_WIDTH;
 125
 126   ey--;
 127   /* Special case for denormal y.  */
 128   if (__glibc_unlikely (ey < 0))
 129     {
 130       my = hy;
 131       ey = 0;
 132       exp_diff--;
 133       lead_zeros_my = __builtin_clz (my);
 134     }
 135
 136   int tail_zeros_my = __builtin_ctz (my);
 137   int sides_zeroes = lead_zeros_my + tail_zeros_my;
 138
 139   int right_shift = exp_diff < tail_zeros_my ? exp_diff : tail_zeros_my;
 140   my >>= right_shift;
 141   exp_diff -= right_shift;
 142   ey += right_shift;
 143
 144   int left_shift = exp_diff < EXPONENT_WIDTH ? exp_diff : EXPONENT_WIDTH;
 145   mx <<= left_shift;
 146   exp_diff -= left_shift;
 147
 148   mx %= my;
 149
 150   if (__glibc_unlikely (mx == 0))
 151     return asfloat (sx);
 152
 153   if (exp_diff == 0)
 154     return make_float (mx, ey, sx);
 155
 156   /* Multiplication with the inverse is faster than repeated modulo.  */
 157   uint32_t inv_hy = UINT32_MAX / my;
 158   while (exp_diff > sides_zeroes) {
 159     exp_diff -= sides_zeroes;
 160     uint32_t hd = (mx * inv_hy) >> (BIT_WIDTH - sides_zeroes);
 161     mx <<= sides_zeroes;
 162     mx -= hd * my;
 163     while (__glibc_unlikely (mx > my))
 164       mx -= my;
 165   }
 166   uint32_t hd = (mx * inv_hy) >> (BIT_WIDTH - exp_diff);
 167   mx <<= exp_diff;
 168   mx -= hd * my;
 169   while (__glibc_unlikely (mx > my))
 170     mx -= my;
 171
 172   return make_float (mx, ey, sx);
 173 }
 174 strong_alias (__fmodf, __ieee754_fmodf)
 175 #if LIBM_SVID_COMPAT
 176 versioned_symbol (libm, __fmodf, fmodf, GLIBC_2_38);
 177 libm_alias_float_other (__fmod, fmod)
 178 #else
 179 libm_alias_float (__fmod, fmod)
 180 #endif
 181 libm_alias_finite (__ieee754_fmodf, __fmodf)