lib/libm/i1.c

   1 /*                                                      i1.c
   2  *
   3  *      Modified Bessel function of order one
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * double x, y, i1();
  10  *
  11  * y = i1( x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Returns modified Bessel function of order one of the
  18  * argument.
  19  *
  20  * The function is defined as i1(x) = -i j1( ix ).
  21  *
  22  * The range is partitioned into the two intervals [0,8] and
  23  * (8, infinity).  Chebyshev polynomial expansions are employed
  24  * in each interval.
  25  *
  26  *
  27  *
  28  * ACCURACY:
  29  *
  30  *                      Relative error:
  31  * arithmetic   domain     # trials      peak         rms
  32  *    DEC       0, 30        3400       1.2e-16     2.3e-17
  33  *    IEEE      0, 30       30000       1.9e-15     2.1e-16
  34  *
  35  *
  36  */
  37 \f/*                                                     i1e.c
  38  *
  39  *      Modified Bessel function of order one,
  40  *      exponentially scaled
  41  *
  42  *
  43  *
  44  * SYNOPSIS:
  45  *
  46  * double x, y, i1e();
  47  *
  48  * y = i1e( x );
  49  *
  50  *
  51  *
  52  * DESCRIPTION:
  53  *
  54  * Returns exponentially scaled modified Bessel function
  55  * of order one of the argument.
  56  *
  57  * The function is defined as i1(x) = -i exp(-|x|) j1( ix ).
  58  *
  59  *
  60  *
  61  * ACCURACY:
  62  *
  63  *                      Relative error:
  64  * arithmetic   domain     # trials      peak         rms
  65  *    IEEE      0, 30       30000       2.0e-15     2.0e-16
  66  * See i1().
  67  *
  68  */
  69 \f
  70 /*                                                      i1.c 2          */
  71
  72
  73 /*
  74 Cephes Math Library Release 2.8:  June, 2000
  75 Copyright 1985, 1987, 2000 by Stephen L. Moshier
  76 */
  77
  78 #include "mconf.h"
  79
  80 /* Chebyshev coefficients for exp(-x) I1(x) / x
  81  * in the interval [0,8].
  82  *
  83  * lim(x->0){ exp(-x) I1(x) / x } = 1/2.
  84  */
  85
  86 #ifdef UNK
  87 static double A[] =
  88 {
  89  2.77791411276104639959E-18,
  90 -2.11142121435816608115E-17,
  91  1.55363195773620046921E-16,
  92 -1.10559694773538630805E-15,
  93  7.60068429473540693410E-15,
  94 -5.04218550472791168711E-14,
  95  3.22379336594557470981E-13,
  96 -1.98397439776494371520E-12,
  97  1.17361862988909016308E-11,
  98 -6.66348972350202774223E-11,
  99  3.62559028155211703701E-10,
 100 -1.88724975172282928790E-9,
 101  9.38153738649577178388E-9,
 102 -4.44505912879632808065E-8,
 103  2.00329475355213526229E-7,
 104 -8.56872026469545474066E-7,
 105  3.47025130813767847674E-6,
 106 -1.32731636560394358279E-5,
 107  4.78156510755005422638E-5,
 108 -1.61760815825896745588E-4,
 109  5.12285956168575772895E-4,
 110 -1.51357245063125314899E-3,
 111  4.15642294431288815669E-3,
 112 -1.05640848946261981558E-2,
 113  2.47264490306265168283E-2,
 114 -5.29459812080949914269E-2,
 115  1.02643658689847095384E-1,
 116 -1.76416518357834055153E-1,
 117  2.52587186443633654823E-1
 118 };
 119 #endif
 120
 121 #ifdef DEC
 122 static unsigned short A[] = {
 123 0021514,0174520,0060742,0000241,
 124 0122302,0137206,0016120,0025663,
 125 0023063,0017437,0026235,0176536,
 126 0123637,0052523,0170150,0125632,
 127 0024410,0165770,0030251,0044134,
 128 0125143,0012160,0162170,0054727,
 129 0025665,0075702,0035716,0145247,
 130 0126413,0116032,0176670,0015462,
 131 0027116,0073425,0110351,0105242,
 132 0127622,0104034,0137530,0037364,
 133 0030307,0050645,0120776,0175535,
 134 0131001,0130331,0043523,0037455,
 135 0031441,0026160,0010712,0100174,
 136 0132076,0164761,0022706,0017500,
 137 0032527,0015045,0115076,0104076,
 138 0133146,0001714,0015434,0144520,
 139 0033550,0161166,0124215,0077050,
 140 0134136,0127715,0143365,0157170,
 141 0034510,0106652,0013070,0064130,
 142 0135051,0117126,0117264,0123761,
 143 0035406,0045355,0133066,0175751,
 144 0135706,0061420,0054746,0122440,
 145 0036210,0031232,0047235,0006640,
 146 0136455,0012373,0144235,0011523,
 147 0036712,0107437,0036731,0015111,
 148 0137130,0156742,0115744,0172743,
 149 0037322,0033326,0124667,0124740,
 150 0137464,0123210,0021510,0144556,
 151 0037601,0051433,0111123,0177721
 152 };
 153 #endif
 154
 155 #ifdef IBMPC
 156 static unsigned short A[] = {
 157 0x4014,0x0c3c,0x9f2a,0x3c49,
 158 0x0576,0xc38a,0x57d0,0xbc78,
 159 0xbfac,0xe593,0x63e3,0x3ca6,
 160 0x1573,0x7e0d,0xeaaa,0xbcd3,
 161 0x290c,0x0615,0x1d7f,0x3d01,
 162 0x0b3b,0x1c8f,0x628e,0xbd2c,
 163 0xd955,0x4779,0xaf78,0x3d56,
 164 0x0366,0x5fb7,0x7383,0xbd81,
 165 0x3154,0xb21d,0xcee2,0x3da9,
 166 0x07de,0x97eb,0x5103,0xbdd2,
 167 0xdf6c,0xb43f,0xea34,0x3df8,
 168 0x67e6,0x28ea,0x361b,0xbe20,
 169 0x5010,0x0239,0x258e,0x3e44,
 170 0xc3e8,0x24b8,0xdd3e,0xbe67,
 171 0xd108,0xb347,0xe344,0x3e8a,
 172 0x992a,0x8363,0xc079,0xbeac,
 173 0xafc5,0xd511,0x1c4e,0x3ecd,
 174 0xbbcf,0xb8de,0xd5f9,0xbeeb,
 175 0x0d0b,0x42c7,0x11b5,0x3f09,
 176 0x94fe,0xd3d6,0x33ca,0xbf25,
 177 0xdf7d,0xb6c6,0xc95d,0x3f40,
 178 0xd4a4,0x0b3c,0xcc62,0xbf58,
 179 0xa1b4,0x49d3,0x0653,0x3f71,
 180 0xa26a,0x7913,0xa29f,0xbf85,
 181 0x2349,0xe7bb,0x51e3,0x3f99,
 182 0x9ebc,0x537c,0x1bbc,0xbfab,
 183 0xf53c,0xd536,0x46da,0x3fba,
 184 0x192e,0x0469,0x94d1,0xbfc6,
 185 0x7ffa,0x724a,0x2a63,0x3fd0
 186 };
 187 #endif
 188
 189 #ifdef MIEEE
 190 static unsigned short A[] = {
 191 0x3c49,0x9f2a,0x0c3c,0x4014,
 192 0xbc78,0x57d0,0xc38a,0x0576,
 193 0x3ca6,0x63e3,0xe593,0xbfac,
 194 0xbcd3,0xeaaa,0x7e0d,0x1573,
 195 0x3d01,0x1d7f,0x0615,0x290c,
 196 0xbd2c,0x628e,0x1c8f,0x0b3b,
 197 0x3d56,0xaf78,0x4779,0xd955,
 198 0xbd81,0x7383,0x5fb7,0x0366,
 199 0x3da9,0xcee2,0xb21d,0x3154,
 200 0xbdd2,0x5103,0x97eb,0x07de,
 201 0x3df8,0xea34,0xb43f,0xdf6c,
 202 0xbe20,0x361b,0x28ea,0x67e6,
 203 0x3e44,0x258e,0x0239,0x5010,
 204 0xbe67,0xdd3e,0x24b8,0xc3e8,
 205 0x3e8a,0xe344,0xb347,0xd108,
 206 0xbeac,0xc079,0x8363,0x992a,
 207 0x3ecd,0x1c4e,0xd511,0xafc5,
 208 0xbeeb,0xd5f9,0xb8de,0xbbcf,
 209 0x3f09,0x11b5,0x42c7,0x0d0b,
 210 0xbf25,0x33ca,0xd3d6,0x94fe,
 211 0x3f40,0xc95d,0xb6c6,0xdf7d,
 212 0xbf58,0xcc62,0x0b3c,0xd4a4,
 213 0x3f71,0x0653,0x49d3,0xa1b4,
 214 0xbf85,0xa29f,0x7913,0xa26a,
 215 0x3f99,0x51e3,0xe7bb,0x2349,
 216 0xbfab,0x1bbc,0x537c,0x9ebc,
 217 0x3fba,0x46da,0xd536,0xf53c,
 218 0xbfc6,0x94d1,0x0469,0x192e,
 219 0x3fd0,0x2a63,0x724a,0x7ffa
 220 };
 221 #endif
 222 \f
 223 /*                                                      i1.c    */
 224
 225 /* Chebyshev coefficients for exp(-x) sqrt(x) I1(x)
 226  * in the inverted interval [8,infinity].
 227  *
 228  * lim(x->inf){ exp(-x) sqrt(x) I1(x) } = 1/sqrt(2pi).
 229  */
 230
 231 #ifdef UNK
 232 static double B[] =
 233 {
 234  7.51729631084210481353E-18,
 235  4.41434832307170791151E-18,
 236 -4.65030536848935832153E-17,
 237 -3.20952592199342395980E-17,
 238  2.96262899764595013876E-16,
 239  3.30820231092092828324E-16,
 240 -1.88035477551078244854E-15,
 241 -3.81440307243700780478E-15,
 242  1.04202769841288027642E-14,
 243  4.27244001671195135429E-14,
 244 -2.10154184277266431302E-14,
 245 -4.08355111109219731823E-13,
 246 -7.19855177624590851209E-13,
 247  2.03562854414708950722E-12,
 248  1.41258074366137813316E-11,
 249  3.25260358301548823856E-11,
 250 -1.89749581235054123450E-11,
 251 -5.58974346219658380687E-10,
 252 -3.83538038596423702205E-9,
 253 -2.63146884688951950684E-8,
 254 -2.51223623787020892529E-7,
 255 -3.88256480887769039346E-6,
 256 -1.10588938762623716291E-4,
 257 -9.76109749136146840777E-3,
 258  7.78576235018280120474E-1
 259 };
 260 #endif
 261
 262 #ifdef DEC
 263 static unsigned short B[] = {
 264 0022012,0125555,0115227,0043456,
 265 0021642,0156127,0052075,0145203,
 266 0122526,0072435,0111231,0011664,
 267 0122424,0001544,0161671,0114403,
 268 0023252,0144257,0163532,0142121,
 269 0023276,0132162,0174045,0013204,
 270 0124007,0077154,0057046,0110517,
 271 0124211,0066650,0116127,0157073,
 272 0024473,0133413,0130551,0107504,
 273 0025100,0064741,0032631,0040364,
 274 0124675,0045101,0071551,0012400,
 275 0125745,0161054,0071637,0011247,
 276 0126112,0117410,0035525,0122231,
 277 0026417,0037237,0131034,0176427,
 278 0027170,0100373,0024742,0025725,
 279 0027417,0006417,0105303,0141446,
 280 0127246,0163716,0121202,0060137,
 281 0130431,0123122,0120436,0166000,
 282 0131203,0144134,0153251,0124500,
 283 0131742,0005234,0122732,0033006,
 284 0132606,0157751,0072362,0121031,
 285 0133602,0043372,0047120,0015626,
 286 0134747,0165774,0001125,0046462,
 287 0136437,0166402,0117746,0155137,
 288 0040107,0050305,0125330,0124241
 289 };
 290 #endif
 291
 292 #ifdef IBMPC
 293 static unsigned short B[] = {
 294 0xe8e6,0xb352,0x556d,0x3c61,
 295 0xb950,0xea87,0x5b8a,0x3c54,
 296 0x2277,0xb253,0xcea3,0xbc8a,
 297 0x3320,0x9c77,0x806c,0xbc82,
 298 0x588a,0xfceb,0x5915,0x3cb5,
 299 0xa2d1,0x5f04,0xd68e,0x3cb7,
 300 0xd22a,0x8bc4,0xefcd,0xbce0,
 301 0xfbc7,0x138a,0x2db5,0xbcf1,
 302 0x31e8,0x762d,0x76e1,0x3d07,
 303 0x281e,0x26b3,0x0d3c,0x3d28,
 304 0x22a0,0x2e6d,0xa948,0xbd17,
 305 0xe255,0x8e73,0xbc45,0xbd5c,
 306 0xb493,0x076a,0x53e1,0xbd69,
 307 0x9fa3,0xf643,0xe7d3,0x3d81,
 308 0x457b,0x653c,0x101f,0x3daf,
 309 0x7865,0xf158,0xe1a1,0x3dc1,
 310 0x4c0c,0xd450,0xdcf9,0xbdb4,
 311 0xdd80,0x5423,0x34ca,0xbe03,
 312 0x3528,0x9ad5,0x790b,0xbe30,
 313 0x46c1,0x94bb,0x4153,0xbe5c,
 314 0x5443,0x2e9e,0xdbfd,0xbe90,
 315 0x0373,0x49ca,0x48df,0xbed0,
 316 0xa9a6,0x804a,0xfd7f,0xbf1c,
 317 0xdb4c,0x53fc,0xfda0,0xbf83,
 318 0x1514,0xb55b,0xea18,0x3fe8
 319 };
 320 #endif
 321
 322 #ifdef MIEEE
 323 static unsigned short B[] = {
 324 0x3c61,0x556d,0xb352,0xe8e6,
 325 0x3c54,0x5b8a,0xea87,0xb950,
 326 0xbc8a,0xcea3,0xb253,0x2277,
 327 0xbc82,0x806c,0x9c77,0x3320,
 328 0x3cb5,0x5915,0xfceb,0x588a,
 329 0x3cb7,0xd68e,0x5f04,0xa2d1,
 330 0xbce0,0xefcd,0x8bc4,0xd22a,
 331 0xbcf1,0x2db5,0x138a,0xfbc7,
 332 0x3d07,0x76e1,0x762d,0x31e8,
 333 0x3d28,0x0d3c,0x26b3,0x281e,
 334 0xbd17,0xa948,0x2e6d,0x22a0,
 335 0xbd5c,0xbc45,0x8e73,0xe255,
 336 0xbd69,0x53e1,0x076a,0xb493,
 337 0x3d81,0xe7d3,0xf643,0x9fa3,
 338 0x3daf,0x101f,0x653c,0x457b,
 339 0x3dc1,0xe1a1,0xf158,0x7865,
 340 0xbdb4,0xdcf9,0xd450,0x4c0c,
 341 0xbe03,0x34ca,0x5423,0xdd80,
 342 0xbe30,0x790b,0x9ad5,0x3528,
 343 0xbe5c,0x4153,0x94bb,0x46c1,
 344 0xbe90,0xdbfd,0x2e9e,0x5443,
 345 0xbed0,0x48df,0x49ca,0x0373,
 346 0xbf1c,0xfd7f,0x804a,0xa9a6,
 347 0xbf83,0xfda0,0x53fc,0xdb4c,
 348 0x3fe8,0xea18,0xb55b,0x1514
 349 };
 350 #endif
 351 \f
 352 /*                                                      i1.c    */
 353 #ifdef ANSIPROT
 354 extern double chbevl ( double, void *, int );
 355 extern double exp ( double );
 356 extern double sqrt ( double );
 357 extern double fabs ( double );
 358 #else
 359 double chbevl(), exp(), sqrt(), fabs();
 360 #endif
 361
 362 double i1(x)
 363 double x;
 364 {
 365 double y, z;
 366
 367 z = fabs(x);
 368 if( z <= 8.0 )
 369         {
 370         y = (z/2.0) - 2.0;
 371         z = chbevl( y, A, 29 ) * z * exp(z);
 372         }
 373 else
 374         {
 375         z = exp(z) * chbevl( 32.0/z - 2.0, B, 25 ) / sqrt(z);
 376         }
 377 if( x < 0.0 )
 378         z = -z;
 379 return( z );
 380 }
 381 \f
 382 /*                                                      i1e()   */
 383
 384 double i1e( x )
 385 double x;
 386 {
 387 double y, z;
 388
 389 z = fabs(x);
 390 if( z <= 8.0 )
 391         {
 392         y = (z/2.0) - 2.0;
 393         z = chbevl( y, A, 29 ) * z;
 394         }
 395 else
 396         {
 397         z = chbevl( 32.0/z - 2.0, B, 25 ) / sqrt(z);
 398         }
 399 if( x < 0.0 )
 400         z = -z;
 401 return( z );
 402 }