[PATCH] msvcrt: Add fallback implementations for Bessel functions
Signed-off-by: Alex Henrie <alexhenrie24(a)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) { -- 2.16.2
Hi Alex, On 03/02/18 05:13, Alex Henrie wrote:
/********************************************************************* @@ -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; } This approximation can only be used for values in specific range to give reasonable results. There are similar problems with j0, j1 fallbacks (I haven't checked other functions).
I think it's better to print FIXME instead of adding very inaccurate fallback. Taking in account the functions are probably mainly used for physics related computations I guess it's better to avoid inaccurate implementation at all. Thanks, Piotr
2018-03-02 7:31 GMT-07:00 Piotr Caban <piotr.caban(a)gmail.com>:
This approximation can only be used for values in specific range to give reasonable results. There are similar problems with j0, j1 fallbacks (I haven't checked other functions).
I think it's better to print FIXME instead of adding very inaccurate fallback. Taking in account the functions are probably mainly used for physics related computations I guess it's better to avoid inaccurate implementation at all.
That's fair for jn, y0, y1, and yn. However, the j0 and j1 fallbacks are exactly correct, as far as I can tell. What values of x did you test? Did you accidentally have your calculator set to use degrees instead of radians? -Alex
On 03/02/18 16:29, Alex Henrie wrote:
That's fair for jn, y0, y1, and yn. However, the j0 and j1 fallbacks are exactly correct, as far as I can tell. What values of x did you test? Did you accidentally have your calculator set to use degrees instead of radians? The j0 you have implemented is the value of spherical Bessel function not Bessel function.
Thanks, Piotr
2018-03-02 9:07 GMT-07:00 Piotr Caban <piotr.caban(a)gmail.com>:
On 03/02/18 16:29, Alex Henrie wrote:
That's fair for jn, y0, y1, and yn. However, the j0 and j1 fallbacks are exactly correct, as far as I can tell. What values of x did you test? Did you accidentally have your calculator set to use degrees instead of radians?
The j0 you have implemented is the value of spherical Bessel function not Bessel function.
Thanks for clarifying that; it was pretty confusing. -Alex
participants (2)
-
Alex Henrie -
Piotr Caban