1
2 /*******************************************************************************
3 MIT License
4 -----------
5
6 Copyright (c) 2002-2019 Advanced Micro Devices, Inc.
7
8 Permission is hereby granted, free of charge, to any person obtaining a copy
9 of this Software and associated documentaon files (the "Software"), to deal
10 in the Software without restriction, including without limitation the rights
11 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12 copies of the Software, and to permit persons to whom the Software is
13 furnished to do so, subject to the following conditions:
14
15 The above copyright notice and this permission notice shall be included in
16 all copies or substantial portions of the Software.
17
18 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
24 THE SOFTWARE.
25 *******************************************************************************/
26
27 #include "libm.h"
28 #include "libm_util.h"
29
30 #define USE_NANF_WITH_FLAGS
31 #define USE_SCALEDOUBLE_1
32 #define USE_GET_FPSW_INLINE
33 #define USE_SET_FPSW_INLINE
34 #define USE_HANDLE_ERRORF
35 #include "libm_inlines.h"
36 #undef USE_NANF_WITH_FLAGS
37 #undef USE_SCALEDOUBLE_1
38 #undef USE_GET_FPSW_INLINE
39 #undef USE_SET_FPSW_INLINE
40 #undef USE_HANDLE_ERRORF
41
42 #if !defined(_CRTBLD_C9X)
43 #define _CRTBLD_C9X
44 #endif
45
46 #include "libm_errno.h"
47
48 // Disable "C4163: not available as intrinsic function" warning that older
49 // compilers may issue here.
50 #pragma warning(disable:4163)
51 #pragma function(remainderf,fmodf)
52
53
54 #undef _FUNCNAME
55 #if defined(COMPILING_FMOD)
fmodf(float x,float y)56 float fmodf(float x, float y)
57 #define _FUNCNAME "fmodf"
58 #define _OPERATION OP_FMOD
59 #else
60 float remainderf(float x, float y)
61 #define _FUNCNAME "remainderf"
62 #define _OPERATION OP_REM
63 #endif
64 {
65 double dx, dy, scale, w, t;
66 int i, ntimes, xexp, yexp;
67 unsigned long long ux, uy, ax, ay;
68
69 unsigned int sw;
70
71 dx = x;
72 dy = y;
73
74
75 GET_BITS_DP64(dx, ux);
76 GET_BITS_DP64(dy, uy);
77 ax = ux & ~SIGNBIT_DP64;
78 ay = uy & ~SIGNBIT_DP64;
79 xexp = (int)((ux & EXPBITS_DP64) >> EXPSHIFTBITS_DP64);
80 yexp = (int)((uy & EXPBITS_DP64) >> EXPSHIFTBITS_DP64);
81
82 if (xexp < 1 || xexp > BIASEDEMAX_DP64 ||
83 yexp < 1 || yexp > BIASEDEMAX_DP64)
84 {
85 /* x or y is zero, NaN or infinity (neither x nor y can be
86 denormalized because we promoted from float to double) */
87 if (xexp > BIASEDEMAX_DP64)
88 {
89 /* x is NaN or infinity */
90 if (ux & MANTBITS_DP64)
91 {
92 /* x is NaN */
93 unsigned int ufx;
94 GET_BITS_SP32(x, ufx);
95 return _handle_errorf(_FUNCNAME, _OPERATION, ufx|0x00400000, _DOMAIN, 0,
96 EDOM, x, y, 2);
97 }
98 else
99 {
100 /* x is infinity; result is NaN */
101 return _handle_errorf(_FUNCNAME, _OPERATION, INDEFBITPATT_SP32, _DOMAIN,
102 AMD_F_INVALID, EDOM, x, y, 2);
103 }
104 }
105 else if (yexp > BIASEDEMAX_DP64)
106 {
107 /* y is NaN or infinity */
108 if (uy & MANTBITS_DP64)
109 {
110 /* y is NaN */
111 unsigned int ufy;
112 GET_BITS_SP32(y, ufy);
113 return _handle_errorf(_FUNCNAME, _OPERATION, ufy|0x00400000, _DOMAIN, 0,
114 EDOM, x, y, 2);
115 }
116 else
117 {
118 #ifdef _CRTBLD_C9X
119 /* C99 return for y = +-inf is x */
120 return x;
121 #else
122 /* y is infinity; result is indefinite */
123 return _handle_errorf(_FUNCNAME, _OPERATION, INDEFBITPATT_SP32, _DOMAIN,
124 AMD_F_INVALID, EDOM, x, y, 2);
125 #endif
126 }
127 }
128 else if (xexp < 1)
129 {
130 /* x must be zero (cannot be denormalized) */
131 if (yexp < 1)
132 {
133 /* y must be zero (cannot be denormalized) */
134 return _handle_errorf(_FUNCNAME, _OPERATION, INDEFBITPATT_SP32, _DOMAIN,
135 AMD_F_INVALID, EDOM, x, y, 2);
136 }
137 else
138 /* C99 return for x = 0 must preserve sign */
139 return x;
140 }
141 else
142 {
143 /* y must be zero */
144 return _handle_errorf(_FUNCNAME, _OPERATION, INDEFBITPATT_SP32, _DOMAIN,
145 AMD_F_INVALID, EDOM, x, y, 2);
146 }
147 }
148 else if (ax == ay)
149 {
150 /* abs(x) == abs(y); return zero with the sign of x */
151 PUT_BITS_DP64(ux & SIGNBIT_DP64, dx);
152 return (float)dx;
153 }
154
155 /* Set dx = abs(x), dy = abs(y) */
156 PUT_BITS_DP64(ax, dx);
157 PUT_BITS_DP64(ay, dy);
158
159 if (ax < ay)
160 {
161 /* abs(x) < abs(y) */
162 #if !defined(COMPILING_FMOD)
163 if (dx > 0.5*dy)
164 dx -= dy;
165 #endif
166 return (float)(x < 0.0? -dx : dx);
167 }
168
169 /* Save the current floating-point status word. We need
170 to do this because the remainder function is always
171 exact for finite arguments, but our algorithm causes
172 the inexact flag to be raised. We therefore need to
173 restore the entry status before exiting. */
174 sw = get_fpsw_inline();
175
176 /* Set ntimes to the number of times we need to do a
177 partial remainder. If the exponent of x is an exact multiple
178 of 24 larger than the exponent of y, and the mantissa of x is
179 less than the mantissa of y, ntimes will be one too large
180 but it doesn't matter - it just means that we'll go round
181 the loop below one extra time. */
182 if (xexp <= yexp)
183 {
184 ntimes = 0;
185 w = dy;
186 scale = 1.0;
187 }
188 else
189 {
190 ntimes = (xexp - yexp) / 24;
191
192 /* Set w = y * 2^(24*ntimes) */
193 PUT_BITS_DP64((unsigned long long)(ntimes * 24 + EXPBIAS_DP64) << EXPSHIFTBITS_DP64,
194 scale);
195 w = scale * dy;
196 /* Set scale = 2^(-24) */
197 PUT_BITS_DP64((unsigned long long)(-24 + EXPBIAS_DP64) << EXPSHIFTBITS_DP64,
198 scale);
199 }
200
201
202 /* Each time round the loop we compute a partial remainder.
203 This is done by subtracting a large multiple of w
204 from x each time, where w is a scaled up version of y.
205 The subtraction can be performed exactly when performed
206 in double precision, and the result at each stage can
207 fit exactly in a single precision number. */
208 for (i = 0; i < ntimes; i++)
209 {
210 /* t is the integer multiple of w that we will subtract.
211 We use a truncated value for t. */
212 t = (double)((int)(dx / w));
213 dx -= w * t;
214 /* Scale w down by 2^(-24) for the next iteration */
215 w *= scale;
216 }
217
218 /* One more time */
219 #if defined(COMPILING_FMOD)
220 t = (double)((int)(dx / w));
221 dx -= w * t;
222 #else
223 {
224 unsigned int todd;
225 /* Variable todd says whether the integer t is odd or not */
226 t = (double)((int)(dx / w));
227 todd = ((int)(dx / w)) & 1;
228 dx -= w * t;
229
230 /* At this point, dx lies in the range [0,dy) */
231 /* For the remainder function, we need to adjust dx
232 so that it lies in the range (-y/2, y/2] by carefully
233 subtracting w (== dy == y) if necessary. */
234 if (dx > 0.5 * w || ((dx == 0.5 * w) && todd))
235 dx -= w;
236 }
237 #endif
238
239 /* **** N.B. for some reason this breaks the 32 bit version
240 of remainder when compiling with optimization. */
241 /* Restore the entry status flags */
242 set_fpsw_inline(sw);
243
244 /* Set the result sign according to input argument x */
245 return (float)(x < 0.0? -dx : dx);
246
247 }
248