Signed-off-by: Alex Henrie alexhenrie24@gmail.com --- Fixes https://bugs.winehq.org/show_bug.cgi?id=37809
The formulas for j0 and j1 are from Wolfram Alpha. The formulas for jn, y0, y1, and yn are from "Fast and Accurate Bessel Function Computation" by John Harrison http://www.cl.cam.ac.uk/~jrh13/papers/bessel.pdf. Recognizing that jn is always between -1 and 1 but this simple formula only roughly approximates that, my implementation clamps the value of jn to [-1, 1].
configure.ac | 8 +++++++- dlls/msvcrt/math.c | 32 +++++++++++++++++++++++++++++++- include/wine/port.h | 4 ++++ 3 files changed, 42 insertions(+), 2 deletions(-)
diff --git a/configure.ac b/configure.ac index 2e99320371..e198018427 100644 --- a/configure.ac +++ b/configure.ac @@ -2698,6 +2698,9 @@ AC_CHECK_FUNCS(\ exp2f \ expm1 \ expm1f \ + j0 \ + j1 \ + jn \ lgamma \ lgammaf \ llrint \ @@ -2722,7 +2725,10 @@ AC_CHECK_FUNCS(\ round \ roundf \ trunc \ - truncf + truncf \ + y0 \ + y1 \ + yn ) LIBS="$ac_save_LIBS"
diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c index 9b480117fc..d1c38b9649 100644 --- a/dlls/msvcrt/math.c +++ b/dlls/msvcrt/math.c @@ -1412,7 +1412,11 @@ INT CDECL MSVCRT__isnan(double num) double CDECL MSVCRT__j0(double num) { /* FIXME: errno handling */ +#ifdef HAVE_J0 return j0(num); +#else + return sin(num) / num; +#endif }
/********************************************************************* @@ -1421,7 +1425,11 @@ double CDECL MSVCRT__j0(double num) double CDECL MSVCRT__j1(double num) { /* FIXME: errno handling */ +#ifdef HAVE_J1 return j1(num); +#else + return sin(num) / (num * num) - cos(num) / num; +#endif }
/********************************************************************* @@ -1429,8 +1437,18 @@ double CDECL MSVCRT__j1(double num) */ double CDECL MSVCRT__jn(int n, double num) { + double retval; /* FIXME: errno handling */ - return jn(n, num); +#ifdef HAVE_JN + retval = jn(n, num); +#else + if (n == 0) return MSVCRT__j0(num); + if (n == 1) return MSVCRT__j1(num); + retval = sqrt(2 / (M_PI * num)) * cos(num - (n / 2.0 + 0.25) * M_PI); + if (retval < -1) return -1; + if (retval > 1) return 1; +#endif + return retval; }
/********************************************************************* @@ -1440,7 +1458,11 @@ double CDECL MSVCRT__y0(double num) { double retval; if (!isfinite(num)) *MSVCRT__errno() = MSVCRT_EDOM; +#ifdef HAVE_Y0 retval = y0(num); +#else + retval = sqrt(2 / (M_PI * num)) * sin(num - M_PI_4); +#endif if (MSVCRT__fpclass(retval) == MSVCRT__FPCLASS_NINF) { *MSVCRT__errno() = MSVCRT_EDOM; @@ -1456,7 +1478,11 @@ double CDECL MSVCRT__y1(double num) { double retval; if (!isfinite(num)) *MSVCRT__errno() = MSVCRT_EDOM; +#ifdef HAVE_Y1 retval = y1(num); +#else + retval = sqrt(2 / (M_PI * num)) * sin(num - 0.75 * M_PI); +#endif if (MSVCRT__fpclass(retval) == MSVCRT__FPCLASS_NINF) { *MSVCRT__errno() = MSVCRT_EDOM; @@ -1472,7 +1498,11 @@ double CDECL MSVCRT__yn(int order, double num) { double retval; if (!isfinite(num)) *MSVCRT__errno() = MSVCRT_EDOM; +#ifdef HAVE_YN retval = yn(order,num); +#else + retval = sqrt(2 / (M_PI * num)) * sin(num - (order / 2.0 + 0.25) * M_PI); +#endif if (MSVCRT__fpclass(retval) == MSVCRT__FPCLASS_NINF) { *MSVCRT__errno() = MSVCRT_EDOM; diff --git a/include/wine/port.h b/include/wine/port.h index fb7251edd6..07366ea498 100644 --- a/include/wine/port.h +++ b/include/wine/port.h @@ -202,6 +202,10 @@ struct statvfs #define M_PI_2 1.570796326794896619 #endif
+#ifndef M_PI_4 +#define M_PI_4 0.785398163397448310 +#endif + #ifndef INFINITY static inline float __port_infinity(void) {