msvcrt: Import jn implementation from musl.

Signed-off-by: Piotr Caban <piotr@codeweavers.com>
Signed-off-by: Alexandre Julliard <julliard@winehq.org>

configure

@@ -19631,7 +19631,6 @@ for ac_func in \
 	expm1f \
 	fma \
 	fmaf \
-	jn \
 	lgamma \
 	lgammaf \
 	llrint \

configure.ac

@@ -2674,7 +2674,6 @@ AC_CHECK_FUNCS(\
 	expm1f \
 	fma \
 	fmaf \
-	jn \
 	lgamma \
 	lgammaf \
 	llrint \

dlls/msvcrt/math.c

@@ -2912,11 +2912,133 @@ double CDECL _j1(double x)
 
 /*********************************************************************
  *		_jn (MSVCRT.@)
+ *
+ * Copied from musl: src/math/jn.c
  */
-double CDECL _jn(int n, double num)
+double CDECL _jn(int n, double x)
 {
-  /* FIXME: errno handling */
+    static const double invsqrtpi = 5.64189583547756279280e-01;
-  return unix_funcs->jn( n, num );
+
+    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;
 }
 
 /*********************************************************************

dlls/msvcrt/unixlib.c

@@ -401,19 +401,6 @@ static float CDECL unix_hypotf(float x, float y)
     return hypotf( x, y );
 }
 
-/*********************************************************************
- *      jn
- */
-static double CDECL unix_jn(int n, double num)
-{
-#ifdef HAVE_JN
-    return jn(n, num);
-#else
-    FIXME("not implemented\n");
-    return 0;
-#endif
-}
-
 /*********************************************************************
  *      ldexp
  */
@@ -1001,7 +988,6 @@ static const struct unix_funcs funcs =
     unix_frexpf,
     unix_hypot,
     unix_hypotf,
-    unix_jn,
     unix_ldexp,
     unix_lgamma,
     unix_lgammaf,

dlls/msvcrt/unixlib.h

@@ -57,7 +57,6 @@ struct unix_funcs
     float           (CDECL *frexpf)(float x, int *exp);
     double          (CDECL *hypot)(double x, double y);
     float           (CDECL *hypotf)(float x, float y);
-    double          (CDECL *jn)(int n, double num);
     double          (CDECL *ldexp)(double x, int exp);
     double          (CDECL *lgamma)(double x);
     float           (CDECL *lgammaf)(float x);

include/config.h.in

@@ -288,9 +288,6 @@
 /* Define to 1 if you have the `isnan' function. */
 #undef HAVE_ISNAN
 
-/* Define to 1 if you have the `jn' function. */
-#undef HAVE_JN
-
 /* Define to 1 if you have the <jpeglib.h> header file. */
 #undef HAVE_JPEGLIB_H