sysdeps/ieee754/dbl-64/e_fmod.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-double.h>
  20 #include <libm-alias-finite.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 double floating point number (which uses
  44    less bits than the storage type), on each step the argument reduction can
  45    be improved by 11 (which is the size of uint64_t minus MANTISSA_WIDTH plus
  46    the signal bit):
  47
  48    mx * 2^ex == 2^11 * mx * 2^(ex - 11)
  49
  50    or
  51
  52    while (ex > ey)
  53      {
  54        mx << 11;
  55        ex -= 11;
  56        mx %= my;
  57      }  */
  58
  59 double
  60 __fmod (double x, double y)
  61 {
  62   uint64_t hx = asuint64 (x);
  63   uint64_t hy = asuint64 (y);
  64
  65   uint64_t sx = hx & SIGN_MASK;
  66   /* Get |x| and |y|.  */
  67   hx ^= sx;
  68   hy &= ~SIGN_MASK;
  69
  70   /* Special cases:
  71      - If x or y is a Nan, NaN is returned.
  72      - If x is an inifinity, a NaN is returned and EDOM is set.
  73      - If y is zero, Nan is returned.
  74      - If x is +0/-0, and y is not zero, +0/-0 is returned.  */
  75   if (__glibc_unlikely (hy == 0
  76                         || hx >= EXPONENT_MASK || hy > EXPONENT_MASK))
  77     {
  78       if (is_nan (hx) || is_nan (hy))
  79         return (x * y) / (x * y);
  80       return __math_edom ((x * y) / (x * y));
  81     }
  82
  83   if (__glibc_unlikely (hx <= hy))
  84     {
  85       if (hx < hy)
  86         return x;
  87       return asdouble (sx);
  88     }
  89
  90   int ex = hx >> MANTISSA_WIDTH;
  91   int ey = hy >> MANTISSA_WIDTH;
  92
  93   /* Common case where exponents are close: ey >= -907 and |x/y| < 2^52,  */
  94   if (__glibc_likely (ey > MANTISSA_WIDTH && ex - ey <= EXPONENT_WIDTH))
  95     {
  96       uint64_t mx = (hx & MANTISSA_MASK) | (MANTISSA_MASK + 1);
  97       uint64_t my = (hy & MANTISSA_MASK) | (MANTISSA_MASK + 1);
  98
  99       uint64_t d = (ex == ey) ? (mx - my) : (mx << (ex - ey)) % my;
 100       return make_double (d, ey - 1, sx);
 101     }
 102
 103   /* Special case, both x and y are subnormal.  */
 104   if (__glibc_unlikely (ex == 0 && ey == 0))
 105     return asdouble (sx | hx % hy);
 106
 107   /* Convert |x| and |y| to 'mx + 2^ex' and 'my + 2^ey'.  Assume that hx is
 108      not subnormal by conditions above.  */
 109   uint64_t mx = get_mantissa (hx) | (MANTISSA_MASK + 1);
 110   ex--;
 111   uint64_t my = get_mantissa (hy) | (MANTISSA_MASK + 1);
 112
 113   int lead_zeros_my = EXPONENT_WIDTH;
 114   if (__glibc_likely (ey > 0))
 115     ey--;
 116   else
 117     {
 118       my = hy;
 119       lead_zeros_my = clz_uint64 (my);
 120     }
 121
 122   /* Assume hy != 0  */
 123   int tail_zeros_my = ctz_uint64 (my);
 124   int sides_zeroes = lead_zeros_my + tail_zeros_my;
 125   int exp_diff = ex - ey;
 126
 127   int right_shift = exp_diff < tail_zeros_my ? exp_diff : tail_zeros_my;
 128   my >>= right_shift;
 129   exp_diff -= right_shift;
 130   ey += right_shift;
 131
 132   int left_shift = exp_diff < EXPONENT_WIDTH ? exp_diff : EXPONENT_WIDTH;
 133   mx <<= left_shift;
 134   exp_diff -= left_shift;
 135
 136   mx %= my;
 137
 138   if (__glibc_unlikely (mx == 0))
 139     return asdouble (sx);
 140
 141   if (exp_diff == 0)
 142     return make_double (mx, ey, sx);
 143
 144   /* Assume modulo/divide operation is slow, so use multiplication with invert
 145      values.  */
 146   uint64_t inv_hy = UINT64_MAX / my;
 147   while (exp_diff > sides_zeroes) {
 148     exp_diff -= sides_zeroes;
 149     uint64_t hd = (mx * inv_hy) >> (BIT_WIDTH - sides_zeroes);
 150     mx <<= sides_zeroes;
 151     mx -= hd * my;
 152     while (__glibc_unlikely (mx > my))
 153       mx -= my;
 154   }
 155   uint64_t hd = (mx * inv_hy) >> (BIT_WIDTH - exp_diff);
 156   mx <<= exp_diff;
 157   mx -= hd * my;
 158   while (__glibc_unlikely (mx > my))
 159     mx -= my;
 160
 161   return make_double (mx, ey, sx);
 162 }
 163 strong_alias (__fmod, __ieee754_fmod)
 164 libm_alias_finite (__ieee754_fmod, __fmod)
 165 #if LIBM_SVID_COMPAT
 166 versioned_symbol (libm, __fmod, fmod, GLIBC_2_38);
 167 libm_alias_double_other (__fmod, fmod)
 168 #else
 169 libm_alias_double (__fmod, fmod)
 170 #endif