123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271 |
- /* 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 "mconf.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);
- }
|