msvcrt: Import j1 implementation from musl.

configure

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

configure.ac

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

dlls/msvcrt/math.c

@@ -139,6 +139,9 @@ static double math_error(int type, const char *name, double arg1, double arg2, d
 
     switch (type)
     {
+    case 0:
+        /* don't set errno */
+        break;
     case _DOMAIN:
         *_errno() = EDOM;
         break;
@@ -2675,13 +2678,236 @@ double CDECL _j0(double x)
     return 1 - x;
 }
 
+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 num)
+double CDECL _j1(double x)
 {
-  /* FIXME: errno handling */
-  return unix_funcs->j1( num );
+    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;
 }
 
 /*********************************************************************

dlls/msvcrt/unixlib.c

@@ -402,19 +402,6 @@ static float CDECL unix_hypotf(float x, float y)
 }
 
 /*********************************************************************
- *      j1
- */
-static double CDECL unix_j1(double num)
-{
-#ifdef HAVE_J1
-    return j1(num);
-#else
-    FIXME("not implemented\n");
-    return 0;
-#endif
-}
-
-/*********************************************************************
  *      jn
  */
 static double CDECL unix_jn(int n, double num)
@@ -1027,7 +1014,6 @@ static const struct unix_funcs funcs =
     unix_frexpf,
     unix_hypot,
     unix_hypotf,
-    unix_j1,
     unix_jn,
     unix_ldexp,
     unix_lgamma,

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 *j1)(double num);
     double          (CDECL *jn)(int n, double num);
     double          (CDECL *ldexp)(double x, int exp);
     double          (CDECL *lgamma)(double 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 `j1' function. */
-#undef HAVE_J1
-
 /* Define to 1 if you have the `jn' function. */
 #undef HAVE_JN