Totally rework the math library, this time based on the MacOs X

math library (which is itself based on the math lib from FreeBSD).
 -Erik

1 files changed, 0 insertions, 551 deletions

diff --git a/libm/ldouble/j1l.c b/libm/ldouble/j1l.c
deleted file mode 100644
index 83428473e..000000000
--- a/libm/ldouble/j1l.c
+++ /dev/null
@@ -1,551 +0,0 @@
-/*							j1l.c
- *
- *	Bessel function of order one
- *
- *
- *
- * SYNOPSIS:
- *
- * long double x, y, j1l();
- *
- * y = j1l( x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns Bessel function of order one of the argument.
- *
- * The domain is divided into the intervals [0, 9] and
- * (9, infinity). In the first interval the rational approximation
- * is (x^2 - r^2) (x^2 - s^2) (x^2 - t^2) x P8(x^2) / Q8(x^2),
- * where r, s, t are the first three zeros of the function.
- * In the second interval the expansion is in terms of the
- * modulus M1(x) = sqrt(J1(x)^2 + Y1(x)^2) and phase  P1(x)
- * = atan(Y1(x)/J1(x)).  M1 is approximated by sqrt(1/x)P7(1/x)/Q8(1/x).
- * The approximation to j1 is M1 * cos(x -  3 pi/4 + 1/x P5(1/x^2)/Q6(1/x^2)).
- *
- *
- * ACCURACY:
- *
- *                      Absolute error:
- * arithmetic   domain      # trials      peak         rms
- *    IEEE      0, 30        40000      1.8e-19      5.0e-20
- *
- *
- */
-/*							y1l.c
- *
- *	Bessel function of the second kind, order zero
- *
- *
- *
- * SYNOPSIS:
- *
- * double x, y, y1l();
- *
- * y = y1l( x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns Bessel function of the second kind, of order
- * zero, of the argument.
- *
- * The domain is divided into the intervals [0, 4.5>, [4.5,9> and
- * [9, infinity). In the first interval a rational approximation
- * R(x) is employed to compute y0(x)  = R(x) + 2/pi * log(x) * j0(x).
- *
- * In the second interval, the approximation is
- *     (x - p)(x - q)(x - r)(x - s)P9(x)/Q10(x)
- * where p, q, r, s are zeros of y1(x).
- *
- * The third interval uses the same approximations to modulus
- * and phase as j1(x), whence y1(x) = modulus * sin(phase).
- *
- * ACCURACY:
- *
- *  Absolute error, when y0(x) < 1; else relative error:
- *
- * arithmetic   domain     # trials      peak         rms
- *    IEEE      0, 30       36000       2.7e-19     5.3e-20
- *
- */
-
-/* Copyright 1994 by Stephen L. Moshier (moshier@world.std.com).  */
-
-#include <math.h>
-
-/*
-j1(x) = (x^2-r0^2)(x^2-r1^2)(x^2-r2^2) x P(x**2)/Q(x**2)
-0 <= x <= 9
-Relative error
-n=8, d=8
-Peak error =  2e-21
-*/
-#if UNK
-static long double j1n[9] = {
--2.63469779622127762897E-4L,
- 9.31329762279632791262E-1L,
--1.46280142797793933909E3L,
- 1.32000129539331214495E6L,
--7.41183271195454042842E8L,
- 2.626500686552841932403E11L,
--5.68263073022183470933E13L,
- 6.80006297997263446982E15L,
--3.41470097444474566748E17L,
-};
-static long double j1d[8] = {
-/* 1.00000000000000000000E0L,*/
- 2.95267951972943745733E3L,
- 4.78723926343829674773E6L,
- 5.37544732957807543920E9L,
- 4.46866213886267829490E12L,
- 2.76959756375961607085E15L,
- 1.23367806884831151194E18L,
- 3.57325874689695599524E20L,
- 5.10779045516141578461E22L,
-};
-#endif
-#if IBMPC
-static short j1n[] = {
-0xf72f,0x18cc,0x50b2,0x8a22,0xbff3, XPD
-0x6dc3,0xc850,0xa096,0xee6b,0x3ffe, XPD
-0x29f3,0x496b,0xa54c,0xb6d9,0xc009, XPD
-0x38f5,0xf72b,0x0a5c,0xa122,0x4013, XPD
-0x1ac8,0xc825,0x3c9c,0xb0b6,0xc01c, XPD
-0x038e,0xbd23,0xa7fa,0xf49c,0x4024, XPD
-0x636c,0x4d29,0x9f71,0xcebb,0xc02c, XPD
-0xd3c2,0xf8f0,0xf852,0xc144,0x4033, XPD
-0xd8d8,0x7311,0xa7d2,0x97a4,0xc039, XPD
-};
-static short j1d[] = {
-/*0x0000,0x0000,0x0000,0x8000,0x3fff, XPD*/
-0xbaf9,0x146e,0xdf50,0xb88a,0x400a, XPD
-0x6a17,0xe162,0x4e86,0x9218,0x4015, XPD
-0x6041,0xc9fe,0x6890,0xa033,0x401f, XPD
-0xb498,0xfdd5,0x209e,0x820e,0x4029, XPD
-0x0122,0x56c0,0xf2ef,0x9d6e,0x4032, XPD
-0xe6c0,0xa725,0x3d56,0x88f7,0x403b, XPD
-0x665d,0xb178,0x242e,0x9af7,0x4043, XPD
-0xdd67,0xf5b3,0x0522,0xad0f,0x404a, XPD
-};
-#endif
-#if MIEEE
-static long j1n[27] = {
-0xbff30000,0x8a2250b2,0x18ccf72f,
-0x3ffe0000,0xee6ba096,0xc8506dc3,
-0xc0090000,0xb6d9a54c,0x496b29f3,
-0x40130000,0xa1220a5c,0xf72b38f5,
-0xc01c0000,0xb0b63c9c,0xc8251ac8,
-0x40240000,0xf49ca7fa,0xbd23038e,
-0xc02c0000,0xcebb9f71,0x4d29636c,
-0x40330000,0xc144f852,0xf8f0d3c2,
-0xc0390000,0x97a4a7d2,0x7311d8d8,
-};
-static long j1d[24] = {
-/*0x3fff0000,0x80000000,0x00000000,*/
-0x400a0000,0xb88adf50,0x146ebaf9,
-0x40150000,0x92184e86,0xe1626a17,
-0x401f0000,0xa0336890,0xc9fe6041,
-0x40290000,0x820e209e,0xfdd5b498,
-0x40320000,0x9d6ef2ef,0x56c00122,
-0x403b0000,0x88f73d56,0xa725e6c0,
-0x40430000,0x9af7242e,0xb178665d,
-0x404a0000,0xad0f0522,0xf5b3dd67,
-};
-#endif
-/*
-sqrt(j0^2(1/x^2) + y0^2(1/x^2)) = z P(z**2)/Q(z**2)
-z(x) = 1/sqrt(x)
-Relative error
-n=7, d=8
-Peak error =  1.35e=20
-Relative error spread =  9.9e0
-*/
-#if UNK
-static long double modulusn[8] = {
--5.041742205078442098874E0L,
- 3.918474430130242177355E-1L,
- 2.527521168680500659056E0L,
- 7.172146812845906480743E0L,
- 2.859499532295180940060E0L,
- 1.014671139779858141347E0L,
- 1.255798064266130869132E-1L,
- 1.596507617085714650238E-2L,
-};
-static long double modulusd[8] = {
-/* 1.000000000000000000000E0L,*/
--6.233092094568239317498E0L,
--9.214128701852838347002E-1L,
- 2.531772200570435289832E0L,
- 8.755081357265851765640E0L,
- 3.554340386955608261463E0L,
- 1.267949948774331531237E0L,
- 1.573909467558180942219E-1L,
- 2.000925566825407466160E-2L,
-};
-#endif
-#if IBMPC
-static short modulusn[] = {
-0x3d53,0xb598,0xf3bf,0xa155,0xc001, XPD
-0x3111,0x863a,0x3a61,0xc8a0,0x3ffd, XPD
-0x7d55,0xdb8c,0xe825,0xa1c2,0x4000, XPD
-0xe5e2,0x6914,0x3a08,0xe582,0x4001, XPD
-0x71e6,0x88a5,0x0a53,0xb702,0x4000, XPD
-0x2cb0,0xc657,0xbe70,0x81e0,0x3fff, XPD
-0x6de4,0x8fae,0xfe26,0x8097,0x3ffc, XPD
-0xa905,0x05fb,0x3101,0x82c9,0x3ff9, XPD
-};
-static short modulusd[] = {
-/*0x0000,0x0000,0x0000,0x8000,0x3fff, XPD*/
-0x2603,0x640e,0x7d8d,0xc775,0xc001, XPD
-0x77b5,0x8f2d,0xb6bf,0xebe1,0xbffe, XPD
-0x6420,0x97ce,0x8e44,0xa208,0x4000, XPD
-0x0260,0x746b,0xd030,0x8c14,0x4002, XPD
-0x77b6,0x34e2,0x501a,0xe37a,0x4000, XPD
-0x37ce,0x79ae,0x2f15,0xa24c,0x3fff, XPD
-0xfc82,0x02c7,0x17a4,0xa12b,0x3ffc, XPD
-0x1237,0xcc6c,0x7356,0xa3ea,0x3ff9, XPD
-};
-#endif
-#if MIEEE
-static long modulusn[24] = {
-0xc0010000,0xa155f3bf,0xb5983d53,
-0x3ffd0000,0xc8a03a61,0x863a3111,
-0x40000000,0xa1c2e825,0xdb8c7d55,
-0x40010000,0xe5823a08,0x6914e5e2,
-0x40000000,0xb7020a53,0x88a571e6,
-0x3fff0000,0x81e0be70,0xc6572cb0,
-0x3ffc0000,0x8097fe26,0x8fae6de4,
-0x3ff90000,0x82c93101,0x05fba905,
-};
-static long modulusd[24] = {
-/*0x3fff0000,0x80000000,0x00000000,*/
-0xc0010000,0xc7757d8d,0x640e2603,
-0xbffe0000,0xebe1b6bf,0x8f2d77b5,
-0x40000000,0xa2088e44,0x97ce6420,
-0x40020000,0x8c14d030,0x746b0260,
-0x40000000,0xe37a501a,0x34e277b6,
-0x3fff0000,0xa24c2f15,0x79ae37ce,
-0x3ffc0000,0xa12b17a4,0x02c7fc82,
-0x3ff90000,0xa3ea7356,0xcc6c1237,
-};
-#endif
-/*
-atan(y1(x)/j1(x))  =  x - 3pi/4 + z P(z**2)/Q(z**2)
-z(x) = 1/x
-Absolute error
-n=5, d=6
-Peak error =  4.83e-21
-Relative error spread =  1.9e0
-*/
-#if UNK
-static long double phasen[6] = {
- 2.010456367705144783933E0L,
- 1.587378144541918176658E0L,
- 2.682837461073751055565E-1L,
- 1.472572645054468815027E-2L,
- 2.884976126715926258586E-4L,
- 1.708502235134706284899E-6L,
-};
-static long double phased[6] = {
-/* 1.000000000000000000000E0L,*/
- 6.809332495854873089362E0L,
- 4.518597941618813112665E0L,
- 7.320149039410806471101E-1L,
- 3.960155028960712309814E-2L,
- 7.713202197319040439861E-4L,
- 4.556005960359216767984E-6L,
-};
-#endif
-#if IBMPC
-static short phasen[] = {
-0xebc0,0x5506,0x512f,0x80ab,0x4000, XPD
-0x6050,0x98aa,0x3500,0xcb2f,0x3fff, XPD
-0xe907,0x28b9,0x7cb7,0x895c,0x3ffd, XPD
-0xa830,0xf4a3,0x2c60,0xf144,0x3ff8, XPD
-0xf74f,0xbe87,0x7e7d,0x9741,0x3ff3, XPD
-0x540c,0xc1d5,0xb096,0xe54f,0x3feb, XPD
-};
-static short phased[] = {
-/*0x0000,0x0000,0x0000,0x8000,0x3fff, XPD*/
-0xefe3,0x292c,0x0d43,0xd9e6,0x4001, XPD
-0xb1f2,0xe0d2,0x5ab5,0x9098,0x4001, XPD
-0xc39e,0x9c8c,0x5428,0xbb65,0x3ffe, XPD
-0x98f8,0xd610,0x3c35,0xa235,0x3ffa, XPD
-0xa853,0x55fb,0x6c79,0xca32,0x3ff4, XPD
-0x8d72,0x2be3,0xcb0f,0x98df,0x3fed, XPD
-};
-#endif
-#if MIEEE
-static long phasen[18] = {
-0x40000000,0x80ab512f,0x5506ebc0,
-0x3fff0000,0xcb2f3500,0x98aa6050,
-0x3ffd0000,0x895c7cb7,0x28b9e907,
-0x3ff80000,0xf1442c60,0xf4a3a830,
-0x3ff30000,0x97417e7d,0xbe87f74f,
-0x3feb0000,0xe54fb096,0xc1d5540c,
-};
-static long phased[18] = {
-/*0x3fff0000,0x80000000,0x00000000,*/
-0x40010000,0xd9e60d43,0x292cefe3,
-0x40010000,0x90985ab5,0xe0d2b1f2,
-0x3ffe0000,0xbb655428,0x9c8cc39e,
-0x3ffa0000,0xa2353c35,0xd61098f8,
-0x3ff40000,0xca326c79,0x55fba853,
-0x3fed0000,0x98dfcb0f,0x2be38d72,
-};
-#endif
-#define JZ1 1.46819706421238932572e1L
-#define JZ2 4.92184563216946036703e1L
-#define JZ3 1.03499453895136580332e2L
-
-#define THPIO4L  2.35619449019234492885L
-#ifdef ANSIPROT
-extern long double sqrtl ( long double );
-extern long double fabsl ( long double );
-extern long double polevll ( long double, void *, int );
-extern long double p1evll ( long double, void *, int );
-extern long double cosl ( long double );
-extern long double sinl ( long double );
-extern long double logl ( long double );
-long double j1l (long double );
-#else
-long double sqrtl(), fabsl(), polevll(), p1evll(), cosl(), sinl(), logl();
-long double j1l();
-#endif
-
-long double j1l(x)
-long double x;
-{
-long double xx, y, z, modulus, phase;
-
-xx = x * x;
-if( xx < 81.0L )
-  {
-    y = (xx - JZ1) * (xx - JZ2) * (xx - JZ3);
-    y *= x * polevll( xx, j1n, 8 ) / p1evll( xx, j1d, 8 );
-    return y;
-  }
-
-y = fabsl(x);
-xx = 1.0/xx;
-phase = polevll( xx, phasen, 5 ) / p1evll( xx, phased, 6 );
-
-z = 1.0/y;
-modulus = polevll( z, modulusn, 7 ) / p1evll( z, modulusd, 8 );
-
-y = modulus * cosl( y -  THPIO4L + z*phase) / sqrtl(y);
-if( x < 0 )
-  y = -y;
-return y;
-}
-
-/*
-y1(x) = 2/pi * (log(x) * j1(x) - 1/x) + R(x^2) z P(z**2)/Q(z**2)
-0 <= x <= 4.5
-z(x) = x
-Absolute error
-n=6, d=7
-Peak error =  7.25e-22
-Relative error spread =  4.5e-2
-*/
-#if UNK
-static long double y1n[7] = {
--1.288901054372751879531E5L,
- 9.914315981558815369372E7L,
--2.906793378120403577274E10L,
- 3.954354656937677136266E12L,
--2.445982226888344140154E14L,
- 5.685362960165615942886E15L,
--2.158855258453711703120E16L,
-};
-static long double y1d[7] = {
-/* 1.000000000000000000000E0L,*/
- 8.926354644853231136073E2L,
- 4.679841933793707979659E5L,
- 1.775133253792677466651E8L,
- 5.089532584184822833416E10L,
- 1.076474894829072923244E13L,
- 1.525917240904692387994E15L,
- 1.101136026928555260168E17L,
-};
-#endif
-#if IBMPC
-static short y1n[] = {
-0x5b16,0xf7f8,0x0d7e,0xfbbd,0xc00f, XPD
-0x53e4,0x194c,0xbefa,0xbd19,0x4019, XPD
-0x7607,0xa687,0xaf0a,0xd892,0xc021, XPD
-0x5633,0xaa6b,0x79e5,0xe62c,0x4028, XPD
-0x69fd,0x1242,0xf62d,0xde75,0xc02e, XPD
-0x7f8b,0x4757,0x75bd,0xa196,0x4033, XPD
-0x3a10,0x0848,0x5930,0x9965,0xc035, XPD
-};
-static short y1d[] = {
-/*0x0000,0x0000,0x0000,0x8000,0x3fff, XPD*/
-0xdd1a,0x3b8e,0xab73,0xdf28,0x4008, XPD
-0x298c,0x29ef,0x0630,0xe482,0x4011, XPD
-0x0e86,0x117b,0x36d6,0xa94a,0x401a, XPD
-0x57e0,0x1d92,0x90a9,0xbd99,0x4022, XPD
-0xaaf0,0x342b,0xd098,0x9ca5,0x402a, XPD
-0x8c6a,0x397e,0x0963,0xad7a,0x4031, XPD
-0x7302,0xb91b,0xde7e,0xc399,0x4037, XPD
-};
-#endif
-#if MIEEE
-static long y1n[21] = {
-0xc00f0000,0xfbbd0d7e,0xf7f85b16,
-0x40190000,0xbd19befa,0x194c53e4,
-0xc0210000,0xd892af0a,0xa6877607,
-0x40280000,0xe62c79e5,0xaa6b5633,
-0xc02e0000,0xde75f62d,0x124269fd,
-0x40330000,0xa19675bd,0x47577f8b,
-0xc0350000,0x99655930,0x08483a10,
-};
-static long y1d[21] = {
-/*0x3fff0000,0x80000000,0x00000000,*/
-0x40080000,0xdf28ab73,0x3b8edd1a,
-0x40110000,0xe4820630,0x29ef298c,
-0x401a0000,0xa94a36d6,0x117b0e86,
-0x40220000,0xbd9990a9,0x1d9257e0,
-0x402a0000,0x9ca5d098,0x342baaf0,
-0x40310000,0xad7a0963,0x397e8c6a,
-0x40370000,0xc399de7e,0xb91b7302,
-};
-#endif
-/*
-y1(x) = (x-YZ1)(x-YZ2)(x-YZ3)(x-YZ4)R(x) P(z)/Q(z)
-z(x) = x
-4.5 <= x <= 9
-Absolute error
-n=9, d=10
-Peak error =  3.27e-22
-Relative error spread =  4.5e-2
-*/
-#if UNK
-static long double y159n[10] = {
--6.806634906054210550896E-1L,
- 4.306669585790359450532E1L,
--9.230477746767243316014E2L,
- 6.171186628598134035237E3L,
- 2.096869860275353982829E4L,
--1.238961670382216747944E5L,
--1.781314136808997406109E6L,
--1.803400156074242435454E6L,
--1.155761550219364178627E6L,
- 3.112221202330688509818E5L,
-};
-static long double y159d[10] = {
-/* 1.000000000000000000000E0L,*/
--6.181482377814679766978E1L,
- 2.238187927382180589099E3L,
--5.225317824142187494326E4L,
- 9.217235006983512475118E5L,
--1.183757638771741974521E7L,
- 1.208072488974110742912E8L,
--8.193431077523942651173E8L,
- 4.282669747880013349981E9L,
--1.171523459555524458808E9L,
- 1.078445545755236785692E8L,
-};
-#endif
-#if IBMPC
-static short y159n[] = {
-0xb5e5,0xbb42,0xf667,0xae3f,0xbffe, XPD
-0xfdf1,0x41e5,0x4beb,0xac44,0x4004, XPD
-0xe917,0x8486,0x0ebd,0xe6c3,0xc008, XPD
-0xdf40,0x226b,0x7e37,0xc0d9,0x400b, XPD
-0xb2bf,0x4296,0x65af,0xa3d1,0x400d, XPD
-0xa33b,0x8229,0x1561,0xf1fc,0xc00f, XPD
-0xcd43,0x2f50,0x1118,0xd972,0xc013, XPD
-0x3811,0xa3da,0x413f,0xdc24,0xc013, XPD
-0xf62f,0xd968,0x8c66,0x8d15,0xc013, XPD
-0x539b,0xf305,0xc3d8,0x97f6,0x4011, XPD
-};
-static short y159d[] = {
-/*0x0000,0x0000,0x0000,0x8000,0x3fff, XPD*/
-0x1a6c,0x1c93,0x612a,0xf742,0xc004, XPD
-0xd0fe,0x2487,0x01c0,0x8be3,0x400a, XPD
-0xbed4,0x3ad5,0x2da1,0xcc1d,0xc00e, XPD
-0x3c4f,0xdc46,0xb802,0xe107,0x4012, XPD
-0xe5e5,0x4172,0x8863,0xb4a0,0xc016, XPD
-0x6de5,0xb797,0xea1c,0xe66b,0x4019, XPD
-0xa46a,0x0273,0xbc0f,0xc358,0xc01c, XPD
-0x8e0e,0xe148,0x5ab3,0xff44,0x401e, XPD
-0xb3ad,0x1c6d,0x0f07,0x8ba8,0xc01d, XPD
-0xa231,0x6ab0,0x7952,0xcdb2,0x4019, XPD
-};
-#endif
-#if MIEEE
-static long y159n[30] = {
-0xbffe0000,0xae3ff667,0xbb42b5e5,
-0x40040000,0xac444beb,0x41e5fdf1,
-0xc0080000,0xe6c30ebd,0x8486e917,
-0x400b0000,0xc0d97e37,0x226bdf40,
-0x400d0000,0xa3d165af,0x4296b2bf,
-0xc00f0000,0xf1fc1561,0x8229a33b,
-0xc0130000,0xd9721118,0x2f50cd43,
-0xc0130000,0xdc24413f,0xa3da3811,
-0xc0130000,0x8d158c66,0xd968f62f,
-0x40110000,0x97f6c3d8,0xf305539b,
-};
-static long y159d[30] = {
-/*0x3fff0000,0x80000000,0x00000000,*/
-0xc0040000,0xf742612a,0x1c931a6c,
-0x400a0000,0x8be301c0,0x2487d0fe,
-0xc00e0000,0xcc1d2da1,0x3ad5bed4,
-0x40120000,0xe107b802,0xdc463c4f,
-0xc0160000,0xb4a08863,0x4172e5e5,
-0x40190000,0xe66bea1c,0xb7976de5,
-0xc01c0000,0xc358bc0f,0x0273a46a,
-0x401e0000,0xff445ab3,0xe1488e0e,
-0xc01d0000,0x8ba80f07,0x1c6db3ad,
-0x40190000,0xcdb27952,0x6ab0a231,
-};
-#endif
-
-extern long double MAXNUML;
-/* #define MAXNUML 1.18973149535723176502e4932L */
-#define TWOOPI 6.36619772367581343075535e-1L
-#define THPIO4 2.35619449019234492885L
-#define Y1Z1 2.19714132603101703515e0L
-#define Y1Z2 5.42968104079413513277e0L
-#define Y1Z3 8.59600586833116892643e0L
-#define Y1Z4 1.17491548308398812434e1L
-
-long double y1l(x)
-long double x;
-{
-long double xx, y, z, modulus, phase;
-
-if( x < 0.0 )
-  {
-    return (-MAXNUML);
-  }
-z = 1.0/x;
-xx = x * x;
-if( xx < 81.0L )
-  {
-    if( xx < 20.25L )
-      {
-	y = TWOOPI * (logl(x) * j1l(x) - z);
-	y += x * polevll( xx, y1n, 6 ) / p1evll( xx, y1d, 7 );
-      }
-    else
-      {
-	y = (x - Y1Z1)*(x - Y1Z2)*(x - Y1Z3)*(x - Y1Z4);
-	y *= polevll( x, y159n, 9 ) / p1evll( x, y159d, 10 );
-      }
-    return y;
-  }
-
-xx = 1.0/xx;
-phase = polevll( xx, phasen, 5 ) / p1evll( xx, phased, 6 );
-
-modulus = polevll( z, modulusn, 7 ) / p1evll( z, modulusd, 8 );
-
-z = modulus * sinl( x -  THPIO4L + z*phase) / sqrtl(x);
-return z;
-}