From 3911ac3f45b8b1372f8677c8ee69c16d9848cd02 Mon Sep 17 00:00:00 2001 From: Piotr Caban Date: Thu, 29 Apr 2021 17:06:15 +0200 Subject: [PATCH] msvcrt: Import j1 implementation from musl. Signed-off-by: Piotr Caban Signed-off-by: Alexandre Julliard --- configure | 1 - configure.ac | 1 - dlls/msvcrt/math.c | 232 +++++++++++++++++++++++++++++++++++++++++- dlls/msvcrt/unixlib.c | 14 --- dlls/msvcrt/unixlib.h | 1 - include/config.h.in | 3 - 6 files changed, 229 insertions(+), 23 deletions(-) diff --git a/configure b/configure index a3d73d0551e..cdc0342b066 100755 --- a/configure +++ b/configure @@ -19631,7 +19631,6 @@ for ac_func in \ expm1f \ fma \ fmaf \ - j1 \ jn \ lgamma \ lgammaf \ diff --git a/configure.ac b/configure.ac index 678d3c8344a..03915d384f2 100644 --- a/configure.ac +++ b/configure.ac @@ -2674,7 +2674,6 @@ AC_CHECK_FUNCS(\ expm1f \ fma \ fmaf \ - j1 \ jn \ lgamma \ lgammaf \ diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c index 1355eb7998f..0ebe7b91cb7 100644 --- a/dlls/msvcrt/math.c +++ b/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; } /********************************************************************* diff --git a/dlls/msvcrt/unixlib.c b/dlls/msvcrt/unixlib.c index 8be8559337c..6b33eed223a 100644 --- a/dlls/msvcrt/unixlib.c +++ b/dlls/msvcrt/unixlib.c @@ -401,19 +401,6 @@ static float CDECL unix_hypotf(float x, float y) return hypotf( x, y ); } -/********************************************************************* - * j1 - */ -static double CDECL unix_j1(double num) -{ -#ifdef HAVE_J1 - return j1(num); -#else - FIXME("not implemented\n"); - return 0; -#endif -} - /********************************************************************* * jn */ @@ -1027,7 +1014,6 @@ static const struct unix_funcs funcs = unix_frexpf, unix_hypot, unix_hypotf, - unix_j1, unix_jn, unix_ldexp, unix_lgamma, diff --git a/dlls/msvcrt/unixlib.h b/dlls/msvcrt/unixlib.h index 6a23637c7e2..d7b1827e08e 100644 --- a/dlls/msvcrt/unixlib.h +++ b/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); diff --git a/include/config.h.in b/include/config.h.in index a37ab5383ff..8985cf6b474 100644 --- a/include/config.h.in +++ b/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