ReactOS 0.4.15-dev-7693-gd0e2924
remainderf.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_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)
56float fmodf(float x, float y)
57#define _FUNCNAME "fmodf"
58#define _OPERATION OP_FMOD
59#else
60float 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 */
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 */
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) */
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 */
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 */
152 return (float)dx;
153 }
154
155 /* Set dx = abs(x), dy = abs(y) */
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}
float __cdecl _handle_errorf(char *fname, int opcode, unsigned long long value, int type, int flags, int error, float arg1, float arg2, int nargs)
Definition: _handle_error.c:56
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
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
GLubyte GLubyte GLubyte GLubyte w
Definition: glext.h:6102
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
#define _DOMAIN
Definition: math.h:39
_Check_return_ float __cdecl fmodf(_In_ float x, _In_ float y)
Definition: math.h:228
#define AMD_F_INVALID
Definition: libm_new.h:86
#define INDEFBITPATT_SP32
Definition: libm_util.h:76
#define GET_BITS_SP32(x, ux)
Definition: libm_util.h:105
#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
float remainderf(float x, float y)
Definition: remainderf.c:60
#define _FUNCNAME
#define _OPERATION
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