]>
Commit | Line | Data |
---|---|---|
cf9cf331 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 | |
220622dd | 19 | #include <libm-alias-finite.h> |
cf9cf331 AZN |
20 | #include <math.h> |
21 | #include "math_config.h" | |
22 | ||
23 | /* With x = mx * 2^ex and y = my * 2^ey (mx, my, ex, ey being integers), the | |
24 | simplest implementation is: | |
25 | ||
26 | mx * 2^ex == 2 * mx * 2^(ex - 1) | |
27 | ||
28 | or | |
f7eac6eb | 29 | |
cf9cf331 AZN |
30 | while (ex > ey) |
31 | { | |
32 | mx *= 2; | |
33 | --ex; | |
34 | mx %= my; | |
35 | } | |
36 | ||
37 | With the mathematical equivalence of: | |
38 | ||
39 | r == x % y == (x % (N * y)) % y | |
40 | ||
41 | And with mx/my being mantissa of double floating point number (which uses | |
42 | less bits than the storage type), on each step the argument reduction can | |
43 | be improved by 8 (which is the size of uint32_t minus MANTISSA_WIDTH plus | |
44 | the signal bit): | |
45 | ||
46 | mx * 2^ex == 2^8 * mx * 2^(ex - 8) | |
47 | ||
48 | or | |
49 | ||
50 | while (ex > ey) | |
51 | { | |
52 | mx << 8; | |
53 | ex -= 8; | |
54 | mx %= my; | |
55 | } */ | |
f7eac6eb | 56 | |
0ac5ae23 UD |
57 | float |
58 | __ieee754_fmodf (float x, float y) | |
f7eac6eb | 59 | { |
cf9cf331 AZN |
60 | uint32_t hx = asuint (x); |
61 | uint32_t hy = asuint (y); | |
62 | ||
63 | uint32_t sx = hx & SIGN_MASK; | |
64 | /* Get |x| and |y|. */ | |
65 | hx ^= sx; | |
66 | hy &= ~SIGN_MASK; | |
67 | ||
68 | /* Special cases: | |
69 | - If x or y is a Nan, NaN is returned. | |
70 | - If x is an inifinity, a NaN is returned. | |
71 | - If y is zero, Nan is returned. | |
72 | - If x is +0/-0, and y is not zero, +0/-0 is returned. */ | |
73 | if (__glibc_unlikely (hy == 0 || hx >= EXPONENT_MASK || hy > EXPONENT_MASK)) | |
74 | return (x * y) / (x * y); | |
75 | ||
76 | if (__glibc_unlikely (hx <= hy)) | |
77 | { | |
78 | if (hx < hy) | |
79 | return x; | |
80 | return asfloat (sx); | |
81 | } | |
82 | ||
83 | int ex = hx >> MANTISSA_WIDTH; | |
84 | int ey = hy >> MANTISSA_WIDTH; | |
85 | ||
86 | /* Common case where exponents are close: ey >= -103 and |x/y| < 2^8, */ | |
87 | if (__glibc_likely (ey > MANTISSA_WIDTH && ex - ey <= EXPONENT_WIDTH)) | |
88 | { | |
89 | uint64_t mx = (hx & MANTISSA_MASK) | (MANTISSA_MASK + 1); | |
90 | uint64_t my = (hy & MANTISSA_MASK) | (MANTISSA_MASK + 1); | |
91 | ||
92 | uint32_t d = (ex == ey) ? (mx - my) : (mx << (ex - ey)) % my; | |
93 | return make_float (d, ey - 1, sx); | |
94 | } | |
95 | ||
96 | /* Special case, both x and y are subnormal. */ | |
97 | if (__glibc_unlikely (ex == 0 && ey == 0)) | |
98 | return asfloat (sx | hx % hy); | |
99 | ||
100 | /* Convert |x| and |y| to 'mx + 2^ex' and 'my + 2^ey'. Assume that hx is | |
101 | not subnormal by conditions above. */ | |
102 | uint32_t mx = get_mantissa (hx) | (MANTISSA_MASK + 1); | |
103 | ex--; | |
104 | ||
105 | uint32_t my = get_mantissa (hy) | (MANTISSA_MASK + 1); | |
106 | int lead_zeros_my = EXPONENT_WIDTH; | |
107 | if (__glibc_likely (ey > 0)) | |
108 | ey--; | |
109 | else | |
110 | { | |
111 | my = hy; | |
112 | lead_zeros_my = __builtin_clz (my); | |
113 | } | |
114 | ||
115 | int tail_zeros_my = __builtin_ctz (my); | |
116 | int sides_zeroes = lead_zeros_my + tail_zeros_my; | |
117 | int exp_diff = ex - ey; | |
118 | ||
119 | int right_shift = exp_diff < tail_zeros_my ? exp_diff : tail_zeros_my; | |
120 | my >>= right_shift; | |
121 | exp_diff -= right_shift; | |
122 | ey += right_shift; | |
123 | ||
124 | int left_shift = exp_diff < EXPONENT_WIDTH ? exp_diff : EXPONENT_WIDTH; | |
125 | mx <<= left_shift; | |
126 | exp_diff -= left_shift; | |
127 | ||
128 | mx %= my; | |
129 | ||
130 | if (__glibc_unlikely (mx == 0)) | |
131 | return asfloat (sx); | |
132 | ||
133 | if (exp_diff == 0) | |
134 | return make_float (mx, ey, sx); | |
135 | ||
136 | /* Assume modulo/divide operation is slow, so use multiplication with invert | |
137 | values. */ | |
138 | uint32_t inv_hy = UINT32_MAX / my; | |
139 | while (exp_diff > sides_zeroes) { | |
140 | exp_diff -= sides_zeroes; | |
141 | uint32_t hd = (mx * inv_hy) >> (BIT_WIDTH - sides_zeroes); | |
142 | mx <<= sides_zeroes; | |
143 | mx -= hd * my; | |
144 | while (__glibc_unlikely (mx > my)) | |
145 | mx -= my; | |
146 | } | |
147 | uint32_t hd = (mx * inv_hy) >> (BIT_WIDTH - exp_diff); | |
148 | mx <<= exp_diff; | |
149 | mx -= hd * my; | |
150 | while (__glibc_unlikely (mx > my)) | |
151 | mx -= my; | |
152 | ||
153 | return make_float (mx, ey, sx); | |
f7eac6eb | 154 | } |
220622dd | 155 | libm_alias_finite (__ieee754_fmodf, __fmodf) |