1 /* Floating-point remainder function.
2 Copyright (C) 2023 Free Software Foundation, Inc.
3 This file is part of the GNU C Library.
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.
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.
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/>. */
19 #include <libm-alias-double.h>
20 #include <libm-alias-finite.h>
21 #include <math-svid-compat.h>
23 #include "math_config.h"
25 /* With x = mx * 2^ex and y = my * 2^ey (mx, my, ex, ey being integers), the
26 simplest implementation is:
28 mx * 2^ex == 2 * mx * 2^(ex - 1)
39 With the mathematical equivalence of:
41 r == x % y == (x % (N * y)) % y
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):
48 mx * 2^ex == 2^11 * mx * 2^(ex - 11)
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. */
65 __fmod (double x
, double y
)
67 uint64_t hx
= asuint64 (x
);
68 uint64_t hy
= asuint64 (y
);
70 uint64_t sx
= hx
& SIGN_MASK
;
71 /* Get |x| and |y|. */
75 /* If x < y, return x (unless y is a NaN). */
76 if (__glibc_likely (hx
< hy
))
78 /* If y is a NaN, return a NaN. */
79 if (__glibc_unlikely (hy
> EXPONENT_MASK
))
84 int ex
= hx
>> MANTISSA_WIDTH
;
85 int ey
= hy
>> MANTISSA_WIDTH
;
86 int exp_diff
= ex
- ey
;
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
))
94 uint64_t mx
= (hx
<< EXPONENT_WIDTH
) | SIGN_MASK
;
95 uint64_t my
= (hy
<< EXPONENT_WIDTH
) | SIGN_MASK
;
97 mx
%= (my
>> exp_diff
);
99 if (__glibc_unlikely (mx
== 0))
100 return asdouble (sx
);
101 int shift
= clz_uint64 (mx
);
104 mx
= sx
| (mx
>> EXPONENT_WIDTH
);
105 return asdouble (mx
+ ((uint64_t)ex
<< MANTISSA_WIDTH
));
108 if (__glibc_unlikely (hy
== 0 || hx
>= EXPONENT_MASK
))
110 /* If x is a NaN, return a NaN. */
111 if (hx
> EXPONENT_MASK
)
114 /* If x is an infinity or y is zero, return a NaN and set EDOM. */
115 return __math_edom ((x
* y
) / (x
* y
));
118 /* Special case, both x and y are denormal. */
119 if (__glibc_unlikely (ex
== 0))
120 return asdouble (sx
| hx
% hy
);
122 /* Extract normalized mantissas - hx is not denormal and hy != 0. */
123 uint64_t mx
= get_mantissa (hx
) | (MANTISSA_MASK
+ 1);
124 uint64_t my
= get_mantissa (hy
) | (MANTISSA_MASK
+ 1);
125 int lead_zeros_my
= EXPONENT_WIDTH
;
128 /* Special case for denormal y. */
129 if (__glibc_unlikely (ey
< 0))
134 lead_zeros_my
= clz_uint64 (my
);
137 int tail_zeros_my
= ctz_uint64 (my
);
138 int sides_zeroes
= lead_zeros_my
+ tail_zeros_my
;
140 int right_shift
= exp_diff
< tail_zeros_my
? exp_diff
: tail_zeros_my
;
142 exp_diff
-= right_shift
;
145 int left_shift
= exp_diff
< EXPONENT_WIDTH
? exp_diff
: EXPONENT_WIDTH
;
147 exp_diff
-= left_shift
;
151 if (__glibc_unlikely (mx
== 0))
152 return asdouble (sx
);
155 return make_double (mx
, ey
, sx
);
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
);
164 while (__glibc_unlikely (mx
> my
))
167 uint64_t hd
= (mx
* inv_hy
) >> (BIT_WIDTH
- exp_diff
);
170 while (__glibc_unlikely (mx
> my
))
173 return make_double (mx
, ey
, sx
);
175 strong_alias (__fmod
, __ieee754_fmod
)
176 libm_alias_finite (__ieee754_fmod
, __fmod
)
178 versioned_symbol (libm
, __fmod
, fmod
, GLIBC_2_38
);
179 libm_alias_double_other (__fmod
, fmod
)
181 libm_alias_double (__fmod
, fmod
)