diff options
author | Eric Andersen <andersen@codepoet.org> | 2001-11-22 14:04:29 +0000 |
---|---|---|
committer | Eric Andersen <andersen@codepoet.org> | 2001-11-22 14:04:29 +0000 |
commit | 7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch) | |
tree | 3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/double/i1.c | |
parent | c117dd5fb183afb1a4790a6f6110d88704be6bf8 (diff) |
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
Diffstat (limited to 'libm/double/i1.c')
-rw-r--r-- | libm/double/i1.c | 402 |
1 files changed, 0 insertions, 402 deletions
diff --git a/libm/double/i1.c b/libm/double/i1.c deleted file mode 100644 index dfde216dc..000000000 --- a/libm/double/i1.c +++ /dev/null @@ -1,402 +0,0 @@ -/* i1.c - * - * Modified Bessel function of order one - * - * - * - * SYNOPSIS: - * - * double x, y, i1(); - * - * y = i1( x ); - * - * - * - * DESCRIPTION: - * - * Returns modified Bessel function of order one of the - * argument. - * - * The function is defined as i1(x) = -i j1( ix ). - * - * The range is partitioned into the two intervals [0,8] and - * (8, infinity). Chebyshev polynomial expansions are employed - * in each interval. - * - * - * - * ACCURACY: - * - * Relative error: - * arithmetic domain # trials peak rms - * DEC 0, 30 3400 1.2e-16 2.3e-17 - * IEEE 0, 30 30000 1.9e-15 2.1e-16 - * - * - */ -/* i1e.c - * - * Modified Bessel function of order one, - * exponentially scaled - * - * - * - * SYNOPSIS: - * - * double x, y, i1e(); - * - * y = i1e( x ); - * - * - * - * DESCRIPTION: - * - * Returns exponentially scaled modified Bessel function - * of order one of the argument. - * - * The function is defined as i1(x) = -i exp(-|x|) j1( ix ). - * - * - * - * ACCURACY: - * - * Relative error: - * arithmetic domain # trials peak rms - * IEEE 0, 30 30000 2.0e-15 2.0e-16 - * See i1(). - * - */ - -/* i1.c 2 */ - - -/* -Cephes Math Library Release 2.8: June, 2000 -Copyright 1985, 1987, 2000 by Stephen L. Moshier -*/ - -#include <math.h> - -/* Chebyshev coefficients for exp(-x) I1(x) / x - * in the interval [0,8]. - * - * lim(x->0){ exp(-x) I1(x) / x } = 1/2. - */ - -#ifdef UNK -static double A[] = -{ - 2.77791411276104639959E-18, --2.11142121435816608115E-17, - 1.55363195773620046921E-16, --1.10559694773538630805E-15, - 7.60068429473540693410E-15, --5.04218550472791168711E-14, - 3.22379336594557470981E-13, --1.98397439776494371520E-12, - 1.17361862988909016308E-11, --6.66348972350202774223E-11, - 3.62559028155211703701E-10, --1.88724975172282928790E-9, - 9.38153738649577178388E-9, --4.44505912879632808065E-8, - 2.00329475355213526229E-7, --8.56872026469545474066E-7, - 3.47025130813767847674E-6, --1.32731636560394358279E-5, - 4.78156510755005422638E-5, --1.61760815825896745588E-4, - 5.12285956168575772895E-4, --1.51357245063125314899E-3, - 4.15642294431288815669E-3, --1.05640848946261981558E-2, - 2.47264490306265168283E-2, --5.29459812080949914269E-2, - 1.02643658689847095384E-1, --1.76416518357834055153E-1, - 2.52587186443633654823E-1 -}; -#endif - -#ifdef DEC -static unsigned short A[] = { -0021514,0174520,0060742,0000241, -0122302,0137206,0016120,0025663, -0023063,0017437,0026235,0176536, -0123637,0052523,0170150,0125632, -0024410,0165770,0030251,0044134, -0125143,0012160,0162170,0054727, -0025665,0075702,0035716,0145247, -0126413,0116032,0176670,0015462, -0027116,0073425,0110351,0105242, -0127622,0104034,0137530,0037364, -0030307,0050645,0120776,0175535, -0131001,0130331,0043523,0037455, -0031441,0026160,0010712,0100174, -0132076,0164761,0022706,0017500, -0032527,0015045,0115076,0104076, -0133146,0001714,0015434,0144520, -0033550,0161166,0124215,0077050, -0134136,0127715,0143365,0157170, -0034510,0106652,0013070,0064130, -0135051,0117126,0117264,0123761, -0035406,0045355,0133066,0175751, -0135706,0061420,0054746,0122440, -0036210,0031232,0047235,0006640, -0136455,0012373,0144235,0011523, -0036712,0107437,0036731,0015111, -0137130,0156742,0115744,0172743, -0037322,0033326,0124667,0124740, -0137464,0123210,0021510,0144556, -0037601,0051433,0111123,0177721 -}; -#endif - -#ifdef IBMPC -static unsigned short A[] = { -0x4014,0x0c3c,0x9f2a,0x3c49, -0x0576,0xc38a,0x57d0,0xbc78, -0xbfac,0xe593,0x63e3,0x3ca6, -0x1573,0x7e0d,0xeaaa,0xbcd3, -0x290c,0x0615,0x1d7f,0x3d01, -0x0b3b,0x1c8f,0x628e,0xbd2c, -0xd955,0x4779,0xaf78,0x3d56, -0x0366,0x5fb7,0x7383,0xbd81, -0x3154,0xb21d,0xcee2,0x3da9, -0x07de,0x97eb,0x5103,0xbdd2, -0xdf6c,0xb43f,0xea34,0x3df8, -0x67e6,0x28ea,0x361b,0xbe20, -0x5010,0x0239,0x258e,0x3e44, -0xc3e8,0x24b8,0xdd3e,0xbe67, -0xd108,0xb347,0xe344,0x3e8a, -0x992a,0x8363,0xc079,0xbeac, -0xafc5,0xd511,0x1c4e,0x3ecd, -0xbbcf,0xb8de,0xd5f9,0xbeeb, -0x0d0b,0x42c7,0x11b5,0x3f09, -0x94fe,0xd3d6,0x33ca,0xbf25, -0xdf7d,0xb6c6,0xc95d,0x3f40, -0xd4a4,0x0b3c,0xcc62,0xbf58, -0xa1b4,0x49d3,0x0653,0x3f71, -0xa26a,0x7913,0xa29f,0xbf85, -0x2349,0xe7bb,0x51e3,0x3f99, -0x9ebc,0x537c,0x1bbc,0xbfab, -0xf53c,0xd536,0x46da,0x3fba, -0x192e,0x0469,0x94d1,0xbfc6, -0x7ffa,0x724a,0x2a63,0x3fd0 -}; -#endif - -#ifdef MIEEE -static unsigned short A[] = { -0x3c49,0x9f2a,0x0c3c,0x4014, -0xbc78,0x57d0,0xc38a,0x0576, -0x3ca6,0x63e3,0xe593,0xbfac, -0xbcd3,0xeaaa,0x7e0d,0x1573, -0x3d01,0x1d7f,0x0615,0x290c, -0xbd2c,0x628e,0x1c8f,0x0b3b, -0x3d56,0xaf78,0x4779,0xd955, -0xbd81,0x7383,0x5fb7,0x0366, -0x3da9,0xcee2,0xb21d,0x3154, -0xbdd2,0x5103,0x97eb,0x07de, -0x3df8,0xea34,0xb43f,0xdf6c, -0xbe20,0x361b,0x28ea,0x67e6, -0x3e44,0x258e,0x0239,0x5010, -0xbe67,0xdd3e,0x24b8,0xc3e8, -0x3e8a,0xe344,0xb347,0xd108, -0xbeac,0xc079,0x8363,0x992a, -0x3ecd,0x1c4e,0xd511,0xafc5, -0xbeeb,0xd5f9,0xb8de,0xbbcf, -0x3f09,0x11b5,0x42c7,0x0d0b, -0xbf25,0x33ca,0xd3d6,0x94fe, -0x3f40,0xc95d,0xb6c6,0xdf7d, -0xbf58,0xcc62,0x0b3c,0xd4a4, -0x3f71,0x0653,0x49d3,0xa1b4, -0xbf85,0xa29f,0x7913,0xa26a, -0x3f99,0x51e3,0xe7bb,0x2349, -0xbfab,0x1bbc,0x537c,0x9ebc, -0x3fba,0x46da,0xd536,0xf53c, -0xbfc6,0x94d1,0x0469,0x192e, -0x3fd0,0x2a63,0x724a,0x7ffa -}; -#endif - -/* i1.c */ - -/* Chebyshev coefficients for exp(-x) sqrt(x) I1(x) - * in the inverted interval [8,infinity]. - * - * lim(x->inf){ exp(-x) sqrt(x) I1(x) } = 1/sqrt(2pi). - */ - -#ifdef UNK -static double B[] = -{ - 7.51729631084210481353E-18, - 4.41434832307170791151E-18, --4.65030536848935832153E-17, --3.20952592199342395980E-17, - 2.96262899764595013876E-16, - 3.30820231092092828324E-16, --1.88035477551078244854E-15, --3.81440307243700780478E-15, - 1.04202769841288027642E-14, - 4.27244001671195135429E-14, --2.10154184277266431302E-14, --4.08355111109219731823E-13, --7.19855177624590851209E-13, - 2.03562854414708950722E-12, - 1.41258074366137813316E-11, - 3.25260358301548823856E-11, --1.89749581235054123450E-11, --5.58974346219658380687E-10, --3.83538038596423702205E-9, --2.63146884688951950684E-8, --2.51223623787020892529E-7, --3.88256480887769039346E-6, --1.10588938762623716291E-4, --9.76109749136146840777E-3, - 7.78576235018280120474E-1 -}; -#endif - -#ifdef DEC -static unsigned short B[] = { -0022012,0125555,0115227,0043456, -0021642,0156127,0052075,0145203, -0122526,0072435,0111231,0011664, -0122424,0001544,0161671,0114403, -0023252,0144257,0163532,0142121, -0023276,0132162,0174045,0013204, -0124007,0077154,0057046,0110517, -0124211,0066650,0116127,0157073, -0024473,0133413,0130551,0107504, -0025100,0064741,0032631,0040364, -0124675,0045101,0071551,0012400, -0125745,0161054,0071637,0011247, -0126112,0117410,0035525,0122231, -0026417,0037237,0131034,0176427, -0027170,0100373,0024742,0025725, -0027417,0006417,0105303,0141446, -0127246,0163716,0121202,0060137, -0130431,0123122,0120436,0166000, -0131203,0144134,0153251,0124500, -0131742,0005234,0122732,0033006, -0132606,0157751,0072362,0121031, -0133602,0043372,0047120,0015626, -0134747,0165774,0001125,0046462, -0136437,0166402,0117746,0155137, -0040107,0050305,0125330,0124241 -}; -#endif - -#ifdef IBMPC -static unsigned short B[] = { -0xe8e6,0xb352,0x556d,0x3c61, -0xb950,0xea87,0x5b8a,0x3c54, -0x2277,0xb253,0xcea3,0xbc8a, -0x3320,0x9c77,0x806c,0xbc82, -0x588a,0xfceb,0x5915,0x3cb5, -0xa2d1,0x5f04,0xd68e,0x3cb7, -0xd22a,0x8bc4,0xefcd,0xbce0, -0xfbc7,0x138a,0x2db5,0xbcf1, -0x31e8,0x762d,0x76e1,0x3d07, -0x281e,0x26b3,0x0d3c,0x3d28, -0x22a0,0x2e6d,0xa948,0xbd17, -0xe255,0x8e73,0xbc45,0xbd5c, -0xb493,0x076a,0x53e1,0xbd69, -0x9fa3,0xf643,0xe7d3,0x3d81, -0x457b,0x653c,0x101f,0x3daf, -0x7865,0xf158,0xe1a1,0x3dc1, -0x4c0c,0xd450,0xdcf9,0xbdb4, -0xdd80,0x5423,0x34ca,0xbe03, -0x3528,0x9ad5,0x790b,0xbe30, -0x46c1,0x94bb,0x4153,0xbe5c, -0x5443,0x2e9e,0xdbfd,0xbe90, -0x0373,0x49ca,0x48df,0xbed0, -0xa9a6,0x804a,0xfd7f,0xbf1c, -0xdb4c,0x53fc,0xfda0,0xbf83, -0x1514,0xb55b,0xea18,0x3fe8 -}; -#endif - -#ifdef MIEEE -static unsigned short B[] = { -0x3c61,0x556d,0xb352,0xe8e6, -0x3c54,0x5b8a,0xea87,0xb950, -0xbc8a,0xcea3,0xb253,0x2277, -0xbc82,0x806c,0x9c77,0x3320, -0x3cb5,0x5915,0xfceb,0x588a, -0x3cb7,0xd68e,0x5f04,0xa2d1, -0xbce0,0xefcd,0x8bc4,0xd22a, -0xbcf1,0x2db5,0x138a,0xfbc7, -0x3d07,0x76e1,0x762d,0x31e8, -0x3d28,0x0d3c,0x26b3,0x281e, -0xbd17,0xa948,0x2e6d,0x22a0, -0xbd5c,0xbc45,0x8e73,0xe255, -0xbd69,0x53e1,0x076a,0xb493, -0x3d81,0xe7d3,0xf643,0x9fa3, -0x3daf,0x101f,0x653c,0x457b, -0x3dc1,0xe1a1,0xf158,0x7865, -0xbdb4,0xdcf9,0xd450,0x4c0c, -0xbe03,0x34ca,0x5423,0xdd80, -0xbe30,0x790b,0x9ad5,0x3528, -0xbe5c,0x4153,0x94bb,0x46c1, -0xbe90,0xdbfd,0x2e9e,0x5443, -0xbed0,0x48df,0x49ca,0x0373, -0xbf1c,0xfd7f,0x804a,0xa9a6, -0xbf83,0xfda0,0x53fc,0xdb4c, -0x3fe8,0xea18,0xb55b,0x1514 -}; -#endif - -/* i1.c */ -#ifdef ANSIPROT -extern double chbevl ( double, void *, int ); -extern double exp ( double ); -extern double sqrt ( double ); -extern double fabs ( double ); -#else -double chbevl(), exp(), sqrt(), fabs(); -#endif - -double i1(x) -double x; -{ -double y, z; - -z = fabs(x); -if( z <= 8.0 ) - { - y = (z/2.0) - 2.0; - z = chbevl( y, A, 29 ) * z * exp(z); - } -else - { - z = exp(z) * chbevl( 32.0/z - 2.0, B, 25 ) / sqrt(z); - } -if( x < 0.0 ) - z = -z; -return( z ); -} - -/* i1e() */ - -double i1e( x ) -double x; -{ -double y, z; - -z = fabs(x); -if( z <= 8.0 ) - { - y = (z/2.0) - 2.0; - z = chbevl( y, A, 29 ) * z; - } -else - { - z = chbevl( 32.0/z - 2.0, B, 25 ) / sqrt(z); - } -if( x < 0.0 ) - z = -z; -return( z ); -} |