?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
msvcrt: Import lround implementation from musl.
[wine.git] / dlls / msvcrt / math.c
blob018a8ae47e425ea7793b35d6bd94da9143dcb74b
1 /*
2 * msvcrt.dll math functions
3 *
4 * Copyright 2000 Jon Griffiths
5 *
6 * This library is free software; you can redistribute it and/or
7 * modify it under the terms of the GNU Lesser General Public
8 * License as published by the Free Software Foundation; either
9 * version 2.1 of the License, or (at your option) any later version.
10 *
11 * This library is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public
17 * License along with this library; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA
19 *
20 *
21 * For functions copied from musl libc (http://musl.libc.org/):
22 * ====================================================
23 * Copyright 2005-2020 Rich Felker, et al.
24 *
25 * Permission is hereby granted, free of charge, to any person obtaining
26 * a copy of this software and associated documentation files (the
27 * "Software"), to deal in the Software without restriction, including
28 * without limitation the rights to use, copy, modify, merge, publish,
29 * distribute, sublicense, and/or sell copies of the Software, and to
30 * permit persons to whom the Software is furnished to do so, subject to
31 * the following conditions:
32 *
33 * The above copyright notice and this permission notice shall be
34 * included in all copies or substantial portions of the Software.
35 * ====================================================
36 */
37
38 #include <complex.h>
39 #include <stdio.h>
40 #include <fenv.h>
41 #include <fpieee.h>
42 #include <limits.h>
43 #include <locale.h>
44 #include <math.h>
45
46 #include "msvcrt.h"
47 #include "winternl.h"
48 #include "unixlib.h"
49
50 #include "wine/asm.h"
51 #include "wine/debug.h"
52
53 WINE_DEFAULT_DEBUG_CHANNEL(msvcrt);
54
55 #undef div
56 #undef ldiv
57
58 #define _DOMAIN 1 /* domain error in argument */
59 #define _SING 2 /* singularity */
60 #define _OVERFLOW 3 /* range overflow */
61 #define _UNDERFLOW 4 /* range underflow */
62
63 typedef int (CDECL *MSVCRT_matherr_func)(struct _exception *);
64
65 static MSVCRT_matherr_func MSVCRT_default_matherr_func = NULL;
66
67 BOOL sse2_supported;
68 static BOOL sse2_enabled;
69
70 static const struct unix_funcs *unix_funcs;
71
72 void msvcrt_init_math( void *module )
73 {
74 sse2_supported = IsProcessorFeaturePresent( PF_XMMI64_INSTRUCTIONS_AVAILABLE );
75 #if _MSVCR_VER <=71
76 sse2_enabled = FALSE;
77 #else
78 sse2_enabled = sse2_supported;
79 #endif
80 __wine_init_unix_lib( module, DLL_PROCESS_ATTACH, NULL, &unix_funcs );
81 }
82
83 /* Copied from musl: src/internal/libm.h */
84 static inline float fp_barrierf(float x)
85 {
86 volatile float y = x;
87 return y;
88 }
89
90 static inline double CDECL ret_nan( BOOL update_sw )
91 {
92 double x = 1.0;
93 if (!update_sw) return -NAN;
94 return (x - x) / (x - x);
95 }
96
97 #define SET_X87_CW(MASK) \
98 "subl $4, %esp\n\t" \
99 __ASM_CFI(".cfi_adjust_cfa_offset 4\n\t") \
100 "fnstcw (%esp)\n\t" \
101 "movw (%esp), %ax\n\t" \
102 "movw %ax, 2(%esp)\n\t" \
103 "testw $" #MASK ", %ax\n\t" \
104 "jz 1f\n\t" \
105 "andw $~" #MASK ", %ax\n\t" \
106 "movw %ax, 2(%esp)\n\t" \
107 "fldcw 2(%esp)\n\t" \
108 "1:\n\t"
109
110 #define RESET_X87_CW \
111 "movw (%esp), %ax\n\t" \
112 "cmpw %ax, 2(%esp)\n\t" \
113 "je 1f\n\t" \
114 "fstpl 8(%esp)\n\t" \
115 "fldcw (%esp)\n\t" \
116 "fldl 8(%esp)\n\t" \
117 "fwait\n\t" \
118 "1:\n\t" \
119 "addl $4, %esp\n\t" \
120 __ASM_CFI(".cfi_adjust_cfa_offset -4\n\t")
121
122 /*********************************************************************
123 * _matherr (CRTDLL.@)
124 */
125 int CDECL _matherr(struct _exception *e)
126 {
127 return 0;
128 }
129
130
131 static double math_error(int type, const char *name, double arg1, double arg2, double retval)
132 {
133 struct _exception exception = {type, (char *)name, arg1, arg2, retval};
134
135 TRACE("(%d, %s, %g, %g, %g)\n", type, debugstr_a(name), arg1, arg2, retval);
136
137 if (MSVCRT_default_matherr_func && MSVCRT_default_matherr_func(&exception))
138 return exception.retval;
139
140 switch (type)
141 {
142 case 0:
143 /* don't set errno */
144 break;
145 case _DOMAIN:
146 *_errno() = EDOM;
147 break;
148 case _SING:
149 case _OVERFLOW:
150 *_errno() = ERANGE;
151 break;
152 case _UNDERFLOW:
153 /* don't set errno */
154 break;
155 default:
156 ERR("Unhandled math error!\n");
157 }
158
159 return exception.retval;
160 }
161
162 /*********************************************************************
163 * __setusermatherr (MSVCRT.@)
164 */
165 void CDECL __setusermatherr(MSVCRT_matherr_func func)
166 {
167 MSVCRT_default_matherr_func = func;
168 TRACE("new matherr handler %p\n", func);
169 }
170
171 /*********************************************************************
172 * _set_SSE2_enable (MSVCRT.@)
173 */
174 int CDECL _set_SSE2_enable(int flag)
175 {
176 sse2_enabled = flag && sse2_supported;
177 return sse2_enabled;
178 }
179
180 #if defined(_WIN64)
181 # if _MSVCR_VER>=140
182 /*********************************************************************
183 * _get_FMA3_enable (UCRTBASE.@)
184 */
185 int CDECL _get_FMA3_enable(void)
186 {
187 FIXME("() stub\n");
188 return 0;
189 }
190 # endif
191
192 # if _MSVCR_VER>=120
193 /*********************************************************************
194 * _set_FMA3_enable (MSVCR120.@)
195 */
196 int CDECL _set_FMA3_enable(int flag)
197 {
198 FIXME("(%x) stub\n", flag);
199 return 0;
200 }
201 # endif
202 #endif
203
204 #if !defined(__i386__) || _MSVCR_VER>=120
205
206 /*********************************************************************
207 * _chgsignf (MSVCRT.@)
208 */
209 float CDECL _chgsignf( float num )
210 {
211 union { float f; UINT32 i; } u = { num };
212 u.i ^= 0x80000000;
213 return u.f;
214 }
215
216 /*********************************************************************
217 * _copysignf (MSVCRT.@)
218 *
219 * Copied from musl: src/math/copysignf.c
220 */
221 float CDECL _copysignf( float x, float y )
222 {
223 union { float f; UINT32 i; } ux = { x }, uy = { y };
224 ux.i &= 0x7fffffff;
225 ux.i |= uy.i & 0x80000000;
226 return ux.f;
227 }
228
229 /*********************************************************************
230 * _nextafterf (MSVCRT.@)
231 */
232 float CDECL _nextafterf( float num, float next )
233 {
234 if (!isfinite(num) || !isfinite(next)) *_errno() = EDOM;
235 return unix_funcs->nextafterf( num, next );
236 }
237
238 /*********************************************************************
239 * _logbf (MSVCRT.@)
240 */
241 float CDECL _logbf( float num )
242 {
243 float ret = unix_funcs->logbf(num);
244 if (isnan(num)) return math_error(_DOMAIN, "_logbf", num, 0, ret);
245 if (!num) return math_error(_SING, "_logbf", num, 0, ret);
246 return ret;
247 }
248
249 #endif
250
251 #ifndef __i386__
252
253 /*********************************************************************
254 * _fpclassf (MSVCRT.@)
255 */
256 int CDECL _fpclassf( float num )
257 {
258 union { float f; UINT32 i; } u = { num };
259 int e = u.i >> 23 & 0xff;
260 int s = u.i >> 31;
261
262 switch (e)
263 {
264 case 0:
265 if (u.i << 1) return s ? _FPCLASS_ND : _FPCLASS_PD;
266 return s ? _FPCLASS_NZ : _FPCLASS_PZ;
267 case 0xff:
268 if (u.i << 9) return ((u.i >> 22) & 1) ? _FPCLASS_QNAN : _FPCLASS_SNAN;
269 return s ? _FPCLASS_NINF : _FPCLASS_PINF;
270 default:
271 return s ? _FPCLASS_NN : _FPCLASS_PN;
272 }
273 }
274
275 /*********************************************************************
276 * _finitef (MSVCRT.@)
277 */
278 int CDECL _finitef( float num )
279 {
280 union { float f; UINT32 i; } u = { num };
281 return (u.i & 0x7fffffff) < 0x7f800000;
282 }
283
284 /*********************************************************************
285 * _isnanf (MSVCRT.@)
286 */
287 int CDECL _isnanf( float num )
288 {
289 union { float f; UINT32 i; } u = { num };
290 return (u.i & 0x7fffffff) > 0x7f800000;
291 }
292
293 static float asinf_R(float z)
294 {
295 /* coefficients for R(x^2) */
296 static const float p1 = 1.66666672e-01,
297 p2 = -5.11644611e-02,
298 p3 = -1.21124933e-02,
299 p4 = -3.58742251e-03,
300 q1 = -7.56982703e-01;
301
302 float p, q;
303 p = z * (p1 + z * (p2 + z * (p3 + z * p4)));
304 q = 1.0f + z * q1;
305 return p / q;
306 }
307
308 /*********************************************************************
309 * acosf (MSVCRT.@)
310 *
311 * Copied from musl: src/math/acosf.c
312 */
313 float CDECL acosf( float x )
314 {
315 static const double pio2_lo = 6.12323399573676603587e-17;
316
317 float z, w, s, c, df;
318 unsigned int hx, ix;
319
320 hx = *(unsigned int*)&x;
321 ix = hx & 0x7fffffff;
322 /* |x| >= 1 or nan */
323 if (ix >= 0x3f800000) {
324 if (ix == 0x3f800000) {
325 if (hx >> 31)
326 return M_PI;
327 return 0;
328 }
329 if (isnan(x)) return x;
330 return math_error(_DOMAIN, "acosf", x, 0, 0 / (x - x));
331 }
332 /* |x| < 0.5 */
333 if (ix < 0x3f000000) {
334 if (ix <= 0x32800000) /* |x| < 2**-26 */
335 return M_PI_2;
336 return M_PI_2 - (x - (pio2_lo - x * asinf_R(x * x)));
337 }
338 /* x < -0.5 */
339 if (hx >> 31) {
340 z = (1 + x) * 0.5f;
341 s = sqrtf(z);
342 return M_PI - 2 * (s + ((double)s * asinf_R(z)));
343 }
344 /* x > 0.5 */
345 z = (1 - x) * 0.5f;
346 s = sqrtf(z);
347 hx = *(unsigned int*)&s & 0xffff0000;
348 df = *(float*)&hx;
349 c = (z - df * df) / (s + df);
350 w = asinf_R(z) * s + c;
351 return 2 * (df + w);
352 }
353
354 /*********************************************************************
355 * asinf (MSVCRT.@)
356 *
357 * Copied from musl: src/math/asinf.c
358 */
359 float CDECL asinf( float x )
360 {
361 static const double pio2 = 1.570796326794896558e+00;
362 static const float pio4_hi = 0.785398125648;
363 static const float pio2_lo = 7.54978941586e-08;
364
365 float s, z, f, c;
366 unsigned int hx, ix;
367
368 hx = *(unsigned int*)&x;
369 ix = hx & 0x7fffffff;
370 if (ix >= 0x3f800000) { /* |x| >= 1 */
371 if (ix == 0x3f800000) /* |x| == 1 */
372 return x * pio2 + 7.5231638453e-37; /* asin(+-1) = +-pi/2 with inexact */
373 if (isnan(x)) return x;
374 return math_error(_DOMAIN, "asinf", x, 0, 0 / (x - x));
375 }
376 if (ix < 0x3f000000) { /* |x| < 0.5 */
377 /* if 0x1p-126 <= |x| < 0x1p-12, avoid raising underflow */
378 if (ix < 0x39800000 && ix >= 0x00800000)
379 return x;
380 return x + x * asinf_R(x * x);
381 }
382 /* 1 > |x| >= 0.5 */
383 z = (1 - fabsf(x)) * 0.5f;
384 s = sqrtf(z);
385 /* f+c = sqrt(z) */
386 *(unsigned int*)&f = *(unsigned int*)&s & 0xffff0000;
387 c = (z - f * f) / (s + f);
388 x = pio4_hi - (2 * s * asinf_R(z) - (pio2_lo - 2 * c) - (pio4_hi - 2 * f));
389 if (hx >> 31)
390 return -x;
391 return x;
392 }
393
394 /*********************************************************************
395 * atanf (MSVCRT.@)
396 *
397 * Copied from musl: src/math/atanf.c
398 */
399 float CDECL atanf( float x )
400 {
401 static const float atanhi[] = {
402 4.6364760399e-01,
403 7.8539812565e-01,
404 9.8279368877e-01,
405 1.5707962513e+00,
406 };
407 static const float atanlo[] = {
408 5.0121582440e-09,
409 3.7748947079e-08,
410 3.4473217170e-08,
411 7.5497894159e-08,
412 };
413 static const float aT[] = {
414 3.3333328366e-01,
415 -1.9999158382e-01,
416 1.4253635705e-01,
417 -1.0648017377e-01,
418 6.1687607318e-02,
419 };
420
421 float w, s1, s2, z;
422 unsigned int ix, sign;
423 int id;
424
425 #if _MSVCR_VER == 0
426 if (isnan(x)) return math_error(_DOMAIN, "atanf", x, 0, x);
427 #endif
428
429 ix = *(unsigned int*)&x;
430 sign = ix >> 31;
431 ix &= 0x7fffffff;
432 if (ix >= 0x4c800000) { /* if |x| >= 2**26 */
433 if (isnan(x))
434 return x;
435 z = atanhi[3] + 7.5231638453e-37;
436 return sign ? -z : z;
437 }
438 if (ix < 0x3ee00000) { /* |x| < 0.4375 */
439 if (ix < 0x39800000) { /* |x| < 2**-12 */
440 if (ix < 0x00800000)
441 /* raise underflow for subnormal x */
442 fp_barrierf(x*x);
443 return x;
444 }
445 id = -1;
446 } else {
447 x = fabsf(x);
448 if (ix < 0x3f980000) { /* |x| < 1.1875 */
449 if (ix < 0x3f300000) { /* 7/16 <= |x| < 11/16 */
450 id = 0;
451 x = (2.0f * x - 1.0f) / (2.0f + x);
452 } else { /* 11/16 <= |x| < 19/16 */
453 id = 1;
454 x = (x - 1.0f) / (x + 1.0f);
455 }
456 } else {
457 if (ix < 0x401c0000) { /* |x| < 2.4375 */
458 id = 2;
459 x = (x - 1.5f) / (1.0f + 1.5f * x);
460 } else { /* 2.4375 <= |x| < 2**26 */
461 id = 3;
462 x = -1.0f / x;
463 }
464 }
465 }
466 /* end of argument reduction */
467 z = x * x;
468 w = z * z;
469 /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
470 s1 = z * (aT[0] + w * (aT[2] + w * aT[4]));
471 s2 = w * (aT[1] + w * aT[3]);
472 if (id < 0)
473 return x - x * (s1 + s2);
474 z = atanhi[id] - ((x * (s1 + s2) - atanlo[id]) - x);
475 return sign ? -z : z;
476 }
477
478 /*********************************************************************
479 * atan2f (MSVCRT.@)
480 *
481 * Copied from musl: src/math/atan2f.c
482 */
483 float CDECL atan2f( float y, float x )
484 {
485 static const float pi = 3.1415927410e+00,
486 pi_lo = -8.7422776573e-08;
487
488 float z;
489 unsigned int m, ix, iy;
490
491 if (isnan(x) || isnan(y))
492 return x + y;
493 ix = *(unsigned int*)&x;
494 iy = *(unsigned int*)&y;
495 if (ix == 0x3f800000) /* x=1.0 */
496 return atanf(y);
497 m = ((iy >> 31) & 1) | ((ix >> 30) & 2); /* 2*sign(x)+sign(y) */
498 ix &= 0x7fffffff;
499 iy &= 0x7fffffff;
500
501 /* when y = 0 */
502 if (iy == 0) {
503 switch (m) {
504 case 0:
505 case 1: return y; /* atan(+-0,+anything)=+-0 */
506 case 2: return pi; /* atan(+0,-anything) = pi */
507 case 3: return -pi; /* atan(-0,-anything) =-pi */
508 }
509 }
510 /* when x = 0 */
511 if (ix == 0)
512 return m & 1 ? -pi / 2 : pi / 2;
513 /* when x is INF */
514 if (ix == 0x7f800000) {
515 if (iy == 0x7f800000) {
516 switch (m) {
517 case 0: return pi / 4; /* atan(+INF,+INF) */
518 case 1: return -pi / 4; /* atan(-INF,+INF) */
519 case 2: return 3 * pi / 4; /*atan(+INF,-INF)*/
520 case 3: return -3 * pi / 4; /*atan(-INF,-INF)*/
521 }
522 } else {
523 switch (m) {
524 case 0: return 0.0f; /* atan(+...,+INF) */
525 case 1: return -0.0f; /* atan(-...,+INF) */
526 case 2: return pi; /* atan(+...,-INF) */
527 case 3: return -pi; /* atan(-...,-INF) */
528 }
529 }
530 }
531 /* |y/x| > 0x1p26 */
532 if (ix + (26 << 23) < iy || iy == 0x7f800000)
533 return m & 1 ? -pi / 2 : pi / 2;
534
535 /* z = atan(|y/x|) with correct underflow */
536 if ((m & 2) && iy + (26 << 23) < ix) /*|y/x| < 0x1p-26, x < 0 */
537 z = 0.0;
538 else
539 z = atanf(fabsf(y / x));
540 switch (m) {
541 case 0: return z; /* atan(+,+) */
542 case 1: return -z; /* atan(-,+) */
543 case 2: return pi - (z - pi_lo); /* atan(+,-) */
544 default: /* case 3 */
545 return (z - pi_lo) - pi; /* atan(-,-) */
546 }
547 }
548
549 /*********************************************************************
550 * cosf (MSVCRT.@)
551 */
552 float CDECL cosf( float x )
553 {
554 float ret = unix_funcs->cosf( x );
555 if (!isfinite(x)) return math_error(_DOMAIN, "cosf", x, 0, ret);
556 return ret;
557 }
558
559 /*********************************************************************
560 * coshf (MSVCRT.@)
561 */
562 float CDECL coshf( float x )
563 {
564 float ret = unix_funcs->coshf( x );
565 if (isnan(x)) return math_error(_DOMAIN, "coshf", x, 0, ret);
566 return ret;
567 }
568
569 /*********************************************************************
570 * expf (MSVCRT.@)
571 */
572 float CDECL expf( float x )
573 {
574 float ret = unix_funcs->expf( x );
575 if (isnan(x)) return math_error(_DOMAIN, "expf", x, 0, ret);
576 if (isfinite(x) && !ret) return math_error(_UNDERFLOW, "expf", x, 0, ret);
577 if (isfinite(x) && !isfinite(ret)) return math_error(_OVERFLOW, "expf", x, 0, ret);
578 return ret;
579 }
580
581 /*********************************************************************
582 * fmodf (MSVCRT.@)
583 */
584 float CDECL fmodf( float x, float y )
585 {
586 float ret = unix_funcs->fmodf( x, y );
587 if (!isfinite(x) || !isfinite(y)) return math_error(_DOMAIN, "fmodf", x, 0, ret);
588 return ret;
589 }
590
591 /*********************************************************************
592 * logf (MSVCRT.@)
593 */
594 float CDECL logf( float x )
595 {
596 float ret = unix_funcs->logf( x );
597 if (x < 0.0) return math_error(_DOMAIN, "logf", x, 0, ret);
598 if (x == 0.0) return math_error(_SING, "logf", x, 0, ret);
599 return ret;
600 }
601
602 /*********************************************************************
603 * log10f (MSVCRT.@)
604 */
605 float CDECL log10f( float x )
606 {
607 float ret = unix_funcs->log10f( x );
608 if (x < 0.0) return math_error(_DOMAIN, "log10f", x, 0, ret);
609 if (x == 0.0) return math_error(_SING, "log10f", x, 0, ret);
610 return ret;
611 }
612
613 /*********************************************************************
614 * powf (MSVCRT.@)
615 */
616 float CDECL powf( float x, float y )
617 {
618 float z = unix_funcs->powf(x,y);
619 if (x < 0 && y != floorf(y)) return math_error(_DOMAIN, "powf", x, y, z);
620 if (!x && isfinite(y) && y < 0) return math_error(_SING, "powf", x, y, z);
621 if (isfinite(x) && isfinite(y) && !isfinite(z)) return math_error(_OVERFLOW, "powf", x, y, z);
622 if (x && isfinite(x) && isfinite(y) && !z) return math_error(_UNDERFLOW, "powf", x, y, z);
623 return z;
624 }
625
626 /*********************************************************************
627 * sinf (MSVCRT.@)
628 */
629 float CDECL sinf( float x )
630 {
631 float ret = unix_funcs->sinf( x );
632 if (!isfinite(x)) return math_error(_DOMAIN, "sinf", x, 0, ret);
633 return ret;
634 }
635
636 /*********************************************************************
637 * sinhf (MSVCRT.@)
638 */
639 float CDECL sinhf( float x )
640 {
641 float ret = unix_funcs->sinhf( x );
642 if (isnan(x)) return math_error(_DOMAIN, "sinhf", x, 0, ret);
643 return ret;
644 }
645
646 static BOOL sqrtf_validate( float *x )
647 {
648 short c = _fdclass(*x);
649
650 if (c == FP_ZERO) return FALSE;
651 if (c == FP_NAN) return FALSE;
652 if (signbit(*x))
653 {
654 *x = math_error(_DOMAIN, "sqrtf", *x, 0, ret_nan(TRUE));
655 return FALSE;
656 }
657 if (c == FP_INFINITE) return FALSE;
658 return TRUE;
659 }
660
661 #if defined(__x86_64__) || defined(__i386__)
662 float CDECL sse2_sqrtf(float);
663 __ASM_GLOBAL_FUNC( sse2_sqrtf,
664 "sqrtss %xmm0, %xmm0\n\t"
665 "ret" )
666 #endif
667
668 /*********************************************************************
669 * sqrtf (MSVCRT.@)
670 *
671 * Copied from musl: src/math/sqrtf.c
672 */
673 float CDECL sqrtf( float x )
674 {
675 #ifdef __x86_64__
676 if (!sqrtf_validate(&x))
677 return x;
678
679 return sse2_sqrtf(x);
680 #else
681 static const float tiny = 1.0e-30;
682
683 float z;
684 int ix,s,q,m,t,i;
685 unsigned int r;
686
687 ix = *(int*)&x;
688
689 if (!sqrtf_validate(&x))
690 return x;
691
692 /* normalize x */
693 m = ix >> 23;
694 if (m == 0) { /* subnormal x */
695 for (i = 0; (ix & 0x00800000) == 0; i++)
696 ix <<= 1;
697 m -= i - 1;
698 }
699 m -= 127; /* unbias exponent */
700 ix = (ix & 0x007fffff) | 0x00800000;
701 if (m & 1) /* odd m, double x to make it even */
702 ix += ix;
703 m >>= 1; /* m = [m/2] */
704
705 /* generate sqrt(x) bit by bit */
706 ix += ix;
707 q = s = 0; /* q = sqrt(x) */
708 r = 0x01000000; /* r = moving bit from right to left */
709
710 while (r != 0) {
711 t = s + r;
712 if (t <= ix) {
713 s = t + r;
714 ix -= t;
715 q += r;
716 }
717 ix += ix;
718 r >>= 1;
719 }
720
721 /* use floating add to find out rounding direction */
722 if (ix != 0) {
723 z = 1.0f - tiny; /* raise inexact flag */
724 if (z >= 1.0f) {
725 z = 1.0f + tiny;
726 if (z > 1.0f)
727 q += 2;
728 else
729 q += q & 1;
730 }
731 }
732 ix = (q >> 1) + 0x3f000000;
733 r = ix + ((unsigned int)m << 23);
734 z = *(float*)&r;
735 return z;
736 #endif
737 }
738
739 /*********************************************************************
740 * tanf (MSVCRT.@)
741 */
742 float CDECL tanf( float x )
743 {
744 float ret = unix_funcs->tanf(x);
745 if (!isfinite(x)) return math_error(_DOMAIN, "tanf", x, 0, ret);
746 return ret;
747 }
748
749 /*********************************************************************
750 * tanhf (MSVCRT.@)
751 */
752 float CDECL tanhf( float x )
753 {
754 float ret = unix_funcs->tanhf(x);
755 if (!isfinite(x)) return math_error(_DOMAIN, "tanhf", x, 0, ret);
756 return ret;
757 }
758
759 /*********************************************************************
760 * ceilf (MSVCRT.@)
761 */
762 float CDECL ceilf( float x )
763 {
764 return unix_funcs->ceilf(x);
765 }
766
767 /*********************************************************************
768 * floorf (MSVCRT.@)
769 */
770 float CDECL floorf( float x )
771 {
772 return unix_funcs->floorf(x);
773 }
774
775 /*********************************************************************
776 * frexpf (MSVCRT.@)
777 */
778 float CDECL frexpf( float x, int *exp )
779 {
780 return unix_funcs->frexpf( x, exp );
781 }
782
783 /*********************************************************************
784 * modff (MSVCRT.@)
785 */
786 float CDECL modff( float x, float *iptr )
787 {
788 return unix_funcs->modff( x, iptr );
789 }
790
791 #endif
792
793 #if !defined(__i386__) && !defined(__x86_64__) && (_MSVCR_VER == 0 || _MSVCR_VER >= 110)
794
795 /*********************************************************************
796 * fabsf (MSVCRT.@)
797 *
798 * Copied from musl: src/math/fabsf.c
799 */
800 float CDECL fabsf( float x )
801 {
802 union { float f; UINT32 i; } u = { x };
803 u.i &= 0x7fffffff;
804 return u.f;
805 }
806
807 #endif
808
809 /*********************************************************************
810 * acos (MSVCRT.@)
811 *
812 * Copied from musl: src/math/acos.c
813 */
814 static double acos_R(double z)
815 {
816 static const double pS0 = 1.66666666666666657415e-01,
817 pS1 = -3.25565818622400915405e-01,
818 pS2 = 2.01212532134862925881e-01,
819 pS3 = -4.00555345006794114027e-02,
820 pS4 = 7.91534994289814532176e-04,
821 pS5 = 3.47933107596021167570e-05,
822 qS1 = -2.40339491173441421878e+00,
823 qS2 = 2.02094576023350569471e+00,
824 qS3 = -6.88283971605453293030e-01,
825 qS4 = 7.70381505559019352791e-02;
826
827 double p, q;
828 p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * pS5)))));
829 q = 1.0 + z * (qS1 + z * (qS2 + z * (qS3 + z * qS4)));
830 return p/q;
831 }
832
833 double CDECL acos( double x )
834 {
835 static const double pio2_hi = 1.57079632679489655800e+00,
836 pio2_lo = 6.12323399573676603587e-17;
837
838 double z, w, s, c, df;
839 unsigned int hx, ix;
840 ULONGLONG llx;
841
842 hx = *(ULONGLONG*)&x >> 32;
843 ix = hx & 0x7fffffff;
844 /* |x| >= 1 or nan */
845 if (ix >= 0x3ff00000) {
846 unsigned int lx;
847
848 lx = *(ULONGLONG*)&x;
849 if (((ix - 0x3ff00000) | lx) == 0) {
850 /* acos(1)=0, acos(-1)=pi */
851 if (hx >> 31)
852 return 2 * pio2_hi + 7.5231638452626401e-37;
853 return 0;
854 }
855 if (isnan(x)) return x;
856 return math_error(_DOMAIN, "acos", x, 0, 0 / (x - x));
857 }
858 /* |x| < 0.5 */
859 if (ix < 0x3fe00000) {
860 if (ix <= 0x3c600000) /* |x| < 2**-57 */
861 return pio2_hi + 7.5231638452626401e-37;
862 return pio2_hi - (x - (pio2_lo - x * acos_R(x * x)));
863 }
864 /* x < -0.5 */
865 if (hx >> 31) {
866 z = (1.0 + x) * 0.5;
867 s = sqrt(z);
868 w = acos_R(z) * s - pio2_lo;
869 return 2 * (pio2_hi - (s + w));
870 }
871 /* x > 0.5 */
872 z = (1.0 - x) * 0.5;
873 s = sqrt(z);
874 df = s;
875 llx = (*(ULONGLONG*)&df >> 32) << 32;
876 df = *(double*)&llx;
877 c = (z - df * df) / (s + df);
878 w = acos_R(z) * s + c;
879 return 2 * (df + w);
880 }
881
882 /*********************************************************************
883 * asin (MSVCRT.@)
884 *
885 * Copied from musl: src/math/asin.c
886 */
887 static double asin_R(double z)
888 {
889 /* coefficients for R(x^2) */
890 static const double pS0 = 1.66666666666666657415e-01,
891 pS1 = -3.25565818622400915405e-01,
892 pS2 = 2.01212532134862925881e-01,
893 pS3 = -4.00555345006794114027e-02,
894 pS4 = 7.91534994289814532176e-04,
895 pS5 = 3.47933107596021167570e-05,
896 qS1 = -2.40339491173441421878e+00,
897 qS2 = 2.02094576023350569471e+00,
898 qS3 = -6.88283971605453293030e-01,
899 qS4 = 7.70381505559019352791e-02;
900
901 double p, q;
902 p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * pS5)))));
903 q = 1.0 + z * (qS1 + z * (qS2 + z * (qS3 + z * qS4)));
904 return p / q;
905 }
906
907 #ifdef __i386__
908 double CDECL x87_asin(double);
909 __ASM_GLOBAL_FUNC( x87_asin,
910 "fldl 4(%esp)\n\t"
911 SET_X87_CW(~0x37f)
912 "fld %st\n\t"
913 "fld1\n\t"
914 "fsubp\n\t"
915 "fld1\n\t"
916 "fadd %st(2)\n\t"
917 "fmulp\n\t"
918 "fsqrt\n\t"
919 "fpatan\n\t"
920 RESET_X87_CW
921 "ret" )
922 #endif
923
924 double CDECL asin( double x )
925 {
926 static const double pio2_hi = 1.57079632679489655800e+00,
927 pio2_lo = 6.12323399573676603587e-17;
928
929 double z, r, s;
930 unsigned int hx, ix;
931 ULONGLONG llx;
932 #ifdef __i386__
933 unsigned int x87_cw, sse2_cw;
934 #endif
935
936 hx = *(ULONGLONG*)&x >> 32;
937 ix = hx & 0x7fffffff;
938 /* |x| >= 1 or nan */
939 if (ix >= 0x3ff00000) {
940 unsigned int lx;
941 lx = *(ULONGLONG*)&x;
942 if (((ix - 0x3ff00000) | lx) == 0)
943 /* asin(1) = +-pi/2 with inexact */
944 return x * pio2_hi + 7.5231638452626401e-37;
945 if (isnan(x))
946 {
947 #ifdef __i386__
948 return math_error(_DOMAIN, "asin", x, 0, x);
949 #else
950 return x;
951 #endif
952 }
953 return math_error(_DOMAIN, "asin", x, 0, 0 / (x - x));
954 }
955
956 #ifdef __i386__
957 __control87_2(0, 0, &x87_cw, &sse2_cw);
958 if (!sse2_enabled || (x87_cw & _MCW_EM) != _MCW_EM
959 || (sse2_cw & (_MCW_EM | _MCW_RC)) != _MCW_EM)
960 return x87_asin(x);
961 #endif
962
963 /* |x| < 0.5 */
964 if (ix < 0x3fe00000) {
965 /* if 0x1p-1022 <= |x| < 0x1p-26, avoid raising underflow */
966 if (ix < 0x3e500000 && ix >= 0x00100000)
967 return x;
968 return x + x * asin_R(x * x);
969 }
970 /* 1 > |x| >= 0.5 */
971 z = (1 - fabs(x)) * 0.5;
972 s = sqrt(z);
973 r = asin_R(z);
974 if (ix >= 0x3fef3333) { /* if |x| > 0.975 */
975 x = pio2_hi - (2 * (s + s * r) - pio2_lo);
976 } else {
977 double f, c;
978 /* f+c = sqrt(z) */
979 f = s;
980 llx = (*(ULONGLONG*)&f >> 32) << 32;
981 f = *(double*)&llx;
982 c = (z - f * f) / (s + f);
983 x = 0.5 * pio2_hi - (2 * s * r - (pio2_lo - 2 * c) - (0.5 * pio2_hi - 2 * f));
984 }
985 if (hx >> 31)
986 return -x;
987 return x;
988 }
989
990 /*********************************************************************
991 * atan (MSVCRT.@)
992 *
993 * Copied from musl: src/math/atan.c
994 */
995 double CDECL atan( double x )
996 {
997 static const double atanhi[] = {
998 4.63647609000806093515e-01,
999 7.85398163397448278999e-01,
1000 9.82793723247329054082e-01,
1001 1.57079632679489655800e+00,
1002 };
1003 static const double atanlo[] = {
1004 2.26987774529616870924e-17,
1005 3.06161699786838301793e-17,
1006 1.39033110312309984516e-17,
1007 6.12323399573676603587e-17,
1008 };
1009 static const double aT[] = {
1010 3.33333333333329318027e-01,
1011 -1.99999999998764832476e-01,
1012 1.42857142725034663711e-01,
1013 -1.11111104054623557880e-01,
1014 9.09088713343650656196e-02,
1015 -7.69187620504482999495e-02,
1016 6.66107313738753120669e-02,
1017 -5.83357013379057348645e-02,
1018 4.97687799461593236017e-02,
1019 -3.65315727442169155270e-02,
1020 1.62858201153657823623e-02,
1021 };
1022
1023 double w, s1, s2, z;
1024 unsigned int ix, sign;
1025 int id;
1026
1027 #if _MSVCR_VER == 0
1028 if (isnan(x)) return math_error(_DOMAIN, "atan", x, 0, x);
1029 #endif
1030
1031 ix = *(ULONGLONG*)&x >> 32;
1032 sign = ix >> 31;
1033 ix &= 0x7fffffff;
1034 if (ix >= 0x44100000) { /* if |x| >= 2^66 */
1035 if (isnan(x))
1036 return x;
1037 z = atanhi[3] + 7.5231638452626401e-37;
1038 return sign ? -z : z;
1039 }
1040 if (ix < 0x3fdc0000) { /* |x| < 0.4375 */
1041 if (ix < 0x3e400000) { /* |x| < 2^-27 */
1042 if (ix < 0x00100000)
1043 /* raise underflow for subnormal x */
1044 fp_barrierf((float)x);
1045 return x;
1046 }
1047 id = -1;
1048 } else {
1049 x = fabs(x);
1050 if (ix < 0x3ff30000) { /* |x| < 1.1875 */
1051 if (ix < 0x3fe60000) { /* 7/16 <= |x| < 11/16 */
1052 id = 0;
1053 x = (2.0 * x - 1.0) / (2.0 + x);
1054 } else { /* 11/16 <= |x| < 19/16 */
1055 id = 1;
1056 x = (x - 1.0) / (x + 1.0);
1057 }
1058 } else {
1059 if (ix < 0x40038000) { /* |x| < 2.4375 */
1060 id = 2;
1061 x = (x - 1.5) / (1.0 + 1.5 * x);
1062 } else { /* 2.4375 <= |x| < 2^66 */
1063 id = 3;
1064 x = -1.0 / x;
1065 }
1066 }
1067 }
1068 /* end of argument reduction */
1069 z = x * x;
1070 w = z * z;
1071 /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
1072 s1 = z * (aT[0] + w * (aT[2] + w * (aT[4] + w * (aT[6] + w * (aT[8] + w * aT[10])))));
1073 s2 = w * (aT[1] + w * (aT[3] + w * (aT[5] + w * (aT[7] + w * aT[9]))));
1074 if (id < 0)
1075 return x - x * (s1 + s2);
1076 z = atanhi[id] - (x * (s1 + s2) - atanlo[id] - x);
1077 return sign ? -z : z;
1078 }
1079
1080 /*********************************************************************
1081 * atan2 (MSVCRT.@)
1082 *
1083 * Copied from musl: src/math/atan2.c
1084 */
1085 double CDECL atan2( double y, double x )
1086 {
1087 static const double pi = 3.1415926535897931160E+00,
1088 pi_lo = 1.2246467991473531772E-16;
1089
1090 double z;
1091 unsigned int m, lx, ly, ix, iy;
1092
1093 if (isnan(x) || isnan(y))
1094 return x+y;
1095 ix = *(ULONGLONG*)&x >> 32;
1096 lx = *(ULONGLONG*)&x;
1097 iy = *(ULONGLONG*)&y >> 32;
1098 ly = *(ULONGLONG*)&y;
1099 if (((ix - 0x3ff00000) | lx) == 0) /* x = 1.0 */
1100 return atan(y);
1101 m = ((iy >> 31) & 1) | ((ix >> 30) & 2); /* 2*sign(x)+sign(y) */
1102 ix = ix & 0x7fffffff;
1103 iy = iy & 0x7fffffff;
1104
1105 /* when y = 0 */
1106 if ((iy | ly) == 0) {
1107 switch(m) {
1108 case 0:
1109 case 1: return y; /* atan(+-0,+anything)=+-0 */
1110 case 2: return pi; /* atan(+0,-anything) = pi */
1111 case 3: return -pi; /* atan(-0,-anything) =-pi */
1112 }
1113 }
1114 /* when x = 0 */
1115 if ((ix | lx) == 0)
1116 return m & 1 ? -pi / 2 : pi / 2;
1117 /* when x is INF */
1118 if (ix == 0x7ff00000) {
1119 if (iy == 0x7ff00000) {
1120 switch(m) {
1121 case 0: return pi / 4; /* atan(+INF,+INF) */
1122 case 1: return -pi / 4; /* atan(-INF,+INF) */
1123 case 2: return 3 * pi / 4; /* atan(+INF,-INF) */
1124 case 3: return -3 * pi / 4; /* atan(-INF,-INF) */
1125 }
1126 } else {
1127 switch(m) {
1128 case 0: return 0.0; /* atan(+...,+INF) */
1129 case 1: return -0.0; /* atan(-...,+INF) */
1130 case 2: return pi; /* atan(+...,-INF) */
1131 case 3: return -pi; /* atan(-...,-INF) */
1132 }
1133 }
1134 }
1135 /* |y/x| > 0x1p64 */
1136 if (ix + (64 << 20) < iy || iy == 0x7ff00000)
1137 return m & 1 ? -pi / 2 : pi / 2;
1138
1139 /* z = atan(|y/x|) without spurious underflow */
1140 if ((m & 2) && iy + (64 << 20) < ix) /* |y/x| < 0x1p-64, x<0 */
1141 z = 0;
1142 else
1143 z = atan(fabs(y / x));
1144 switch (m) {
1145 case 0: return z; /* atan(+,+) */
1146 case 1: return -z; /* atan(-,+) */
1147 case 2: return pi - (z - pi_lo); /* atan(+,-) */
1148 default: /* case 3 */
1149 return (z - pi_lo) - pi; /* atan(-,-) */
1150 }
1151 }
1152
1153 /*********************************************************************
1154 * cos (MSVCRT.@)
1155 */
1156 double CDECL cos( double x )
1157 {
1158 double ret = unix_funcs->cos( x );
1159 if (!isfinite(x)) return math_error(_DOMAIN, "cos", x, 0, ret);
1160 return ret;
1161 }
1162
1163 /*********************************************************************
1164 * cosh (MSVCRT.@)
1165 */
1166 double CDECL cosh( double x )
1167 {
1168 double ret = unix_funcs->cosh( x );
1169 if (isnan(x)) return math_error(_DOMAIN, "cosh", x, 0, ret);
1170 return ret;
1171 }
1172
1173 /*********************************************************************
1174 * exp (MSVCRT.@)
1175 */
1176 double CDECL exp( double x )
1177 {
1178 double ret = unix_funcs->exp( x );
1179 if (isnan(x)) return math_error(_DOMAIN, "exp", x, 0, ret);
1180 if (isfinite(x) && !ret) return math_error(_UNDERFLOW, "exp", x, 0, ret);
1181 if (isfinite(x) && !isfinite(ret)) return math_error(_OVERFLOW, "exp", x, 0, ret);
1182 return ret;
1183 }
1184
1185 /*********************************************************************
1186 * fmod (MSVCRT.@)
1187 */
1188 double CDECL fmod( double x, double y )
1189 {
1190 double ret = unix_funcs->fmod( x, y );
1191 if (!isfinite(x) || !isfinite(y)) return math_error(_DOMAIN, "fmod", x, y, ret);
1192 return ret;
1193 }
1194
1195 /*********************************************************************
1196 * log (MSVCRT.@)
1197 */
1198 double CDECL log( double x )
1199 {
1200 double ret = unix_funcs->log( x );
1201 if (x < 0.0) return math_error(_DOMAIN, "log", x, 0, ret);
1202 if (x == 0.0) return math_error(_SING, "log", x, 0, ret);
1203 return ret;
1204 }
1205
1206 /*********************************************************************
1207 * log10 (MSVCRT.@)
1208 */
1209 double CDECL log10( double x )
1210 {
1211 double ret = unix_funcs->log10( x );
1212 if (x < 0.0) return math_error(_DOMAIN, "log10", x, 0, ret);
1213 if (x == 0.0) return math_error(_SING, "log10", x, 0, ret);
1214 return ret;
1215 }
1216
1217 /*********************************************************************
1218 * pow (MSVCRT.@)
1219 */
1220 double CDECL pow( double x, double y )
1221 {
1222 double z = unix_funcs->pow(x,y);
1223 if (x < 0 && y != floor(y))
1224 return math_error(_DOMAIN, "pow", x, y, z);
1225 if (!x && isfinite(y) && y < 0)
1226 return math_error(_SING, "pow", x, y, z);
1227 if (isfinite(x) && isfinite(y) && !isfinite(z))
1228 return math_error(_OVERFLOW, "pow", x, y, z);
1229 if (x && isfinite(x) && isfinite(y) && !z)
1230 return math_error(_UNDERFLOW, "pow", x, y, z);
1231 return z;
1232 }
1233
1234 /*********************************************************************
1235 * sin (MSVCRT.@)
1236 */
1237 double CDECL sin( double x )
1238 {
1239 double ret = unix_funcs->sin( x );
1240 if (!isfinite(x)) return math_error(_DOMAIN, "sin", x, 0, ret);
1241 return ret;
1242 }
1243
1244 /*********************************************************************
1245 * sinh (MSVCRT.@)
1246 */
1247 double CDECL sinh( double x )
1248 {
1249 double ret = unix_funcs->sinh( x );
1250 if (isnan(x)) return math_error(_DOMAIN, "sinh", x, 0, ret);
1251 return ret;
1252 }
1253
1254 static BOOL sqrt_validate( double *x, BOOL update_sw )
1255 {
1256 short c = _dclass(*x);
1257
1258 if (c == FP_ZERO) return FALSE;
1259 if (c == FP_NAN)
1260 {
1261 #ifdef __i386__
1262 if (update_sw)
1263 *x = math_error(_DOMAIN, "sqrt", *x, 0, *x);
1264 #else
1265 /* set signaling bit */
1266 *(ULONGLONG*)x |= 0x8000000000000ULL;
1267 #endif
1268 return FALSE;
1269 }
1270 if (signbit(*x))
1271 {
1272 *x = math_error(_DOMAIN, "sqrt", *x, 0, ret_nan(update_sw));
1273 return FALSE;
1274 }
1275 if (c == FP_INFINITE) return FALSE;
1276 return TRUE;
1277 }
1278
1279 #if defined(__x86_64__) || defined(__i386__)
1280 double CDECL sse2_sqrt(double);
1281 __ASM_GLOBAL_FUNC( sse2_sqrt,
1282 "sqrtsd %xmm0, %xmm0\n\t"
1283 "ret" )
1284 #endif
1285
1286 #ifdef __i386__
1287 double CDECL x87_sqrt(double);
1288 __ASM_GLOBAL_FUNC( x87_sqrt,
1289 "fldl 4(%esp)\n\t"
1290 SET_X87_CW(0xc00)
1291 "fsqrt\n\t"
1292 RESET_X87_CW
1293 "ret" )
1294 #endif
1295
1296 /*********************************************************************
1297 * sqrt (MSVCRT.@)
1298 *
1299 * Copied from musl: src/math/sqrt.c
1300 */
1301 double CDECL sqrt( double x )
1302 {
1303 #ifdef __x86_64__
1304 if (!sqrt_validate(&x, TRUE))
1305 return x;
1306
1307 return sse2_sqrt(x);
1308 #elif defined( __i386__ )
1309 if (!sqrt_validate(&x, TRUE))
1310 return x;
1311
1312 return x87_sqrt(x);
1313 #else
1314 static const double tiny = 1.0e-300;
1315
1316 double z;
1317 int sign = 0x80000000;
1318 int ix0,s0,q,m,t,i;
1319 unsigned int r,t1,s1,ix1,q1;
1320 ULONGLONG ix;
1321
1322 if (!sqrt_validate(&x, TRUE))
1323 return x;
1324
1325 ix = *(ULONGLONG*)&x;
1326 ix0 = ix >> 32;
1327 ix1 = ix;
1328
1329 /* normalize x */
1330 m = ix0 >> 20;
1331 if (m == 0) { /* subnormal x */
1332 while (ix0 == 0) {
1333 m -= 21;
1334 ix0 |= (ix1 >> 11);
1335 ix1 <<= 21;
1336 }
1337 for (i=0; (ix0 & 0x00100000) == 0; i++)
1338 ix0 <<= 1;
1339 m -= i - 1;
1340 ix0 |= ix1 >> (32 - i);
1341 ix1 <<= i;
1342 }
1343 m -= 1023; /* unbias exponent */
1344 ix0 = (ix0 & 0x000fffff) | 0x00100000;
1345 if (m & 1) { /* odd m, double x to make it even */
1346 ix0 += ix0 + ((ix1 & sign) >> 31);
1347 ix1 += ix1;
1348 }
1349 m >>= 1; /* m = [m/2] */
1350
1351 /* generate sqrt(x) bit by bit */
1352 ix0 += ix0 + ((ix1 & sign) >> 31);
1353 ix1 += ix1;
1354 q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
1355 r = 0x00200000; /* r = moving bit from right to left */
1356
1357 while (r != 0) {
1358 t = s0 + r;
1359 if (t <= ix0) {
1360 s0 = t + r;
1361 ix0 -= t;
1362 q += r;
1363 }
1364 ix0 += ix0 + ((ix1 & sign) >> 31);
1365 ix1 += ix1;
1366 r >>= 1;
1367 }
1368
1369 r = sign;
1370 while (r != 0) {
1371 t1 = s1 + r;
1372 t = s0;
1373 if (t < ix0 || (t == ix0 && t1 <= ix1)) {
1374 s1 = t1 + r;
1375 if ((t1&sign) == sign && (s1 & sign) == 0)
1376 s0++;
1377 ix0 -= t;
1378 if (ix1 < t1)
1379 ix0--;
1380 ix1 -= t1;
1381 q1 += r;
1382 }
1383 ix0 += ix0 + ((ix1 & sign) >> 31);
1384 ix1 += ix1;
1385 r >>= 1;
1386 }
1387
1388 /* use floating add to find out rounding direction */
1389 if ((ix0 | ix1) != 0) {
1390 z = 1.0 - tiny; /* raise inexact flag */
1391 if (z >= 1.0) {
1392 z = 1.0 + tiny;
1393 if (q1 == (unsigned int)0xffffffff) {
1394 q1 = 0;
1395 q++;
1396 } else if (z > 1.0) {
1397 if (q1 == (unsigned int)0xfffffffe)
1398 q++;
1399 q1 += 2;
1400 } else
1401 q1 += q1 & 1;
1402 }
1403 }
1404 ix0 = (q >> 1) + 0x3fe00000;
1405 ix1 = q1 >> 1;
1406 if (q & 1)
1407 ix1 |= sign;
1408 ix = ix0 + ((unsigned int)m << 20);
1409 ix <<= 32;
1410 ix |= ix1;
1411 return *(double*)&ix;
1412 #endif
1413 }
1414
1415 /*********************************************************************
1416 * tan (MSVCRT.@)
1417 */
1418 double CDECL tan( double x )
1419 {
1420 double ret = unix_funcs->tan(x);
1421 if (!isfinite(x)) return math_error(_DOMAIN, "tan", x, 0, ret);
1422 return ret;
1423 }
1424
1425 /*********************************************************************
1426 * tanh (MSVCRT.@)
1427 */
1428 double CDECL tanh( double x )
1429 {
1430 double ret = unix_funcs->tanh(x);
1431 if (isnan(x)) return math_error(_DOMAIN, "tanh", x, 0, ret);
1432 return ret;
1433 }
1434
1435
1436 #if (defined(__GNUC__) || defined(__clang__)) && defined(__i386__)
1437
1438 #define CREATE_FPU_FUNC1(name, call) \
1439 __ASM_GLOBAL_FUNC(name, \
1440 "pushl %ebp\n\t" \
1441 __ASM_CFI(".cfi_adjust_cfa_offset 4\n\t") \
1442 __ASM_CFI(".cfi_rel_offset %ebp,0\n\t") \
1443 "movl %esp, %ebp\n\t" \
1444 __ASM_CFI(".cfi_def_cfa_register %ebp\n\t") \
1445 "subl $68, %esp\n\t" /* sizeof(double)*8 + sizeof(int) */ \
1446 "fstpl (%esp)\n\t" /* store function argument */ \
1447 "fwait\n\t" \
1448 "movl $1, %ecx\n\t" /* empty FPU stack */ \
1449 "1:\n\t" \
1450 "fxam\n\t" \
1451 "fstsw %ax\n\t" \
1452 "and $0x4500, %ax\n\t" \
1453 "cmp $0x4100, %ax\n\t" \
1454 "je 2f\n\t" \
1455 "fstpl (%esp,%ecx,8)\n\t" \
1456 "fwait\n\t" \
1457 "incl %ecx\n\t" \
1458 "jmp 1b\n\t" \
1459 "2:\n\t" \
1460 "movl %ecx, -4(%ebp)\n\t" \
1461 "call " __ASM_NAME( #call ) "\n\t" \
1462 "movl -4(%ebp), %ecx\n\t" \
1463 "fstpl (%esp)\n\t" /* save result */ \
1464 "3:\n\t" /* restore FPU stack */ \
1465 "decl %ecx\n\t" \
1466 "fldl (%esp,%ecx,8)\n\t" \
1467 "cmpl $0, %ecx\n\t" \
1468 "jne 3b\n\t" \
1469 "leave\n\t" \
1470 __ASM_CFI(".cfi_def_cfa %esp,4\n\t") \
1471 __ASM_CFI(".cfi_same_value %ebp\n\t") \
1472 "ret")
1473
1474 #define CREATE_FPU_FUNC2(name, call) \
1475 __ASM_GLOBAL_FUNC(name, \
1476 "pushl %ebp\n\t" \
1477 __ASM_CFI(".cfi_adjust_cfa_offset 4\n\t") \
1478 __ASM_CFI(".cfi_rel_offset %ebp,0\n\t") \
1479 "movl %esp, %ebp\n\t" \
1480 __ASM_CFI(".cfi_def_cfa_register %ebp\n\t") \
1481 "subl $68, %esp\n\t" /* sizeof(double)*8 + sizeof(int) */ \
1482 "fstpl 8(%esp)\n\t" /* store function argument */ \
1483 "fwait\n\t" \
1484 "fstpl (%esp)\n\t" \
1485 "fwait\n\t" \
1486 "movl $2, %ecx\n\t" /* empty FPU stack */ \
1487 "1:\n\t" \
1488 "fxam\n\t" \
1489 "fstsw %ax\n\t" \
1490 "and $0x4500, %ax\n\t" \
1491 "cmp $0x4100, %ax\n\t" \
1492 "je 2f\n\t" \
1493 "fstpl (%esp,%ecx,8)\n\t" \
1494 "fwait\n\t" \
1495 "incl %ecx\n\t" \
1496 "jmp 1b\n\t" \
1497 "2:\n\t" \
1498 "movl %ecx, -4(%ebp)\n\t" \
1499 "call " __ASM_NAME( #call ) "\n\t" \
1500 "movl -4(%ebp), %ecx\n\t" \
1501 "fstpl 8(%esp)\n\t" /* save result */ \
1502 "3:\n\t" /* restore FPU stack */ \
1503 "decl %ecx\n\t" \
1504 "fldl (%esp,%ecx,8)\n\t" \
1505 "cmpl $1, %ecx\n\t" \
1506 "jne 3b\n\t" \
1507 "leave\n\t" \
1508 __ASM_CFI(".cfi_def_cfa %esp,4\n\t") \
1509 __ASM_CFI(".cfi_same_value %ebp\n\t") \
1510 "ret")
1511
1512 CREATE_FPU_FUNC1(_CIacos, acos)
1513 CREATE_FPU_FUNC1(_CIasin, asin)
1514 CREATE_FPU_FUNC1(_CIatan, atan)
1515 CREATE_FPU_FUNC2(_CIatan2, atan2)
1516 CREATE_FPU_FUNC1(_CIcos, cos)
1517 CREATE_FPU_FUNC1(_CIcosh, cosh)
1518 CREATE_FPU_FUNC1(_CIexp, exp)
1519 CREATE_FPU_FUNC2(_CIfmod, fmod)
1520 CREATE_FPU_FUNC1(_CIlog, log)
1521 CREATE_FPU_FUNC1(_CIlog10, log10)
1522 CREATE_FPU_FUNC2(_CIpow, pow)
1523 CREATE_FPU_FUNC1(_CIsin, sin)
1524 CREATE_FPU_FUNC1(_CIsinh, sinh)
1525 CREATE_FPU_FUNC1(_CIsqrt, sqrt)
1526 CREATE_FPU_FUNC1(_CItan, tan)
1527 CREATE_FPU_FUNC1(_CItanh, tanh)
1528
1529 __ASM_GLOBAL_FUNC(_ftol,
1530 "pushl %ebp\n\t"
1531 __ASM_CFI(".cfi_adjust_cfa_offset 4\n\t")
1532 __ASM_CFI(".cfi_rel_offset %ebp,0\n\t")
1533 "movl %esp, %ebp\n\t"
1534 __ASM_CFI(".cfi_def_cfa_register %ebp\n\t")
1535 "subl $12, %esp\n\t" /* sizeof(LONGLONG) + 2*sizeof(WORD) */
1536 "fnstcw (%esp)\n\t"
1537 "mov (%esp), %ax\n\t"
1538 "or $0xc00, %ax\n\t"
1539 "mov %ax, 2(%esp)\n\t"
1540 "fldcw 2(%esp)\n\t"
1541 "fistpq 4(%esp)\n\t"
1542 "fldcw (%esp)\n\t"
1543 "movl 4(%esp), %eax\n\t"
1544 "movl 8(%esp), %edx\n\t"
1545 "leave\n\t"
1546 __ASM_CFI(".cfi_def_cfa %esp,4\n\t")
1547 __ASM_CFI(".cfi_same_value %ebp\n\t")
1548 "ret")
1549
1550 #endif /* (defined(__GNUC__) || defined(__clang__)) && defined(__i386__) */
1551
1552 /*********************************************************************
1553 * _fpclass (MSVCRT.@)
1554 */
1555 int CDECL _fpclass(double num)
1556 {
1557 union { double f; UINT64 i; } u = { num };
1558 int e = u.i >> 52 & 0x7ff;
1559 int s = u.i >> 63;
1560
1561 switch (e)
1562 {
1563 case 0:
1564 if (u.i << 1) return s ? _FPCLASS_ND : _FPCLASS_PD;
1565 return s ? _FPCLASS_NZ : _FPCLASS_PZ;
1566 case 0x7ff:
1567 if (u.i << 12) return ((u.i >> 51) & 1) ? _FPCLASS_QNAN : _FPCLASS_SNAN;
1568 return s ? _FPCLASS_NINF : _FPCLASS_PINF;
1569 default:
1570 return s ? _FPCLASS_NN : _FPCLASS_PN;
1571 }
1572 }
1573
1574 /*********************************************************************
1575 * _rotl (MSVCRT.@)
1576 */
1577 unsigned int CDECL MSVCRT__rotl(unsigned int num, int shift)
1578 {
1579 shift &= 31;
1580 return (num << shift) | (num >> (32-shift));
1581 }
1582
1583 /*********************************************************************
1584 * _lrotl (MSVCRT.@)
1585 */
1586 __msvcrt_ulong CDECL MSVCRT__lrotl(__msvcrt_ulong num, int shift)
1587 {
1588 shift &= 0x1f;
1589 return (num << shift) | (num >> (32-shift));
1590 }
1591
1592 /*********************************************************************
1593 * _lrotr (MSVCRT.@)
1594 */
1595 __msvcrt_ulong CDECL MSVCRT__lrotr(__msvcrt_ulong num, int shift)
1596 {
1597 shift &= 0x1f;
1598 return (num >> shift) | (num << (32-shift));
1599 }
1600
1601 /*********************************************************************
1602 * _rotr (MSVCRT.@)
1603 */
1604 unsigned int CDECL MSVCRT__rotr(unsigned int num, int shift)
1605 {
1606 shift &= 0x1f;
1607 return (num >> shift) | (num << (32-shift));
1608 }
1609
1610 /*********************************************************************
1611 * _rotl64 (MSVCRT.@)
1612 */
1613 unsigned __int64 CDECL MSVCRT__rotl64(unsigned __int64 num, int shift)
1614 {
1615 shift &= 63;
1616 return (num << shift) | (num >> (64-shift));
1617 }
1618
1619 /*********************************************************************
1620 * _rotr64 (MSVCRT.@)
1621 */
1622 unsigned __int64 CDECL MSVCRT__rotr64(unsigned __int64 num, int shift)
1623 {
1624 shift &= 63;
1625 return (num >> shift) | (num << (64-shift));
1626 }
1627
1628 /*********************************************************************
1629 * abs (MSVCRT.@)
1630 */
1631 int CDECL abs( int n )
1632 {
1633 return n >= 0 ? n : -n;
1634 }
1635
1636 /*********************************************************************
1637 * labs (MSVCRT.@)
1638 */
1639 __msvcrt_long CDECL labs( __msvcrt_long n )
1640 {
1641 return n >= 0 ? n : -n;
1642 }
1643
1644 #if _MSVCR_VER>=100
1645 /*********************************************************************
1646 * llabs (MSVCR100.@)
1647 */
1648 __int64 CDECL llabs( __int64 n )
1649 {
1650 return n >= 0 ? n : -n;
1651 }
1652 #endif
1653
1654 #if _MSVCR_VER>=120
1655 /*********************************************************************
1656 * imaxabs (MSVCR120.@)
1657 */
1658 intmax_t CDECL imaxabs( intmax_t n )
1659 {
1660 return n >= 0 ? n : -n;
1661 }
1662 #endif
1663
1664 /*********************************************************************
1665 * _abs64 (MSVCRT.@)
1666 */
1667 __int64 CDECL _abs64( __int64 n )
1668 {
1669 return n >= 0 ? n : -n;
1670 }
1671
1672 /*********************************************************************
1673 * _logb (MSVCRT.@)
1674 */
1675 double CDECL _logb(double num)
1676 {
1677 double ret = unix_funcs->logb(num);
1678 if (isnan(num)) return math_error(_DOMAIN, "_logb", num, 0, ret);
1679 if (!num) return math_error(_SING, "_logb", num, 0, ret);
1680 return ret;
1681 }
1682
1683 /*********************************************************************
1684 * _hypot (MSVCRT.@)
1685 */
1686 double CDECL _hypot(double x, double y)
1687 {
1688 /* FIXME: errno handling */
1689 return unix_funcs->hypot( x, y );
1690 }
1691
1692 /*********************************************************************
1693 * _hypotf (MSVCRT.@)
1694 */
1695 float CDECL _hypotf(float x, float y)
1696 {
1697 /* FIXME: errno handling */
1698 return unix_funcs->hypotf( x, y );
1699 }
1700
1701 /*********************************************************************
1702 * ceil (MSVCRT.@)
1703 */
1704 double CDECL ceil( double x )
1705 {
1706 return unix_funcs->ceil(x);
1707 }
1708
1709 /*********************************************************************
1710 * floor (MSVCRT.@)
1711 */
1712 double CDECL floor( double x )
1713 {
1714 return unix_funcs->floor(x);
1715 }
1716
1717 /*********************************************************************
1718 * fma (MSVCRT.@)
1719 */
1720 double CDECL fma( double x, double y, double z )
1721 {
1722 double w = unix_funcs->fma(x, y, z);
1723 if ((isinf(x) && y == 0) || (x == 0 && isinf(y))) *_errno() = EDOM;
1724 else if (isinf(x) && isinf(z) && x != z) *_errno() = EDOM;
1725 else if (isinf(y) && isinf(z) && y != z) *_errno() = EDOM;
1726 return w;
1727 }
1728
1729 /*********************************************************************
1730 * fmaf (MSVCRT.@)
1731 */
1732 float CDECL fmaf( float x, float y, float z )
1733 {
1734 float w = unix_funcs->fmaf(x, y, z);
1735 if ((isinf(x) && y == 0) || (x == 0 && isinf(y))) *_errno() = EDOM;
1736 else if (isinf(x) && isinf(z) && x != z) *_errno() = EDOM;
1737 else if (isinf(y) && isinf(z) && y != z) *_errno() = EDOM;
1738 return w;
1739 }
1740
1741 /*********************************************************************
1742 * fabs (MSVCRT.@)
1743 *
1744 * Copied from musl: src/math/fabsf.c
1745 */
1746 double CDECL fabs( double x )
1747 {
1748 union { double f; UINT64 i; } u = { x };
1749 u.i &= ~0ull >> 1;
1750 return u.f;
1751 }
1752
1753 /*********************************************************************
1754 * frexp (MSVCRT.@)
1755 */
1756 double CDECL frexp( double x, int *exp )
1757 {
1758 return unix_funcs->frexp( x, exp );
1759 }
1760
1761 /*********************************************************************
1762 * modf (MSVCRT.@)
1763 */
1764 double CDECL modf( double x, double *iptr )
1765 {
1766 return unix_funcs->modf( x, iptr );
1767 }
1768
1769 /**********************************************************************
1770 * _statusfp2 (MSVCRT.@)
1771 *
1772 * Not exported by native msvcrt, added in msvcr80.
1773 */
1774 #if defined(__i386__) || defined(__x86_64__)
1775 void CDECL _statusfp2( unsigned int *x86_sw, unsigned int *sse2_sw )
1776 {
1777 #if defined(__GNUC__) || defined(__clang__)
1778 unsigned int flags;
1779 unsigned long fpword;
1780
1781 if (x86_sw)
1782 {
1783 __asm__ __volatile__( "fstsw %0" : "=m" (fpword) );
1784 flags = 0;
1785 if (fpword & 0x1) flags |= _SW_INVALID;
1786 if (fpword & 0x2) flags |= _SW_DENORMAL;
1787 if (fpword & 0x4) flags |= _SW_ZERODIVIDE;
1788 if (fpword & 0x8) flags |= _SW_OVERFLOW;
1789 if (fpword & 0x10) flags |= _SW_UNDERFLOW;
1790 if (fpword & 0x20) flags |= _SW_INEXACT;
1791 *x86_sw = flags;
1792 }
1793
1794 if (!sse2_sw) return;
1795
1796 if (sse2_supported)
1797 {
1798 __asm__ __volatile__( "stmxcsr %0" : "=m" (fpword) );
1799 flags = 0;
1800 if (fpword & 0x1) flags |= _SW_INVALID;
1801 if (fpword & 0x2) flags |= _SW_DENORMAL;
1802 if (fpword & 0x4) flags |= _SW_ZERODIVIDE;
1803 if (fpword & 0x8) flags |= _SW_OVERFLOW;
1804 if (fpword & 0x10) flags |= _SW_UNDERFLOW;
1805 if (fpword & 0x20) flags |= _SW_INEXACT;
1806 *sse2_sw = flags;
1807 }
1808 else *sse2_sw = 0;
1809 #else
1810 FIXME( "not implemented\n" );
1811 #endif
1812 }
1813 #endif
1814
1815 /**********************************************************************
1816 * _statusfp (MSVCRT.@)
1817 */
1818 unsigned int CDECL _statusfp(void)
1819 {
1820 unsigned int flags = 0;
1821 #if defined(__i386__) || defined(__x86_64__)
1822 unsigned int x86_sw, sse2_sw;
1823
1824 _statusfp2( &x86_sw, &sse2_sw );
1825 /* FIXME: there's no definition for ambiguous status, just return all status bits for now */
1826 flags = x86_sw | sse2_sw;
1827 #elif defined(__aarch64__)
1828 ULONG_PTR fpsr;
1829
1830 __asm__ __volatile__( "mrs %0, fpsr" : "=r" (fpsr) );
1831 if (fpsr & 0x1) flags |= _SW_INVALID;
1832 if (fpsr & 0x2) flags |= _SW_ZERODIVIDE;
1833 if (fpsr & 0x4) flags |= _SW_OVERFLOW;
1834 if (fpsr & 0x8) flags |= _SW_UNDERFLOW;
1835 if (fpsr & 0x10) flags |= _SW_INEXACT;
1836 if (fpsr & 0x80) flags |= _SW_DENORMAL;
1837 #else
1838 FIXME( "not implemented\n" );
1839 #endif
1840 return flags;
1841 }
1842
1843 /*********************************************************************
1844 * _clearfp (MSVCRT.@)
1845 */
1846 unsigned int CDECL _clearfp(void)
1847 {
1848 unsigned int flags = 0;
1849 #if (defined(__GNUC__) || defined(__clang__)) && (defined(__i386__) || defined(__x86_64__))
1850 unsigned long fpword;
1851
1852 __asm__ __volatile__( "fnstsw %0; fnclex" : "=m" (fpword) );
1853 if (fpword & 0x1) flags |= _SW_INVALID;
1854 if (fpword & 0x2) flags |= _SW_DENORMAL;
1855 if (fpword & 0x4) flags |= _SW_ZERODIVIDE;
1856 if (fpword & 0x8) flags |= _SW_OVERFLOW;
1857 if (fpword & 0x10) flags |= _SW_UNDERFLOW;
1858 if (fpword & 0x20) flags |= _SW_INEXACT;
1859
1860 if (sse2_supported)
1861 {
1862 __asm__ __volatile__( "stmxcsr %0" : "=m" (fpword) );
1863 if (fpword & 0x1) flags |= _SW_INVALID;
1864 if (fpword & 0x2) flags |= _SW_DENORMAL;
1865 if (fpword & 0x4) flags |= _SW_ZERODIVIDE;
1866 if (fpword & 0x8) flags |= _SW_OVERFLOW;
1867 if (fpword & 0x10) flags |= _SW_UNDERFLOW;
1868 if (fpword & 0x20) flags |= _SW_INEXACT;
1869 fpword &= ~0x3f;
1870 __asm__ __volatile__( "ldmxcsr %0" : : "m" (fpword) );
1871 }
1872 #elif defined(__aarch64__)
1873 ULONG_PTR fpsr;
1874
1875 __asm__ __volatile__( "mrs %0, fpsr" : "=r" (fpsr) );
1876 if (fpsr & 0x1) flags |= _SW_INVALID;
1877 if (fpsr & 0x2) flags |= _SW_ZERODIVIDE;
1878 if (fpsr & 0x4) flags |= _SW_OVERFLOW;
1879 if (fpsr & 0x8) flags |= _SW_UNDERFLOW;
1880 if (fpsr & 0x10) flags |= _SW_INEXACT;
1881 if (fpsr & 0x80) flags |= _SW_DENORMAL;
1882 fpsr &= ~0x9f;
1883 __asm__ __volatile__( "msr fpsr, %0" :: "r" (fpsr) );
1884 #else
1885 FIXME( "not implemented\n" );
1886 #endif
1887 return flags;
1888 }
1889
1890 /*********************************************************************
1891 * __fpecode (MSVCRT.@)
1892 */
1893 int * CDECL __fpecode(void)
1894 {
1895 return &msvcrt_get_thread_data()->fpecode;
1896 }
1897
1898 /*********************************************************************
1899 * ldexp (MSVCRT.@)
1900 */
1901 double CDECL ldexp(double num, int exp)
1902 {
1903 double z = unix_funcs->ldexp(num,exp);
1904
1905 if (isfinite(num) && !isfinite(z))
1906 return math_error(_OVERFLOW, "ldexp", num, exp, z);
1907 if (num && isfinite(num) && !z)
1908 return math_error(_UNDERFLOW, "ldexp", num, exp, z);
1909 if (z == 0 && signbit(z))
1910 z = 0.0; /* Convert -0 -> +0 */
1911 return z;
1912 }
1913
1914 /*********************************************************************
1915 * _cabs (MSVCRT.@)
1916 */
1917 double CDECL _cabs(struct _complex num)
1918 {
1919 return sqrt(num.x * num.x + num.y * num.y);
1920 }
1921
1922 /*********************************************************************
1923 * _chgsign (MSVCRT.@)
1924 */
1925 double CDECL _chgsign(double num)
1926 {
1927 union { double f; UINT64 i; } u = { num };
1928 u.i ^= 1ull << 63;
1929 return u.f;
1930 }
1931
1932 /*********************************************************************
1933 * __control87_2 (MSVCR80.@)
1934 *
1935 * Not exported by native msvcrt, added in msvcr80.
1936 */
1937 #ifdef __i386__
1938 int CDECL __control87_2( unsigned int newval, unsigned int mask,
1939 unsigned int *x86_cw, unsigned int *sse2_cw )
1940 {
1941 #if defined(__GNUC__) || defined(__clang__)
1942 unsigned long fpword;
1943 unsigned int flags;
1944 unsigned int old_flags;
1945
1946 if (x86_cw)
1947 {
1948 __asm__ __volatile__( "fstcw %0" : "=m" (fpword) );
1949
1950 /* Convert into mask constants */
1951 flags = 0;
1952 if (fpword & 0x1) flags |= _EM_INVALID;
1953 if (fpword & 0x2) flags |= _EM_DENORMAL;
1954 if (fpword & 0x4) flags |= _EM_ZERODIVIDE;
1955 if (fpword & 0x8) flags |= _EM_OVERFLOW;
1956 if (fpword & 0x10) flags |= _EM_UNDERFLOW;
1957 if (fpword & 0x20) flags |= _EM_INEXACT;
1958 switch (fpword & 0xc00)
1959 {
1960 case 0xc00: flags |= _RC_UP|_RC_DOWN; break;
1961 case 0x800: flags |= _RC_UP; break;
1962 case 0x400: flags |= _RC_DOWN; break;
1963 }
1964 switch (fpword & 0x300)
1965 {
1966 case 0x0: flags |= _PC_24; break;
1967 case 0x200: flags |= _PC_53; break;
1968 case 0x300: flags |= _PC_64; break;
1969 }
1970 if (fpword & 0x1000) flags |= _IC_AFFINE;
1971
1972 TRACE( "x86 flags=%08x newval=%08x mask=%08x\n", flags, newval, mask );
1973 if (mask)
1974 {
1975 flags = (flags & ~mask) | (newval & mask);
1976
1977 /* Convert (masked) value back to fp word */
1978 fpword = 0;
1979 if (flags & _EM_INVALID) fpword |= 0x1;
1980 if (flags & _EM_DENORMAL) fpword |= 0x2;
1981 if (flags & _EM_ZERODIVIDE) fpword |= 0x4;
1982 if (flags & _EM_OVERFLOW) fpword |= 0x8;
1983 if (flags & _EM_UNDERFLOW) fpword |= 0x10;
1984 if (flags & _EM_INEXACT) fpword |= 0x20;
1985 switch (flags & _MCW_RC)
1986 {
1987 case _RC_UP|_RC_DOWN: fpword |= 0xc00; break;
1988 case _RC_UP: fpword |= 0x800; break;
1989 case _RC_DOWN: fpword |= 0x400; break;
1990 }
1991 switch (flags & _MCW_PC)
1992 {
1993 case _PC_64: fpword |= 0x300; break;
1994 case _PC_53: fpword |= 0x200; break;
1995 case _PC_24: fpword |= 0x0; break;
1996 }
1997 if (flags & _IC_AFFINE) fpword |= 0x1000;
1998
1999 __asm__ __volatile__( "fldcw %0" : : "m" (fpword) );
2000 }
2001 *x86_cw = flags;
2002 }
2003
2004 if (!sse2_cw) return 1;
2005
2006 if (sse2_supported)
2007 {
2008 __asm__ __volatile__( "stmxcsr %0" : "=m" (fpword) );
2009
2010 /* Convert into mask constants */
2011 flags = 0;
2012 if (fpword & 0x80) flags |= _EM_INVALID;
2013 if (fpword & 0x100) flags |= _EM_DENORMAL;
2014 if (fpword & 0x200) flags |= _EM_ZERODIVIDE;
2015 if (fpword & 0x400) flags |= _EM_OVERFLOW;
2016 if (fpword & 0x800) flags |= _EM_UNDERFLOW;
2017 if (fpword & 0x1000) flags |= _EM_INEXACT;
2018 switch (fpword & 0x6000)
2019 {
2020 case 0x6000: flags |= _RC_UP|_RC_DOWN; break;
2021 case 0x4000: flags |= _RC_UP; break;
2022 case 0x2000: flags |= _RC_DOWN; break;
2023 }
2024 switch (fpword & 0x8040)
2025 {
2026 case 0x0040: flags |= _DN_FLUSH_OPERANDS_SAVE_RESULTS; break;
2027 case 0x8000: flags |= _DN_SAVE_OPERANDS_FLUSH_RESULTS; break;
2028 case 0x8040: flags |= _DN_FLUSH; break;
2029 }
2030
2031 TRACE( "sse2 flags=%08x newval=%08x mask=%08x\n", flags, newval, mask );
2032 if (mask)
2033 {
2034 old_flags = flags;
2035 mask &= _MCW_EM | _MCW_RC | _MCW_DN;
2036 flags = (flags & ~mask) | (newval & mask);
2037
2038 if (flags != old_flags)
2039 {
2040 /* Convert (masked) value back to fp word */
2041 fpword = 0;
2042 if (flags & _EM_INVALID) fpword |= 0x80;
2043 if (flags & _EM_DENORMAL) fpword |= 0x100;
2044 if (flags & _EM_ZERODIVIDE) fpword |= 0x200;
2045 if (flags & _EM_OVERFLOW) fpword |= 0x400;
2046 if (flags & _EM_UNDERFLOW) fpword |= 0x800;
2047 if (flags & _EM_INEXACT) fpword |= 0x1000;
2048 switch (flags & _MCW_RC)
2049 {
2050 case _RC_UP|_RC_DOWN: fpword |= 0x6000; break;
2051 case _RC_UP: fpword |= 0x4000; break;
2052 case _RC_DOWN: fpword |= 0x2000; break;
2053 }
2054 switch (flags & _MCW_DN)
2055 {
2056 case _DN_FLUSH_OPERANDS_SAVE_RESULTS: fpword |= 0x0040; break;
2057 case _DN_SAVE_OPERANDS_FLUSH_RESULTS: fpword |= 0x8000; break;
2058 case _DN_FLUSH: fpword |= 0x8040; break;
2059 }
2060 __asm__ __volatile__( "ldmxcsr %0" : : "m" (fpword) );
2061 }
2062 }
2063 *sse2_cw = flags;
2064 }
2065 else *sse2_cw = 0;
2066
2067 return 1;
2068 #else
2069 FIXME( "not implemented\n" );
2070 return 0;
2071 #endif
2072 }
2073 #endif
2074
2075 /*********************************************************************
2076 * _control87 (MSVCRT.@)
2077 */
2078 unsigned int CDECL _control87(unsigned int newval, unsigned int mask)
2079 {
2080 unsigned int flags = 0;
2081 #ifdef __i386__
2082 unsigned int sse2_cw;
2083
2084 __control87_2( newval, mask, &flags, &sse2_cw );
2085
2086 if ((flags ^ sse2_cw) & (_MCW_EM | _MCW_RC)) flags |= _EM_AMBIGUOUS;
2087 flags |= sse2_cw;
2088 #elif defined(__x86_64__)
2089 unsigned long fpword;
2090 unsigned int old_flags;
2091
2092 __asm__ __volatile__( "stmxcsr %0" : "=m" (fpword) );
2093 if (fpword & 0x80) flags |= _EM_INVALID;
2094 if (fpword & 0x100) flags |= _EM_DENORMAL;
2095 if (fpword & 0x200) flags |= _EM_ZERODIVIDE;
2096 if (fpword & 0x400) flags |= _EM_OVERFLOW;
2097 if (fpword & 0x800) flags |= _EM_UNDERFLOW;
2098 if (fpword & 0x1000) flags |= _EM_INEXACT;
2099 switch (fpword & 0x6000)
2100 {
2101 case 0x6000: flags |= _RC_CHOP; break;
2102 case 0x4000: flags |= _RC_UP; break;
2103 case 0x2000: flags |= _RC_DOWN; break;
2104 }
2105 switch (fpword & 0x8040)
2106 {
2107 case 0x0040: flags |= _DN_FLUSH_OPERANDS_SAVE_RESULTS; break;
2108 case 0x8000: flags |= _DN_SAVE_OPERANDS_FLUSH_RESULTS; break;
2109 case 0x8040: flags |= _DN_FLUSH; break;
2110 }
2111 old_flags = flags;
2112 mask &= _MCW_EM | _MCW_RC | _MCW_DN;
2113 flags = (flags & ~mask) | (newval & mask);
2114 if (flags != old_flags)
2115 {
2116 fpword = 0;
2117 if (flags & _EM_INVALID) fpword |= 0x80;
2118 if (flags & _EM_DENORMAL) fpword |= 0x100;
2119 if (flags & _EM_ZERODIVIDE) fpword |= 0x200;
2120 if (flags & _EM_OVERFLOW) fpword |= 0x400;
2121 if (flags & _EM_UNDERFLOW) fpword |= 0x800;
2122 if (flags & _EM_INEXACT) fpword |= 0x1000;
2123 switch (flags & _MCW_RC)
2124 {
2125 case _RC_CHOP: fpword |= 0x6000; break;
2126 case _RC_UP: fpword |= 0x4000; break;
2127 case _RC_DOWN: fpword |= 0x2000; break;
2128 }
2129 switch (flags & _MCW_DN)
2130 {
2131 case _DN_FLUSH_OPERANDS_SAVE_RESULTS: fpword |= 0x0040; break;
2132 case _DN_SAVE_OPERANDS_FLUSH_RESULTS: fpword |= 0x8000; break;
2133 case _DN_FLUSH: fpword |= 0x8040; break;
2134 }
2135 __asm__ __volatile__( "ldmxcsr %0" :: "m" (fpword) );
2136 }
2137 #elif defined(__aarch64__)
2138 ULONG_PTR fpcr;
2139
2140 __asm__ __volatile__( "mrs %0, fpcr" : "=r" (fpcr) );
2141 if (!(fpcr & 0x100)) flags |= _EM_INVALID;
2142 if (!(fpcr & 0x200)) flags |= _EM_ZERODIVIDE;
2143 if (!(fpcr & 0x400)) flags |= _EM_OVERFLOW;
2144 if (!(fpcr & 0x800)) flags |= _EM_UNDERFLOW;
2145 if (!(fpcr & 0x1000)) flags |= _EM_INEXACT;
2146 if (!(fpcr & 0x8000)) flags |= _EM_DENORMAL;
2147 switch (fpcr & 0xc00000)
2148 {
2149 case 0x400000: flags |= _RC_UP; break;
2150 case 0x800000: flags |= _RC_DOWN; break;
2151 case 0xc00000: flags |= _RC_CHOP; break;
2152 }
2153 flags = (flags & ~mask) | (newval & mask);
2154 fpcr &= ~0xc09f00ul;
2155 if (!(flags & _EM_INVALID)) fpcr |= 0x100;
2156 if (!(flags & _EM_ZERODIVIDE)) fpcr |= 0x200;
2157 if (!(flags & _EM_OVERFLOW)) fpcr |= 0x400;
2158 if (!(flags & _EM_UNDERFLOW)) fpcr |= 0x800;
2159 if (!(flags & _EM_INEXACT)) fpcr |= 0x1000;
2160 if (!(flags & _EM_DENORMAL)) fpcr |= 0x8000;
2161 switch (flags & _MCW_RC)
2162 {
2163 case _RC_CHOP: fpcr |= 0xc00000; break;
2164 case _RC_UP: fpcr |= 0x400000; break;
2165 case _RC_DOWN: fpcr |= 0x800000; break;
2166 }
2167 __asm__ __volatile__( "msr fpcr, %0" :: "r" (fpcr) );
2168 #else
2169 FIXME( "not implemented\n" );
2170 #endif
2171 return flags;
2172 }
2173
2174 /*********************************************************************
2175 * _controlfp (MSVCRT.@)
2176 */
2177 unsigned int CDECL _controlfp(unsigned int newval, unsigned int mask)
2178 {
2179 return _control87( newval, mask & ~_EM_DENORMAL );
2180 }
2181
2182 /*********************************************************************
2183 * _set_controlfp (MSVCRT.@)
2184 */
2185 void CDECL _set_controlfp( unsigned int newval, unsigned int mask )
2186 {
2187 _controlfp( newval, mask );
2188 }
2189
2190 /*********************************************************************
2191 * _controlfp_s (MSVCRT.@)
2192 */
2193 int CDECL _controlfp_s(unsigned int *cur, unsigned int newval, unsigned int mask)
2194 {
2195 static const unsigned int all_flags = (_MCW_EM | _MCW_IC | _MCW_RC |
2196 _MCW_PC | _MCW_DN);
2197 unsigned int val;
2198
2199 if (!MSVCRT_CHECK_PMT( !(newval & mask & ~all_flags) ))
2200 {
2201 if (cur) *cur = _controlfp( 0, 0 ); /* retrieve it anyway */
2202 return EINVAL;
2203 }
2204 val = _controlfp( newval, mask );
2205 if (cur) *cur = val;
2206 return 0;
2207 }
2208
2209 #if _MSVCR_VER>=120
2210 /*********************************************************************
2211 * fegetenv (MSVCR120.@)
2212 */
2213 int CDECL fegetenv(fenv_t *env)
2214 {
2215 env->_Fe_ctl = _controlfp(0, 0) & (_EM_INEXACT | _EM_UNDERFLOW |
2216 _EM_OVERFLOW | _EM_ZERODIVIDE | _EM_INVALID | _RC_CHOP);
2217 env->_Fe_stat = _statusfp();
2218 return 0;
2219 }
2220
2221 /*********************************************************************
2222 * fetestexcept (MSVCR120.@)
2223 */
2224 int CDECL fetestexcept(int flags)
2225 {
2226 return _statusfp() & flags;
2227 }
2228
2229 /*********************************************************************
2230 * fesetexceptflag (MSVCR120.@)
2231 */
2232 int CDECL fesetexceptflag(const fexcept_t *status, int excepts)
2233 {
2234 fenv_t env;
2235
2236 if(!excepts)
2237 return 0;
2238
2239 fegetenv(&env);
2240 excepts &= FE_ALL_EXCEPT;
2241 env._Fe_stat &= ~excepts;
2242 env._Fe_stat |= (*status & excepts);
2243 return fesetenv(&env);
2244 }
2245
2246 /*********************************************************************
2247 * feclearexcept (MSVCR120.@)
2248 */
2249 int CDECL feclearexcept(int flags)
2250 {
2251 fenv_t env;
2252
2253 fegetenv(&env);
2254 env._Fe_stat &= ~(flags & FE_ALL_EXCEPT);
2255 return fesetenv(&env);
2256 }
2257
2258 /*********************************************************************
2259 * fegetexceptflag (MSVCR120.@)
2260 */
2261 int CDECL fegetexceptflag(fexcept_t *status, int excepts)
2262 {
2263 *status = _statusfp() & excepts;
2264 return 0;
2265 }
2266 #endif
2267
2268 #if _MSVCR_VER>=140
2269 /*********************************************************************
2270 * __fpe_flt_rounds (UCRTBASE.@)
2271 */
2272 int CDECL __fpe_flt_rounds(void)
2273 {
2274 unsigned int fpc = _controlfp(0, 0) & _RC_CHOP;
2275
2276 TRACE("()\n");
2277
2278 switch(fpc) {
2279 case _RC_CHOP: return 0;
2280 case _RC_NEAR: return 1;
2281 case _RC_UP: return 2;
2282 default: return 3;
2283 }
2284 }
2285 #endif
2286
2287 #if _MSVCR_VER>=120
2288
2289 /*********************************************************************
2290 * fegetround (MSVCR120.@)
2291 */
2292 int CDECL fegetround(void)
2293 {
2294 return _controlfp(0, 0) & _RC_CHOP;
2295 }
2296
2297 /*********************************************************************
2298 * fesetround (MSVCR120.@)
2299 */
2300 int CDECL fesetround(int round_mode)
2301 {
2302 if (round_mode & (~_RC_CHOP))
2303 return 1;
2304 _controlfp(round_mode, _RC_CHOP);
2305 return 0;
2306 }
2307
2308 #endif /* _MSVCR_VER>=120 */
2309
2310 /*********************************************************************
2311 * _copysign (MSVCRT.@)
2312 *
2313 * Copied from musl: src/math/copysign.c
2314 */
2315 double CDECL _copysign( double x, double y )
2316 {
2317 union { double f; UINT64 i; } ux = { x }, uy = { y };
2318 ux.i &= ~0ull >> 1;
2319 ux.i |= uy.i & 1ull << 63;
2320 return ux.f;
2321 }
2322
2323 /*********************************************************************
2324 * _finite (MSVCRT.@)
2325 */
2326 int CDECL _finite(double num)
2327 {
2328 union { double f; UINT64 i; } u = { num };
2329 return (u.i & ~0ull >> 1) < 0x7ffull << 52;
2330 }
2331
2332 /*********************************************************************
2333 * _fpreset (MSVCRT.@)
2334 */
2335 void CDECL _fpreset(void)
2336 {
2337 #if (defined(__GNUC__) || defined(__clang__)) && (defined(__i386__) || defined(__x86_64__))
2338 const unsigned int x86_cw = 0x27f;
2339 __asm__ __volatile__( "fninit; fldcw %0" : : "m" (x86_cw) );
2340 if (sse2_supported)
2341 {
2342 const unsigned long sse2_cw = 0x1f80;
2343 __asm__ __volatile__( "ldmxcsr %0" : : "m" (sse2_cw) );
2344 }
2345 #else
2346 FIXME( "not implemented\n" );
2347 #endif
2348 }
2349
2350 #if _MSVCR_VER>=120
2351 /*********************************************************************
2352 * fesetenv (MSVCR120.@)
2353 */
2354 int CDECL fesetenv(const fenv_t *env)
2355 {
2356 #if (defined(__GNUC__) || defined(__clang__)) && (defined(__i386__) || defined(__x86_64__))
2357 struct {
2358 WORD control_word;
2359 WORD unused1;
2360 WORD status_word;
2361 WORD unused2;
2362 WORD tag_word;
2363 WORD unused3;
2364 DWORD instruction_pointer;
2365 WORD code_segment;
2366 WORD unused4;
2367 DWORD operand_addr;
2368 WORD data_segment;
2369 WORD unused5;
2370 } fenv;
2371
2372 TRACE( "(%p)\n", env );
2373
2374 if (!env->_Fe_ctl && !env->_Fe_stat) {
2375 _fpreset();
2376 return 0;
2377 }
2378
2379 __asm__ __volatile__( "fnstenv %0" : "=m" (fenv) );
2380
2381 fenv.control_word &= ~0xc3d;
2382 if (env->_Fe_ctl & _EM_INVALID) fenv.control_word |= 0x1;
2383 if (env->_Fe_ctl & _EM_ZERODIVIDE) fenv.control_word |= 0x4;
2384 if (env->_Fe_ctl & _EM_OVERFLOW) fenv.control_word |= 0x8;
2385 if (env->_Fe_ctl & _EM_UNDERFLOW) fenv.control_word |= 0x10;
2386 if (env->_Fe_ctl & _EM_INEXACT) fenv.control_word |= 0x20;
2387 switch (env->_Fe_ctl & _MCW_RC)
2388 {
2389 case _RC_UP|_RC_DOWN: fenv.control_word |= 0xc00; break;
2390 case _RC_UP: fenv.control_word |= 0x800; break;
2391 case _RC_DOWN: fenv.control_word |= 0x400; break;
2392 }
2393
2394 fenv.status_word &= ~0x3d;
2395 if (env->_Fe_stat & FE_INVALID) fenv.status_word |= 0x1;
2396 if (env->_Fe_stat & FE_DIVBYZERO) fenv.status_word |= 0x4;
2397 if (env->_Fe_stat & FE_OVERFLOW) fenv.status_word |= 0x8;
2398 if (env->_Fe_stat & FE_UNDERFLOW) fenv.status_word |= 0x10;
2399 if (env->_Fe_stat & FE_INEXACT) fenv.status_word |= 0x20;
2400
2401 __asm__ __volatile__( "fldenv %0" : : "m" (fenv) : "st", "st(1)",
2402 "st(2)", "st(3)", "st(4)", "st(5)", "st(6)", "st(7)" );
2403
2404 if (sse2_supported)
2405 {
2406 DWORD fpword;
2407 __asm__ __volatile__( "stmxcsr %0" : "=m" (fpword) );
2408 fpword &= ~0x7e80;
2409 if (env->_Fe_ctl & _EM_INVALID) fpword |= 0x80;
2410 if (env->_Fe_ctl & _EM_ZERODIVIDE) fpword |= 0x200;
2411 if (env->_Fe_ctl & _EM_OVERFLOW) fpword |= 0x400;
2412 if (env->_Fe_ctl & _EM_UNDERFLOW) fpword |= 0x800;
2413 if (env->_Fe_ctl & _EM_INEXACT) fpword |= 0x1000;
2414 switch (env->_Fe_ctl & _MCW_RC)
2415 {
2416 case _RC_CHOP: fpword |= 0x6000; break;
2417 case _RC_UP: fpword |= 0x4000; break;
2418 case _RC_DOWN: fpword |= 0x2000; break;
2419 }
2420 __asm__ __volatile__( "ldmxcsr %0" : : "m" (fpword) );
2421 }
2422
2423 return 0;
2424 #else
2425 FIXME( "not implemented\n" );
2426 #endif
2427 return 1;
2428 }
2429 #endif
2430
2431 /*********************************************************************
2432 * _isnan (MSVCRT.@)
2433 */
2434 int CDECL _isnan(double num)
2435 {
2436 union { double f; UINT64 i; } u = { num };
2437 return (u.i & ~0ull >> 1) > 0x7ffull << 52;
2438 }
2439
2440 static double pzero(double x)
2441 {
2442 static const double pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
2443 0.00000000000000000000e+00,
2444 -7.03124999999900357484e-02,
2445 -8.08167041275349795626e+00,
2446 -2.57063105679704847262e+02,
2447 -2.48521641009428822144e+03,
2448 -5.25304380490729545272e+03,
2449 }, pS8[5] = {
2450 1.16534364619668181717e+02,
2451 3.83374475364121826715e+03,
2452 4.05978572648472545552e+04,
2453 1.16752972564375915681e+05,
2454 4.76277284146730962675e+04,
2455 }, pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
2456 -1.14125464691894502584e-11,
2457 -7.03124940873599280078e-02,
2458 -4.15961064470587782438e+00,
2459 -6.76747652265167261021e+01,
2460 -3.31231299649172967747e+02,
2461 -3.46433388365604912451e+02,
2462 }, pS5[5] = {
2463 6.07539382692300335975e+01,
2464 1.05125230595704579173e+03,
2465 5.97897094333855784498e+03,
2466 9.62544514357774460223e+03,
2467 2.40605815922939109441e+03,
2468 }, pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
2469 -2.54704601771951915620e-09,
2470 -7.03119616381481654654e-02,
2471 -2.40903221549529611423e+00,
2472 -2.19659774734883086467e+01,
2473 -5.80791704701737572236e+01,
2474 -3.14479470594888503854e+01,
2475 }, pS3[5] = {
2476 3.58560338055209726349e+01,
2477 3.61513983050303863820e+02,
2478 1.19360783792111533330e+03,
2479 1.12799679856907414432e+03,
2480 1.73580930813335754692e+02,
2481 }, pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
2482 -8.87534333032526411254e-08,
2483 -7.03030995483624743247e-02,
2484 -1.45073846780952986357e+00,
2485 -7.63569613823527770791e+00,
2486 -1.11931668860356747786e+01,
2487 -3.23364579351335335033e+00,
2488 }, pS2[5] = {
2489 2.22202997532088808441e+01,
2490 1.36206794218215208048e+02,
2491 2.70470278658083486789e+02,
2492 1.53875394208320329881e+02,
2493 1.46576176948256193810e+01,
2494 };
2495
2496 const double *p, *q;
2497 double z, r, s;
2498 uint32_t ix;
2499
2500 ix = *(ULONGLONG*)&x >> 32;
2501 ix &= 0x7fffffff;
2502 if (ix >= 0x40200000) {
2503 p = pR8;
2504 q = pS8;
2505 } else if (ix >= 0x40122E8B) {
2506 p = pR5;
2507 q = pS5;
2508 } else if (ix >= 0x4006DB6D) {
2509 p = pR3;
2510 q = pS3;
2511 } else /*ix >= 0x40000000*/ {
2512 p = pR2;
2513 q = pS2;
2514 }
2515
2516 z = 1.0 / (x * x);
2517 r = p[0] + z * (p[1] + z * (p[2] + z * (p[3] + z * (p[4] + z * p[5]))));
2518 s = 1.0 + z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * q[4]))));
2519 return 1.0 + r / s;
2520 }
2521
2522 static double qzero(double x)
2523 {
2524 static const double qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
2525 0.00000000000000000000e+00,
2526 7.32421874999935051953e-02,
2527 1.17682064682252693899e+01,
2528 5.57673380256401856059e+02,
2529 8.85919720756468632317e+03,
2530 3.70146267776887834771e+04,
2531 }, qS8[6] = {
2532 1.63776026895689824414e+02,
2533 8.09834494656449805916e+03,
2534 1.42538291419120476348e+05,
2535 8.03309257119514397345e+05,
2536 8.40501579819060512818e+05,
2537 -3.43899293537866615225e+05,
2538 }, qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
2539 1.84085963594515531381e-11,
2540 7.32421766612684765896e-02,
2541 5.83563508962056953777e+00,
2542 1.35111577286449829671e+02,
2543 1.02724376596164097464e+03,
2544 1.98997785864605384631e+03,
2545 }, qS5[6] = {
2546 8.27766102236537761883e+01,
2547 2.07781416421392987104e+03,
2548 1.88472887785718085070e+04,
2549 5.67511122894947329769e+04,
2550 3.59767538425114471465e+04,
2551 -5.35434275601944773371e+03,
2552 }, qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
2553 4.37741014089738620906e-09,
2554 7.32411180042911447163e-02,
2555 3.34423137516170720929e+00,
2556 4.26218440745412650017e+01,
2557 1.70808091340565596283e+02,
2558 1.66733948696651168575e+02,
2559 }, qS3[6] = {
2560 4.87588729724587182091e+01,
2561 7.09689221056606015736e+02,
2562 3.70414822620111362994e+03,
2563 6.46042516752568917582e+03,
2564 2.51633368920368957333e+03,
2565 -1.49247451836156386662e+02,
2566 }, qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
2567 1.50444444886983272379e-07,
2568 7.32234265963079278272e-02,
2569 1.99819174093815998816e+00,
2570 1.44956029347885735348e+01,
2571 3.16662317504781540833e+01,
2572 1.62527075710929267416e+01,
2573 }, qS2[6] = {
2574 3.03655848355219184498e+01,
2575 2.69348118608049844624e+02,
2576 8.44783757595320139444e+02,
2577 8.82935845112488550512e+02,
2578 2.12666388511798828631e+02,
2579 -5.31095493882666946917e+00,
2580 };
2581
2582 const double *p, *q;
2583 double s, r, z;
2584 unsigned int ix;
2585
2586 ix = *(ULONGLONG*)&x >> 32;
2587 ix &= 0x7fffffff;
2588 if (ix >= 0x40200000) {
2589 p = qR8;
2590 q = qS8;
2591 } else if (ix >= 0x40122E8B) {
2592 p = qR5;
2593 q = qS5;
2594 } else if (ix >= 0x4006DB6D) {
2595 p = qR3;
2596 q = qS3;
2597 } else /*ix >= 0x40000000*/ {
2598 p = qR2;
2599 q = qS2;
2600 }
2601
2602 z = 1.0 / (x * x);
2603 r = p[0] + z * (p[1] + z * (p[2] + z * (p[3] + z * (p[4] + z * p[5]))));
2604 s = 1.0 + z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * q[5])))));
2605 return (-0.125 + r / s) / x;
2606 }
2607
2608 /* j0 and y0 approximation for |x|>=2 */
2609 static double j0_y0_approx(unsigned int ix, double x, BOOL y0)
2610 {
2611 static const double invsqrtpi = 5.64189583547756279280e-01;
2612
2613 double s, c, ss, cc, z;
2614
2615 s = sin(x);
2616 c = cos(x);
2617 if (y0) c = -c;
2618 cc = s + c;
2619 /* avoid overflow in 2*x, big ulp error when x>=0x1p1023 */
2620 if (ix < 0x7fe00000) {
2621 ss = s - c;
2622 z = -cos(2 * x);
2623 if (s * c < 0) cc = z / ss;
2624 else ss = z / cc;
2625 if (ix < 0x48000000) {
2626 if (y0) ss = -ss;
2627 cc = pzero(x) * cc - qzero(x) * ss;
2628 }
2629 }
2630 return invsqrtpi * cc / sqrt(x);
2631 }
2632
2633 /*********************************************************************
2634 * _j0 (MSVCRT.@)
2635 *
2636 * Copied from musl: src/math/j0.c
2637 */
2638 double CDECL _j0(double x)
2639 {
2640 static const double R02 = 1.56249999999999947958e-02,
2641 R03 = -1.89979294238854721751e-04,
2642 R04 = 1.82954049532700665670e-06,
2643 R05 = -4.61832688532103189199e-09,
2644 S01 = 1.56191029464890010492e-02,
2645 S02 = 1.16926784663337450260e-04,
2646 S03 = 5.13546550207318111446e-07,
2647 S04 = 1.16614003333790000205e-09;
2648
2649 double z, r, s;
2650 unsigned int ix;
2651
2652 ix = *(ULONGLONG*)&x >> 32;
2653 ix &= 0x7fffffff;
2654
2655 /* j0(+-inf)=0, j0(nan)=nan */
2656 if (ix >= 0x7ff00000)
2657 return math_error(_DOMAIN, "_j0", x, 0, 1 / (x * x));
2658 x = fabs(x);
2659
2660 if (ix >= 0x40000000) { /* |x| >= 2 */
2661 /* large ulp error near zeros: 2.4, 5.52, 8.6537,.. */
2662 return j0_y0_approx(ix, x, FALSE);
2663 }
2664
2665 if (ix >= 0x3f200000) { /* |x| >= 2**-13 */
2666 /* up to 4ulp error close to 2 */
2667 z = x * x;
2668 r = z * (R02 + z * (R03 + z * (R04 + z * R05)));
2669 s = 1 + z * (S01 + z * (S02 + z * (S03 + z * S04)));
2670 return (1 + x / 2) * (1 - x / 2) + z * (r / s);
2671 }
2672
2673 /* 1 - x*x/4 */
2674 /* prevent underflow */
2675 /* inexact should be raised when x!=0, this is not done correctly */
2676 if (ix >= 0x38000000) /* |x| >= 2**-127 */
2677 x = 0.25 * x * x;
2678 return 1 - x;
2679 }
2680
2681 static double pone(double x)
2682 {
2683 static const double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
2684 0.00000000000000000000e+00,
2685 1.17187499999988647970e-01,
2686 1.32394806593073575129e+01,
2687 4.12051854307378562225e+02,
2688 3.87474538913960532227e+03,
2689 7.91447954031891731574e+03,
2690 }, ps8[5] = {
2691 1.14207370375678408436e+02,
2692 3.65093083420853463394e+03,
2693 3.69562060269033463555e+04,
2694 9.76027935934950801311e+04,
2695 3.08042720627888811578e+04,
2696 }, pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
2697 1.31990519556243522749e-11,
2698 1.17187493190614097638e-01,
2699 6.80275127868432871736e+00,
2700 1.08308182990189109773e+02,
2701 5.17636139533199752805e+02,
2702 5.28715201363337541807e+02,
2703 }, ps5[5] = {
2704 5.92805987221131331921e+01,
2705 9.91401418733614377743e+02,
2706 5.35326695291487976647e+03,
2707 7.84469031749551231769e+03,
2708 1.50404688810361062679e+03,
2709 }, pr3[6] = {
2710 3.02503916137373618024e-09,
2711 1.17186865567253592491e-01,
2712 3.93297750033315640650e+00,
2713 3.51194035591636932736e+01,
2714 9.10550110750781271918e+01,
2715 4.85590685197364919645e+01,
2716 }, ps3[5] = {
2717 3.47913095001251519989e+01,
2718 3.36762458747825746741e+02,
2719 1.04687139975775130551e+03,
2720 8.90811346398256432622e+02,
2721 1.03787932439639277504e+02,
2722 }, pr2[6] = { /* for x in [2.8570,2]=1/[0.3499,0.5] */
2723 1.07710830106873743082e-07,
2724 1.17176219462683348094e-01,
2725 2.36851496667608785174e+00,
2726 1.22426109148261232917e+01,
2727 1.76939711271687727390e+01,
2728 5.07352312588818499250e+00,
2729 }, ps2[5] = {
2730 2.14364859363821409488e+01,
2731 1.25290227168402751090e+02,
2732 2.32276469057162813669e+02,
2733 1.17679373287147100768e+02,
2734 8.36463893371618283368e+00,
2735 };
2736
2737 const double *p, *q;
2738 double z, r, s;
2739 unsigned int ix;
2740
2741 ix = *(ULONGLONG*)&x >> 32;
2742 ix &= 0x7fffffff;
2743 if (ix >= 0x40200000) {
2744 p = pr8;
2745 q = ps8;
2746 } else if (ix >= 0x40122E8B) {
2747 p = pr5;
2748 q = ps5;
2749 } else if (ix >= 0x4006DB6D) {
2750 p = pr3;
2751 q = ps3;
2752 } else /*ix >= 0x40000000*/ {
2753 p = pr2;
2754 q = ps2;
2755 }
2756 z = 1.0 / (x * x);
2757 r = p[0] + z * (p[1] + z * (p[2] + z * (p[3] + z * (p[4] + z * p[5]))));
2758 s = 1.0 + z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * q[4]))));
2759 return 1.0 + r / s;
2760 }
2761
2762 static double qone(double x)
2763 {
2764 static const double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
2765 0.00000000000000000000e+00,
2766 -1.02539062499992714161e-01,
2767 -1.62717534544589987888e+01,
2768 -7.59601722513950107896e+02,
2769 -1.18498066702429587167e+04,
2770 -4.84385124285750353010e+04,
2771 }, qs8[6] = {
2772 1.61395369700722909556e+02,
2773 7.82538599923348465381e+03,
2774 1.33875336287249578163e+05,
2775 7.19657723683240939863e+05,
2776 6.66601232617776375264e+05,
2777 -2.94490264303834643215e+05,
2778 }, qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
2779 -2.08979931141764104297e-11,
2780 -1.02539050241375426231e-01,
2781 -8.05644828123936029840e+00,
2782 -1.83669607474888380239e+02,
2783 -1.37319376065508163265e+03,
2784 -2.61244440453215656817e+03,
2785 }, qs5[6] = {
2786 8.12765501384335777857e+01,
2787 1.99179873460485964642e+03,
2788 1.74684851924908907677e+04,
2789 4.98514270910352279316e+04,
2790 2.79480751638918118260e+04,
2791 -4.71918354795128470869e+03,
2792 }, qr3[6] = {
2793 -5.07831226461766561369e-09,
2794 -1.02537829820837089745e-01,
2795 -4.61011581139473403113e+00,
2796 -5.78472216562783643212e+01,
2797 -2.28244540737631695038e+02,
2798 -2.19210128478909325622e+02,
2799 }, qs3[6] = {
2800 4.76651550323729509273e+01,
2801 6.73865112676699709482e+02,
2802 3.38015286679526343505e+03,
2803 5.54772909720722782367e+03,
2804 1.90311919338810798763e+03,
2805 -1.35201191444307340817e+02,
2806 }, qr2[6] = { /* for x in [2.8570,2]=1/[0.3499,0.5] */
2807 -1.78381727510958865572e-07,
2808 -1.02517042607985553460e-01,
2809 -2.75220568278187460720e+00,
2810 -1.96636162643703720221e+01,
2811 -4.23253133372830490089e+01,
2812 -2.13719211703704061733e+01,
2813 }, qs2[6] = {
2814 2.95333629060523854548e+01,
2815 2.52981549982190529136e+02,
2816 7.57502834868645436472e+02,
2817 7.39393205320467245656e+02,
2818 1.55949003336666123687e+02,
2819 -4.95949898822628210127e+00,
2820 };
2821
2822 const double *p, *q;
2823 double s, r, z;
2824 unsigned int ix;
2825
2826 ix = *(ULONGLONG*)&x >> 32;
2827 ix &= 0x7fffffff;
2828 if (ix >= 0x40200000) {
2829 p = qr8;
2830 q = qs8;
2831 } else if (ix >= 0x40122E8B) {
2832 p = qr5;
2833 q = qs5;
2834 } else if (ix >= 0x4006DB6D) {
2835 p = qr3;
2836 q = qs3;
2837 } else /*ix >= 0x40000000*/ {
2838 p = qr2;
2839 q = qs2;
2840 }
2841 z = 1.0 / (x * x);
2842 r = p[0] + z * (p[1] + z * (p[2] + z * (p[3] + z * (p[4] + z * p[5]))));
2843 s = 1.0 + z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * q[5])))));
2844 return (0.375 + r / s) / x;
2845 }
2846
2847 static double j1_y1_approx(unsigned int ix, double x, BOOL y1, int sign)
2848 {
2849 static const double invsqrtpi = 5.64189583547756279280e-01;
2850
2851 double z, s, c, ss, cc;
2852
2853 s = sin(x);
2854 if (y1) s = -s;
2855 c = cos(x);
2856 cc = s - c;
2857 if (ix < 0x7fe00000) {
2858 ss = -s - c;
2859 z = cos(2 * x);
2860 if (s * c > 0) cc = z / ss;
2861 else ss = z / cc;
2862 if (ix < 0x48000000) {
2863 if (y1)
2864 ss = -ss;
2865 cc = pone(x) * cc - qone(x) * ss;
2866 }
2867 }
2868 if (sign)
2869 cc = -cc;
2870 return invsqrtpi * cc / sqrt(x);
2871 }
2872
2873 /*********************************************************************
2874 * _j1 (MSVCRT.@)
2875 *
2876 * Copied from musl: src/math/j1.c
2877 */
2878 double CDECL _j1(double x)
2879 {
2880 static const double r00 = -6.25000000000000000000e-02,
2881 r01 = 1.40705666955189706048e-03,
2882 r02 = -1.59955631084035597520e-05,
2883 r03 = 4.96727999609584448412e-08,
2884 s01 = 1.91537599538363460805e-02,
2885 s02 = 1.85946785588630915560e-04,
2886 s03 = 1.17718464042623683263e-06,
2887 s04 = 5.04636257076217042715e-09,
2888 s05 = 1.23542274426137913908e-11;
2889
2890 double z, r, s;
2891 unsigned int ix;
2892 int sign;
2893
2894 ix = *(ULONGLONG*)&x >> 32;
2895 sign = ix >> 31;
2896 ix &= 0x7fffffff;
2897 if (ix >= 0x7ff00000)
2898 return math_error(isnan(x) ? 0 : _DOMAIN, "_j1", x, 0, 1 / (x * x));
2899 if (ix >= 0x40000000) /* |x| >= 2 */
2900 return j1_y1_approx(ix, fabs(x), FALSE, sign);
2901 if (ix >= 0x38000000) { /* |x| >= 2**-127 */
2902 z = x * x;
2903 r = z * (r00 + z * (r01 + z * (r02 + z * r03)));
2904 s = 1 + z * (s01 + z * (s02 + z * (s03 + z * (s04 + z * s05))));
2905 z = r / s;
2906 } else {
2907 /* avoid underflow, raise inexact if x!=0 */
2908 z = x;
2909 }
2910 return (0.5 + z) * x;
2911 }
2912
2913 /*********************************************************************
2914 * _jn (MSVCRT.@)
2915 *
2916 * Copied from musl: src/math/jn.c
2917 */
2918 double CDECL _jn(int n, double x)
2919 {
2920 static const double invsqrtpi = 5.64189583547756279280e-01;
2921
2922 unsigned int ix, lx;
2923 int nm1, i, sign;
2924 double a, b, temp;
2925
2926 ix = *(ULONGLONG*)&x >> 32;
2927 lx = *(ULONGLONG*)&x;
2928 sign = ix >> 31;
2929 ix &= 0x7fffffff;
2930
2931 if ((ix | (lx | -lx) >> 31) > 0x7ff00000) /* nan */
2932 return x;
2933
2934 if (n == 0)
2935 return _j0(x);
2936 if (n < 0) {
2937 nm1 = -(n + 1);
2938 x = -x;
2939 sign ^= 1;
2940 } else {
2941 nm1 = n-1;
2942 }
2943 if (nm1 == 0)
2944 return j1(x);
2945
2946 sign &= n; /* even n: 0, odd n: signbit(x) */
2947 x = fabs(x);
2948 if ((ix | lx) == 0 || ix == 0x7ff00000) /* if x is 0 or inf */
2949 b = 0.0;
2950 else if (nm1 < x) {
2951 if (ix >= 0x52d00000) { /* x > 2**302 */
2952 switch(nm1 & 3) {
2953 case 0:
2954 temp = -cos(x) + sin(x);
2955 break;
2956 case 1:
2957 temp = -cos(x) - sin(x);
2958 break;
2959 case 2:
2960 temp = cos(x) - sin(x);
2961 break;
2962 default:
2963 temp = cos(x) + sin(x);
2964 break;
2965 }
2966 b = invsqrtpi * temp / sqrt(x);
2967 } else {
2968 a = _j0(x);
2969 b = _j1(x);
2970 for (i = 0; i < nm1; ) {
2971 i++;
2972 temp = b;
2973 b = b * (2.0 * i / x) - a; /* avoid underflow */
2974 a = temp;
2975 }
2976 }
2977 } else {
2978 if (ix < 0x3e100000) { /* x < 2**-29 */
2979 if (nm1 > 32) /* underflow */
2980 b = 0.0;
2981 else {
2982 temp = x * 0.5;
2983 b = temp;
2984 a = 1.0;
2985 for (i = 2; i <= nm1 + 1; i++) {
2986 a *= (double)i; /* a = n! */
2987 b *= temp; /* b = (x/2)^n */
2988 }
2989 b = b / a;
2990 }
2991 } else {
2992 double t, q0, q1, w, h, z, tmp, nf;
2993 int k;
2994
2995 nf = nm1 + 1.0;
2996 w = 2 * nf / x;
2997 h = 2 / x;
2998 z = w + h;
2999 q0 = w;
3000 q1 = w * z - 1.0;
3001 k = 1;
3002 while (q1 < 1.0e9) {
3003 k += 1;
3004 z += h;
3005 tmp = z * q1 - q0;
3006 q0 = q1;
3007 q1 = tmp;
3008 }
3009 for (t = 0.0, i = k; i >= 0; i--)
3010 t = 1 / (2 * (i + nf) / x - t);
3011 a = t;
3012 b = 1.0;
3013 tmp = nf * log(fabs(w));
3014 if (tmp < 7.09782712893383973096e+02) {
3015 for (i = nm1; i > 0; i--) {
3016 temp = b;
3017 b = b * (2.0 * i) / x - a;
3018 a = temp;
3019 }
3020 } else {
3021 for (i = nm1; i > 0; i--) {
3022 temp = b;
3023 b = b * (2.0 * i) / x - a;
3024 a = temp;
3025 /* scale b to avoid spurious overflow */
3026 if (b > 0x1p500) {
3027 a /= b;
3028 t /= b;
3029 b = 1.0;
3030 }
3031 }
3032 }
3033 z = j0(x);
3034 w = j1(x);
3035 if (fabs(z) >= fabs(w))
3036 b = t * z / b;
3037 else
3038 b = t * w / a;
3039 }
3040 }
3041 return sign ? -b : b;
3042 }
3043
3044 /*********************************************************************
3045 * _y0 (MSVCRT.@)
3046 */
3047 double CDECL _y0(double x)
3048 {
3049 static const double tpi = 6.36619772367581382433e-01,
3050 u00 = -7.38042951086872317523e-02,
3051 u01 = 1.76666452509181115538e-01,
3052 u02 = -1.38185671945596898896e-02,
3053 u03 = 3.47453432093683650238e-04,
3054 u04 = -3.81407053724364161125e-06,
3055 u05 = 1.95590137035022920206e-08,
3056 u06 = -3.98205194132103398453e-11,
3057 v01 = 1.27304834834123699328e-02,
3058 v02 = 7.60068627350353253702e-05,
3059 v03 = 2.59150851840457805467e-07,
3060 v04 = 4.41110311332675467403e-10;
3061
3062 double z, u, v;
3063 unsigned int ix, lx;
3064
3065 ix = *(ULONGLONG*)&x >> 32;
3066 lx = *(ULONGLONG*)&x;
3067
3068 /* y0(nan)=nan, y0(<0)=nan, y0(0)=-inf, y0(inf)=0 */
3069 if ((ix << 1 | lx) == 0)
3070 return math_error(_OVERFLOW, "_y0", x, 0, -INFINITY);
3071 if (isnan(x))
3072 return x;
3073 if (ix >> 31)
3074 return math_error(_DOMAIN, "_y0", x, 0, 0 / (x - x));
3075 if (ix >= 0x7ff00000)
3076 return 1 / x;
3077
3078 if (ix >= 0x40000000) { /* x >= 2 */
3079 /* large ulp errors near zeros: 3.958, 7.086,.. */
3080 return j0_y0_approx(ix, x, TRUE);
3081 }
3082
3083 if (ix >= 0x3e400000) { /* x >= 2**-27 */
3084 /* large ulp error near the first zero, x ~= 0.89 */
3085 z = x * x;
3086 u = u00 + z * (u01 + z * (u02 + z * (u03 + z * (u04 + z * (u05 + z * u06)))));
3087 v = 1.0 + z * (v01 + z * (v02 + z * (v03 + z * v04)));
3088 return u / v + tpi * (j0(x) * log(x));
3089 }
3090 return u00 + tpi * log(x);
3091 }
3092
3093 /*********************************************************************
3094 * _y1 (MSVCRT.@)
3095 */
3096 double CDECL _y1(double x)
3097 {
3098 static const double tpi = 6.36619772367581382433e-01,
3099 u00 = -1.96057090646238940668e-01,
3100 u01 = 5.04438716639811282616e-02,
3101 u02 = -1.91256895875763547298e-03,
3102 u03 = 2.35252600561610495928e-05,
3103 u04 = -9.19099158039878874504e-08,
3104 v00 = 1.99167318236649903973e-02,
3105 v01 = 2.02552581025135171496e-04,
3106 v02 = 1.35608801097516229404e-06,
3107 v03 = 6.22741452364621501295e-09,
3108 v04 = 1.66559246207992079114e-11;
3109
3110 double z, u, v;
3111 unsigned int ix, lx;
3112
3113 ix = *(ULONGLONG*)&x >> 32;
3114 lx = *(ULONGLONG*)&x;
3115
3116 /* y1(nan)=nan, y1(<0)=nan, y1(0)=-inf, y1(inf)=0 */
3117 if ((ix << 1 | lx) == 0)
3118 return math_error(_OVERFLOW, "_y1", x, 0, -INFINITY);
3119 if (isnan(x))
3120 return x;
3121 if (ix >> 31)
3122 return math_error(_DOMAIN, "_y1", x, 0, 0 / (x - x));
3123 if (ix >= 0x7ff00000)
3124 return 1 / x;
3125
3126 if (ix >= 0x40000000) /* x >= 2 */
3127 return j1_y1_approx(ix, x, TRUE, 0);
3128 if (ix < 0x3c900000) /* x < 2**-54 */
3129 return -tpi / x;
3130 z = x * x;
3131 u = u00 + z * (u01 + z * (u02 + z * (u03 + z * u04)));
3132 v = 1 + z * (v00 + z * (v01 + z * (v02 + z * (v03 + z * v04))));
3133 return x * (u / v) + tpi * (j1(x) * log(x) - 1 / x);
3134 }
3135
3136 /*********************************************************************
3137 * _yn (MSVCRT.@)
3138 *
3139 * Copied from musl: src/math/jn.c
3140 */
3141 double CDECL _yn(int n, double x)
3142 {
3143 static const double invsqrtpi = 5.64189583547756279280e-01;
3144
3145 unsigned int ix, lx, ib;
3146 int nm1, sign, i;
3147 double a, b, temp;
3148
3149 ix = *(ULONGLONG*)&x >> 32;
3150 lx = *(ULONGLONG*)&x;
3151 sign = ix >> 31;
3152 ix &= 0x7fffffff;
3153
3154 if ((ix | (lx | -lx) >> 31) > 0x7ff00000) /* nan */
3155 return x;
3156 if (sign && (ix | lx) != 0) /* x < 0 */
3157 return math_error(_DOMAIN, "_y1", x, 0, 0 / (x - x));
3158 if (ix == 0x7ff00000)
3159 return 0.0;
3160
3161 if (n == 0)
3162 return y0(x);
3163 if (n < 0) {
3164 nm1 = -(n + 1);
3165 sign = n & 1;
3166 } else {
3167 nm1 = n - 1;
3168 sign = 0;
3169 }
3170 if (nm1 == 0)
3171 return sign ? -y1(x) : y1(x);
3172
3173 if (ix >= 0x52d00000) { /* x > 2**302 */
3174 switch(nm1 & 3) {
3175 case 0:
3176 temp = -sin(x) - cos(x);
3177 break;
3178 case 1:
3179 temp = -sin(x) + cos(x);
3180 break;
3181 case 2:
3182 temp = sin(x) + cos(x);
3183 break;
3184 default:
3185 temp = sin(x) - cos(x);
3186 break;
3187 }
3188 b = invsqrtpi * temp / sqrt(x);
3189 } else {
3190 a = y0(x);
3191 b = y1(x);
3192 /* quit if b is -inf */
3193 ib = *(ULONGLONG*)&b >> 32;
3194 for (i = 0; i < nm1 && ib != 0xfff00000;) {
3195 i++;
3196 temp = b;
3197 b = (2.0 * i / x) * b - a;
3198 ib = *(ULONGLONG*)&b >> 32;
3199 a = temp;
3200 }
3201 }
3202 return sign ? -b : b;
3203 }
3204
3205 #if _MSVCR_VER>=120
3206
3207 /*********************************************************************
3208 * _nearbyint (MSVCR120.@)
3209 */
3210 double CDECL nearbyint(double num)
3211 {
3212 return unix_funcs->nearbyint( num );
3213 }
3214
3215 /*********************************************************************
3216 * _nearbyintf (MSVCR120.@)
3217 */
3218 float CDECL nearbyintf(float num)
3219 {
3220 return unix_funcs->nearbyintf( num );
3221 }
3222
3223 /*********************************************************************
3224 * nexttoward (MSVCR120.@)
3225 */
3226 double CDECL MSVCRT_nexttoward(double num, double next)
3227 {
3228 double ret = unix_funcs->nexttoward(num, next);
3229 if (!(_fpclass(ret) & (_FPCLASS_PN | _FPCLASS_NN
3230 | _FPCLASS_SNAN | _FPCLASS_QNAN)) && !isinf(num))
3231 {
3232 *_errno() = ERANGE;
3233 }
3234 return ret;
3235 }
3236
3237 /*********************************************************************
3238 * nexttowardf (MSVCR120.@)
3239 */
3240 float CDECL MSVCRT_nexttowardf(float num, double next)
3241 {
3242 float ret = unix_funcs->nexttowardf( num, next );
3243 if (!(_fpclass(ret) & (_FPCLASS_PN | _FPCLASS_NN
3244 | _FPCLASS_SNAN | _FPCLASS_QNAN)) && !isinf(num))
3245 {
3246 *_errno() = ERANGE;
3247 }
3248 return ret;
3249 }
3250
3251 #endif /* _MSVCR_VER>=120 */
3252
3253 /*********************************************************************
3254 * _nextafter (MSVCRT.@)
3255 */
3256 double CDECL _nextafter(double num, double next)
3257 {
3258 double retval;
3259 if (!isfinite(num) || !isfinite(next)) *_errno() = EDOM;
3260 retval = unix_funcs->nextafter(num,next);
3261 return retval;
3262 }
3263
3264 /*********************************************************************
3265 * _ecvt (MSVCRT.@)
3266 */
3267 char * CDECL _ecvt( double number, int ndigits, int *decpt, int *sign )
3268 {
3269 int prec, len;
3270 thread_data_t *data = msvcrt_get_thread_data();
3271 /* FIXME: check better for overflow (native supports over 300 chars) */
3272 ndigits = min( ndigits, 80 - 8); /* 8 : space for sign, dec point, "e",
3273 * 4 for exponent and one for
3274 * terminating '\0' */
3275 if (!data->efcvt_buffer)
3276 data->efcvt_buffer = malloc( 80 ); /* ought to be enough */
3277
3278 /* handle cases with zero ndigits or less */
3279 prec = ndigits;
3280 if( prec < 1) prec = 2;
3281 len = _snprintf(data->efcvt_buffer, 80, "%.*le", prec - 1, number);
3282
3283 if (data->efcvt_buffer[0] == '-') {
3284 memmove( data->efcvt_buffer, data->efcvt_buffer + 1, len-- );
3285 *sign = 1;
3286 } else *sign = 0;
3287
3288 /* take the decimal "point away */
3289 if( prec != 1)
3290 memmove( data->efcvt_buffer + 1, data->efcvt_buffer + 2, len - 1 );
3291 /* take the exponential "e" out */
3292 data->efcvt_buffer[ prec] = '\0';
3293 /* read the exponent */
3294 sscanf( data->efcvt_buffer + prec + 1, "%d", decpt);
3295 (*decpt)++;
3296 /* adjust for some border cases */
3297 if( data->efcvt_buffer[0] == '0')/* value is zero */
3298 *decpt = 0;
3299 /* handle cases with zero ndigits or less */
3300 if( ndigits < 1){
3301 if( data->efcvt_buffer[ 0] >= '5')
3302 (*decpt)++;
3303 data->efcvt_buffer[ 0] = '\0';
3304 }
3305 TRACE("out=\"%s\"\n",data->efcvt_buffer);
3306 return data->efcvt_buffer;
3307 }
3308
3309 /*********************************************************************
3310 * _ecvt_s (MSVCRT.@)
3311 */
3312 int CDECL _ecvt_s( char *buffer, size_t length, double number, int ndigits, int *decpt, int *sign )
3313 {
3314 int prec, len;
3315 char *result;
3316
3317 if (!MSVCRT_CHECK_PMT(buffer != NULL)) return EINVAL;
3318 if (!MSVCRT_CHECK_PMT(decpt != NULL)) return EINVAL;
3319 if (!MSVCRT_CHECK_PMT(sign != NULL)) return EINVAL;
3320 if (!MSVCRT_CHECK_PMT_ERR( length > 2, ERANGE )) return ERANGE;
3321 if (!MSVCRT_CHECK_PMT_ERR(ndigits < (int)length - 1, ERANGE )) return ERANGE;
3322
3323 /* handle cases with zero ndigits or less */
3324 prec = ndigits;
3325 if( prec < 1) prec = 2;
3326 result = malloc(prec + 8);
3327
3328 len = _snprintf(result, prec + 8, "%.*le", prec - 1, number);
3329 if (result[0] == '-') {
3330 memmove( result, result + 1, len-- );
3331 *sign = 1;
3332 } else *sign = 0;
3333
3334 /* take the decimal "point away */
3335 if( prec != 1)
3336 memmove( result + 1, result + 2, len - 1 );
3337 /* take the exponential "e" out */
3338 result[ prec] = '\0';
3339 /* read the exponent */
3340 sscanf( result + prec + 1, "%d", decpt);
3341 (*decpt)++;
3342 /* adjust for some border cases */
3343 if( result[0] == '0')/* value is zero */
3344 *decpt = 0;
3345 /* handle cases with zero ndigits or less */
3346 if( ndigits < 1){
3347 if( result[ 0] >= '5')
3348 (*decpt)++;
3349 result[ 0] = '\0';
3350 }
3351 memcpy( buffer, result, max(ndigits + 1, 1) );
3352 free( result );
3353 return 0;
3354 }
3355
3356 /***********************************************************************
3357 * _fcvt (MSVCRT.@)
3358 */
3359 char * CDECL _fcvt( double number, int ndigits, int *decpt, int *sign )
3360 {
3361 thread_data_t *data = msvcrt_get_thread_data();
3362 int stop, dec1, dec2;
3363 char *ptr1, *ptr2, *first;
3364 char buf[80]; /* ought to be enough */
3365 char decimal_separator = get_locinfo()->lconv->decimal_point[0];
3366
3367 if (!data->efcvt_buffer)
3368 data->efcvt_buffer = malloc( 80 ); /* ought to be enough */
3369
3370 stop = _snprintf(buf, 80, "%.*f", ndigits < 0 ? 0 : ndigits, number);
3371 ptr1 = buf;
3372 ptr2 = data->efcvt_buffer;
3373 first = NULL;
3374 dec1 = 0;
3375 dec2 = 0;
3376
3377 if (*ptr1 == '-') {
3378 *sign = 1;
3379 ptr1++;
3380 } else *sign = 0;
3381
3382 /* For numbers below the requested resolution, work out where
3383 the decimal point will be rather than finding it in the string */
3384 if (number < 1.0 && number > 0.0) {
3385 dec2 = log10(number + 1e-10);
3386 if (-dec2 <= ndigits) dec2 = 0;
3387 }
3388
3389 /* If requested digits is zero or less, we will need to truncate
3390 * the returned string */
3391 if (ndigits < 1) {
3392 stop += ndigits;
3393 }
3394
3395 while (*ptr1 == '0') ptr1++; /* Skip leading zeroes */
3396 while (*ptr1 != '\0' && *ptr1 != decimal_separator) {
3397 if (!first) first = ptr2;
3398 if ((ptr1 - buf) < stop) {
3399 *ptr2++ = *ptr1++;
3400 } else {
3401 ptr1++;
3402 }
3403 dec1++;
3404 }
3405
3406 if (ndigits > 0) {
3407 ptr1++;
3408 if (!first) {
3409 while (*ptr1 == '0') { /* Process leading zeroes */
3410 *ptr2++ = *ptr1++;
3411 dec1--;
3412 }
3413 }
3414 while (*ptr1 != '\0') {
3415 if (!first) first = ptr2;
3416 *ptr2++ = *ptr1++;
3417 }
3418 }
3419
3420 *ptr2 = '\0';
3421
3422 /* We never found a non-zero digit, then our number is either
3423 * smaller than the requested precision, or 0.0 */
3424 if (!first) {
3425 if (number > 0.0) {
3426 first = ptr2;
3427 } else {
3428 first = data->efcvt_buffer;
3429 dec1 = 0;
3430 }
3431 }
3432
3433 *decpt = dec2 ? dec2 : dec1;
3434 return first;
3435 }
3436
3437 /***********************************************************************
3438 * _fcvt_s (MSVCRT.@)
3439 */
3440 int CDECL _fcvt_s(char* outbuffer, size_t size, double number, int ndigits, int *decpt, int *sign)
3441 {
3442 int stop, dec1, dec2;
3443 char *ptr1, *ptr2, *first;
3444 char buf[80]; /* ought to be enough */
3445 char decimal_separator = get_locinfo()->lconv->decimal_point[0];
3446
3447 if (!outbuffer || !decpt || !sign || size == 0)
3448 {
3449 *_errno() = EINVAL;
3450 return EINVAL;
3451 }
3452
3453 stop = _snprintf(buf, 80, "%.*f", ndigits < 0 ? 0 : ndigits, number);
3454 ptr1 = buf;
3455 ptr2 = outbuffer;
3456 first = NULL;
3457 dec1 = 0;
3458 dec2 = 0;
3459
3460 if (*ptr1 == '-') {
3461 *sign = 1;
3462 ptr1++;
3463 } else *sign = 0;
3464
3465 /* For numbers below the requested resolution, work out where
3466 the decimal point will be rather than finding it in the string */
3467 if (number < 1.0 && number > 0.0) {
3468 dec2 = log10(number + 1e-10);
3469 if (-dec2 <= ndigits) dec2 = 0;
3470 }
3471
3472 /* If requested digits is zero or less, we will need to truncate
3473 * the returned string */
3474 if (ndigits < 1) {
3475 stop += ndigits;
3476 }
3477
3478 while (*ptr1 == '0') ptr1++; /* Skip leading zeroes */
3479 while (*ptr1 != '\0' && *ptr1 != decimal_separator) {
3480 if (!first) first = ptr2;
3481 if ((ptr1 - buf) < stop) {
3482 if (size > 1) {
3483 *ptr2++ = *ptr1++;
3484 size--;
3485 }
3486 } else {
3487 ptr1++;
3488 }
3489 dec1++;
3490 }
3491
3492 if (ndigits > 0) {
3493 ptr1++;
3494 if (!first) {
3495 while (*ptr1 == '0') { /* Process leading zeroes */
3496 if (number == 0.0 && size > 1) {
3497 *ptr2++ = '0';
3498 size--;
3499 }
3500 ptr1++;
3501 dec1--;
3502 }
3503 }
3504 while (*ptr1 != '\0') {
3505 if (!first) first = ptr2;
3506 if (size > 1) {
3507 *ptr2++ = *ptr1++;
3508 size--;
3509 }
3510 }
3511 }
3512
3513 *ptr2 = '\0';
3514
3515 /* We never found a non-zero digit, then our number is either
3516 * smaller than the requested precision, or 0.0 */
3517 if (!first && (number <= 0.0))
3518 dec1 = 0;
3519
3520 *decpt = dec2 ? dec2 : dec1;
3521 return 0;
3522 }
3523
3524 /***********************************************************************
3525 * _gcvt (MSVCRT.@)
3526 */
3527 char * CDECL _gcvt( double number, int ndigit, char *buff )
3528 {
3529 if(!buff) {
3530 *_errno() = EINVAL;
3531 return NULL;
3532 }
3533
3534 if(ndigit < 0) {
3535 *_errno() = ERANGE;
3536 return NULL;
3537 }
3538
3539 sprintf(buff, "%.*g", ndigit, number);
3540 return buff;
3541 }
3542
3543 /***********************************************************************
3544 * _gcvt_s (MSVCRT.@)
3545 */
3546 int CDECL _gcvt_s(char *buff, size_t size, double number, int digits)
3547 {
3548 int len;
3549
3550 if(!buff) {
3551 *_errno() = EINVAL;
3552 return EINVAL;
3553 }
3554
3555 if( digits<0 || digits>=size) {
3556 if(size)
3557 buff[0] = '\0';
3558
3559 *_errno() = ERANGE;
3560 return ERANGE;
3561 }
3562
3563 len = _scprintf("%.*g", digits, number);
3564 if(len > size) {
3565 buff[0] = '\0';
3566 *_errno() = ERANGE;
3567 return ERANGE;
3568 }
3569
3570 sprintf(buff, "%.*g", digits, number);
3571 return 0;
3572 }
3573
3574 #include <stdlib.h> /* div_t, ldiv_t */
3575
3576 /*********************************************************************
3577 * div (MSVCRT.@)
3578 * VERSION
3579 * [i386] Windows binary compatible - returns the struct in eax/edx.
3580 */
3581 #ifdef __i386__
3582 unsigned __int64 CDECL div(int num, int denom)
3583 {
3584 union {
3585 div_t div;
3586 unsigned __int64 uint64;
3587 } ret;
3588
3589 ret.div.quot = num / denom;
3590 ret.div.rem = num % denom;
3591 return ret.uint64;
3592 }
3593 #else
3594 /*********************************************************************
3595 * div (MSVCRT.@)
3596 * VERSION
3597 * [!i386] Non-x86 can't run win32 apps so we don't need binary compatibility
3598 */
3599 div_t CDECL div(int num, int denom)
3600 {
3601 div_t ret;
3602
3603 ret.quot = num / denom;
3604 ret.rem = num % denom;
3605 return ret;
3606 }
3607 #endif /* ifdef __i386__ */
3608
3609
3610 /*********************************************************************
3611 * ldiv (MSVCRT.@)
3612 * VERSION
3613 * [i386] Windows binary compatible - returns the struct in eax/edx.
3614 */
3615 #ifdef __i386__
3616 unsigned __int64 CDECL ldiv(__msvcrt_long num, __msvcrt_long denom)
3617 {
3618 union {
3619 ldiv_t ldiv;
3620 unsigned __int64 uint64;
3621 } ret;
3622
3623 ret.ldiv.quot = num / denom;
3624 ret.ldiv.rem = num % denom;
3625 return ret.uint64;
3626 }
3627 #else
3628 /*********************************************************************
3629 * ldiv (MSVCRT.@)
3630 * VERSION
3631 * [!i386] Non-x86 can't run win32 apps so we don't need binary compatibility
3632 */
3633 ldiv_t CDECL ldiv(__msvcrt_long num, __msvcrt_long denom)
3634 {
3635 ldiv_t ret;
3636
3637 ret.quot = num / denom;
3638 ret.rem = num % denom;
3639 return ret;
3640 }
3641 #endif /* ifdef __i386__ */
3642
3643 #if _MSVCR_VER>=100
3644 /*********************************************************************
3645 * lldiv (MSVCR100.@)
3646 */
3647 lldiv_t CDECL lldiv(__int64 num, __int64 denom)
3648 {
3649 lldiv_t ret;
3650
3651 ret.quot = num / denom;
3652 ret.rem = num % denom;
3653
3654 return ret;
3655 }
3656 #endif
3657
3658 #ifdef __i386__
3659
3660 /*********************************************************************
3661 * _adjust_fdiv (MSVCRT.@)
3662 * Used by the MSVC compiler to work around the Pentium FDIV bug.
3663 */
3664 int MSVCRT__adjust_fdiv = 0;
3665
3666 /***********************************************************************
3667 * _adj_fdiv_m16i (MSVCRT.@)
3668 *
3669 * NOTE
3670 * I _think_ this function is intended to work around the Pentium
3671 * fdiv bug.
3672 */
3673 void __stdcall _adj_fdiv_m16i( short arg )
3674 {
3675 TRACE("(): stub\n");
3676 }
3677
3678 /***********************************************************************
3679 * _adj_fdiv_m32 (MSVCRT.@)
3680 *
3681 * NOTE
3682 * I _think_ this function is intended to work around the Pentium
3683 * fdiv bug.
3684 */
3685 void __stdcall _adj_fdiv_m32( unsigned int arg )
3686 {
3687 TRACE("(): stub\n");
3688 }
3689
3690 /***********************************************************************
3691 * _adj_fdiv_m32i (MSVCRT.@)
3692 *
3693 * NOTE
3694 * I _think_ this function is intended to work around the Pentium
3695 * fdiv bug.
3696 */
3697 void __stdcall _adj_fdiv_m32i( int arg )
3698 {
3699 TRACE("(): stub\n");
3700 }
3701
3702 /***********************************************************************
3703 * _adj_fdiv_m64 (MSVCRT.@)
3704 *
3705 * NOTE
3706 * I _think_ this function is intended to work around the Pentium
3707 * fdiv bug.
3708 */
3709 void __stdcall _adj_fdiv_m64( unsigned __int64 arg )
3710 {
3711 TRACE("(): stub\n");
3712 }
3713
3714 /***********************************************************************
3715 * _adj_fdiv_r (MSVCRT.@)
3716 * FIXME
3717 * This function is likely to have the wrong number of arguments.
3718 *
3719 * NOTE
3720 * I _think_ this function is intended to work around the Pentium
3721 * fdiv bug.
3722 */
3723 void _adj_fdiv_r(void)
3724 {
3725 TRACE("(): stub\n");
3726 }
3727
3728 /***********************************************************************
3729 * _adj_fdivr_m16i (MSVCRT.@)
3730 *
3731 * NOTE
3732 * I _think_ this function is intended to work around the Pentium
3733 * fdiv bug.
3734 */
3735 void __stdcall _adj_fdivr_m16i( short arg )
3736 {
3737 TRACE("(): stub\n");
3738 }
3739
3740 /***********************************************************************
3741 * _adj_fdivr_m32 (MSVCRT.@)
3742 *
3743 * NOTE
3744 * I _think_ this function is intended to work around the Pentium
3745 * fdiv bug.
3746 */
3747 void __stdcall _adj_fdivr_m32( unsigned int arg )
3748 {
3749 TRACE("(): stub\n");
3750 }
3751
3752 /***********************************************************************
3753 * _adj_fdivr_m32i (MSVCRT.@)
3754 *
3755 * NOTE
3756 * I _think_ this function is intended to work around the Pentium
3757 * fdiv bug.
3758 */
3759 void __stdcall _adj_fdivr_m32i( int arg )
3760 {
3761 TRACE("(): stub\n");
3762 }
3763
3764 /***********************************************************************
3765 * _adj_fdivr_m64 (MSVCRT.@)
3766 *
3767 * NOTE
3768 * I _think_ this function is intended to work around the Pentium
3769 * fdiv bug.
3770 */
3771 void __stdcall _adj_fdivr_m64( unsigned __int64 arg )
3772 {
3773 TRACE("(): stub\n");
3774 }
3775
3776 /***********************************************************************
3777 * _adj_fpatan (MSVCRT.@)
3778 * FIXME
3779 * This function is likely to have the wrong number of arguments.
3780 *
3781 * NOTE
3782 * I _think_ this function is intended to work around the Pentium
3783 * fdiv bug.
3784 */
3785 void _adj_fpatan(void)
3786 {
3787 TRACE("(): stub\n");
3788 }
3789
3790 /***********************************************************************
3791 * _adj_fprem (MSVCRT.@)
3792 * FIXME
3793 * This function is likely to have the wrong number of arguments.
3794 *
3795 * NOTE
3796 * I _think_ this function is intended to work around the Pentium
3797 * fdiv bug.
3798 */
3799 void _adj_fprem(void)
3800 {
3801 TRACE("(): stub\n");
3802 }
3803
3804 /***********************************************************************
3805 * _adj_fprem1 (MSVCRT.@)
3806 * FIXME
3807 * This function is likely to have the wrong number of arguments.
3808 *
3809 * NOTE
3810 * I _think_ this function is intended to work around the Pentium
3811 * fdiv bug.
3812 */
3813 void _adj_fprem1(void)
3814 {
3815 TRACE("(): stub\n");
3816 }
3817
3818 /***********************************************************************
3819 * _adj_fptan (MSVCRT.@)
3820 * FIXME
3821 * This function is likely to have the wrong number of arguments.
3822 *
3823 * NOTE
3824 * I _think_ this function is intended to work around the Pentium
3825 * fdiv bug.
3826 */
3827 void _adj_fptan(void)
3828 {
3829 TRACE("(): stub\n");
3830 }
3831
3832 /***********************************************************************
3833 * _safe_fdiv (MSVCRT.@)
3834 * FIXME
3835 * This function is likely to have the wrong number of arguments.
3836 *
3837 * NOTE
3838 * I _think_ this function is intended to work around the Pentium
3839 * fdiv bug.
3840 */
3841 void _safe_fdiv(void)
3842 {
3843 TRACE("(): stub\n");
3844 }
3845
3846 /***********************************************************************
3847 * _safe_fdivr (MSVCRT.@)
3848 * FIXME
3849 * This function is likely to have the wrong number of arguments.
3850 *
3851 * NOTE
3852 * I _think_ this function is intended to work around the Pentium
3853 * fdiv bug.
3854 */
3855 void _safe_fdivr(void)
3856 {
3857 TRACE("(): stub\n");
3858 }
3859
3860 /***********************************************************************
3861 * _safe_fprem (MSVCRT.@)
3862 * FIXME
3863 * This function is likely to have the wrong number of arguments.
3864 *
3865 * NOTE
3866 * I _think_ this function is intended to work around the Pentium
3867 * fdiv bug.
3868 */
3869 void _safe_fprem(void)
3870 {
3871 TRACE("(): stub\n");
3872 }
3873
3874 /***********************************************************************
3875 * _safe_fprem1 (MSVCRT.@)
3876 *
3877 * FIXME
3878 * This function is likely to have the wrong number of arguments.
3879 *
3880 * NOTE
3881 * I _think_ this function is intended to work around the Pentium
3882 * fdiv bug.
3883 */
3884 void _safe_fprem1(void)
3885 {
3886 TRACE("(): stub\n");
3887 }
3888
3889 /***********************************************************************
3890 * __libm_sse2_acos (MSVCRT.@)
3891 */
3892 void __cdecl __libm_sse2_acos(void)
3893 {
3894 double d;
3895 __asm__ __volatile__( "movq %%xmm0,%0" : "=m" (d) );
3896 d = acos( d );
3897 __asm__ __volatile__( "movq %0,%%xmm0" : : "m" (d) );
3898 }
3899
3900 /***********************************************************************
3901 * __libm_sse2_acosf (MSVCRT.@)
3902 */
3903 void __cdecl __libm_sse2_acosf(void)
3904 {
3905 float f;
3906 __asm__ __volatile__( "movd %%xmm0,%0" : "=g" (f) );
3907 f = acosf( f );
3908 __asm__ __volatile__( "movd %0,%%xmm0" : : "g" (f) );
3909 }
3910
3911 /***********************************************************************
3912 * __libm_sse2_asin (MSVCRT.@)
3913 */
3914 void __cdecl __libm_sse2_asin(void)
3915 {
3916 double d;
3917 __asm__ __volatile__( "movq %%xmm0,%0" : "=m" (d) );
3918 d = asin( d );
3919 __asm__ __volatile__( "movq %0,%%xmm0" : : "m" (d) );
3920 }
3921
3922 /***********************************************************************
3923 * __libm_sse2_asinf (MSVCRT.@)
3924 */
3925 void __cdecl __libm_sse2_asinf(void)
3926 {
3927 float f;
3928 __asm__ __volatile__( "movd %%xmm0,%0" : "=g" (f) );
3929 f = asinf( f );
3930 __asm__ __volatile__( "movd %0,%%xmm0" : : "g" (f) );
3931 }
3932
3933 /***********************************************************************
3934 * __libm_sse2_atan (MSVCRT.@)
3935 */
3936 void __cdecl __libm_sse2_atan(void)
3937 {
3938 double d;
3939 __asm__ __volatile__( "movq %%xmm0,%0" : "=m" (d) );
3940 d = atan( d );
3941 __asm__ __volatile__( "movq %0,%%xmm0" : : "m" (d) );
3942 }
3943
3944 /***********************************************************************
3945 * __libm_sse2_atan2 (MSVCRT.@)
3946 */
3947 void __cdecl __libm_sse2_atan2(void)
3948 {
3949 double d1, d2;
3950 __asm__ __volatile__( "movq %%xmm0,%0; movq %%xmm1,%1 " : "=m" (d1), "=m" (d2) );
3951 d1 = atan2( d1, d2 );
3952 __asm__ __volatile__( "movq %0,%%xmm0" : : "m" (d1) );
3953 }
3954
3955 /***********************************************************************
3956 * __libm_sse2_atanf (MSVCRT.@)
3957 */
3958 void __cdecl __libm_sse2_atanf(void)
3959 {
3960 float f;
3961 __asm__ __volatile__( "movd %%xmm0,%0" : "=g" (f) );
3962 f = atanf( f );
3963 __asm__ __volatile__( "movd %0,%%xmm0" : : "g" (f) );
3964 }
3965
3966 /***********************************************************************
3967 * __libm_sse2_cos (MSVCRT.@)
3968 */
3969 void __cdecl __libm_sse2_cos(void)
3970 {
3971 double d;
3972 __asm__ __volatile__( "movq %%xmm0,%0" : "=m" (d) );
3973 d = cos( d );
3974 __asm__ __volatile__( "movq %0,%%xmm0" : : "m" (d) );
3975 }
3976
3977 /***********************************************************************
3978 * __libm_sse2_cosf (MSVCRT.@)
3979 */
3980 void __cdecl __libm_sse2_cosf(void)
3981 {
3982 float f;
3983 __asm__ __volatile__( "movd %%xmm0,%0" : "=g" (f) );
3984 f = cosf( f );
3985 __asm__ __volatile__( "movd %0,%%xmm0" : : "g" (f) );
3986 }
3987
3988 /***********************************************************************
3989 * __libm_sse2_exp (MSVCRT.@)
3990 */
3991 void __cdecl __libm_sse2_exp(void)
3992 {
3993 double d;
3994 __asm__ __volatile__( "movq %%xmm0,%0" : "=m" (d) );
3995 d = exp( d );
3996 __asm__ __volatile__( "movq %0,%%xmm0" : : "m" (d) );
3997 }
3998
3999 /***********************************************************************
4000 * __libm_sse2_expf (MSVCRT.@)
4001 */
4002 void __cdecl __libm_sse2_expf(void)
4003 {
4004 float f;
4005 __asm__ __volatile__( "movd %%xmm0,%0" : "=g" (f) );
4006 f = expf( f );
4007 __asm__ __volatile__( "movd %0,%%xmm0" : : "g" (f) );
4008 }
4009
4010 /***********************************************************************
4011 * __libm_sse2_log (MSVCRT.@)
4012 */
4013 void __cdecl __libm_sse2_log(void)
4014 {
4015 double d;
4016 __asm__ __volatile__( "movq %%xmm0,%0" : "=m" (d) );
4017 d = log( d );
4018 __asm__ __volatile__( "movq %0,%%xmm0" : : "m" (d) );
4019 }
4020
4021 /***********************************************************************
4022 * __libm_sse2_log10 (MSVCRT.@)
4023 */
4024 void __cdecl __libm_sse2_log10(void)
4025 {
4026 double d;
4027 __asm__ __volatile__( "movq %%xmm0,%0" : "=m" (d) );
4028 d = log10( d );
4029 __asm__ __volatile__( "movq %0,%%xmm0" : : "m" (d) );
4030 }
4031
4032 /***********************************************************************
4033 * __libm_sse2_log10f (MSVCRT.@)
4034 */
4035 void __cdecl __libm_sse2_log10f(void)
4036 {
4037 float f;
4038 __asm__ __volatile__( "movd %%xmm0,%0" : "=g" (f) );
4039 f = log10f( f );
4040 __asm__ __volatile__( "movd %0,%%xmm0" : : "g" (f) );
4041 }
4042
4043 /***********************************************************************
4044 * __libm_sse2_logf (MSVCRT.@)
4045 */
4046 void __cdecl __libm_sse2_logf(void)
4047 {
4048 float f;
4049 __asm__ __volatile__( "movd %%xmm0,%0" : "=g" (f) );
4050 f = logf( f );
4051 __asm__ __volatile__( "movd %0,%%xmm0" : : "g" (f) );
4052 }
4053
4054 /***********************************************************************
4055 * __libm_sse2_pow (MSVCRT.@)
4056 */
4057 void __cdecl __libm_sse2_pow(void)
4058 {
4059 double d1, d2;
4060 __asm__ __volatile__( "movq %%xmm0,%0; movq %%xmm1,%1 " : "=m" (d1), "=m" (d2) );
4061 d1 = pow( d1, d2 );
4062 __asm__ __volatile__( "movq %0,%%xmm0" : : "m" (d1) );
4063 }
4064
4065 /***********************************************************************
4066 * __libm_sse2_powf (MSVCRT.@)
4067 */
4068 void __cdecl __libm_sse2_powf(void)
4069 {
4070 float f1, f2;
4071 __asm__ __volatile__( "movd %%xmm0,%0; movd %%xmm1,%1" : "=g" (f1), "=g" (f2) );
4072 f1 = powf( f1, f2 );
4073 __asm__ __volatile__( "movd %0,%%xmm0" : : "g" (f1) );
4074 }
4075
4076 /***********************************************************************
4077 * __libm_sse2_sin (MSVCRT.@)
4078 */
4079 void __cdecl __libm_sse2_sin(void)
4080 {
4081 double d;
4082 __asm__ __volatile__( "movq %%xmm0,%0" : "=m" (d) );
4083 d = sin( d );
4084 __asm__ __volatile__( "movq %0,%%xmm0" : : "m" (d) );
4085 }
4086
4087 /***********************************************************************
4088 * __libm_sse2_sinf (MSVCRT.@)
4089 */
4090 void __cdecl __libm_sse2_sinf(void)
4091 {
4092 float f;
4093 __asm__ __volatile__( "movd %%xmm0,%0" : "=g" (f) );
4094 f = sinf( f );
4095 __asm__ __volatile__( "movd %0,%%xmm0" : : "g" (f) );
4096 }
4097
4098 /***********************************************************************
4099 * __libm_sse2_tan (MSVCRT.@)
4100 */
4101 void __cdecl __libm_sse2_tan(void)
4102 {
4103 double d;
4104 __asm__ __volatile__( "movq %%xmm0,%0" : "=m" (d) );
4105 d = tan( d );
4106 __asm__ __volatile__( "movq %0,%%xmm0" : : "m" (d) );
4107 }
4108
4109 /***********************************************************************
4110 * __libm_sse2_tanf (MSVCRT.@)
4111 */
4112 void __cdecl __libm_sse2_tanf(void)
4113 {
4114 float f;
4115 __asm__ __volatile__( "movd %%xmm0,%0" : "=g" (f) );
4116 f = tanf( f );
4117 __asm__ __volatile__( "movd %0,%%xmm0" : : "g" (f) );
4118 }
4119
4120 /***********************************************************************
4121 * __libm_sse2_sqrt_precise (MSVCR110.@)
4122 */
4123 void __cdecl __libm_sse2_sqrt_precise(void)
4124 {
4125 unsigned int cw;
4126 double d;
4127
4128 __asm__ __volatile__( "movq %%xmm0,%0" : "=m" (d) );
4129 __control87_2(0, 0, NULL, &cw);
4130 if (cw & _MCW_RC)
4131 {
4132 d = sqrt(d);
4133 __asm__ __volatile__( "movq %0,%%xmm0" : : "m" (d) );
4134 return;
4135 }
4136
4137 if (!sqrt_validate(&d, FALSE))
4138 {
4139 __asm__ __volatile__( "movq %0,%%xmm0" : : "m" (d) );
4140 return;
4141 }
4142 __asm__ __volatile__( "call " __ASM_NAME( "sse2_sqrt" ) );
4143 }
4144 #endif /* __i386__ */
4145
4146 /*********************************************************************
4147 * cbrt (MSVCR120.@)
4148 */
4149 double CDECL cbrt(double x)
4150 {
4151 return unix_funcs->cbrt( x );
4152 }
4153
4154 /*********************************************************************
4155 * cbrtf (MSVCR120.@)
4156 */
4157 float CDECL cbrtf(float x)
4158 {
4159 return unix_funcs->cbrtf( x );
4160 }
4161
4162 /*********************************************************************
4163 * exp2 (MSVCR120.@)
4164 */
4165 double CDECL exp2(double x)
4166 {
4167 double ret = unix_funcs->exp2( x );
4168 if (isfinite(x) && !isfinite(ret)) *_errno() = ERANGE;
4169 return ret;
4170 }
4171
4172 /*********************************************************************
4173 * exp2f (MSVCR120.@)
4174 */
4175 float CDECL exp2f(float x)
4176 {
4177 float ret = unix_funcs->exp2f( x );
4178 if (isfinite(x) && !isfinite(ret)) *_errno() = ERANGE;
4179 return ret;
4180 }
4181
4182 /*********************************************************************
4183 * expm1 (MSVCR120.@)
4184 */
4185 double CDECL expm1(double x)
4186 {
4187 double ret = unix_funcs->expm1( x );
4188 if (isfinite(x) && !isfinite(ret)) *_errno() = ERANGE;
4189 return ret;
4190 }
4191
4192 /*********************************************************************
4193 * expm1f (MSVCR120.@)
4194 */
4195 float CDECL expm1f(float x)
4196 {
4197 float ret = unix_funcs->expm1f( x );
4198 if (isfinite(x) && !isfinite(ret)) *_errno() = ERANGE;
4199 return ret;
4200 }
4201
4202 /*********************************************************************
4203 * log1p (MSVCR120.@)
4204 */
4205 double CDECL log1p(double x)
4206 {
4207 if (x < -1) *_errno() = EDOM;
4208 else if (x == -1) *_errno() = ERANGE;
4209 return unix_funcs->log1p( x );
4210 }
4211
4212 /*********************************************************************
4213 * log1pf (MSVCR120.@)
4214 */
4215 float CDECL log1pf(float x)
4216 {
4217 if (x < -1) *_errno() = EDOM;
4218 else if (x == -1) *_errno() = ERANGE;
4219 return unix_funcs->log1pf( x );
4220 }
4221
4222 /*********************************************************************
4223 * log2 (MSVCR120.@)
4224 */
4225 double CDECL log2(double x)
4226 {
4227 if (x < 0) *_errno() = EDOM;
4228 else if (x == 0) *_errno() = ERANGE;
4229 return unix_funcs->log2( x );
4230 }
4231
4232 /*********************************************************************
4233 * log2f (MSVCR120.@)
4234 */
4235 float CDECL log2f(float x)
4236 {
4237 if (x < 0) *_errno() = EDOM;
4238 else if (x == 0) *_errno() = ERANGE;
4239 return unix_funcs->log2f( x );
4240 }
4241
4242 /*********************************************************************
4243 * rint (MSVCR120.@)
4244 */
4245 double CDECL rint(double x)
4246 {
4247 return unix_funcs->rint(x);
4248 }
4249
4250 /*********************************************************************
4251 * rintf (MSVCR120.@)
4252 */
4253 float CDECL rintf(float x)
4254 {
4255 return unix_funcs->rintf(x);
4256 }
4257
4258 /*********************************************************************
4259 * lrint (MSVCR120.@)
4260 */
4261 __msvcrt_long CDECL lrint(double x)
4262 {
4263 return unix_funcs->lrint( x );
4264 }
4265
4266 /*********************************************************************
4267 * lrintf (MSVCR120.@)
4268 */
4269 __msvcrt_long CDECL lrintf(float x)
4270 {
4271 return unix_funcs->lrintf( x );
4272 }
4273
4274 /*********************************************************************
4275 * llrint (MSVCR120.@)
4276 */
4277 __int64 CDECL llrint(double x)
4278 {
4279 return unix_funcs->llrint( x );
4280 }
4281
4282 /*********************************************************************
4283 * llrintf (MSVCR120.@)
4284 */
4285 __int64 CDECL llrintf(float x)
4286 {
4287 return unix_funcs->llrintf( x );
4288 }
4289
4290 /*********************************************************************
4291 * _fdclass (MSVCR120.@)
4292 *
4293 * Copied from musl: src/math/__fpclassifyf.c
4294 */
4295 short CDECL _fdclass(float x)
4296 {
4297 union { float f; UINT32 i; } u = { x };
4298 int e = u.i >> 23 & 0xff;
4299
4300 if (!e) return u.i << 1 ? FP_SUBNORMAL : FP_ZERO;
4301 if (e == 0xff) return u.i << 9 ? FP_NAN : FP_INFINITE;
4302 return FP_NORMAL;
4303 }
4304
4305 /*********************************************************************
4306 * _dclass (MSVCR120.@)
4307 *
4308 * Copied from musl: src/math/__fpclassify.c
4309 */
4310 short CDECL _dclass(double x)
4311 {
4312 union { double f; UINT64 i; } u = { x };
4313 int e = u.i >> 52 & 0x7ff;
4314
4315 if (!e) return u.i << 1 ? FP_SUBNORMAL : FP_ZERO;
4316 if (e == 0x7ff) return (u.i << 12) ? FP_NAN : FP_INFINITE;
4317 return FP_NORMAL;
4318 }
4319
4320 #if _MSVCR_VER>=120
4321
4322 /*********************************************************************
4323 * round (MSVCR120.@)
4324 */
4325 double CDECL round(double x)
4326 {
4327 return unix_funcs->round(x);
4328 }
4329
4330 /*********************************************************************
4331 * roundf (MSVCR120.@)
4332 *
4333 * Copied from musl: src/math/roundf.c
4334 */
4335 float CDECL roundf(float x)
4336 {
4337 static const float toint = 1 / FLT_EPSILON;
4338
4339 unsigned int ix = *(unsigned int*)&x;
4340 int e = ix >> 23 & 0xff;
4341 float y;
4342
4343 if (e >= 0x7f + 23)
4344 return x;
4345 if (ix >> 31)
4346 x = -x;
4347 if (e < 0x7f - 1)
4348 return 0 * *(float*)&ix;
4349 y = fp_barrierf(x + toint) - toint - x;
4350 if (y > 0.5f)
4351 y = y + x - 1;
4352 else if (y <= -0.5f)
4353 y = y + x + 1;
4354 else
4355 y = y + x;
4356 if (ix >> 31)
4357 y = -y;
4358 return y;
4359 }
4360
4361 /*********************************************************************
4362 * lround (MSVCR120.@)
4363 *
4364 * Copied from musl: src/math/lround.c
4365 */
4366 __msvcrt_long CDECL lround(double x)
4367 {
4368 double d = round(x);
4369 if (d != (double)(__msvcrt_long)d) {
4370 *_errno() = EDOM;
4371 return 0;
4372 }
4373 return d;
4374 }
4375
4376 /*********************************************************************
4377 * lroundf (MSVCR120.@)
4378 *
4379 * Copied from musl: src/math/lroundf.c
4380 */
4381 __msvcrt_long CDECL lroundf(float x)
4382 {
4383 float f = roundf(x);
4384 if (f != (float)(__msvcrt_long)f) {
4385 *_errno() = EDOM;
4386 return 0;
4387 }
4388 return f;
4389 }
4390
4391 /*********************************************************************
4392 * llround (MSVCR120.@)
4393 *
4394 * Copied from musl: src/math/llround.c
4395 */
4396 __int64 CDECL llround(double x)
4397 {
4398 double d = round(x);
4399 if (d != (double)(__int64)d) {
4400 *_errno() = EDOM;
4401 return 0;
4402 }
4403 return d;
4404 }
4405
4406 /*********************************************************************
4407 * llroundf (MSVCR120.@)
4408 *
4409 * Copied from musl: src/math/llroundf.c
4410 */
4411 __int64 CDECL llroundf(float x)
4412 {
4413 float f = roundf(x);
4414 if (f != (float)(__int64)f) {
4415 *_errno() = EDOM;
4416 return 0;
4417 }
4418 return f;
4419 }
4420
4421 /*********************************************************************
4422 * trunc (MSVCR120.@)
4423 */
4424 double CDECL trunc(double x)
4425 {
4426 return unix_funcs->trunc(x);
4427 }
4428
4429 /*********************************************************************
4430 * truncf (MSVCR120.@)
4431 */
4432 float CDECL truncf(float x)
4433 {
4434 return unix_funcs->truncf(x);
4435 }
4436
4437 /*********************************************************************
4438 * _dtest (MSVCR120.@)
4439 */
4440 short CDECL _dtest(double *x)
4441 {
4442 return _dclass(*x);
4443 }
4444
4445 /*********************************************************************
4446 * _fdtest (MSVCR120.@)
4447 */
4448 short CDECL _fdtest(float *x)
4449 {
4450 return _fdclass(*x);
4451 }
4452
4453 /*********************************************************************
4454 * erf (MSVCR120.@)
4455 */
4456 double CDECL erf(double x)
4457 {
4458 return unix_funcs->erf( x );
4459 }
4460
4461 /*********************************************************************
4462 * erff (MSVCR120.@)
4463 */
4464 float CDECL erff(float x)
4465 {
4466 return unix_funcs->erff( x );
4467 }
4468
4469 /*********************************************************************
4470 * erfc (MSVCR120.@)
4471 */
4472 double CDECL erfc(double x)
4473 {
4474 return unix_funcs->erfc( x );
4475 }
4476
4477 /*********************************************************************
4478 * erfcf (MSVCR120.@)
4479 */
4480 float CDECL erfcf(float x)
4481 {
4482 return unix_funcs->erfcf( x );
4483 }
4484
4485 /*********************************************************************
4486 * fmaxf (MSVCR120.@)
4487 */
4488 float CDECL fmaxf(float x, float y)
4489 {
4490 if(isnan(x))
4491 return y;
4492 if(isnan(y))
4493 return x;
4494 if(x==0 && y==0)
4495 return signbit(x) ? y : x;
4496 return x<y ? y : x;
4497 }
4498
4499 /*********************************************************************
4500 * fmax (MSVCR120.@)
4501 */
4502 double CDECL fmax(double x, double y)
4503 {
4504 if(isnan(x))
4505 return y;
4506 if(isnan(y))
4507 return x;
4508 if(x==0 && y==0)
4509 return signbit(x) ? y : x;
4510 return x<y ? y : x;
4511 }
4512
4513 /*********************************************************************
4514 * fdimf (MSVCR120.@)
4515 */
4516 float CDECL fdimf(float x, float y)
4517 {
4518 if(isnan(x))
4519 return x;
4520 if(isnan(y))
4521 return y;
4522 return x>y ? x-y : 0;
4523 }
4524
4525 /*********************************************************************
4526 * fdim (MSVCR120.@)
4527 */
4528 double CDECL fdim(double x, double y)
4529 {
4530 if(isnan(x))
4531 return x;
4532 if(isnan(y))
4533 return y;
4534 return x>y ? x-y : 0;
4535 }
4536
4537 /*********************************************************************
4538 * _fdsign (MSVCR120.@)
4539 */
4540 int CDECL _fdsign(float x)
4541 {
4542 union { float f; UINT32 i; } u = { x };
4543 return (u.i >> 16) & 0x8000;
4544 }
4545
4546 /*********************************************************************
4547 * _dsign (MSVCR120.@)
4548 */
4549 int CDECL _dsign(double x)
4550 {
4551 union { double f; UINT64 i; } u = { x };
4552 return (u.i >> 48) & 0x8000;
4553 }
4554
4555
4556 /*********************************************************************
4557 * _dpcomp (MSVCR120.@)
4558 */
4559 int CDECL _dpcomp(double x, double y)
4560 {
4561 if(isnan(x) || isnan(y))
4562 return 0;
4563
4564 if(x == y) return 2;
4565 return x < y ? 1 : 4;
4566 }
4567
4568 /*********************************************************************
4569 * _fdpcomp (MSVCR120.@)
4570 */
4571 int CDECL _fdpcomp(float x, float y)
4572 {
4573 return _dpcomp(x, y);
4574 }
4575
4576 /*********************************************************************
4577 * fminf (MSVCR120.@)
4578 */
4579 float CDECL fminf(float x, float y)
4580 {
4581 if(isnan(x))
4582 return y;
4583 if(isnan(y))
4584 return x;
4585 if(x==0 && y==0)
4586 return signbit(x) ? x : y;
4587 return x<y ? x : y;
4588 }
4589
4590 /*********************************************************************
4591 * fmin (MSVCR120.@)
4592 */
4593 double CDECL fmin(double x, double y)
4594 {
4595 if(isnan(x))
4596 return y;
4597 if(isnan(y))
4598 return x;
4599 if(x==0 && y==0)
4600 return signbit(x) ? x : y;
4601 return x<y ? x : y;
4602 }
4603
4604 /*********************************************************************
4605 * asinh (MSVCR120.@)
4606 */
4607 double CDECL asinh(double x)
4608 {
4609 return unix_funcs->asinh( x );
4610 }
4611
4612 /*********************************************************************
4613 * asinhf (MSVCR120.@)
4614 */
4615 float CDECL asinhf(float x)
4616 {
4617 return unix_funcs->asinhf( x );
4618 }
4619
4620 /*********************************************************************
4621 * acosh (MSVCR120.@)
4622 */
4623 double CDECL acosh(double x)
4624 {
4625 if (x < 1)
4626 {
4627 fenv_t env;
4628
4629 *_errno() = EDOM;
4630 fegetenv(&env);
4631 env._Fe_stat |= FE_INVALID;
4632 fesetenv(&env);
4633 return NAN;
4634 }
4635 return unix_funcs->acosh( x );
4636 }
4637
4638 /*********************************************************************
4639 * acoshf (MSVCR120.@)
4640 */
4641 float CDECL acoshf(float x)
4642 {
4643 if (x < 1)
4644 {
4645 fenv_t env;
4646
4647 *_errno() = EDOM;
4648 fegetenv(&env);
4649 env._Fe_stat |= FE_INVALID;
4650 fesetenv(&env);
4651 return NAN;
4652 }
4653 return unix_funcs->acoshf( x );
4654 }
4655
4656 /*********************************************************************
4657 * atanh (MSVCR120.@)
4658 */
4659 double CDECL atanh(double x)
4660 {
4661 double ret;
4662
4663 if (x > 1 || x < -1) {
4664 fenv_t env;
4665
4666 *_errno() = EDOM;
4667
4668 /* on Linux atanh returns -NAN in this case */
4669 fegetenv(&env);
4670 env._Fe_stat |= FE_INVALID;
4671 fesetenv(&env);
4672 return NAN;
4673 }
4674 ret = unix_funcs->atanh( x );
4675
4676 if (!isfinite(ret)) *_errno() = ERANGE;
4677 return ret;
4678 }
4679
4680 /*********************************************************************
4681 * atanhf (MSVCR120.@)
4682 */
4683 float CDECL atanhf(float x)
4684 {
4685 float ret;
4686
4687 if (x > 1 || x < -1) {
4688 fenv_t env;
4689
4690 *_errno() = EDOM;
4691
4692 fegetenv(&env);
4693 env._Fe_stat |= FE_INVALID;
4694 fesetenv(&env);
4695 return NAN;
4696 }
4697
4698 ret = unix_funcs->atanh( x );
4699
4700 if (!isfinite(ret)) *_errno() = ERANGE;
4701 return ret;
4702 }
4703
4704 #endif /* _MSVCR_VER>=120 */
4705
4706 /*********************************************************************
4707 * _scalb (MSVCRT.@)
4708 * scalbn (MSVCR120.@)
4709 * scalbln (MSVCR120.@)
4710 */
4711 double CDECL _scalb(double num, __msvcrt_long power)
4712 {
4713 return ldexp(num, power);
4714 }
4715
4716 /*********************************************************************
4717 * _scalbf (MSVCRT.@)
4718 * scalbnf (MSVCR120.@)
4719 * scalblnf (MSVCR120.@)
4720 */
4721 float CDECL _scalbf(float num, __msvcrt_long power)
4722 {
4723 return ldexp(num, power);
4724 }
4725
4726 #if _MSVCR_VER>=120
4727
4728 /*********************************************************************
4729 * remainder (MSVCR120.@)
4730 */
4731 double CDECL remainder(double x, double y)
4732 {
4733 /* this matches 64-bit Windows. 32-bit Windows is slightly different */
4734 if(!isfinite(x)) *_errno() = EDOM;
4735 if(isnan(y) || y==0.0) *_errno() = EDOM;
4736 return unix_funcs->remainder( x, y );
4737 }
4738
4739 /*********************************************************************
4740 * remainderf (MSVCR120.@)
4741 */
4742 float CDECL remainderf(float x, float y)
4743 {
4744 /* this matches 64-bit Windows. 32-bit Windows is slightly different */
4745 if(!isfinite(x)) *_errno() = EDOM;
4746 if(isnan(y) || y==0.0f) *_errno() = EDOM;
4747 return unix_funcs->remainderf( x, y );
4748 }
4749
4750 /*********************************************************************
4751 * remquo (MSVCR120.@)
4752 */
4753 double CDECL remquo(double x, double y, int *quo)
4754 {
4755 if(!isfinite(x)) *_errno() = EDOM;
4756 if(isnan(y) || y==0.0) *_errno() = EDOM;
4757 return unix_funcs->remquo( x, y, quo );
4758 }
4759
4760 /*********************************************************************
4761 * remquof (MSVCR120.@)
4762 */
4763 float CDECL remquof(float x, float y, int *quo)
4764 {
4765 if(!isfinite(x)) *_errno() = EDOM;
4766 if(isnan(y) || y==0.0f) *_errno() = EDOM;
4767 return unix_funcs->remquof( x, y, quo );
4768 }
4769
4770 /*********************************************************************
4771 * lgamma (MSVCR120.@)
4772 */
4773 double CDECL lgamma(double x)
4774 {
4775 return unix_funcs->lgamma( x );
4776 }
4777
4778 /*********************************************************************
4779 * lgammaf (MSVCR120.@)
4780 */
4781 float CDECL lgammaf(float x)
4782 {
4783 return unix_funcs->lgammaf( x );
4784 }
4785
4786 /*********************************************************************
4787 * tgamma (MSVCR120.@)
4788 */
4789 double CDECL tgamma(double x)
4790 {
4791 return unix_funcs->tgamma( x );
4792 }
4793
4794 /*********************************************************************
4795 * tgammaf (MSVCR120.@)
4796 */
4797 float CDECL tgammaf(float x)
4798 {
4799 return unix_funcs->tgammaf( x );
4800 }
4801
4802 /*********************************************************************
4803 * nan (MSVCR120.@)
4804 */
4805 double CDECL nan(const char *tagp)
4806 {
4807 /* Windows ignores input (MSDN) */
4808 return NAN;
4809 }
4810
4811 /*********************************************************************
4812 * nanf (MSVCR120.@)
4813 */
4814 float CDECL nanf(const char *tagp)
4815 {
4816 return NAN;
4817 }
4818
4819 /*********************************************************************
4820 * _except1 (MSVCR120.@)
4821 * TODO:
4822 * - find meaning of ignored cw and operation bits
4823 * - unk parameter
4824 */
4825 double CDECL _except1(DWORD fpe, _FP_OPERATION_CODE op, double arg, double res, DWORD cw, void *unk)
4826 {
4827 ULONG_PTR exception_arg;
4828 DWORD exception = 0;
4829 fenv_t env;
4830 DWORD fpword = 0;
4831 WORD operation;
4832
4833 TRACE("(%x %x %lf %lf %x %p)\n", fpe, op, arg, res, cw, unk);
4834
4835 #ifdef _WIN64
4836 cw = ((cw >> 7) & 0x3f) | ((cw >> 3) & 0xc00);
4837 #endif
4838 operation = op << 5;
4839 exception_arg = (ULONG_PTR)&operation;
4840
4841 fegetenv(&env);
4842
4843 if (fpe & 0x1) { /* overflow */
4844 if ((fpe == 0x1 && (cw & 0x8)) || (fpe==0x11 && (cw & 0x28))) {
4845 /* 32-bit version also sets SW_INEXACT here */
4846 env._Fe_stat |= FE_OVERFLOW;
4847 if (fpe & 0x10) env._Fe_stat |= FE_INEXACT;
4848 res = signbit(res) ? -INFINITY : INFINITY;
4849 } else {
4850 exception = EXCEPTION_FLT_OVERFLOW;
4851 }
4852 } else if (fpe & 0x2) { /* underflow */
4853 if ((fpe == 0x2 && (cw & 0x10)) || (fpe==0x12 && (cw & 0x30))) {
4854 env._Fe_stat |= FE_UNDERFLOW;
4855 if (fpe & 0x10) env._Fe_stat |= FE_INEXACT;
4856 res = signbit(res) ? -0.0 : 0.0;
4857 } else {
4858 exception = EXCEPTION_FLT_UNDERFLOW;
4859 }
4860 } else if (fpe & 0x4) { /* zerodivide */
4861 if ((fpe == 0x4 && (cw & 0x4)) || (fpe==0x14 && (cw & 0x24))) {
4862 env._Fe_stat |= FE_DIVBYZERO;
4863 if (fpe & 0x10) env._Fe_stat |= FE_INEXACT;
4864 } else {
4865 exception = EXCEPTION_FLT_DIVIDE_BY_ZERO;
4866 }
4867 } else if (fpe & 0x8) { /* invalid */
4868 if (fpe == 0x8 && (cw & 0x1)) {
4869 env._Fe_stat |= FE_INVALID;
4870 } else {
4871 exception = EXCEPTION_FLT_INVALID_OPERATION;
4872 }
4873 } else if (fpe & 0x10) { /* inexact */
4874 if (fpe == 0x10 && (cw & 0x20)) {
4875 env._Fe_stat |= FE_INEXACT;
4876 } else {
4877 exception = EXCEPTION_FLT_INEXACT_RESULT;
4878 }
4879 }
4880
4881 if (exception)
4882 env._Fe_stat = 0;
4883 fesetenv(&env);
4884 if (exception)
4885 RaiseException(exception, 0, 1, &exception_arg);
4886
4887 if (cw & 0x1) fpword |= _EM_INVALID;
4888 if (cw & 0x2) fpword |= _EM_DENORMAL;
4889 if (cw & 0x4) fpword |= _EM_ZERODIVIDE;
4890 if (cw & 0x8) fpword |= _EM_OVERFLOW;
4891 if (cw & 0x10) fpword |= _EM_UNDERFLOW;
4892 if (cw & 0x20) fpword |= _EM_INEXACT;
4893 switch (cw & 0xc00)
4894 {
4895 case 0xc00: fpword |= _RC_UP|_RC_DOWN; break;
4896 case 0x800: fpword |= _RC_UP; break;
4897 case 0x400: fpword |= _RC_DOWN; break;
4898 }
4899 switch (cw & 0x300)
4900 {
4901 case 0x0: fpword |= _PC_24; break;
4902 case 0x200: fpword |= _PC_53; break;
4903 case 0x300: fpword |= _PC_64; break;
4904 }
4905 if (cw & 0x1000) fpword |= _IC_AFFINE;
4906 _control87(fpword, 0xffffffff);
4907
4908 return res;
4909 }
4910
4911 _Dcomplex* CDECL _Cbuild(_Dcomplex *ret, double r, double i)
4912 {
4913 ret->_Val[0] = r;
4914 ret->_Val[1] = i;
4915 return ret;
4916 }
4917
4918 double CDECL MSVCR120_creal(_Dcomplex z)
4919 {
4920 return z._Val[0];
4921 }
4922
4923 /*********************************************************************
4924 * ilogb (MSVCR120.@)
4925 *
4926 * Copied from musl: src/math/ilogb.c
4927 */
4928 int CDECL ilogb(double x)
4929 {
4930 union { double f; UINT64 i; } u = { x };
4931 int e = u.i >> 52 & 0x7ff;
4932
4933 if (!e)
4934 {
4935 u.i <<= 12;
4936 if (u.i == 0) return FP_ILOGB0;
4937 /* subnormal x */
4938 for (e = -0x3ff; u.i >> 63 == 0; e--, u.i <<= 1);
4939 return e;
4940 }
4941 if (e == 0x7ff) return u.i << 12 ? FP_ILOGBNAN : INT_MAX;
4942 return e - 0x3ff;
4943 }
4944
4945 /*********************************************************************
4946 * ilogbf (MSVCR120.@)
4947 *
4948 * Copied from musl: src/math/ilogbf.c
4949 */
4950 int CDECL ilogbf(float x)
4951 {
4952 union { float f; UINT32 i; } u = { x };
4953 int e = u.i >> 23 & 0xff;
4954
4955 if (!e)
4956 {
4957 u.i <<= 9;
4958 if (u.i == 0) return FP_ILOGB0;
4959 /* subnormal x */
4960 for (e = -0x7f; u.i >> 31 == 0; e--, u.i <<= 1);
4961 return e;
4962 }
4963 if (e == 0xff) return u.i << 9 ? FP_ILOGBNAN : INT_MAX;
4964 return e - 0x7f;
4965 }
4966 #endif /* _MSVCR_VER>=120 */
Windows API implementation