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