1/* Floating-point remainder function.
2 Copyright (C) 2023-2024 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 a 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 implicit one 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 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
64double
65__fmod (double x, double y)
66{
67 uint64_t hx = asuint64 (f: x);
68 uint64_t hy = asuint64 (f: y);
69
70 uint64_t sx = hx & SIGN_MASK;
71 /* Get |x| and |y|. */
72 hx ^= sx;
73 hy &= ~SIGN_MASK;
74
75 /* If x < y, return x (unless y is a NaN). */
76 if (__glibc_likely (hx < hy))
77 {
78 /* If y is a NaN, return a NaN. */
79 if (__glibc_unlikely (hy > EXPONENT_MASK))
80 return x * y;
81 return x;
82 }
83
84 int ex = hx >> MANTISSA_WIDTH;
85 int ey = hy >> MANTISSA_WIDTH;
86 int exp_diff = ex - ey;
87
88 /* Common case where exponents are close: |x/y| < 2^12, x not inf/NaN
89 and |x%y| not denormal. */
90 if (__glibc_likely (ey < (EXPONENT_MASK >> MANTISSA_WIDTH) - EXPONENT_WIDTH
91 && ey > MANTISSA_WIDTH
92 && exp_diff <= EXPONENT_WIDTH))
93 {
94 uint64_t mx = (hx << EXPONENT_WIDTH) | SIGN_MASK;
95 uint64_t my = (hy << EXPONENT_WIDTH) | SIGN_MASK;
96
97 mx %= (my >> exp_diff);
98
99 if (__glibc_unlikely (mx == 0))
100 return asdouble (i: sx);
101 int shift = clz_uint64 (x: mx);
102 ex -= shift + 1;
103 mx <<= shift;
104 mx = sx | (mx >> EXPONENT_WIDTH);
105 return asdouble (i: mx + ((uint64_t)ex << MANTISSA_WIDTH));
106 }
107
108 if (__glibc_unlikely (hy == 0 || hx >= EXPONENT_MASK))
109 {
110 /* If x is a NaN, return a NaN. */
111 if (hx > EXPONENT_MASK)
112 return x * y;
113
114 /* If x is an infinity or y is zero, return a NaN and set EDOM. */
115 return __math_edom (x: (x * y) / (x * y));
116 }
117
118 /* Special case, both x and y are denormal. */
119 if (__glibc_unlikely (ex == 0))
120 return asdouble (i: sx | hx % hy);
121
122 /* Extract normalized mantissas - hx is not denormal and hy != 0. */
123 uint64_t mx = get_mantissa (x: hx) | (MANTISSA_MASK + 1);
124 uint64_t my = get_mantissa (x: hy) | (MANTISSA_MASK + 1);
125 int lead_zeros_my = EXPONENT_WIDTH;
126
127 ey--;
128 /* Special case for denormal y. */
129 if (__glibc_unlikely (ey < 0))
130 {
131 my = hy;
132 ey = 0;
133 exp_diff--;
134 lead_zeros_my = clz_uint64 (x: my);
135 }
136
137 int tail_zeros_my = ctz_uint64 (x: my);
138 int sides_zeroes = lead_zeros_my + tail_zeros_my;
139
140 int right_shift = exp_diff < tail_zeros_my ? exp_diff : tail_zeros_my;
141 my >>= right_shift;
142 exp_diff -= right_shift;
143 ey += right_shift;
144
145 int left_shift = exp_diff < EXPONENT_WIDTH ? exp_diff : EXPONENT_WIDTH;
146 mx <<= left_shift;
147 exp_diff -= left_shift;
148
149 mx %= my;
150
151 if (__glibc_unlikely (mx == 0))
152 return asdouble (i: sx);
153
154 if (exp_diff == 0)
155 return make_double (x: mx, ep: ey, s: sx);
156
157 /* Multiplication with the inverse is faster than repeated modulo. */
158 uint64_t inv_hy = UINT64_MAX / my;
159 while (exp_diff > sides_zeroes) {
160 exp_diff -= sides_zeroes;
161 uint64_t hd = (mx * inv_hy) >> (BIT_WIDTH - sides_zeroes);
162 mx <<= sides_zeroes;
163 mx -= hd * my;
164 while (__glibc_unlikely (mx > my))
165 mx -= my;
166 }
167 uint64_t hd = (mx * inv_hy) >> (BIT_WIDTH - exp_diff);
168 mx <<= exp_diff;
169 mx -= hd * my;
170 while (__glibc_unlikely (mx > my))
171 mx -= my;
172
173 return make_double (x: mx, ep: ey, s: sx);
174}
175strong_alias (__fmod, __ieee754_fmod)
176libm_alias_finite (__ieee754_fmod, __fmod)
177#if LIBM_SVID_COMPAT
178versioned_symbol (libm, __fmod, fmod, GLIBC_2_38);
179libm_alias_double_other (__fmod, fmod)
180#else
181libm_alias_double (__fmod, fmod)
182#endif
183

source code of glibc/sysdeps/ieee754/dbl-64/e_fmod.c