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