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