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