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