]>
Commit | Line | Data |
---|---|---|
34b9f8bc AZN |
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/>. */ | |
f7eac6eb | 18 | |
16439f41 | 19 | #include <libm-alias-double.h> |
220622dd | 20 | #include <libm-alias-finite.h> |
16439f41 | 21 | #include <math-svid-compat.h> |
34b9f8bc AZN |
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 | |
f7eac6eb | 31 | |
34b9f8bc AZN |
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 | } */ | |
f7eac6eb | 58 | |
0ac5ae23 | 59 | double |
16439f41 | 60 | __fmod (double x, double y) |
f7eac6eb | 61 | { |
34b9f8bc AZN |
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. | |
16439f41 | 72 | - If x is an inifinity, a NaN is returned and EDOM is set. |
34b9f8bc AZN |
73 | - If y is zero, Nan is returned. |
74 | - If x is +0/-0, and y is not zero, +0/-0 is returned. */ | |
16439f41 AZN |
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 | } | |
34b9f8bc AZN |
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); | |
f7eac6eb | 162 | } |
16439f41 | 163 | strong_alias (__fmod, __ieee754_fmod) |
220622dd | 164 | libm_alias_finite (__ieee754_fmod, __fmod) |
16439f41 AZN |
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 |