[PATCH] msvcrt: Add fallback implementations for Bessel functions

[Alex Henrie]

Signed-off-by: Alex Henrie <alexhenrie24 at 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