From: Alexandre Julliard julliard@winehq.org
--- dlls/msvcrt/math.c | 17 ----------------- 1 file changed, 17 deletions(-)
diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c index c0d9baab361..e87441d1427 100644 --- a/dlls/msvcrt/math.c +++ b/dlls/msvcrt/math.c @@ -9796,23 +9796,6 @@ float CDECL tgammaf(float x) return tgamma(x); }
-/********************************************************************* - * nan (MSVCR120.@) - */ -double CDECL nan(const char *tagp) -{ - /* Windows ignores input (MSDN) */ - return NAN; -} - -/********************************************************************* - * nanf (MSVCR120.@) - */ -float CDECL nanf(const char *tagp) -{ - return NAN; -} - /********************************************************************* * _except1 (MSVCR120.@) * TODO:
From: Alexandre Julliard julliard@winehq.org
--- dlls/msvcrt/math.c | 44 ++++++-------------------------------------- 1 file changed, 6 insertions(+), 38 deletions(-)
diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c index e87441d1427..61c95512448 100644 --- a/dlls/msvcrt/math.c +++ b/dlls/msvcrt/math.c @@ -295,38 +295,6 @@ float CDECL _logbf(float x)
#endif
-/* Copied from musl: src/math/scalbn.c */ -static double __scalbn(double x, int n) -{ - union {double f; UINT64 i;} u; - double y = x; - - if (n > 1023) { - y *= 0x1p1023; - n -= 1023; - if (n > 1023) { - y *= 0x1p1023; - n -= 1023; - if (n > 1023) - n = 1023; - } - } else if (n < -1022) { - /* make sure final n < -53 to avoid double - rounding in the subnormal range */ - y *= 0x1p-1022 * 0x1p53; - n += 1022 - 53; - if (n < -1022) { - y *= 0x1p-1022 * 0x1p53; - n += 1022 - 53; - if (n < -1022) - n = -1022; - } - } - u.i = (UINT64)(0x3ff + n) << 52; - x = y * u.f; - return x; -} - /* Copied from musl: src/math/__rem_pio2_large.c */ static int __rem_pio2_large(double *x, double *y, int e0, int nx, int prec) { @@ -391,7 +359,7 @@ recompute: }
/* compute n */ - z = __scalbn(z, q0); /* actual value of z */ + z = scalbn(z, q0); /* actual value of z */ z -= 8.0 * floor(z * 0.125); /* trim off integer >= 8 */ n = (INT32)z; z -= (double)n; @@ -431,7 +399,7 @@ recompute: if (ih == 2) { z = 1.0 - z; if (carry != 0) - z -= __scalbn(1.0, q0); + z -= scalbn(1.0, q0); } }
@@ -462,7 +430,7 @@ recompute: q0 -= 24; } } else { /* break z into 24-bit if necessary */ - z = __scalbn(z, -q0); + z = scalbn(z, -q0); if (z >= 0x1p24) { fw = (double)(INT32)(0x1p-24 * z); iq[jz] = (INT32)(z - 0x1p24 * fw); @@ -474,7 +442,7 @@ recompute: }
/* convert integer "bit" chunk to floating-point value */ - fw = __scalbn(1.0, q0); + fw = scalbn(1.0, q0); for (i = jz; i >= 0; i--) { q[i] = fw * (double)iq[i]; fw *= 0x1p-24; @@ -4853,7 +4821,7 @@ double CDECL fma( double x, double y, double z ) r = i; } } - return __scalbn(r, e); + return scalbn(r, e); }
/********************************************************************* @@ -5338,7 +5306,7 @@ int * CDECL __fpecode(void) */ double CDECL ldexp(double num, int exp) { - double z = __scalbn(num, exp); + double z = scalbn(num, exp);
if (isfinite(num) && !isfinite(z)) return math_error(_OVERFLOW, "ldexp", num, exp, z);
From: Alexandre Julliard julliard@winehq.org
--- dlls/msvcrt/math.c | 55 ++----------------------------------- libs/musl/src/math/ilogb.c | 3 -- libs/musl/src/math/ilogbf.c | 3 -- 3 files changed, 2 insertions(+), 59 deletions(-)
diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c index 61c95512448..a000d9196df 100644 --- a/dlls/msvcrt/math.c +++ b/dlls/msvcrt/math.c @@ -259,24 +259,6 @@ float CDECL _nextafterf( float x, float y ) return y; }
-/* Copied from musl: src/math/ilogbf.c */ -static int __ilogbf(float x) -{ - union { float f; UINT32 i; } u = { x }; - int e = u.i >> 23 & 0xff; - - if (!e) - { - u.i <<= 9; - if (u.i == 0) return FP_ILOGB0; - /* subnormal x */ - for (e = -0x7f; u.i >> 31 == 0; e--, u.i <<= 1); - return e; - } - if (e == 0xff) return u.i << 9 ? FP_ILOGBNAN : INT_MAX; - return e - 0x7f; -} - /********************************************************************* * _logbf (MSVCRT.@) * @@ -290,7 +272,7 @@ float CDECL _logbf(float x) *_errno() = ERANGE; return -1 / (x * x); } - return __ilogbf(x); + return ilogbf(x); }
#endif @@ -4522,24 +4504,6 @@ __int64 CDECL _abs64( __int64 n ) return n >= 0 ? n : -n; }
-/* Copied from musl: src/math/ilogb.c */ -static int __ilogb(double x) -{ - union { double f; UINT64 i; } u = { x }; - int e = u.i >> 52 & 0x7ff; - - if (!e) - { - u.i <<= 12; - if (u.i == 0) return FP_ILOGB0; - /* subnormal x */ - for (e = -0x3ff; u.i >> 63 == 0; e--, u.i <<= 1); - return e; - } - if (e == 0x7ff) return u.i << 12 ? FP_ILOGBNAN : INT_MAX; - return e - 0x3ff; -} - /********************************************************************* * _logb (MSVCRT.@) * @@ -4551,7 +4515,7 @@ double CDECL _logb(double x) return x * x; if (x == 0) return math_error(_SING, "_logb", x, 0, -1 / (x * x)); - return __ilogb(x); + return ilogb(x); }
/********************************************************************* @@ -9866,19 +9830,4 @@ double CDECL MSVCR120_creal(_Dcomplex z) return z._Val[0]; }
-/********************************************************************* - * ilogb (MSVCR120.@) - */ -int CDECL ilogb(double x) -{ - return __ilogb(x); -} - -/********************************************************************* - * ilogbf (MSVCR120.@) - */ -int CDECL ilogbf(float x) -{ - return __ilogbf(x); -} #endif /* _MSVCR_VER>=120 */ diff --git a/libs/musl/src/math/ilogb.c b/libs/musl/src/math/ilogb.c index 9351411a5aa..12738446ef0 100644 --- a/libs/musl/src/math/ilogb.c +++ b/libs/musl/src/math/ilogb.c @@ -3,7 +3,6 @@
int __cdecl ilogb(double x) { - #pragma STDC FENV_ACCESS ON union {double f; uint64_t i;} u = {x}; uint64_t i = u.i; int e = i>>52 & 0x7ff; @@ -11,7 +10,6 @@ int __cdecl ilogb(double x) if (!e) { i <<= 12; if (i == 0) { - FORCE_EVAL(0/0.0f); return FP_ILOGB0; } /* subnormal x */ @@ -19,7 +17,6 @@ int __cdecl ilogb(double x) return e; } if (e == 0x7ff) { - FORCE_EVAL(0/0.0f); return i<<12 ? FP_ILOGBNAN : INT_MAX; } return e - 0x3ff; diff --git a/libs/musl/src/math/ilogbf.c b/libs/musl/src/math/ilogbf.c index 65ee300531b..2eab40cd1a4 100644 --- a/libs/musl/src/math/ilogbf.c +++ b/libs/musl/src/math/ilogbf.c @@ -3,7 +3,6 @@
int __cdecl ilogbf(float x) { - #pragma STDC FENV_ACCESS ON union {float f; uint32_t i;} u = {x}; uint32_t i = u.i; int e = i>>23 & 0xff; @@ -11,7 +10,6 @@ int __cdecl ilogbf(float x) if (!e) { i <<= 9; if (i == 0) { - FORCE_EVAL(0/0.0f); return FP_ILOGB0; } /* subnormal x */ @@ -19,7 +17,6 @@ int __cdecl ilogbf(float x) return e; } if (e == 0xff) { - FORCE_EVAL(0/0.0f); return i<<9 ? FP_ILOGBNAN : INT_MAX; } return e - 0x7f;
From: Alexandre Julliard julliard@winehq.org
--- dlls/crtdll/crtdll.spec | 2 +- dlls/msvcr100/msvcr100.spec | 4 ++-- dlls/msvcr110/msvcr110.spec | 4 ++-- dlls/msvcr120/msvcr120.spec | 10 +++++----- dlls/msvcr70/msvcr70.spec | 2 +- dlls/msvcr71/msvcr71.spec | 2 +- dlls/msvcr80/msvcr80.spec | 4 ++-- dlls/msvcr90/msvcr90.spec | 4 ++-- dlls/msvcrt/math.c | 32 +------------------------------- dlls/msvcrt/msvcrt.spec | 4 ++-- dlls/msvcrtd/msvcrtd.spec | 2 +- dlls/ucrtbase/ucrtbase.spec | 20 ++++++++++---------- libs/musl/src/internal/libm.h | 2 ++ libs/musl/src/math/logb.c | 2 +- libs/musl/src/math/logbf.c | 5 ++++- 15 files changed, 37 insertions(+), 62 deletions(-)
diff --git a/dlls/crtdll/crtdll.spec b/dlls/crtdll/crtdll.spec index c66cea5bade..3a90ad7ef74 100644 --- a/dlls/crtdll/crtdll.spec +++ b/dlls/crtdll/crtdll.spec @@ -179,7 +179,7 @@ @ cdecl _loaddll(str) @ cdecl -arch=i386 _local_unwind2(ptr long) @ cdecl _locking(long long long) -@ cdecl _logb(double) +@ cdecl _logb(double) logb @ cdecl _lrotl(long long) MSVCRT__lrotl @ cdecl _lrotr(long long) MSVCRT__lrotr @ cdecl _lsearch(ptr ptr ptr long ptr) diff --git a/dlls/msvcr100/msvcr100.spec b/dlls/msvcr100/msvcr100.spec index a42fe349cf5..1b9abbfdec2 100644 --- a/dlls/msvcr100/msvcr100.spec +++ b/dlls/msvcr100/msvcr100.spec @@ -1052,8 +1052,8 @@ @ cdecl _lock(long) @ cdecl _lock_file(ptr) @ cdecl _locking(long long long) -@ cdecl _logb(double) -@ cdecl -arch=!i386 _logbf(float) +@ cdecl _logb(double) logb +@ cdecl -arch=!i386 _logbf(float) logbf @ cdecl -arch=i386 _longjmpex(ptr long) MSVCRT_longjmp @ cdecl _lrotl(long long) MSVCRT__lrotl @ cdecl _lrotr(long long) MSVCRT__lrotr diff --git a/dlls/msvcr110/msvcr110.spec b/dlls/msvcr110/msvcr110.spec index 2bccd8c6fe9..bba1222a44e 100644 --- a/dlls/msvcr110/msvcr110.spec +++ b/dlls/msvcr110/msvcr110.spec @@ -1409,8 +1409,8 @@ @ cdecl _lock(long) @ cdecl _lock_file(ptr) @ cdecl _locking(long long long) -@ cdecl _logb(double) -@ cdecl -arch=!i386 _logbf(float) +@ cdecl _logb(double) logb +@ cdecl -arch=!i386 _logbf(float) logbf @ cdecl -arch=i386 _longjmpex(ptr long) MSVCRT_longjmp @ cdecl _lrotl(long long) MSVCRT__lrotl @ cdecl _lrotr(long long) MSVCRT__lrotr diff --git a/dlls/msvcr120/msvcr120.spec b/dlls/msvcr120/msvcr120.spec index ed0f69b49bc..654fb4071f8 100644 --- a/dlls/msvcr120/msvcr120.spec +++ b/dlls/msvcr120/msvcr120.spec @@ -1420,8 +1420,8 @@ @ cdecl _lock(long) @ cdecl _lock_file(ptr) @ cdecl _locking(long long long) -@ cdecl _logb(double) -@ cdecl -arch=!i386 _logbf(float) +@ cdecl _logb(double) logb +@ cdecl -arch=!i386 _logbf(float) logbf @ cdecl -arch=i386 _longjmpex(ptr long) MSVCRT_longjmp @ cdecl _lrotl(long long) MSVCRT__lrotl @ cdecl _lrotr(long long) MSVCRT__lrotr @@ -2266,9 +2266,9 @@ @ cdecl log2(double) @ cdecl log2f(float) @ cdecl log2l(double) log2 -@ cdecl logb(double) _logb -@ cdecl logbf(float) _logbf -@ cdecl logbl(double) _logb +@ cdecl logb(double) +@ cdecl logbf(float) +@ cdecl logbl(double) logb @ cdecl -arch=i386,x86_64,arm,arm64 longjmp(ptr long) MSVCRT_longjmp @ cdecl lrint(double) @ cdecl lrintf(float) diff --git a/dlls/msvcr70/msvcr70.spec b/dlls/msvcr70/msvcr70.spec index cb3a92cf395..1e02dc95d2d 100644 --- a/dlls/msvcr70/msvcr70.spec +++ b/dlls/msvcr70/msvcr70.spec @@ -408,7 +408,7 @@ @ cdecl _localtime64(ptr) @ cdecl _lock(long) @ cdecl _locking(long long long) -@ cdecl _logb(double) +@ cdecl _logb(double) logb @ cdecl -arch=i386 _longjmpex(ptr long) MSVCRT_longjmp @ cdecl _lrotl(long long) MSVCRT__lrotl @ cdecl _lrotr(long long) MSVCRT__lrotr diff --git a/dlls/msvcr71/msvcr71.spec b/dlls/msvcr71/msvcr71.spec index a559b5cc376..1a075986ddd 100644 --- a/dlls/msvcr71/msvcr71.spec +++ b/dlls/msvcr71/msvcr71.spec @@ -403,7 +403,7 @@ @ cdecl _localtime64(ptr) @ cdecl _lock(long) @ cdecl _locking(long long long) -@ cdecl _logb(double) +@ cdecl _logb(double) logb @ cdecl -arch=i386 _longjmpex(ptr long) MSVCRT_longjmp @ cdecl _lrotl(long long) MSVCRT__lrotl @ cdecl _lrotr(long long) MSVCRT__lrotr diff --git a/dlls/msvcr80/msvcr80.spec b/dlls/msvcr80/msvcr80.spec index 06d0bdf82cd..2ea91c0ce8a 100644 --- a/dlls/msvcr80/msvcr80.spec +++ b/dlls/msvcr80/msvcr80.spec @@ -724,8 +724,8 @@ @ cdecl _lock(long) @ cdecl _lock_file(ptr) @ cdecl _locking(long long long) -@ cdecl _logb(double) -@ cdecl -arch=!i386 _logbf(float) +@ cdecl _logb(double) logb +@ cdecl -arch=!i386 _logbf(float) logbf @ cdecl -arch=i386 _longjmpex(ptr long) MSVCRT_longjmp @ cdecl _lrotl(long long) MSVCRT__lrotl @ cdecl _lrotr(long long) MSVCRT__lrotr diff --git a/dlls/msvcr90/msvcr90.spec b/dlls/msvcr90/msvcr90.spec index d020de98a02..fb274878d39 100644 --- a/dlls/msvcr90/msvcr90.spec +++ b/dlls/msvcr90/msvcr90.spec @@ -702,8 +702,8 @@ @ cdecl _lock(long) @ cdecl _lock_file(ptr) @ cdecl _locking(long long long) -@ cdecl _logb(double) -@ cdecl -arch=!i386 _logbf(float) +@ cdecl _logb(double) logb +@ cdecl -arch=!i386 _logbf(float) logbf @ cdecl -arch=i386 _longjmpex(ptr long) MSVCRT_longjmp @ cdecl _lrotl(long long) MSVCRT__lrotl @ cdecl _lrotr(long long) MSVCRT__lrotr diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c index a000d9196df..94da43ec3f6 100644 --- a/dlls/msvcrt/math.c +++ b/dlls/msvcrt/math.c @@ -131,7 +131,7 @@ int CDECL _matherr(struct _exception *e) }
-static double math_error(int type, const char *name, double arg1, double arg2, double retval) +double math_error(int type, const char *name, double arg1, double arg2, double retval) { struct _exception exception = {type, (char *)name, arg1, arg2, retval};
@@ -259,22 +259,6 @@ float CDECL _nextafterf( float x, float y ) return y; }
-/********************************************************************* - * _logbf (MSVCRT.@) - * - * Copied from musl: src/math/logbf.c - */ -float CDECL _logbf(float x) -{ - if (!isfinite(x)) - return x * x; - if (x == 0) { - *_errno() = ERANGE; - return -1 / (x * x); - } - return ilogbf(x); -} - #endif
/* Copied from musl: src/math/__rem_pio2_large.c */ @@ -4504,20 +4488,6 @@ __int64 CDECL _abs64( __int64 n ) return n >= 0 ? n : -n; }
-/********************************************************************* - * _logb (MSVCRT.@) - * - * Copied from musl: src/math/logb.c - */ -double CDECL _logb(double x) -{ - if (!isfinite(x)) - return x * x; - if (x == 0) - return math_error(_SING, "_logb", x, 0, -1 / (x * x)); - return ilogb(x); -} - /********************************************************************* * ceil (MSVCRT.@) * diff --git a/dlls/msvcrt/msvcrt.spec b/dlls/msvcrt/msvcrt.spec index 43259d7abad..f339757438f 100644 --- a/dlls/msvcrt/msvcrt.spec +++ b/dlls/msvcrt/msvcrt.spec @@ -670,8 +670,8 @@ @ cdecl _lock(long) @ cdecl _lock_file(ptr) @ cdecl _locking(long long long) -@ cdecl _logb(double) -@ cdecl -arch=!i386 _logbf(float) +@ cdecl _logb(double) logb +@ cdecl -arch=!i386 _logbf(float) logbf @ cdecl -arch=i386 _longjmpex(ptr long) MSVCRT_longjmp @ cdecl _lrotl(long long) MSVCRT__lrotl @ cdecl _lrotr(long long) MSVCRT__lrotr diff --git a/dlls/msvcrtd/msvcrtd.spec b/dlls/msvcrtd/msvcrtd.spec index 5cb09d92c74..68fa80d532f 100644 --- a/dlls/msvcrtd/msvcrtd.spec +++ b/dlls/msvcrtd/msvcrtd.spec @@ -385,7 +385,7 @@ @ cdecl -arch=i386 _local_unwind2(ptr long) @ cdecl _lock(long) @ cdecl _locking(long long long) -@ cdecl _logb(double) +@ cdecl _logb(double) logb @ cdecl -arch=i386 _longjmpex(ptr long) MSVCRT_longjmp @ cdecl _lrotl(long long) MSVCRT__lrotl @ cdecl _lrotr(long long) MSVCRT__lrotr diff --git a/dlls/ucrtbase/ucrtbase.spec b/dlls/ucrtbase/ucrtbase.spec index 87a266198bd..51137681f14 100644 --- a/dlls/ucrtbase/ucrtbase.spec +++ b/dlls/ucrtbase/ucrtbase.spec @@ -568,8 +568,8 @@ @ cdecl _lock_file(ptr) @ cdecl _lock_locales() @ cdecl _locking(long long long) -@ cdecl _logb(double) -@ cdecl -arch=!i386 _logbf(float) +@ cdecl _logb(double) logb +@ cdecl -arch=!i386 _logbf(float) logbf @ cdecl -arch=i386 _longjmpex(ptr long) MSVCRT_longjmp @ cdecl _lrotl(long long) MSVCRT__lrotl @ cdecl _lrotr(long long) MSVCRT__lrotr @@ -1142,8 +1142,8 @@ @ cdecl _o__localtime64_s(ptr ptr) _localtime64_s @ cdecl _o__lock_file(ptr) _lock_file @ cdecl _o__locking(long long long) _locking -@ cdecl _o__logb(double) _logb -@ cdecl -arch=!i386 _o__logbf(float) _logbf +@ cdecl _o__logb(double) logb +@ cdecl -arch=!i386 _o__logbf(float) logbf @ cdecl _o__lsearch(ptr ptr ptr long ptr) _lsearch @ stub _o__lsearch_s @ cdecl _o__lseek(long long long) _lseek @@ -1699,9 +1699,9 @@ @ cdecl _o_log2(double) log2 @ cdecl _o_log2f(float) log2f @ cdecl _o_log2l(double) log2 -@ cdecl _o_logb(double) _logb -@ cdecl _o_logbf(float) _logbf -@ cdecl _o_logbl(double) _logb +@ cdecl _o_logb(double) logb +@ cdecl _o_logbf(float) logbf +@ cdecl _o_logbl(double) logb @ cdecl -arch=!i386 _o_logf(float) logf @ cdecl _o_lrint(double) lrint @ cdecl _o_lrintf(float) lrintf @@ -2401,9 +2401,9 @@ @ cdecl log2(double) @ cdecl log2f(float) @ cdecl log2l(double) log2 -@ cdecl logb(double) _logb -@ cdecl logbf(float) _logbf -@ cdecl logbl(double) _logb +@ cdecl logb(double) +@ cdecl logbf(float) +@ cdecl logbl(double) logb @ cdecl -arch=!i386 logf(float) @ cdecl -arch=i386,x86_64,arm,arm64 longjmp(ptr long) MSVCRT_longjmp @ cdecl lrint(double) diff --git a/libs/musl/src/internal/libm.h b/libs/musl/src/internal/libm.h index 39233282750..a1e9bc08716 100644 --- a/libs/musl/src/internal/libm.h +++ b/libs/musl/src/internal/libm.h @@ -10,6 +10,8 @@ typedef float float_t; typedef double double_t;
+hidden double math_error(int type, const char *name, double arg1, double arg2, double retval); + #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __LITTLE_ENDIAN union ldshape { diff --git a/libs/musl/src/math/logb.c b/libs/musl/src/math/logb.c index 7c0423fb6c0..ce496e6bf18 100644 --- a/libs/musl/src/math/logb.c +++ b/libs/musl/src/math/logb.c @@ -13,6 +13,6 @@ double __cdecl logb(double x) if (!isfinite(x)) return x * x; if (x == 0) - return -1/(x*x); + return math_error(_SING, "_logb", x, 0, -1 / (x * x)); return ilogb(x); } diff --git a/libs/musl/src/math/logbf.c b/libs/musl/src/math/logbf.c index 743af679ab0..9f1728fc118 100644 --- a/libs/musl/src/math/logbf.c +++ b/libs/musl/src/math/logbf.c @@ -1,10 +1,13 @@ #include <math.h> +#include "libm.h"
float __cdecl logbf(float x) { if (!isfinite(x)) return x * x; - if (x == 0) + if (x == 0) { + errno = ERANGE; return -1/(x*x); + } return ilogbf(x); }
From: Alexandre Julliard julliard@winehq.org
--- dlls/msvcrt/math.c | 169 ----------------------------------- libs/musl/src/math/remquo.c | 2 + libs/musl/src/math/remquof.c | 2 + 3 files changed, 4 insertions(+), 169 deletions(-)
diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c index 94da43ec3f6..0c84dc23c8d 100644 --- a/dlls/msvcrt/math.c +++ b/dlls/msvcrt/math.c @@ -8994,175 +8994,6 @@ float CDECL remainderf(float x, float y) return remquof(x, y, &q); }
-/********************************************************************* - * remquo (MSVCR120.@) - * - * Copied from musl: src/math/remquo.c - */ -double CDECL remquo(double x, double y, int *quo) -{ - UINT64 uxi = *(UINT64*)&x; - UINT64 uyi = *(UINT64*)&y; - int ex = uxi >> 52 & 0x7ff; - int ey = uyi >> 52 & 0x7ff; - int sx = uxi >> 63; - int sy = uyi >> 63; - UINT32 q; - UINT64 i; - - *quo = 0; - if (y == 0 || isinf(x)) *_errno() = EDOM; - if (uyi << 1 == 0 || isnan(y) || ex == 0x7ff) - return (x * y) / (x * y); - if (uxi << 1 == 0) - return x; - - /* normalize x and y */ - if (!ex) { - for (i = uxi << 12; i >> 63 == 0; ex--, i <<= 1); - uxi <<= -ex + 1; - } else { - uxi &= -1ULL >> 12; - uxi |= 1ULL << 52; - } - if (!ey) { - for (i = uyi << 12; i >> 63 == 0; ey--, i <<= 1); - uyi <<= -ey + 1; - } else { - uyi &= -1ULL >> 12; - uyi |= 1ULL << 52; - } - - q = 0; - if (ex < ey) { - if (ex+1 == ey) - goto end; - return x; - } - - /* x mod y */ - for (; ex > ey; ex--) { - i = uxi - uyi; - if (i >> 63 == 0) { - uxi = i; - q++; - } - uxi <<= 1; - q <<= 1; - } - i = uxi - uyi; - if (i >> 63 == 0) { - uxi = i; - q++; - } - if (uxi == 0) - ex = -60; - else - for (; uxi >> 52 == 0; uxi <<= 1, ex--); -end: - /* scale result and decide between |x| and |x|-|y| */ - if (ex > 0) { - uxi -= 1ULL << 52; - uxi |= (UINT64)ex << 52; - } else { - uxi >>= -ex + 1; - } - x = *(double*)&uxi; - if (sy) - y = -y; - if (ex == ey || (ex + 1 == ey && (2 * x > y || (2 * x == y && q % 2)))) { - x -= y; - q++; - } - q &= 0x7fffffff; - *quo = sx ^ sy ? -(int)q : (int)q; - return sx ? -x : x; -} - -/********************************************************************* - * remquof (MSVCR120.@) - * - * Copied from musl: src/math/remquof.c - */ -float CDECL remquof(float x, float y, int *quo) -{ - UINT32 uxi = *(UINT32*)&x; - UINT32 uyi = *(UINT32*)&y; - int ex = uxi >> 23 & 0xff; - int ey = uyi >> 23 & 0xff; - int sx = uxi >> 31; - int sy = uyi>> 31; - UINT32 q, i; - - *quo = 0; - if (y == 0 || isinf(x)) *_errno() = EDOM; - if (uyi << 1 == 0 || isnan(y) || ex == 0xff) - return (x * y) / (x * y); - if (uxi << 1 == 0) - return x; - - /* normalize x and y */ - if (!ex) { - for (i = uxi << 9; i >> 31 == 0; ex--, i <<= 1); - uxi <<= -ex + 1; - } else { - uxi &= -1U >> 9; - uxi |= 1U << 23; - } - if (!ey) { - for (i = uyi << 9; i >> 31 == 0; ey--, i <<= 1); - uyi <<= -ey + 1; - } else { - uyi &= -1U >> 9; - uyi |= 1U << 23; - } - - q = 0; - if (ex < ey) { - if (ex + 1 == ey) - goto end; - return x; - } - - /* x mod y */ - for (; ex > ey; ex--) { - i = uxi - uyi; - if (i >> 31 == 0) { - uxi = i; - q++; - } - uxi <<= 1; - q <<= 1; - } - i = uxi - uyi; - if (i >> 31 == 0) { - uxi = i; - q++; - } - if (uxi == 0) - ex = -30; - else - for (; uxi >> 23 == 0; uxi <<= 1, ex--); -end: - /* scale result and decide between |x| and |x|-|y| */ - if (ex > 0) { - uxi -= 1U << 23; - uxi |= (UINT32)ex << 23; - } else { - uxi >>= -ex + 1; - } - x = *(float*)&uxi; - if (sy) - y = -y; - if (ex == ey || (ex + 1 == ey && (2 * x > y || (2 * x == y && q % 2)))) { - x -= y; - q++; - } - q &= 0x7fffffff; - *quo = sx ^ sy ? -(int)q : (int)q; - return sx ? -x : x; -} - /* sin(pi*x) assuming x > 2^-100, if sin(pi*x)==0 the sign is arbitrary */ static double sin_pi(double x) { diff --git a/libs/musl/src/math/remquo.c b/libs/musl/src/math/remquo.c index 6ba88688623..04bdce5027f 100644 --- a/libs/musl/src/math/remquo.c +++ b/libs/musl/src/math/remquo.c @@ -1,5 +1,6 @@ #include <math.h> #include <stdint.h> +#include "libm.h"
double __cdecl remquo(double x, double y, int *quo) { @@ -13,6 +14,7 @@ double __cdecl remquo(double x, double y, int *quo) uint64_t uxi = ux.i;
*quo = 0; + if (y == 0 || isinf(x)) errno = EDOM; if (uy.i<<1 == 0 || isnan(y) || ex == 0x7ff) return (x*y)/(x*y); if (ux.i<<1 == 0) diff --git a/libs/musl/src/math/remquof.c b/libs/musl/src/math/remquof.c index cd3bde759a8..26d48091db0 100644 --- a/libs/musl/src/math/remquof.c +++ b/libs/musl/src/math/remquof.c @@ -1,5 +1,6 @@ #include <math.h> #include <stdint.h> +#include "libm.h"
float __cdecl remquof(float x, float y, int *quo) { @@ -13,6 +14,7 @@ float __cdecl remquof(float x, float y, int *quo) uint32_t uxi = ux.i;
*quo = 0; + if (y == 0 || isinf(x)) errno = EDOM; if (uy.i<<1 == 0 || isnan(y) || ex == 0xff) return (x*y)/(x*y); if (ux.i<<1 == 0)
From: Alexandre Julliard julliard@winehq.org
--- dlls/msvcrt/math.c | 382 +-------------------------------- libs/musl/src/math/lgamma_r.c | 4 +- libs/musl/src/math/lgammaf_r.c | 4 +- 3 files changed, 8 insertions(+), 382 deletions(-)
diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c index 0c84dc23c8d..20ce8bf9599 100644 --- a/dlls/msvcrt/math.c +++ b/dlls/msvcrt/math.c @@ -576,6 +576,7 @@ static float __expm1f(float x) return y; }
+#ifndef __i386__ /* Copied from musl: src/math/__sindf.c */ static float __sindf(double x) { @@ -611,6 +612,7 @@ static float __cosdf(double x) return 1 + C0 * z; return 1.0 + z * (C0 + z * (C1 + z * (C2 + z * (C3 + z * C4)))); } +#endif
static const UINT64 exp2f_T[] = { 0x3ff0000000000000ULL, 0x3fefd9b0d3158574ULL, 0x3fefb5586cf9890fULL, 0x3fef9301d0125b51ULL, @@ -9016,386 +9018,6 @@ static double sin_pi(double x) } }
-/********************************************************************* - * lgamma (MSVCR120.@) - * - * Copied from musl: src/math/lgamma_r.c - */ -double CDECL lgamma(double x) -{ - static const double pi = 3.14159265358979311600e+00, - a0 = 7.72156649015328655494e-02, - a1 = 3.22467033424113591611e-01, - a2 = 6.73523010531292681824e-02, - a3 = 2.05808084325167332806e-02, - a4 = 7.38555086081402883957e-03, - a5 = 2.89051383673415629091e-03, - a6 = 1.19270763183362067845e-03, - a7 = 5.10069792153511336608e-04, - a8 = 2.20862790713908385557e-04, - a9 = 1.08011567247583939954e-04, - a10 = 2.52144565451257326939e-05, - a11 = 4.48640949618915160150e-05, - tc = 1.46163214496836224576e+00, - tf = -1.21486290535849611461e-01, - tt = -3.63867699703950536541e-18, - t0 = 4.83836122723810047042e-01, - t1 = -1.47587722994593911752e-01, - t2 = 6.46249402391333854778e-02, - t3 = -3.27885410759859649565e-02, - t4 = 1.79706750811820387126e-02, - t5 = -1.03142241298341437450e-02, - t6 = 6.10053870246291332635e-03, - t7 = -3.68452016781138256760e-03, - t8 = 2.25964780900612472250e-03, - t9 = -1.40346469989232843813e-03, - t10 = 8.81081882437654011382e-04, - t11 = -5.38595305356740546715e-04, - t12 = 3.15632070903625950361e-04, - t13 = -3.12754168375120860518e-04, - t14 = 3.35529192635519073543e-04, - u0 = -7.72156649015328655494e-02, - u1 = 6.32827064025093366517e-01, - u2 = 1.45492250137234768737e+00, - u3 = 9.77717527963372745603e-01, - u4 = 2.28963728064692451092e-01, - u5 = 1.33810918536787660377e-02, - v1 = 2.45597793713041134822e+00, - v2 = 2.12848976379893395361e+00, - v3 = 7.69285150456672783825e-01, - v4 = 1.04222645593369134254e-01, - v5 = 3.21709242282423911810e-03, - s0 = -7.72156649015328655494e-02, - s1 = 2.14982415960608852501e-01, - s2 = 3.25778796408930981787e-01, - s3 = 1.46350472652464452805e-01, - s4 = 2.66422703033638609560e-02, - s5 = 1.84028451407337715652e-03, - s6 = 3.19475326584100867617e-05, - r1 = 1.39200533467621045958e+00, - r2 = 7.21935547567138069525e-01, - r3 = 1.71933865632803078993e-01, - r4 = 1.86459191715652901344e-02, - r5 = 7.77942496381893596434e-04, - r6 = 7.32668430744625636189e-06, - w0 = 4.18938533204672725052e-01, - w1 = 8.33333333333329678849e-02, - w2 = -2.77777777728775536470e-03, - w3 = 7.93650558643019558500e-04, - w4 = -5.95187557450339963135e-04, - w5 = 8.36339918996282139126e-04, - w6 = -1.63092934096575273989e-03; - - union {double f; UINT64 i;} u = {x}; - double t, y, z, nadj, p, p1, p2, p3, q, r, w; - UINT32 ix; - int sign,i; - - /* purge off +-inf, NaN, +-0, tiny and negative arguments */ - sign = u.i >> 63; - ix = u.i >> 32 & 0x7fffffff; - if (ix >= 0x7ff00000) - return x * x; - if (ix < (0x3ff - 70) << 20) { /* |x|<2**-70, return -log(|x|) */ - if(sign) - x = -x; - return -log(x); - } - if (sign) { - x = -x; - t = sin_pi(x); - if (t == 0.0) { /* -integer */ - *_errno() = ERANGE; - return 1.0 / (x - x); - } - if (t <= 0.0) - t = -t; - nadj = log(pi / (t * x)); - } - - /* purge off 1 and 2 */ - if ((ix == 0x3ff00000 || ix == 0x40000000) && (UINT32)u.i == 0) - r = 0; - /* for x < 2.0 */ - else if (ix < 0x40000000) { - if (ix <= 0x3feccccc) { /* lgamma(x) = lgamma(x+1)-log(x) */ - r = -log(x); - if (ix >= 0x3FE76944) { - y = 1.0 - x; - i = 0; - } else if (ix >= 0x3FCDA661) { - y = x - (tc - 1.0); - i = 1; - } else { - y = x; - i = 2; - } - } else { - r = 0.0; - if (ix >= 0x3FFBB4C3) { /* [1.7316,2] */ - y = 2.0 - x; - i = 0; - } else if(ix >= 0x3FF3B4C4) { /* [1.23,1.73] */ - y = x - tc; - i = 1; - } else { - y = x - 1.0; - i = 2; - } - } - switch (i) { - case 0: - z = y * y; - p1 = a0 + z * (a2 + z * (a4 + z * (a6 + z * (a8 + z * a10)))); - p2 = z * (a1 + z * (a3 + z * (a5 + z * (a7 + z * (a9 + z * a11))))); - p = y * p1 + p2; - r += (p - 0.5 * y); - break; - case 1: - z = y * y; - w = z * y; - p1 = t0 + w * (t3 + w * (t6 + w * (t9 + w * t12))); /* parallel comp */ - p2 = t1 + w * (t4 + w * (t7 + w * (t10 + w * t13))); - p3 = t2 + w * (t5 + w * (t8 + w * (t11 + w * t14))); - p = z * p1 - (tt - w * (p2 + y * p3)); - r += tf + p; - break; - case 2: - p1 = y * (u0 + y * (u1 + y * (u2 + y * (u3 + y * (u4 + y * u5))))); - p2 = 1.0 + y * (v1 + y * (v2 + y * (v3 + y * (v4 + y * v5)))); - r += -0.5 * y + p1 / p2; - } - } else if (ix < 0x40200000) { /* x < 8.0 */ - i = (int)x; - y = x - (double)i; - p = y * (s0 + y * (s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6)))))); - q = 1.0 + y * (r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * r6))))); - r = 0.5 * y + p / q; - z = 1.0; /* lgamma(1+s) = log(s) + lgamma(s) */ - switch (i) { - case 7: z *= y + 6.0; /* fall through */ - case 6: z *= y + 5.0; /* fall through */ - case 5: z *= y + 4.0; /* fall through */ - case 4: z *= y + 3.0; /* fall through */ - case 3: - z *= y + 2.0; - r += log(z); - break; - } - } else if (ix < 0x43900000) { /* 8.0 <= x < 2**58 */ - t = log(x); - z = 1.0 / x; - y = z * z; - w = w0 + z * (w1 + y * (w2 + y * (w3 + y * (w4 + y * (w5 + y * w6))))); - r = (x - 0.5) * (t - 1.0) + w; - } else /* 2**58 <= x <= inf */ - r = x * (log(x) - 1.0); - if (sign) - r = nadj - r; - return r; -} - -/* sin(pi*x) assuming x > 2^-100, if sin(pi*x)==0 the sign is arbitrary */ -static float sinf_pi(float x) -{ - double y; - int n; - - /* spurious inexact if odd int */ - x = 2 * (x * 0.5f - floorf(x * 0.5f)); /* x mod 2.0 */ - - n = (int)(x * 4); - n = (n + 1) / 2; - y = x - n * 0.5f; - y *= M_PI; - switch (n) { - default: /* case 4: */ - case 0: return __sindf(y); - case 1: return __cosdf(y); - case 2: return __sindf(-y); - case 3: return -__cosdf(y); - } -} - -/********************************************************************* - * lgammaf (MSVCR120.@) - * - * Copied from musl: src/math/lgammaf_r.c - */ -float CDECL lgammaf(float x) -{ - static const float pi = 3.1415927410e+00, - a0 = 7.7215664089e-02, - a1 = 3.2246702909e-01, - a2 = 6.7352302372e-02, - a3 = 2.0580807701e-02, - a4 = 7.3855509982e-03, - a5 = 2.8905137442e-03, - a6 = 1.1927076848e-03, - a7 = 5.1006977446e-04, - a8 = 2.2086278477e-04, - a9 = 1.0801156895e-04, - a10 = 2.5214456400e-05, - a11 = 4.4864096708e-05, - tc = 1.4616321325e+00, - tf = -1.2148628384e-01, - tt = 6.6971006518e-09, - t0 = 4.8383611441e-01, - t1 = -1.4758771658e-01, - t2 = 6.4624942839e-02, - t3 = -3.2788541168e-02, - t4 = 1.7970675603e-02, - t5 = -1.0314224288e-02, - t6 = 6.1005386524e-03, - t7 = -3.6845202558e-03, - t8 = 2.2596477065e-03, - t9 = -1.4034647029e-03, - t10 = 8.8108185446e-04, - t11 = -5.3859531181e-04, - t12 = 3.1563205994e-04, - t13 = -3.1275415677e-04, - t14 = 3.3552918467e-04, - u0 = -7.7215664089e-02, - u1 = 6.3282704353e-01, - u2 = 1.4549225569e+00, - u3 = 9.7771751881e-01, - u4 = 2.2896373272e-01, - u5 = 1.3381091878e-02, - v1 = 2.4559779167e+00, - v2 = 2.1284897327e+00, - v3 = 7.6928514242e-01, - v4 = 1.0422264785e-01, - v5 = 3.2170924824e-03, - s0 = -7.7215664089e-02, - s1 = 2.1498242021e-01, - s2 = 3.2577878237e-01, - s3 = 1.4635047317e-01, - s4 = 2.6642270386e-02, - s5 = 1.8402845599e-03, - s6 = 3.1947532989e-05, - r1 = 1.3920053244e+00, - r2 = 7.2193557024e-01, - r3 = 1.7193385959e-01, - r4 = 1.8645919859e-02, - r5 = 7.7794247773e-04, - r6 = 7.3266842264e-06, - w0 = 4.1893854737e-01, - w1 = 8.3333335817e-02, - w2 = -2.7777778450e-03, - w3 = 7.9365057172e-04, - w4 = -5.9518753551e-04, - w5 = 8.3633989561e-04, - w6 = -1.6309292987e-03; - - union {float f; UINT32 i;} u = {x}; - float t, y, z, nadj, p, p1, p2, p3, q, r, w; - UINT32 ix; - int i, sign; - - /* purge off +-inf, NaN, +-0, tiny and negative arguments */ - sign = u.i >> 31; - ix = u.i & 0x7fffffff; - if (ix >= 0x7f800000) - return x * x; - if (ix < 0x35000000) { /* |x| < 2**-21, return -log(|x|) */ - if (sign) - x = -x; - return -logf(x); - } - if (sign) { - x = -x; - t = sinf_pi(x); - if (t == 0.0f) { /* -integer */ - *_errno() = ERANGE; - return 1.0f / (x - x); - } - if (t <= 0.0f) - t = -t; - nadj = logf(pi / (t * x)); - } - - /* purge off 1 and 2 */ - if (ix == 0x3f800000 || ix == 0x40000000) - r = 0; - /* for x < 2.0 */ - else if (ix < 0x40000000) { - if (ix <= 0x3f666666) { /* lgamma(x) = lgamma(x+1)-log(x) */ - r = -logf(x); - if (ix >= 0x3f3b4a20) { - y = 1.0f - x; - i = 0; - } else if (ix >= 0x3e6d3308) { - y = x - (tc - 1.0f); - i = 1; - } else { - y = x; - i = 2; - } - } else { - r = 0.0f; - if (ix >= 0x3fdda618) { /* [1.7316,2] */ - y = 2.0f - x; - i = 0; - } else if (ix >= 0x3F9da620) { /* [1.23,1.73] */ - y = x - tc; - i = 1; - } else { - y = x - 1.0f; - i = 2; - } - } - switch(i) { - case 0: - z = y * y; - p1 = a0 + z * (a2 + z * (a4 + z * (a6 + z * (a8 + z * a10)))); - p2 = z * (a1 + z * (a3 + z * (a5 + z * (a7 + z * (a9 + z * a11))))); - p = y * p1 + p2; - r += p - 0.5f * y; - break; - case 1: - z = y * y; - w = z * y; - p1 = t0 + w * (t3 + w * (t6 + w * (t9 + w * t12))); /* parallel comp */ - p2 = t1 + w * (t4 + w * (t7 + w * (t10 + w * t13))); - p3 = t2 + w * (t5 + w * (t8 + w * (t11 + w * t14))); - p = z * p1 - (tt - w * (p2 + y * p3)); - r += (tf + p); - break; - case 2: - p1 = y * (u0 + y * (u1 + y * (u2 + y * (u3 + y * (u4 + y * u5))))); - p2 = 1.0f + y * (v1 + y * (v2 + y * (v3 + y * (v4 + y * v5)))); - r += -0.5f * y + p1 / p2; - } - } else if (ix < 0x41000000) { /* x < 8.0 */ - i = (int)x; - y = x - (float)i; - p = y * (s0 + y * (s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6)))))); - q = 1.0f + y * (r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * r6))))); - r = 0.5f * y + p / q; - z = 1.0f; /* lgamma(1+s) = log(s) + lgamma(s) */ - switch (i) { - case 7: z *= y + 6.0f; /* fall through */ - case 6: z *= y + 5.0f; /* fall through */ - case 5: z *= y + 4.0f; /* fall through */ - case 4: z *= y + 3.0f; /* fall through */ - case 3: - z *= y + 2.0f; - r += logf(z); - break; - } - } else if (ix < 0x5c800000) { /* 8.0 <= x < 2**58 */ - t = logf(x); - z = 1.0f / x; - y = z * z; - w = w0 + z * (w1 + y * (w2 + y * (w3 + y * (w4 + y * (w5 + y * w6))))); - r = (x - 0.5f) * (t - 1.0f) + w; - } else /* 2**58 <= x <= inf */ - r = x * (logf(x) - 1.0f); - if (sign) - r = nadj - r; - return r; -} - static double tgamma_S(double x) { static const double Snum[] = { diff --git a/libs/musl/src/math/lgamma_r.c b/libs/musl/src/math/lgamma_r.c index f9984cd0c66..b265bdf290d 100644 --- a/libs/musl/src/math/lgamma_r.c +++ b/libs/musl/src/math/lgamma_r.c @@ -190,8 +190,10 @@ double __lgamma_r(double x, int *signgamp) if (sign) { x = -x; t = sin_pi(x); - if (t == 0.0) /* -integer */ + if (t == 0.0) { /* -integer */ + errno = ERANGE; return 1.0/(x-x); + } if (t > 0.0) *signgamp = -1; else diff --git a/libs/musl/src/math/lgammaf_r.c b/libs/musl/src/math/lgammaf_r.c index 3f353f19b33..51964ad38f7 100644 --- a/libs/musl/src/math/lgammaf_r.c +++ b/libs/musl/src/math/lgammaf_r.c @@ -125,8 +125,10 @@ float __lgammaf_r(float x, int *signgamp) if (sign) { x = -x; t = sin_pi(x); - if (t == 0.0f) /* -integer */ + if (t == 0.0f) { /* -integer */ + errno = ERANGE; return 1.0f/(x-x); + } if (t > 0.0f) *signgamp = -1; else
From: Alexandre Julliard julliard@winehq.org
--- dlls/msvcrt/math.c | 155 ------------------------------------ libs/musl/src/math/tgamma.c | 13 ++- 2 files changed, 10 insertions(+), 158 deletions(-)
diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c index 20ce8bf9599..0c0e9df41c7 100644 --- a/dlls/msvcrt/math.c +++ b/dlls/msvcrt/math.c @@ -8996,161 +8996,6 @@ float CDECL remainderf(float x, float y) return remquof(x, y, &q); }
-/* sin(pi*x) assuming x > 2^-100, if sin(pi*x)==0 the sign is arbitrary */ -static double sin_pi(double x) -{ - int n; - - /* spurious inexact if odd int */ - x = 2.0 * (x * 0.5 - floor(x * 0.5)); /* x mod 2.0 */ - - n = x * 4.0; - n = (n + 1) / 2; - x -= n * 0.5f; - x *= M_PI; - - switch (n) { - default: /* case 4: */ - case 0: return __sin(x, 0.0, 0); - case 1: return __cos(x, 0.0); - case 2: return __sin(-x, 0.0, 0); - case 3: return -__cos(x, 0.0); - } -} - -static double tgamma_S(double x) -{ - static const double Snum[] = { - 23531376880.410759688572007674451636754734846804940, - 42919803642.649098768957899047001988850926355848959, - 35711959237.355668049440185451547166705960488635843, - 17921034426.037209699919755754458931112671403265390, - 6039542586.3520280050642916443072979210699388420708, - 1439720407.3117216736632230727949123939715485786772, - 248874557.86205415651146038641322942321632125127801, - 31426415.585400194380614231628318205362874684987640, - 2876370.6289353724412254090516208496135991145378768, - 186056.26539522349504029498971604569928220784236328, - 8071.6720023658162106380029022722506138218516325024, - 210.82427775157934587250973392071336271166969580291, - 2.5066282746310002701649081771338373386264310793408, - }; - static const double Sden[] = { - 0, 39916800, 120543840, 150917976, 105258076, 45995730, 13339535, - 2637558, 357423, 32670, 1925, 66, 1, - }; - - double num = 0, den = 0; - int i; - - /* to avoid overflow handle large x differently */ - if (x < 8) - for (i = ARRAY_SIZE(Snum) - 1; i >= 0; i--) { - num = num * x + Snum[i]; - den = den * x + Sden[i]; - } - else - for (i = 0; i < ARRAY_SIZE(Snum); i++) { - num = num / x + Snum[i]; - den = den / x + Sden[i]; - } - return num / den; -} - -/********************************************************************* - * tgamma (MSVCR120.@) - * - * Copied from musl: src/math/tgamma.c - */ -double CDECL tgamma(double x) -{ - static const double gmhalf = 5.524680040776729583740234375; - static const double fact[] = { - 1, 1, 2, 6, 24, 120, 720, 5040.0, 40320.0, 362880.0, 3628800.0, 39916800.0, - 479001600.0, 6227020800.0, 87178291200.0, 1307674368000.0, 20922789888000.0, - 355687428096000.0, 6402373705728000.0, 121645100408832000.0, - 2432902008176640000.0, 51090942171709440000.0, 1124000727777607680000.0, - }; - - union {double f; UINT64 i;} u = {x}; - double absx, y, dy, z, r; - UINT32 ix = u.i >> 32 & 0x7fffffff; - int sign = u.i >> 63; - - /* special cases */ - if (ix >= 0x7ff00000) { - /* tgamma(nan)=nan, tgamma(inf)=inf, tgamma(-inf)=nan with invalid */ - if (u.i == 0xfff0000000000000ULL) - *_errno() = EDOM; - return x + INFINITY; - } - if (ix < (0x3ff - 54) << 20) { - /* |x| < 2^-54: tgamma(x) ~ 1/x, +-0 raises div-by-zero */ - if (x == 0.0) - *_errno() = ERANGE; - return 1 / x; - } - - /* integer arguments */ - /* raise inexact when non-integer */ - if (x == floor(x)) { - if (sign) { - *_errno() = EDOM; - return 0 / (x - x); - } - if (x <= ARRAY_SIZE(fact)) - return fact[(int)x - 1]; - } - - /* x >= 172: tgamma(x)=inf with overflow */ - /* x =< -184: tgamma(x)=+-0 with underflow */ - if (ix >= 0x40670000) { /* |x| >= 184 */ - *_errno() = ERANGE; - if (sign) { - fp_barrierf(0x1p-126 / x); - return 0; - } - x *= 0x1p1023; - return x; - } - - absx = sign ? -x : x; - - /* handle the error of x + g - 0.5 */ - y = absx + gmhalf; - if (absx > gmhalf) { - dy = y - absx; - dy -= gmhalf; - } else { - dy = y - gmhalf; - dy -= absx; - } - - z = absx - 0.5; - r = tgamma_S(absx) * exp(-y); - if (x < 0) { - /* reflection formula for negative x */ - /* sinpi(absx) is not 0, integers are already handled */ - r = -M_PI / (sin_pi(absx) * absx * r); - dy = -dy; - z = -z; - } - r += dy * (gmhalf + 0.5) * r / y; - z = pow(y, 0.5 * z); - y = r * z * z; - return y; -} - -/********************************************************************* - * tgammaf (MSVCR120.@) - * - * Copied from musl: src/math/tgammaf.c - */ -float CDECL tgammaf(float x) -{ - return tgamma(x); -} - /********************************************************************* * _except1 (MSVCR120.@) * TODO: diff --git a/libs/musl/src/math/tgamma.c b/libs/musl/src/math/tgamma.c index 94c8d14f000..8bea688565b 100644 --- a/libs/musl/src/math/tgamma.c +++ b/libs/musl/src/math/tgamma.c @@ -114,18 +114,24 @@ double __cdecl tgamma(double x) int sign = u.i>>63;
/* special cases */ - if (ix >= 0x7ff00000) + if (ix >= 0x7ff00000) { /* tgamma(nan)=nan, tgamma(inf)=inf, tgamma(-inf)=nan with invalid */ + if (u.i == 0xfff0000000000000ULL) errno = EDOM; return x + INFINITY; - if (ix < (0x3ff-54)<<20) + } + if (ix < (0x3ff-54)<<20) { /* |x| < 2^-54: tgamma(x) ~ 1/x, +-0 raises div-by-zero */ + if (x == 0.0) errno = ERANGE; return 1/x; + }
/* integer arguments */ /* raise inexact when non-integer */ if (x == floor(x)) { - if (sign) + if (sign) { + errno = EDOM; return 0/0.0; + } if (x <= sizeof fact/sizeof *fact) return fact[(int)x - 1]; } @@ -133,6 +139,7 @@ double __cdecl tgamma(double x) /* x >= 172: tgamma(x)=inf with overflow */ /* x =< -184: tgamma(x)=+-0 with underflow */ if (ix >= 0x40670000) { /* |x| >= 184 */ + errno = ERANGE; if (sign) { FORCE_EVAL((float)(0x1p-126/x)); if (floor(x) * 0.5 == floor(x * 0.5))
From: Alexandre Julliard julliard@winehq.org
--- dlls/msvcrt/math.c | 217 ++---------------------------------- libs/musl/src/math/expm1.c | 7 +- libs/musl/src/math/expm1f.c | 7 +- 3 files changed, 17 insertions(+), 214 deletions(-)
diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c index 0c0e9df41c7..be24a75a14f 100644 --- a/dlls/msvcrt/math.c +++ b/dlls/msvcrt/math.c @@ -489,93 +489,6 @@ static double __round(double x) }
#if !defined(__i386__) || _MSVCR_VER >= 120 -/* Copied from musl: src/math/expm1f.c */ -static float __expm1f(float x) -{ - static const float ln2_hi = 6.9313812256e-01, - ln2_lo = 9.0580006145e-06, - invln2 = 1.4426950216e+00, - Q1 = -3.3333212137e-2, - Q2 = 1.5807170421e-3; - - float y, hi, lo, c, t, e, hxs, hfx, r1, twopk; - union {float f; UINT32 i;} u = {x}; - UINT32 hx = u.i & 0x7fffffff; - int k, sign = u.i >> 31; - - /* filter out huge and non-finite argument */ - if (hx >= 0x4195b844) { /* if |x|>=27*ln2 */ - if (hx >= 0x7f800000) /* NaN */ - return u.i == 0xff800000 ? -1 : x; - if (sign) - return math_error(_UNDERFLOW, "exp", x, 0, -1); - if (hx > 0x42b17217) /* x > log(FLT_MAX) */ - return math_error(_OVERFLOW, "exp", x, 0, fp_barrierf(x * FLT_MAX)); - } - - /* argument reduction */ - if (hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */ - if (hx < 0x3F851592) { /* and |x| < 1.5 ln2 */ - if (!sign) { - hi = x - ln2_hi; - lo = ln2_lo; - k = 1; - } else { - hi = x + ln2_hi; - lo = -ln2_lo; - k = -1; - } - } else { - k = invln2 * x + (sign ? -0.5f : 0.5f); - t = k; - hi = x - t * ln2_hi; /* t*ln2_hi is exact here */ - lo = t * ln2_lo; - } - x = hi - lo; - c = (hi - x) - lo; - } else if (hx < 0x33000000) { /* when |x|<2**-25, return x */ - if (hx < 0x00800000) - fp_barrierf(x * x); - return x; - } else - k = 0; - - /* x is now in primary range */ - hfx = 0.5f * x; - hxs = x * hfx; - r1 = 1.0f + hxs * (Q1 + hxs * Q2); - t = 3.0f - r1 * hfx; - e = hxs * ((r1 - t) / (6.0f - x * t)); - if (k == 0) /* c is 0 */ - return x - (x * e - hxs); - e = x * (e - c) - c; - e -= hxs; - /* exp(x) ~ 2^k (x_reduced - e + 1) */ - if (k == -1) - return 0.5f * (x - e) - 0.5f; - if (k == 1) { - if (x < -0.25f) - return -2.0f * (e - (x + 0.5f)); - return 1.0f + 2.0f * (x - e); - } - u.i = (0x7f + k) << 23; /* 2^k */ - twopk = u.f; - if (k < 0 || k > 56) { /* suffice to return exp(x)-1 */ - y = x - e + 1.0f; - if (k == 128) - y = y * 2.0f * 0x1p127f; - else - y = y * twopk; - return y - 1.0f; - } - u.i = (0x7f-k) << 23; /* 2^-k */ - if (k < 23) - y = (x - e + (1 - u.f)) * twopk; - else - y = (x - (e + u.f) + 1) * twopk; - return y; -} - #ifndef __i386__ /* Copied from musl: src/math/__sindf.c */ static float __sindf(double x) @@ -1102,7 +1015,7 @@ float CDECL coshf( float x ) fp_barrierf(x + 0x1p120f); return 1; } - t = __expm1f(x); + t = expm1f(x); return 1 + t * t / (2 * (1 + t)); }
@@ -1641,7 +1554,7 @@ float CDECL sinhf( float x )
/* |x| < log(FLT_MAX) */ if (ui < 0x42b17217) { - t = __expm1f(absx); + t = expm1f(absx); if (ui < 0x3f800000) { if (ui < 0x3f800000 - (12 << 23)) return x; @@ -1855,16 +1768,16 @@ float CDECL tanhf( float x ) fp_barrierf(x + 0x1p120f); t = 1 + 0 / x; } else { - t = __expm1f(2 * x); + t = expm1f(2 * x); t = 1 - 2 / (t + 2); } } else if (ui > 0x3e82c578) { /* |x| > log(5/3)/2 ~= 0.2554 */ - t = __expm1f(2 * x); + t = expm1f(2 * x); t = t / (t + 2); } else if (ui >= 0x00800000) { /* |x| >= 0x1p-126 */ - t = __expm1f(-2 * x); + t = expm1f(-2 * x); t = -t / (t + 2); } else { /* |x| is subnormal */ @@ -2536,100 +2449,6 @@ double CDECL cos( double x ) } }
-/* Copied from musl: src/math/expm1.c */ -static double __expm1(double x) -{ - static const double o_threshold = 7.09782712893383973096e+02, - ln2_hi = 6.93147180369123816490e-01, - ln2_lo = 1.90821492927058770002e-10, - invln2 = 1.44269504088896338700e+00, - Q1 = -3.33333333333331316428e-02, - Q2 = 1.58730158725481460165e-03, - Q3 = -7.93650757867487942473e-05, - Q4 = 4.00821782732936239552e-06, - Q5 = -2.01099218183624371326e-07; - - double y, hi, lo, c, t, e, hxs, hfx, r1, twopk; - union {double f; UINT64 i;} u = {x}; - UINT32 hx = u.i >> 32 & 0x7fffffff; - int k, sign = u.i >> 63; - - /* filter out huge and non-finite argument */ - if (hx >= 0x4043687A) { /* if |x|>=56*ln2 */ - if (isnan(x)) - return x; - if (isinf(x)) - return sign ? -1 : x; - if (sign) - return math_error(_UNDERFLOW, "exp", x, 0, -1); - if (x > o_threshold) - return math_error(_OVERFLOW, "exp", x, 0, x * 0x1p1023); - } - - /* argument reduction */ - if (hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ - if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ - if (!sign) { - hi = x - ln2_hi; - lo = ln2_lo; - k = 1; - } else { - hi = x + ln2_hi; - lo = -ln2_lo; - k = -1; - } - } else { - k = invln2 * x + (sign ? -0.5 : 0.5); - t = k; - hi = x - t * ln2_hi; /* t*ln2_hi is exact here */ - lo = t * ln2_lo; - } - x = hi - lo; - c = (hi - x) - lo; - } else if (hx < 0x3c900000) { /* |x| < 2**-54, return x */ - fp_barrier(x + 0x1p120f); - if (hx < 0x00100000) - fp_barrier((float)x); - return x; - } else - k = 0; - - /* x is now in primary range */ - hfx = 0.5 * x; - hxs = x * hfx; - r1 = 1.0 + hxs * (Q1 + hxs * (Q2 + hxs * (Q3 + hxs * (Q4 + hxs * Q5)))); - t = 3.0 - r1 * hfx; - e = hxs * ((r1 - t) / (6.0 - x * t)); - if (k == 0) /* c is 0 */ - return x - (x * e - hxs); - e = x * (e - c) - c; - e -= hxs; - /* exp(x) ~ 2^k (x_reduced - e + 1) */ - if (k == -1) - return 0.5 * (x - e) - 0.5; - if (k == 1) { - if (x < -0.25) - return -2.0 * (e - (x + 0.5)); - return 1.0 + 2.0 * (x - e); - } - u.i = (UINT64)(0x3ff + k) << 52; /* 2^k */ - twopk = u.f; - if (k < 0 || k > 56) { /* suffice to return exp(x)-1 */ - y = x - e + 1.0; - if (k == 1024) - y = y * 2.0 * 0x1p1023; - else - y = y * twopk; - return y - 1.0; - } - u.i = (UINT64)(0x3ff - k) << 52; /* 2^-k */ - if (k < 20) - y = (x - e + (1 - u.f)) * twopk; - else - y = (x - (e + u.f) + 1) * twopk; - return y; -} - static double __expo2(double x, double sign) { static const int k = 2043; @@ -2663,7 +2482,7 @@ double CDECL cosh( double x ) fp_barrier(x + 0x1p120f); return 1; } - t = __expm1(x); + t = expm1(x); return 1 + t * t / (2 * (1 + t)); }
@@ -3932,7 +3751,7 @@ double CDECL sinh( double x )
/* |x| < log(DBL_MAX) */ if (w < 0x40862e42) { - t = __expm1(absx); + t = expm1(absx); if (w < 0x3ff00000) { if (w < 0x3ff00000 - (26 << 20)) return x; @@ -4233,16 +4052,16 @@ double CDECL tanh( double x ) fp_barrier(x + 0x1p120f); t = 1 - 0 / x; } else { - t = __expm1(2 * x); + t = expm1(2 * x); t = 1 - 2 / (t + 2); } } else if (w > 0x3fd058ae) { /* |x| > log(5/3)/2 ~= 0.2554 */ - t = __expm1(2 * x); + t = expm1(2 * x); t = t / (t + 2); } else if (w >= 0x00100000) { /* |x| >= 0x1p-1022, up to 2ulp error in [0.1,0.2554] */ - t = __expm1(-2 * x); + t = expm1(-2 * x); t = -t / (t + 2); } else { /* |x| is subnormal */ @@ -7666,22 +7485,6 @@ float CDECL exp2f(float x) return y; }
-/********************************************************************* - * expm1 (MSVCR120.@) - */ -double CDECL expm1(double x) -{ - return __expm1(x); -} - -/********************************************************************* - * expm1f (MSVCR120.@) - */ -float CDECL expm1f(float x) -{ - return __expm1f(x); -} - /********************************************************************* * log1p (MSVCR120.@) * diff --git a/libs/musl/src/math/expm1.c b/libs/musl/src/math/expm1.c index 01c303c3e7d..d29f73fac69 100644 --- a/libs/musl/src/math/expm1.c +++ b/libs/musl/src/math/expm1.c @@ -129,11 +129,12 @@ double __cdecl expm1(double x) if (hx >= 0x4043687A) { /* if |x|>=56*ln2 */ if (isnan(x)) return x; + if (isinf(x)) + return sign ? -1 : x; if (sign) - return -1; + return math_error(_UNDERFLOW, "exp", x, 0, -1); if (x > o_threshold) { - x *= 0x1p1023; - return x; + return math_error(_OVERFLOW, "exp", x, 0, x * 0x1p1023); } }
diff --git a/libs/musl/src/math/expm1f.c b/libs/musl/src/math/expm1f.c index 4fb10625c0c..3cd81903cc9 100644 --- a/libs/musl/src/math/expm1f.c +++ b/libs/musl/src/math/expm1f.c @@ -37,12 +37,11 @@ float __cdecl expm1f(float x) /* filter out huge and non-finite argument */ if (hx >= 0x4195b844) { /* if |x|>=27*ln2 */ if (hx > 0x7f800000) /* NaN */ - return x; + return u.i == 0xff800000 ? -1 : x; if (sign) - return -1; + return math_error(_UNDERFLOW, "exp", x, 0, -1); if (hx > 0x42b17217) { /* x > log(FLT_MAX) */ - x *= 0x1p127f; - return x; + return math_error(_OVERFLOW, "exp", x, 0, fp_barrierf(x * FLT_MAX)); } }
From: Alexandre Julliard julliard@winehq.org
--- dlls/msvcrt/math.c | 290 ---------------------------------------- libs/musl/src/math/j0.c | 8 +- 2 files changed, 5 insertions(+), 293 deletions(-)
diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c index be24a75a14f..13fe11ff10f 100644 --- a/dlls/msvcrt/math.c +++ b/dlls/msvcrt/math.c @@ -5510,247 +5510,6 @@ int CDECL _isnan(double num) return (u.i & ~0ull >> 1) > 0x7ffull << 52; }
-static double pzero(double x) -{ - static const double pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ - 0.00000000000000000000e+00, - -7.03124999999900357484e-02, - -8.08167041275349795626e+00, - -2.57063105679704847262e+02, - -2.48521641009428822144e+03, - -5.25304380490729545272e+03, - }, pS8[5] = { - 1.16534364619668181717e+02, - 3.83374475364121826715e+03, - 4.05978572648472545552e+04, - 1.16752972564375915681e+05, - 4.76277284146730962675e+04, - }, pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ - -1.14125464691894502584e-11, - -7.03124940873599280078e-02, - -4.15961064470587782438e+00, - -6.76747652265167261021e+01, - -3.31231299649172967747e+02, - -3.46433388365604912451e+02, - }, pS5[5] = { - 6.07539382692300335975e+01, - 1.05125230595704579173e+03, - 5.97897094333855784498e+03, - 9.62544514357774460223e+03, - 2.40605815922939109441e+03, - }, pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ - -2.54704601771951915620e-09, - -7.03119616381481654654e-02, - -2.40903221549529611423e+00, - -2.19659774734883086467e+01, - -5.80791704701737572236e+01, - -3.14479470594888503854e+01, - }, pS3[5] = { - 3.58560338055209726349e+01, - 3.61513983050303863820e+02, - 1.19360783792111533330e+03, - 1.12799679856907414432e+03, - 1.73580930813335754692e+02, - }, pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ - -8.87534333032526411254e-08, - -7.03030995483624743247e-02, - -1.45073846780952986357e+00, - -7.63569613823527770791e+00, - -1.11931668860356747786e+01, - -3.23364579351335335033e+00, - }, pS2[5] = { - 2.22202997532088808441e+01, - 1.36206794218215208048e+02, - 2.70470278658083486789e+02, - 1.53875394208320329881e+02, - 1.46576176948256193810e+01, - }; - - const double *p, *q; - double z, r, s; - UINT32 ix; - - ix = *(ULONGLONG*)&x >> 32; - ix &= 0x7fffffff; - if (ix >= 0x40200000) { - p = pR8; - q = pS8; - } else if (ix >= 0x40122E8B) { - p = pR5; - q = pS5; - } else if (ix >= 0x4006DB6D) { - p = pR3; - q = pS3; - } else /*ix >= 0x40000000*/ { - p = pR2; - q = pS2; - } - - z = 1.0 / (x * x); - r = p[0] + z * (p[1] + z * (p[2] + z * (p[3] + z * (p[4] + z * p[5])))); - s = 1.0 + z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * q[4])))); - return 1.0 + r / s; -} - -static double qzero(double x) -{ - static const double qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ - 0.00000000000000000000e+00, - 7.32421874999935051953e-02, - 1.17682064682252693899e+01, - 5.57673380256401856059e+02, - 8.85919720756468632317e+03, - 3.70146267776887834771e+04, - }, qS8[6] = { - 1.63776026895689824414e+02, - 8.09834494656449805916e+03, - 1.42538291419120476348e+05, - 8.03309257119514397345e+05, - 8.40501579819060512818e+05, - -3.43899293537866615225e+05, - }, qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ - 1.84085963594515531381e-11, - 7.32421766612684765896e-02, - 5.83563508962056953777e+00, - 1.35111577286449829671e+02, - 1.02724376596164097464e+03, - 1.98997785864605384631e+03, - }, qS5[6] = { - 8.27766102236537761883e+01, - 2.07781416421392987104e+03, - 1.88472887785718085070e+04, - 5.67511122894947329769e+04, - 3.59767538425114471465e+04, - -5.35434275601944773371e+03, - }, qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ - 4.37741014089738620906e-09, - 7.32411180042911447163e-02, - 3.34423137516170720929e+00, - 4.26218440745412650017e+01, - 1.70808091340565596283e+02, - 1.66733948696651168575e+02, - }, qS3[6] = { - 4.87588729724587182091e+01, - 7.09689221056606015736e+02, - 3.70414822620111362994e+03, - 6.46042516752568917582e+03, - 2.51633368920368957333e+03, - -1.49247451836156386662e+02, - }, qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ - 1.50444444886983272379e-07, - 7.32234265963079278272e-02, - 1.99819174093815998816e+00, - 1.44956029347885735348e+01, - 3.16662317504781540833e+01, - 1.62527075710929267416e+01, - }, qS2[6] = { - 3.03655848355219184498e+01, - 2.69348118608049844624e+02, - 8.44783757595320139444e+02, - 8.82935845112488550512e+02, - 2.12666388511798828631e+02, - -5.31095493882666946917e+00, - }; - - const double *p, *q; - double s, r, z; - unsigned int ix; - - ix = *(ULONGLONG*)&x >> 32; - ix &= 0x7fffffff; - if (ix >= 0x40200000) { - p = qR8; - q = qS8; - } else if (ix >= 0x40122E8B) { - p = qR5; - q = qS5; - } else if (ix >= 0x4006DB6D) { - p = qR3; - q = qS3; - } else /*ix >= 0x40000000*/ { - p = qR2; - q = qS2; - } - - z = 1.0 / (x * x); - r = p[0] + z * (p[1] + z * (p[2] + z * (p[3] + z * (p[4] + z * p[5])))); - s = 1.0 + z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * q[5]))))); - return (-0.125 + r / s) / x; -} - -/* j0 and y0 approximation for |x|>=2 */ -static double j0_y0_approx(unsigned int ix, double x, BOOL y0) -{ - static const double invsqrtpi = 5.64189583547756279280e-01; - - double s, c, ss, cc, z; - - s = sin(x); - c = cos(x); - if (y0) c = -c; - cc = s + c; - /* avoid overflow in 2*x, big ulp error when x>=0x1p1023 */ - if (ix < 0x7fe00000) { - ss = s - c; - z = -cos(2 * x); - if (s * c < 0) cc = z / ss; - else ss = z / cc; - if (ix < 0x48000000) { - if (y0) ss = -ss; - cc = pzero(x) * cc - qzero(x) * ss; - } - } - return invsqrtpi * cc / sqrt(x); -} - -/********************************************************************* - * _j0 (MSVCRT.@) - * - * Copied from musl: src/math/j0.c - */ -double CDECL _j0(double x) -{ - static const double R02 = 1.56249999999999947958e-02, - R03 = -1.89979294238854721751e-04, - R04 = 1.82954049532700665670e-06, - R05 = -4.61832688532103189199e-09, - S01 = 1.56191029464890010492e-02, - S02 = 1.16926784663337450260e-04, - S03 = 5.13546550207318111446e-07, - S04 = 1.16614003333790000205e-09; - - double z, r, s; - unsigned int ix; - - ix = *(ULONGLONG*)&x >> 32; - ix &= 0x7fffffff; - - /* j0(+-inf)=0, j0(nan)=nan */ - if (ix >= 0x7ff00000) - return math_error(_DOMAIN, "_j0", x, 0, 1 / (x * x)); - x = fabs(x); - - if (ix >= 0x40000000) { /* |x| >= 2 */ - /* large ulp error near zeros: 2.4, 5.52, 8.6537,.. */ - return j0_y0_approx(ix, x, FALSE); - } - - if (ix >= 0x3f200000) { /* |x| >= 2**-13 */ - /* up to 4ulp error close to 2 */ - z = x * x; - r = z * (R02 + z * (R03 + z * (R04 + z * R05))); - s = 1 + z * (S01 + z * (S02 + z * (S03 + z * S04))); - return (1 + x / 2) * (1 - x / 2) + z * (r / s); - } - - /* 1 - x*x/4 */ - /* prevent underflow */ - /* inexact should be raised when x!=0, this is not done correctly */ - if (ix >= 0x38000000) /* |x| >= 2**-127 */ - x = 0.25 * x * x; - return 1 - x; -} - static double pone(double x) { static const double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ @@ -6114,55 +5873,6 @@ double CDECL _jn(int n, double x) return sign ? -b : b; }
-/********************************************************************* - * _y0 (MSVCRT.@) - */ -double CDECL _y0(double x) -{ - static const double tpi = 6.36619772367581382433e-01, - u00 = -7.38042951086872317523e-02, - u01 = 1.76666452509181115538e-01, - u02 = -1.38185671945596898896e-02, - u03 = 3.47453432093683650238e-04, - u04 = -3.81407053724364161125e-06, - u05 = 1.95590137035022920206e-08, - u06 = -3.98205194132103398453e-11, - v01 = 1.27304834834123699328e-02, - v02 = 7.60068627350353253702e-05, - v03 = 2.59150851840457805467e-07, - v04 = 4.41110311332675467403e-10; - - double z, u, v; - unsigned int ix, lx; - - ix = *(ULONGLONG*)&x >> 32; - lx = *(ULONGLONG*)&x; - - /* y0(nan)=nan, y0(<0)=nan, y0(0)=-inf, y0(inf)=0 */ - if ((ix << 1 | lx) == 0) - return math_error(_OVERFLOW, "_y0", x, 0, -INFINITY); - if (isnan(x)) - return x; - if (ix >> 31) - return math_error(_DOMAIN, "_y0", x, 0, 0 / (x - x)); - if (ix >= 0x7ff00000) - return 1 / x; - - if (ix >= 0x40000000) { /* x >= 2 */ - /* large ulp errors near zeros: 3.958, 7.086,.. */ - return j0_y0_approx(ix, x, TRUE); - } - - if (ix >= 0x3e400000) { /* x >= 2**-27 */ - /* large ulp error near the first zero, x ~= 0.89 */ - z = x * x; - u = u00 + z * (u01 + z * (u02 + z * (u03 + z * (u04 + z * (u05 + z * u06))))); - v = 1.0 + z * (v01 + z * (v02 + z * (v03 + z * v04))); - return u / v + tpi * (j0(x) * log(x)); - } - return u00 + tpi * log(x); -} - /********************************************************************* * _y1 (MSVCRT.@) */ diff --git a/libs/musl/src/math/j0.c b/libs/musl/src/math/j0.c index 9b5e6c0c914..3462c803acf 100644 --- a/libs/musl/src/math/j0.c +++ b/libs/musl/src/math/j0.c @@ -118,7 +118,7 @@ double __cdecl _j0(double x)
/* j0(+-inf)=0, j0(nan)=nan */ if (ix >= 0x7ff00000) - return 1/(x*x); + return math_error(_DOMAIN, "_j0", x, 0, 1 / (x * x)); x = fabs(x);
if (ix >= 0x40000000) { /* |x| >= 2 */ @@ -165,9 +165,11 @@ double __cdecl _y0(double x)
/* y0(nan)=nan, y0(<0)=nan, y0(0)=-inf, y0(inf)=0 */ if ((ix<<1 | lx) == 0) - return -1/0.0; + return math_error(_OVERFLOW, "_y0", x, 0, -INFINITY); + if (isnan(x)) + return x; if (ix>>31) - return 0/0.0; + return math_error(_DOMAIN, "_y0", x, 0, 0 / (x - x)); if (ix >= 0x7ff00000) return 1/x;
From: Alexandre Julliard julliard@winehq.org
--- dlls/msvcrt/math.c | 275 ---------------------------------------- libs/musl/src/math/j1.c | 8 +- 2 files changed, 5 insertions(+), 278 deletions(-)
diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c index 13fe11ff10f..4437cbbf48e 100644 --- a/dlls/msvcrt/math.c +++ b/dlls/msvcrt/math.c @@ -5510,238 +5510,6 @@ int CDECL _isnan(double num) return (u.i & ~0ull >> 1) > 0x7ffull << 52; }
-static double pone(double x) -{ - static const double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ - 0.00000000000000000000e+00, - 1.17187499999988647970e-01, - 1.32394806593073575129e+01, - 4.12051854307378562225e+02, - 3.87474538913960532227e+03, - 7.91447954031891731574e+03, - }, ps8[5] = { - 1.14207370375678408436e+02, - 3.65093083420853463394e+03, - 3.69562060269033463555e+04, - 9.76027935934950801311e+04, - 3.08042720627888811578e+04, - }, pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ - 1.31990519556243522749e-11, - 1.17187493190614097638e-01, - 6.80275127868432871736e+00, - 1.08308182990189109773e+02, - 5.17636139533199752805e+02, - 5.28715201363337541807e+02, - }, ps5[5] = { - 5.92805987221131331921e+01, - 9.91401418733614377743e+02, - 5.35326695291487976647e+03, - 7.84469031749551231769e+03, - 1.50404688810361062679e+03, - }, pr3[6] = { - 3.02503916137373618024e-09, - 1.17186865567253592491e-01, - 3.93297750033315640650e+00, - 3.51194035591636932736e+01, - 9.10550110750781271918e+01, - 4.85590685197364919645e+01, - }, ps3[5] = { - 3.47913095001251519989e+01, - 3.36762458747825746741e+02, - 1.04687139975775130551e+03, - 8.90811346398256432622e+02, - 1.03787932439639277504e+02, - }, pr2[6] = { /* for x in [2.8570,2]=1/[0.3499,0.5] */ - 1.07710830106873743082e-07, - 1.17176219462683348094e-01, - 2.36851496667608785174e+00, - 1.22426109148261232917e+01, - 1.76939711271687727390e+01, - 5.07352312588818499250e+00, - }, ps2[5] = { - 2.14364859363821409488e+01, - 1.25290227168402751090e+02, - 2.32276469057162813669e+02, - 1.17679373287147100768e+02, - 8.36463893371618283368e+00, - }; - - const double *p, *q; - double z, r, s; - unsigned int ix; - - ix = *(ULONGLONG*)&x >> 32; - ix &= 0x7fffffff; - if (ix >= 0x40200000) { - p = pr8; - q = ps8; - } else if (ix >= 0x40122E8B) { - p = pr5; - q = ps5; - } else if (ix >= 0x4006DB6D) { - p = pr3; - q = ps3; - } else /*ix >= 0x40000000*/ { - p = pr2; - q = ps2; - } - z = 1.0 / (x * x); - r = p[0] + z * (p[1] + z * (p[2] + z * (p[3] + z * (p[4] + z * p[5])))); - s = 1.0 + z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * q[4])))); - return 1.0 + r / s; -} - -static double qone(double x) -{ - static const double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ - 0.00000000000000000000e+00, - -1.02539062499992714161e-01, - -1.62717534544589987888e+01, - -7.59601722513950107896e+02, - -1.18498066702429587167e+04, - -4.84385124285750353010e+04, - }, qs8[6] = { - 1.61395369700722909556e+02, - 7.82538599923348465381e+03, - 1.33875336287249578163e+05, - 7.19657723683240939863e+05, - 6.66601232617776375264e+05, - -2.94490264303834643215e+05, - }, qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ - -2.08979931141764104297e-11, - -1.02539050241375426231e-01, - -8.05644828123936029840e+00, - -1.83669607474888380239e+02, - -1.37319376065508163265e+03, - -2.61244440453215656817e+03, - }, qs5[6] = { - 8.12765501384335777857e+01, - 1.99179873460485964642e+03, - 1.74684851924908907677e+04, - 4.98514270910352279316e+04, - 2.79480751638918118260e+04, - -4.71918354795128470869e+03, - }, qr3[6] = { - -5.07831226461766561369e-09, - -1.02537829820837089745e-01, - -4.61011581139473403113e+00, - -5.78472216562783643212e+01, - -2.28244540737631695038e+02, - -2.19210128478909325622e+02, - }, qs3[6] = { - 4.76651550323729509273e+01, - 6.73865112676699709482e+02, - 3.38015286679526343505e+03, - 5.54772909720722782367e+03, - 1.90311919338810798763e+03, - -1.35201191444307340817e+02, - }, qr2[6] = { /* for x in [2.8570,2]=1/[0.3499,0.5] */ - -1.78381727510958865572e-07, - -1.02517042607985553460e-01, - -2.75220568278187460720e+00, - -1.96636162643703720221e+01, - -4.23253133372830490089e+01, - -2.13719211703704061733e+01, - }, qs2[6] = { - 2.95333629060523854548e+01, - 2.52981549982190529136e+02, - 7.57502834868645436472e+02, - 7.39393205320467245656e+02, - 1.55949003336666123687e+02, - -4.95949898822628210127e+00, - }; - - const double *p, *q; - double s, r, z; - unsigned int ix; - - ix = *(ULONGLONG*)&x >> 32; - ix &= 0x7fffffff; - if (ix >= 0x40200000) { - p = qr8; - q = qs8; - } else if (ix >= 0x40122E8B) { - p = qr5; - q = qs5; - } else if (ix >= 0x4006DB6D) { - p = qr3; - q = qs3; - } else /*ix >= 0x40000000*/ { - p = qr2; - q = qs2; - } - z = 1.0 / (x * x); - r = p[0] + z * (p[1] + z * (p[2] + z * (p[3] + z * (p[4] + z * p[5])))); - s = 1.0 + z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * q[5]))))); - return (0.375 + r / s) / x; -} - -static double j1_y1_approx(unsigned int ix, double x, BOOL y1, int sign) -{ - static const double invsqrtpi = 5.64189583547756279280e-01; - - double z, s, c, ss, cc; - - s = sin(x); - if (y1) s = -s; - c = cos(x); - cc = s - c; - if (ix < 0x7fe00000) { - ss = -s - c; - z = cos(2 * x); - if (s * c > 0) cc = z / ss; - else ss = z / cc; - if (ix < 0x48000000) { - if (y1) - ss = -ss; - cc = pone(x) * cc - qone(x) * ss; - } - } - if (sign) - cc = -cc; - return invsqrtpi * cc / sqrt(x); -} - -/********************************************************************* - * _j1 (MSVCRT.@) - * - * Copied from musl: src/math/j1.c - */ -double CDECL _j1(double x) -{ - static const double r00 = -6.25000000000000000000e-02, - r01 = 1.40705666955189706048e-03, - r02 = -1.59955631084035597520e-05, - r03 = 4.96727999609584448412e-08, - s01 = 1.91537599538363460805e-02, - s02 = 1.85946785588630915560e-04, - s03 = 1.17718464042623683263e-06, - s04 = 5.04636257076217042715e-09, - s05 = 1.23542274426137913908e-11; - - double z, r, s; - unsigned int ix; - int sign; - - ix = *(ULONGLONG*)&x >> 32; - sign = ix >> 31; - ix &= 0x7fffffff; - if (ix >= 0x7ff00000) - return math_error(isnan(x) ? 0 : _DOMAIN, "_j1", x, 0, 1 / (x * x)); - if (ix >= 0x40000000) /* |x| >= 2 */ - return j1_y1_approx(ix, fabs(x), FALSE, sign); - if (ix >= 0x38000000) { /* |x| >= 2**-127 */ - z = x * x; - r = z * (r00 + z * (r01 + z * (r02 + z * r03))); - s = 1 + z * (s01 + z * (s02 + z * (s03 + z * (s04 + z * s05)))); - z = r / s; - } else { - /* avoid underflow, raise inexact if x!=0 */ - z = x; - } - return (0.5 + z) * x; -} - /********************************************************************* * _jn (MSVCRT.@) * @@ -5873,49 +5641,6 @@ double CDECL _jn(int n, double x) return sign ? -b : b; }
-/********************************************************************* - * _y1 (MSVCRT.@) - */ -double CDECL _y1(double x) -{ - static const double tpi = 6.36619772367581382433e-01, - u00 = -1.96057090646238940668e-01, - u01 = 5.04438716639811282616e-02, - u02 = -1.91256895875763547298e-03, - u03 = 2.35252600561610495928e-05, - u04 = -9.19099158039878874504e-08, - v00 = 1.99167318236649903973e-02, - v01 = 2.02552581025135171496e-04, - v02 = 1.35608801097516229404e-06, - v03 = 6.22741452364621501295e-09, - v04 = 1.66559246207992079114e-11; - - double z, u, v; - unsigned int ix, lx; - - ix = *(ULONGLONG*)&x >> 32; - lx = *(ULONGLONG*)&x; - - /* y1(nan)=nan, y1(<0)=nan, y1(0)=-inf, y1(inf)=0 */ - if ((ix << 1 | lx) == 0) - return math_error(_OVERFLOW, "_y1", x, 0, -INFINITY); - if (isnan(x)) - return x; - if (ix >> 31) - return math_error(_DOMAIN, "_y1", x, 0, 0 / (x - x)); - if (ix >= 0x7ff00000) - return 1 / x; - - if (ix >= 0x40000000) /* x >= 2 */ - return j1_y1_approx(ix, x, TRUE, 0); - if (ix < 0x3c900000) /* x < 2**-54 */ - return -tpi / x; - z = x * x; - u = u00 + z * (u01 + z * (u02 + z * (u03 + z * u04))); - v = 1 + z * (v00 + z * (v01 + z * (v02 + z * (v03 + z * v04)))); - return x * (u / v) + tpi * (j1(x) * log(x) - 1 / x); -} - /********************************************************************* * _yn (MSVCRT.@) * diff --git a/libs/musl/src/math/j1.c b/libs/musl/src/math/j1.c index 594711b88bb..178c1912e26 100644 --- a/libs/musl/src/math/j1.c +++ b/libs/musl/src/math/j1.c @@ -120,7 +120,7 @@ double __cdecl _j1(double x) sign = ix>>31; ix &= 0x7fffffff; if (ix >= 0x7ff00000) - return 1/(x*x); + return math_error(isnan(x) ? 0 : _DOMAIN, "_j1", x, 0, 1 / (x * x)); if (ix >= 0x40000000) /* |x| >= 2 */ return common(ix, fabs(x), 0, sign); if (ix >= 0x38000000) { /* |x| >= 2**-127 */ @@ -157,9 +157,11 @@ double __cdecl _y1(double x) EXTRACT_WORDS(ix, lx, x); /* y1(nan)=nan, y1(<0)=nan, y1(0)=-inf, y1(inf)=0 */ if ((ix<<1 | lx) == 0) - return -1/0.0; + return math_error(_OVERFLOW, "_y1", x, 0, -INFINITY); + if (isnan(x)) + return x; if (ix>>31) - return 0/0.0; + return math_error(_DOMAIN, "_y1", x, 0, 0 / (x - x)); if (ix >= 0x7ff00000) return 1/x;
From: Alexandre Julliard julliard@winehq.org
--- dlls/msvcrt/math.c | 200 ---------------------------------------- libs/musl/src/math/jn.c | 2 +- 2 files changed, 1 insertion(+), 201 deletions(-)
diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c index 4437cbbf48e..5ff75274dbc 100644 --- a/dlls/msvcrt/math.c +++ b/dlls/msvcrt/math.c @@ -5510,206 +5510,6 @@ int CDECL _isnan(double num) return (u.i & ~0ull >> 1) > 0x7ffull << 52; }
-/********************************************************************* - * _jn (MSVCRT.@) - * - * Copied from musl: src/math/jn.c - */ -double CDECL _jn(int n, double x) -{ - static const double invsqrtpi = 5.64189583547756279280e-01; - - unsigned int ix, lx; - int nm1, i, sign; - double a, b, temp; - - ix = *(ULONGLONG*)&x >> 32; - lx = *(ULONGLONG*)&x; - sign = ix >> 31; - ix &= 0x7fffffff; - - if ((ix | (lx | -lx) >> 31) > 0x7ff00000) /* nan */ - return x; - - if (n == 0) - return _j0(x); - if (n < 0) { - nm1 = -(n + 1); - x = -x; - sign ^= 1; - } else { - nm1 = n-1; - } - if (nm1 == 0) - return j1(x); - - sign &= n; /* even n: 0, odd n: signbit(x) */ - x = fabs(x); - if ((ix | lx) == 0 || ix == 0x7ff00000) /* if x is 0 or inf */ - b = 0.0; - else if (nm1 < x) { - if (ix >= 0x52d00000) { /* x > 2**302 */ - switch(nm1 & 3) { - case 0: - temp = -cos(x) + sin(x); - break; - case 1: - temp = -cos(x) - sin(x); - break; - case 2: - temp = cos(x) - sin(x); - break; - default: - temp = cos(x) + sin(x); - break; - } - b = invsqrtpi * temp / sqrt(x); - } else { - a = _j0(x); - b = _j1(x); - for (i = 0; i < nm1; ) { - i++; - temp = b; - b = b * (2.0 * i / x) - a; /* avoid underflow */ - a = temp; - } - } - } else { - if (ix < 0x3e100000) { /* x < 2**-29 */ - if (nm1 > 32) /* underflow */ - b = 0.0; - else { - temp = x * 0.5; - b = temp; - a = 1.0; - for (i = 2; i <= nm1 + 1; i++) { - a *= (double)i; /* a = n! */ - b *= temp; /* b = (x/2)^n */ - } - b = b / a; - } - } else { - double t, q0, q1, w, h, z, tmp, nf; - int k; - - nf = nm1 + 1.0; - w = 2 * nf / x; - h = 2 / x; - z = w + h; - q0 = w; - q1 = w * z - 1.0; - k = 1; - while (q1 < 1.0e9) { - k += 1; - z += h; - tmp = z * q1 - q0; - q0 = q1; - q1 = tmp; - } - for (t = 0.0, i = k; i >= 0; i--) - t = 1 / (2 * (i + nf) / x - t); - a = t; - b = 1.0; - tmp = nf * log(fabs(w)); - if (tmp < 7.09782712893383973096e+02) { - for (i = nm1; i > 0; i--) { - temp = b; - b = b * (2.0 * i) / x - a; - a = temp; - } - } else { - for (i = nm1; i > 0; i--) { - temp = b; - b = b * (2.0 * i) / x - a; - a = temp; - /* scale b to avoid spurious overflow */ - if (b > 0x1p500) { - a /= b; - t /= b; - b = 1.0; - } - } - } - z = j0(x); - w = j1(x); - if (fabs(z) >= fabs(w)) - b = t * z / b; - else - b = t * w / a; - } - } - return sign ? -b : b; -} - -/********************************************************************* - * _yn (MSVCRT.@) - * - * Copied from musl: src/math/jn.c - */ -double CDECL _yn(int n, double x) -{ - static const double invsqrtpi = 5.64189583547756279280e-01; - - unsigned int ix, lx, ib; - int nm1, sign, i; - double a, b, temp; - - ix = *(ULONGLONG*)&x >> 32; - lx = *(ULONGLONG*)&x; - sign = ix >> 31; - ix &= 0x7fffffff; - - if ((ix | (lx | -lx) >> 31) > 0x7ff00000) /* nan */ - return x; - if (sign && (ix | lx) != 0) /* x < 0 */ - return math_error(_DOMAIN, "_y1", x, 0, 0 / (x - x)); - if (ix == 0x7ff00000) - return 0.0; - - if (n == 0) - return y0(x); - if (n < 0) { - nm1 = -(n + 1); - sign = n & 1; - } else { - nm1 = n - 1; - sign = 0; - } - if (nm1 == 0) - return sign ? -y1(x) : y1(x); - - if (ix >= 0x52d00000) { /* x > 2**302 */ - switch(nm1 & 3) { - case 0: - temp = -sin(x) - cos(x); - break; - case 1: - temp = -sin(x) + cos(x); - break; - case 2: - temp = sin(x) + cos(x); - break; - default: - temp = sin(x) - cos(x); - break; - } - b = invsqrtpi * temp / sqrt(x); - } else { - a = y0(x); - b = y1(x); - /* quit if b is -inf */ - ib = *(ULONGLONG*)&b >> 32; - for (i = 0; i < nm1 && ib != 0xfff00000;) { - i++; - temp = b; - b = (2.0 * i / x) * b - a; - ib = *(ULONGLONG*)&b >> 32; - a = temp; - } - } - return sign ? -b : b; -} - #if _MSVCR_VER>=120
/********************************************************************* diff --git a/libs/musl/src/math/jn.c b/libs/musl/src/math/jn.c index 55c05157be0..f7e679e25ea 100644 --- a/libs/musl/src/math/jn.c +++ b/libs/musl/src/math/jn.c @@ -225,7 +225,7 @@ double __cdecl _yn(int n, double x) if ((ix | (lx|-lx)>>31) > 0x7ff00000) /* nan */ return x; if (sign && (ix|lx)!=0) /* x < 0 */ - return 0/0.0; + return math_error(_DOMAIN, "_y1", x, 0, 0 / (x - x)); if (ix == 0x7ff00000) return 0.0;
From: Alexandre Julliard julliard@winehq.org
--- dlls/msvcrt/math.c | 136 ------------------------------------- libs/musl/src/math/fmod.c | 3 + libs/musl/src/math/fmodf.c | 2 + 3 files changed, 5 insertions(+), 136 deletions(-)
diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c index 5ff75274dbc..1b3d380bc53 100644 --- a/dlls/msvcrt/math.c +++ b/dlls/msvcrt/math.c @@ -1084,74 +1084,6 @@ float CDECL expf( float x ) return y; }
-/********************************************************************* - * fmodf (MSVCRT.@) - * - * Copied from musl: src/math/fmodf.c - */ -float CDECL fmodf( float x, float y ) -{ - UINT32 xi = *(UINT32*)&x; - UINT32 yi = *(UINT32*)&y; - int ex = xi>>23 & 0xff; - int ey = yi>>23 & 0xff; - UINT32 sx = xi & 0x80000000; - UINT32 i; - - if (isinf(x)) return math_error(_DOMAIN, "fmodf", x, y, (x * y) / (x * y)); - if (yi << 1 == 0 || isnan(y) || ex == 0xff) - return (x * y) / (x * y); - if (xi << 1 <= yi << 1) { - if (xi << 1 == yi << 1) - return 0 * x; - return x; - } - - /* normalize x and y */ - if (!ex) { - for (i = xi << 9; i >> 31 == 0; ex--, i <<= 1); - xi <<= -ex + 1; - } else { - xi &= -1U >> 9; - xi |= 1U << 23; - } - if (!ey) { - for (i = yi << 9; i >> 31 == 0; ey--, i <<= 1); - yi <<= -ey + 1; - } else { - yi &= -1U >> 9; - yi |= 1U << 23; - } - - /* x mod y */ - for (; ex > ey; ex--) { - i = xi - yi; - if (i >> 31 == 0) { - if (i == 0) - return 0 * x; - xi = i; - } - xi <<= 1; - } - i = xi - yi; - if (i >> 31 == 0) { - if (i == 0) - return 0 * x; - xi = i; - } - for (; xi>>23 == 0; xi <<= 1, ex--); - - /* scale result up */ - if (ex > 0) { - xi -= 1U << 23; - xi |= (UINT32)ex << 23; - } else { - xi >>= -ex + 1; - } - xi |= sx; - return *(float*)ξ -} - /********************************************************************* * logf (MSVCRT.@) * @@ -2748,74 +2680,6 @@ double CDECL exp( double x ) return scale + scale * tmp; }
-/********************************************************************* - * fmod (MSVCRT.@) - * - * Copied from musl: src/math/fmod.c - */ -double CDECL fmod( double x, double y ) -{ - UINT64 xi = *(UINT64*)&x; - UINT64 yi = *(UINT64*)&y; - int ex = xi >> 52 & 0x7ff; - int ey = yi >> 52 & 0x7ff; - int sx = xi >> 63; - UINT64 i; - - if (isinf(x)) return math_error(_DOMAIN, "fmod", x, y, (x * y) / (x * y)); - if (yi << 1 == 0 || isnan(y) || ex == 0x7ff) - return (x * y) / (x * y); - if (xi << 1 <= yi << 1) { - if (xi << 1 == yi << 1) - return 0 * x; - return x; - } - - /* normalize x and y */ - if (!ex) { - for (i = xi << 12; i >> 63 == 0; ex--, i <<= 1); - xi <<= -ex + 1; - } else { - xi &= -1ULL >> 12; - xi |= 1ULL << 52; - } - if (!ey) { - for (i = yi << 12; i >> 63 == 0; ey--, i <<= 1); - yi <<= -ey + 1; - } else { - yi &= -1ULL >> 12; - yi |= 1ULL << 52; - } - - /* x mod y */ - for (; ex > ey; ex--) { - i = xi - yi; - if (i >> 63 == 0) { - if (i == 0) - return 0 * x; - xi = i; - } - xi <<= 1; - } - i = xi - yi; - if (i >> 63 == 0) { - if (i == 0) - return 0 * x; - xi = i; - } - for (; xi >> 52 == 0; xi <<= 1, ex--); - - /* scale result */ - if (ex > 0) { - xi -= 1ULL << 52; - xi |= (UINT64)ex << 52; - } else { - xi >>= -ex + 1; - } - xi |= (UINT64)sx << 63; - return *(double*)ξ -} - /********************************************************************* * log (MSVCRT.@) * diff --git a/libs/musl/src/math/fmod.c b/libs/musl/src/math/fmod.c index 89d882e2554..484f65a6b5f 100644 --- a/libs/musl/src/math/fmod.c +++ b/libs/musl/src/math/fmod.c @@ -1,5 +1,6 @@ #include <math.h> #include <stdint.h> +#include "libm.h"
double __cdecl fmod(double x, double y) { @@ -13,6 +14,8 @@ double __cdecl fmod(double x, double y) /* float load/store to inner loops ruining performance and code size */ uint64_t uxi = ux.i;
+ if (isinf(x)) + return math_error(_DOMAIN, "fmod", x, y, (x * y) / (x * y)); if (uy.i<<1 == 0 || isnan(y) || ex == 0x7ff) return (x*y)/(x*y); if (uxi<<1 <= uy.i<<1) { diff --git a/libs/musl/src/math/fmodf.c b/libs/musl/src/math/fmodf.c index 2b79521902e..12431b06da0 100644 --- a/libs/musl/src/math/fmodf.c +++ b/libs/musl/src/math/fmodf.c @@ -11,6 +11,8 @@ float __cdecl fmodf(float x, float y) uint32_t i; uint32_t uxi = ux.i;
+ if (isinf(x)) + return math_error(_DOMAIN, "fmodf", x, y, (x * y) / (x * y)); if (uy.i<<1 == 0 || isnan(y) || ex == 0xff) return (x*y)/(x*y); if (uxi<<1 <= uy.i<<1) {
This merge request was approved by Piotr Caban.
-
Alexandre Julliard -
Alexandre Julliard (＠julliard)
-
Piotr Caban (＠piotr)