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)