source: sasview/src/sas/models/c_extension/cephes/i1.c @ b9dbd6b

/*                                                      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 );
}