]>
git.ipfire.org Git - thirdparty/glibc.git/blob - 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.
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-finite.h>
20 #include <libm-alias-float.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 double 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
48 mx * 2^ex == 2^8 * mx * 2^(ex - 8)
60 __fmodf (float x
, float y
)
62 uint32_t hx
= asuint (x
);
63 uint32_t hy
= asuint (y
);
65 uint32_t sx
= hx
& SIGN_MASK
;
66 /* Get |x| and |y|. */
71 - If x or y is a Nan, NaN is returned.
72 - If x is an inifinity, a NaN is returned.
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
))
78 if (is_nan (hx
) || is_nan (hy
))
79 return (x
* y
) / (x
* y
);
80 return __math_edomf ((x
* y
) / (x
* y
));
83 if (__glibc_unlikely (hx
<= hy
))
90 int ex
= hx
>> MANTISSA_WIDTH
;
91 int ey
= hy
>> MANTISSA_WIDTH
;
93 /* Common case where exponents are close: ey >= -103 and |x/y| < 2^8, */
94 if (__glibc_likely (ey
> MANTISSA_WIDTH
&& ex
- ey
<= EXPONENT_WIDTH
))
96 uint64_t mx
= (hx
& MANTISSA_MASK
) | (MANTISSA_MASK
+ 1);
97 uint64_t my
= (hy
& MANTISSA_MASK
) | (MANTISSA_MASK
+ 1);
99 uint32_t d
= (ex
== ey
) ? (mx
- my
) : (mx
<< (ex
- ey
)) % my
;
100 return make_float (d
, ey
- 1, sx
);
103 /* Special case, both x and y are subnormal. */
104 if (__glibc_unlikely (ex
== 0 && ey
== 0))
105 return asfloat (sx
| hx
% hy
);
107 /* Convert |x| and |y| to 'mx + 2^ex' and 'my + 2^ey'. Assume that hx is
108 not subnormal by conditions above. */
109 uint32_t mx
= get_mantissa (hx
) | (MANTISSA_MASK
+ 1);
112 uint32_t my
= get_mantissa (hy
) | (MANTISSA_MASK
+ 1);
113 int lead_zeros_my
= EXPONENT_WIDTH
;
114 if (__glibc_likely (ey
> 0))
119 lead_zeros_my
= __builtin_clz (my
);
122 int tail_zeros_my
= __builtin_ctz (my
);
123 int sides_zeroes
= lead_zeros_my
+ tail_zeros_my
;
124 int exp_diff
= ex
- ey
;
126 int right_shift
= exp_diff
< tail_zeros_my
? exp_diff
: tail_zeros_my
;
128 exp_diff
-= right_shift
;
131 int left_shift
= exp_diff
< EXPONENT_WIDTH
? exp_diff
: EXPONENT_WIDTH
;
133 exp_diff
-= left_shift
;
137 if (__glibc_unlikely (mx
== 0))
141 return make_float (mx
, ey
, sx
);
143 /* Assume modulo/divide operation is slow, so use multiplication with invert
145 uint32_t inv_hy
= UINT32_MAX
/ my
;
146 while (exp_diff
> sides_zeroes
) {
147 exp_diff
-= sides_zeroes
;
148 uint32_t hd
= (mx
* inv_hy
) >> (BIT_WIDTH
- sides_zeroes
);
151 while (__glibc_unlikely (mx
> my
))
154 uint32_t hd
= (mx
* inv_hy
) >> (BIT_WIDTH
- exp_diff
);
157 while (__glibc_unlikely (mx
> my
))
160 return make_float (mx
, ey
, sx
);
162 strong_alias (__fmodf
, __ieee754_fmodf
)
164 versioned_symbol (libm
, __fmodf
, fmodf
, GLIBC_2_38
);
165 libm_alias_float_other (__fmod
, fmod
)
167 libm_alias_float (__fmod
, fmod
)
169 libm_alias_finite (__ieee754_fmodf
, __fmodf
)