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 <math.h>
  21 #include "math_config.h"
  22
  23 /* With x = mx * 2^ex and y = my * 2^ey (mx, my, ex, ey being integers), the
  24    simplest implementation is:
  25
  26    mx * 2^ex == 2 * mx * 2^(ex - 1)
  27
  28    or
  29
  30    while (ex > ey)
  31      {
  32        mx *= 2;
  33        --ex;
  34        mx %= my;
  35      }
  36
  37    With the mathematical equivalence of:
  38
  39    r == x % y == (x % (N * y)) % y
  40
  41    And with mx/my being mantissa of double floating point number (which uses
  42    less bits than the storage type), on each step the argument reduction can
  43    be improved by 8 (which is the size of uint32_t minus MANTISSA_WIDTH plus
  44    the signal bit):
  45
  46    mx * 2^ex == 2^8 * mx * 2^(ex - 8)
  47
  48    or
  49
  50    while (ex > ey)
  51      {
  52        mx << 8;
  53        ex -= 8;
  54        mx %= my;
  55      }  */
  56
  57 float
  58 __ieee754_fmodf (float x, float y)
  59 {
  60   uint32_t hx = asuint (x);
  61   uint32_t hy = asuint (y);
  62
  63   uint32_t sx = hx & SIGN_MASK;
  64   /* Get |x| and |y|.  */
  65   hx ^= sx;
  66   hy &= ~SIGN_MASK;
  67
  68   /* Special cases:
  69      - If x or y is a Nan, NaN is returned.
  70      - If x is an inifinity, a NaN is returned.
  71      - If y is zero, Nan is returned.
  72      - If x is +0/-0, and y is not zero, +0/-0 is returned.  */
  73   if (__glibc_unlikely (hy == 0 || hx >= EXPONENT_MASK || hy > EXPONENT_MASK))
  74     return (x * y) / (x * y);
  75
  76   if (__glibc_unlikely (hx <= hy))
  77     {
  78       if (hx < hy)
  79         return x;
  80       return asfloat (sx);
  81     }
  82
  83   int ex = hx >> MANTISSA_WIDTH;
  84   int ey = hy >> MANTISSA_WIDTH;
  85
  86   /* Common case where exponents are close: ey >= -103 and |x/y| < 2^8,  */
  87   if (__glibc_likely (ey > MANTISSA_WIDTH && ex - ey <= EXPONENT_WIDTH))
  88     {
  89       uint64_t mx = (hx & MANTISSA_MASK) | (MANTISSA_MASK + 1);
  90       uint64_t my = (hy & MANTISSA_MASK) | (MANTISSA_MASK + 1);
  91
  92       uint32_t d = (ex == ey) ? (mx - my) : (mx << (ex - ey)) % my;
  93       return make_float (d, ey - 1, sx);
  94     }
  95
  96   /* Special case, both x and y are subnormal.  */
  97   if (__glibc_unlikely (ex == 0 && ey == 0))
  98     return asfloat (sx | hx % hy);
  99
 100   /* Convert |x| and |y| to 'mx + 2^ex' and 'my + 2^ey'.  Assume that hx is
 101      not subnormal by conditions above.  */
 102   uint32_t mx = get_mantissa (hx) | (MANTISSA_MASK + 1);
 103   ex--;
 104
 105   uint32_t my = get_mantissa (hy) | (MANTISSA_MASK + 1);
 106   int lead_zeros_my = EXPONENT_WIDTH;
 107   if (__glibc_likely (ey > 0))
 108     ey--;
 109   else
 110     {
 111       my = hy;
 112       lead_zeros_my = __builtin_clz (my);
 113     }
 114
 115   int tail_zeros_my = __builtin_ctz (my);
 116   int sides_zeroes = lead_zeros_my + tail_zeros_my;
 117   int exp_diff = ex - ey;
 118
 119   int right_shift = exp_diff < tail_zeros_my ? exp_diff : tail_zeros_my;
 120   my >>= right_shift;
 121   exp_diff -= right_shift;
 122   ey += right_shift;
 123
 124   int left_shift = exp_diff < EXPONENT_WIDTH ? exp_diff : EXPONENT_WIDTH;
 125   mx <<= left_shift;
 126   exp_diff -= left_shift;
 127
 128   mx %= my;
 129
 130   if (__glibc_unlikely (mx == 0))
 131     return asfloat (sx);
 132
 133   if (exp_diff == 0)
 134     return make_float (mx, ey, sx);
 135
 136   /* Assume modulo/divide operation is slow, so use multiplication with invert
 137      values.  */
 138   uint32_t inv_hy = UINT32_MAX / my;
 139   while (exp_diff > sides_zeroes) {
 140     exp_diff -= sides_zeroes;
 141     uint32_t hd = (mx * inv_hy) >> (BIT_WIDTH - sides_zeroes);
 142     mx <<= sides_zeroes;
 143     mx -= hd * my;
 144     while (__glibc_unlikely (mx > my))
 145       mx -= my;
 146   }
 147   uint32_t hd = (mx * inv_hy) >> (BIT_WIDTH - exp_diff);
 148   mx <<= exp_diff;
 149   mx -= hd * my;
 150   while (__glibc_unlikely (mx > my))
 151     mx -= my;
 152
 153   return make_float (mx, ey, sx);
 154 }
 155 libm_alias_finite (__ieee754_fmodf, __fmodf)