ReactOS 0.4.15-dev-7846-g8ba6c66
remainder.c
Go to the documentation of this file.
1
2/*******************************************************************************
3MIT License
4-----------
5
6Copyright (c) 2002-2019 Advanced Micro Devices, Inc.
7
8Permission is hereby granted, free of charge, to any person obtaining a copy
9of this Software and associated documentaon files (the "Software"), to deal
10in the Software without restriction, including without limitation the rights
11to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12copies of the Software, and to permit persons to whom the Software is
13furnished to do so, subject to the following conditions:
14
15The above copyright notice and this permission notice shall be included in
16all copies or substantial portions of the Software.
17
18THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
24THE SOFTWARE.
25*******************************************************************************/
26
27#include "libm.h"
28#include "libm_util.h"
29
30#define USE_NAN_WITH_FLAGS
31#define USE_SCALEDOUBLE_3
32#define USE_GET_FPSW_INLINE
33#define USE_SET_FPSW_INLINE
34#define USE_HANDLE_ERROR
35#include "libm_inlines.h"
36#undef USE_NAN_WITH_FLAGS
37#undef USE_SCALEDOUBLE_3
38#undef USE_GET_FPSW_INLINE
39#undef USE_SET_FPSW_INLINE
40#undef USE_HANDLE_ERROR
41
42#if !defined(_CRTBLD_C9X)
43#define _CRTBLD_C9X
44#endif
45
46#include "libm_errno.h"
47
48/* Computes the exact product of x and y, the result being the
49 nearly doublelength number (z,zz) */
50static inline void dekker_mul12(double x, double y,
51 double *z, double *zz)
52{
53 double hx, tx, hy, ty;
54 /* Split x into hx (head) and tx (tail). Do the same for y. */
55 unsigned long long u;
57 u &= 0xfffffffff8000000;
58 PUT_BITS_DP64(u, hx);
59 tx = x - hx;
61 u &= 0xfffffffff8000000;
62 PUT_BITS_DP64(u, hy);
63 ty = y - hy;
64 *z = x * y;
65 *zz = (((hx * hy - *z) + hx * ty) + tx * hy) + tx * ty;
66}
67
68#pragma function(fmod)
69#undef _FUNCNAME
70#if defined(COMPILING_FMOD)
71double fmod(double x, double y)
72#define _FUNCNAME "fmod"
73#define _OPERATION OP_FMOD
74#else
75double remainder(double x, double y)
76#define _FUNCNAME "remainder"
77#define _OPERATION OP_REM
78#endif
79{
80 double dx, dy, scale, w, t, v, c, cc;
81 int i, ntimes, xexp, yexp;
82 unsigned long long u, ux, uy, ax, ay, todd;
83 unsigned int sw;
84
85 dx = x;
86 dy = y;
87
88
89 GET_BITS_DP64(dx, ux);
90 GET_BITS_DP64(dy, uy);
91 ax = ux & ~SIGNBIT_DP64;
92 ay = uy & ~SIGNBIT_DP64;
93 xexp = (int)((ux & EXPBITS_DP64) >> EXPSHIFTBITS_DP64);
94 yexp = (int)((uy & EXPBITS_DP64) >> EXPSHIFTBITS_DP64);
95
96 if (xexp < 1 || xexp > BIASEDEMAX_DP64 ||
97 yexp < 1 || yexp > BIASEDEMAX_DP64)
98 {
99 /* x or y is zero, denormalized, NaN or infinity */
100 if (xexp > BIASEDEMAX_DP64)
101 {
102 /* x is NaN or infinity */
103 if (ux & MANTBITS_DP64)
104 {
105 /* x is NaN */
106 return _handle_error(_FUNCNAME, _OPERATION, ux|0x0008000000000000, _DOMAIN, 0,
107 EDOM, x, y, 2);
108 }
109 else
110 {
111 /* x is infinity; result is NaN */
113 AMD_F_INVALID, EDOM, x, y, 2);
114 }
115 }
116 else if (yexp > BIASEDEMAX_DP64)
117 {
118 /* y is NaN or infinity */
119 if (uy & MANTBITS_DP64)
120 {
121 /* y is NaN */
122 return _handle_error(_FUNCNAME, _OPERATION, uy|0x0008000000000000, _DOMAIN, 0,
123 EDOM, x, y, 2);
124 }
125 else
126 {
127#ifdef _CRTBLD_C9X
128 /* C99 return for y = +-inf is x */
129 return x;
130#else
131 /* y is infinity; result is indefinite */
133 AMD_F_INVALID, EDOM, x, y, 2);
134#endif
135 }
136 }
137 else if (ax == 0x0000000000000000)
138 {
139 /* x is zero */
140 if (ay == 0x0000000000000000)
141 {
142 /* y is zero */
144 AMD_F_INVALID, EDOM, x, y, 2);
145 }
146 else
147 /* C99 return for x = 0 must preserve sign */
148 return x;
149 }
150 else if (ay == 0x0000000000000000)
151 {
152 /* y is zero */
154 AMD_F_INVALID, EDOM, x, y, 2);
155 }
156
157 /* We've exhausted all other possibilities. One or both of x and
158 y must be denormalized */
159 if (xexp < 1)
160 {
161 /* x is denormalized. Figure out its exponent. */
162 u = ax;
163 while (u < IMPBIT_DP64)
164 {
165 xexp--;
166 u <<= 1;
167 }
168 }
169 if (yexp < 1)
170 {
171 /* y is denormalized. Figure out its exponent. */
172 u = ay;
173 while (u < IMPBIT_DP64)
174 {
175 yexp--;
176 u <<= 1;
177 }
178 }
179 }
180 else if (ax == ay)
181 {
182 /* abs(x) == abs(y); return zero with the sign of x */
184 return dx;
185 }
186
187 /* Set x = abs(x), y = abs(y) */
189 PUT_BITS_DP64(ay, dy);
190
191 if (ax < ay)
192 {
193 /* abs(x) < abs(y) */
194#if !defined(COMPILING_FMOD)
195 if (dx > 0.5*dy)
196 dx -= dy;
197#endif
198 return x < 0.0? -dx : dx;
199 }
200
201 /* Save the current floating-point status word. We need
202 to do this because the remainder function is always
203 exact for finite arguments, but our algorithm causes
204 the inexact flag to be raised. We therefore need to
205 restore the entry status before exiting. */
206 sw = get_fpsw_inline();
207
208 /* Set ntimes to the number of times we need to do a
209 partial remainder. If the exponent of x is an exact multiple
210 of 52 larger than the exponent of y, and the mantissa of x is
211 less than the mantissa of y, ntimes will be one too large
212 but it doesn't matter - it just means that we'll go round
213 the loop below one extra time. */
214 if (xexp <= yexp)
215 ntimes = 0;
216 else
217 ntimes = (xexp - yexp) / 52;
218
219 if (ntimes == 0)
220 {
221 w = dy;
222 scale = 1.0;
223 }
224 else
225 {
226 /* Set w = y * 2^(52*ntimes) */
227 w = scaleDouble_3(dy, ntimes * 52);
228
229 /* Set scale = 2^(-52) */
230 PUT_BITS_DP64((unsigned long long)(-52 + EXPBIAS_DP64) << EXPSHIFTBITS_DP64,
231 scale);
232 }
233
234
235 /* Each time round the loop we compute a partial remainder.
236 This is done by subtracting a large multiple of w
237 from x each time, where w is a scaled up version of y.
238 The subtraction must be performed exactly in quad
239 precision, though the result at each stage can
240 fit exactly in a double precision number. */
241 for (i = 0; i < ntimes; i++)
242 {
243 /* t is the integer multiple of w that we will subtract.
244 We use a truncated value for t.
245
246 N.B. w has been chosen so that the integer t will have
247 at most 52 significant bits. This is the amount by
248 which the exponent of the partial remainder dx gets reduced
249 every time around the loop. In theory we could use
250 53 bits in t, but the quad precision multiplication
251 routine dekker_mul12 does not allow us to do that because
252 it loses the last (106th) bit of its quad precision result. */
253
254 /* Set dx = dx - w * t, where t is equal to trunc(dx/w). */
255 t = (double)(long long)(dx / w);
256 /* At this point, t may be one too large due to
257 rounding of dx/w */
258
259 /* Compute w * t in quad precision */
260 dekker_mul12(w, t, &c, &cc);
261
262 /* Subtract w * t from dx */
263 v = dx - c;
264 dx = v + (((dx - v) - c) - cc);
265
266 /* If t was one too large, dx will be negative. Add back
267 one w */
268 /* It might be possible to speed up this loop by finding
269 a way to compute correctly truncated t directly from dx and w.
270 We would then avoid the need for this check on negative dx. */
271 if (dx < 0.0)
272 dx += w;
273
274 /* Scale w down by 2^(-52) for the next iteration */
275 w *= scale;
276 }
277
278 /* One more time */
279 /* Variable todd says whether the integer t is odd or not */
280 t = (double)(long long)(dx / w);
281 todd = ((long long)(dx / w)) & 1;
282 dekker_mul12(w, t, &c, &cc);
283 v = dx - c;
284 dx = v + (((dx - v) - c) - cc);
285 if (dx < 0.0)
286 {
287 todd = !todd;
288 dx += w;
289 }
290
291 /* At this point, dx lies in the range [0,dy) */
292#if !defined(COMPILING_FMOD)
293 /* For the fmod function, we're done apart from setting
294 the correct sign. */
295 /* For the remainder function, we need to adjust dx
296 so that it lies in the range (-y/2, y/2] by carefully
297 subtracting w (== dy == y) if necessary. The rigmarole
298 with todd is to get the correct sign of the result
299 when x/y lies exactly half way between two integers,
300 when we need to choose the even integer. */
301 if (ay < 0x7fd0000000000000)
302 {
303 if (dx + dx > w || (todd && (dx + dx == w)))
304 dx -= w;
305 }
306 else if (dx > 0.5 * w || (todd && (dx == 0.5 * w)))
307 dx -= w;
308
309#endif
310
311 /* **** N.B. for some reason this breaks the 32 bit version
312 of remainder when compiling with optimization. */
313 /* Restore the entry status flags */
314 set_fpsw_inline(sw);
315
316 /* Set the result sign according to input argument x */
317 return x < 0.0? -dx : dx;
318
319}
double __cdecl _handle_error(char *fname, int opcode, unsigned long long value, int type, int flags, int error, double arg1, double arg2, int nargs)
Handles an error condition.
Definition: _handle_error.c:34
unsigned int(__cdecl typeof(jpeg_read_scanlines))(struct jpeg_decompress_struct *
Definition: typeof.h:31
#define EDOM
Definition: errno.h:39
GLint GLint GLint GLint GLint x
Definition: gl.h:1548
const GLdouble * v
Definition: gl.h:2040
GLint GLint GLint GLint GLint GLint y
Definition: gl.h:1548
GLdouble GLdouble t
Definition: gl.h:2047
GLenum GLenum GLenum GLenum GLenum scale
Definition: glext.h:9032
const GLubyte * c
Definition: glext.h:8905
GLubyte GLubyte GLubyte GLubyte w
Definition: glext.h:6102
GLbyte ty
Definition: glext.h:8756
GLdouble GLdouble z
Definition: glext.h:5874
GLsizei GLenum const GLvoid GLsizei GLenum GLbyte GLbyte GLbyte GLdouble GLdouble GLdouble GLfloat GLfloat GLfloat GLint GLint GLint GLshort GLshort GLshort GLubyte GLubyte GLubyte GLuint GLuint GLuint GLushort GLushort GLushort GLbyte GLbyte GLbyte GLbyte GLdouble GLdouble GLdouble GLdouble GLfloat GLfloat GLfloat GLfloat GLint GLint GLint GLint GLshort GLshort GLshort GLshort GLubyte GLubyte GLubyte GLubyte GLuint GLuint GLuint GLuint GLushort GLushort GLushort GLushort GLboolean const GLdouble const GLfloat const GLint const GLshort const GLbyte const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLdouble const GLfloat const GLfloat const GLint const GLint const GLshort const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort GLenum GLenum GLenum GLfloat GLenum GLint GLenum GLenum GLenum GLfloat GLenum GLenum GLint GLenum GLfloat GLenum GLint GLint GLushort GLenum GLenum GLfloat GLenum GLenum GLint GLfloat const GLubyte GLenum GLenum GLenum const GLfloat GLenum GLenum const GLint GLenum GLint GLint GLsizei GLsizei GLint GLenum GLenum const GLvoid GLenum GLenum const GLfloat GLenum GLenum const GLint GLenum GLenum const GLdouble GLenum GLenum const GLfloat GLenum GLenum const GLint GLsizei GLuint GLfloat GLuint GLbitfield GLfloat GLint GLuint GLboolean GLenum GLfloat GLenum GLbitfield GLenum GLfloat GLfloat GLint GLint const GLfloat GLenum GLfloat GLfloat GLint GLint GLfloat GLfloat GLint GLint const GLfloat GLint GLfloat GLfloat GLint GLfloat GLfloat GLint GLfloat GLfloat const GLdouble const GLfloat const GLdouble const GLfloat GLint i
Definition: glfuncs.h:248
GLsizei GLenum const GLvoid GLsizei GLenum GLbyte GLbyte GLbyte GLdouble GLdouble GLdouble GLfloat GLfloat GLfloat GLint GLint GLint GLshort GLshort GLshort GLubyte GLubyte GLubyte GLuint GLuint GLuint GLushort GLushort GLushort GLbyte GLbyte GLbyte GLbyte GLdouble GLdouble GLdouble GLdouble GLfloat GLfloat GLfloat GLfloat GLint GLint GLint GLint GLshort GLshort GLshort GLshort GLubyte GLubyte GLubyte GLubyte GLuint GLuint GLuint GLuint GLushort GLushort GLushort GLushort GLboolean const GLdouble const GLfloat const GLint const GLshort const GLbyte const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLdouble const GLfloat const GLfloat const GLint const GLint const GLshort const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort const GLdouble const GLfloat const GLint const GLshort GLenum GLenum GLenum GLfloat GLenum GLint GLenum GLenum GLenum GLfloat GLenum GLenum GLint GLenum GLfloat GLenum GLint GLint GLushort GLenum GLenum GLfloat GLenum GLenum GLint GLfloat const GLubyte GLenum GLenum GLenum const GLfloat GLenum GLenum const GLint GLenum GLint GLint GLsizei GLsizei GLint GLenum GLenum const GLvoid GLenum GLenum const GLfloat GLenum GLenum const GLint GLenum GLenum const GLdouble GLenum GLenum const GLfloat GLenum GLenum const GLint GLsizei GLuint GLfloat GLuint GLbitfield GLfloat GLint GLuint GLboolean GLenum GLfloat GLenum GLbitfield GLenum GLfloat GLfloat GLint GLint const GLfloat GLenum GLfloat GLfloat GLint GLint GLfloat GLfloat GLint GLint const GLfloat GLint GLfloat GLfloat GLint GLfloat GLfloat GLint GLfloat GLfloat const GLdouble * u
Definition: glfuncs.h:240
#define _DOMAIN
Definition: math.h:39
_Check_return_ double __cdecl fmod(_In_ double x, _In_ double y)
uint32_t cc
Definition: isohybrid.c:75
#define c
Definition: ke_i.h:80
#define AMD_F_INVALID
Definition: libm_new.h:86
#define IMPBIT_DP64
Definition: libm_util.h:50
#define INDEFBITPATT_DP64
Definition: libm_util.h:52
#define EXPSHIFTBITS_DP64
Definition: libm_util.h:56
#define GET_BITS_DP64(x, ux)
Definition: libm_util.h:118
#define BIASEDEMAX_DP64
Definition: libm_util.h:59
#define SIGNBIT_DP64
Definition: libm_util.h:44
#define EXPBITS_DP64
Definition: libm_util.h:45
#define EXPBIAS_DP64
Definition: libm_util.h:55
#define MANTBITS_DP64
Definition: libm_util.h:46
#define PUT_BITS_DP64(ux, x)
Definition: libm_util.h:124
GLint dy
Definition: linetemp.h:97
GLint dx
Definition: linetemp.h:97
static const char mbstate_t *static wchar_t const char mbstate_t *static const wchar_t int *static double
Definition: string.c:80
#define long
Definition: qsort.c:33
#define _FUNCNAME
#define _OPERATION
static void dekker_mul12(double x, double y, double *z, double *zz)
Definition: remainder.c:50
double remainder(double x, double y)
Definition: remainder.c:75
ecx edi movl ebx edx edi decl ecx esi eax jecxz decl eax andl eax esi movl edx movl TEMP incl eax andl eax ecx incl ebx testl eax jnz xchgl ecx incl TEMP esp ecx subl ebx pushl ecx ecx edx ecx shrl ecx mm0 mm4 mm0 mm4 mm1 mm5 mm1 mm5 mm2 mm6 mm2 mm6 mm3 mm7 mm3 mm7 paddd mm0 paddd mm4 paddd mm0 paddd mm4 paddd mm0 paddd mm4 movq mm1 movq mm5 psrlq mm1 psrlq mm5 paddd mm0 paddd mm4 psrad mm0 psrad mm4 packssdw mm0 packssdw mm4 mm1 punpckldq mm0 pand mm1 pand mm0 por mm1 movq edi esi edx edi decl ecx jnz popl ecx andl ecx jecxz mm0 mm0 mm1 mm1 mm2 mm2 mm3 mm3 paddd mm0 paddd mm0 paddd mm0 movq mm1 psrlq mm1 paddd mm0 psrad mm0 packssdw mm0 movd eax movw ax
Definition: synth_sse3d.h:180