lib/libm/j1.c

   1 /*                                                      j1.c
   2  *
   3  *      Bessel function of order one
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * double x, y, j1();
  10  *
  11  * y = j1( x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Returns Bessel function of order one of the argument.
  18  *
  19  * The domain is divided into the intervals [0, 8] and
  20  * (8, infinity). In the first interval a 24 term Chebyshev
  21  * expansion is used. In the second, the asymptotic
  22  * trigonometric representation is employed using two
  23  * rational functions of degree 5/5.
  24  *
  25  *
  26  *
  27  * ACCURACY:
  28  *
  29  *                      Absolute error:
  30  * arithmetic   domain      # trials      peak         rms
  31  *    DEC       0, 30       10000       4.0e-17     1.1e-17
  32  *    IEEE      0, 30       30000       2.6e-16     1.1e-16
  33  *
  34  *
  35  */
  36 \f/*                                                     y1.c
  37  *
  38  *      Bessel function of second kind of order one
  39  *
  40  *
  41  *
  42  * SYNOPSIS:
  43  *
  44  * double x, y, y1();
  45  *
  46  * y = y1( x );
  47  *
  48  *
  49  *
  50  * DESCRIPTION:
  51  *
  52  * Returns Bessel function of the second kind of order one
  53  * of the argument.
  54  *
  55  * The domain is divided into the intervals [0, 8] and
  56  * (8, infinity). In the first interval a 25 term Chebyshev
  57  * expansion is used, and a call to j1() is required.
  58  * In the second, the asymptotic trigonometric representation
  59  * is employed using two rational functions of degree 5/5.
  60  *
  61  *
  62  *
  63  * ACCURACY:
  64  *
  65  *                      Absolute error:
  66  * arithmetic   domain      # trials      peak         rms
  67  *    DEC       0, 30       10000       8.6e-17     1.3e-17
  68  *    IEEE      0, 30       30000       1.0e-15     1.3e-16
  69  *
  70  * (error criterion relative when |y1| > 1).
  71  *
  72  */
  73 \f
  74
  75 /*
  76 Cephes Math Library Release 2.8:  June, 2000
  77 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
  78 */
  79
  80 /*
  81 #define PIO4 .78539816339744830962
  82 #define THPIO4 2.35619449019234492885
  83 #define SQ2OPI .79788456080286535588
  84 */
  85
  86 #include "mconf.h"
  87
  88 #ifdef UNK
  89 static double RP[4] = {
  90 -8.99971225705559398224E8,
  91  4.52228297998194034323E11,
  92 -7.27494245221818276015E13,
  93  3.68295732863852883286E15,
  94 };
  95 static double RQ[8] = {
  96 /* 1.00000000000000000000E0,*/
  97  6.20836478118054335476E2,
  98  2.56987256757748830383E5,
  99  8.35146791431949253037E7,
 100  2.21511595479792499675E10,
 101  4.74914122079991414898E12,
 102  7.84369607876235854894E14,
 103  8.95222336184627338078E16,
 104  5.32278620332680085395E18,
 105 };
 106 #endif
 107 #ifdef DEC
 108 static unsigned short RP[16] = {
 109 0147526,0110742,0063322,0077052,
 110 0051722,0112720,0065034,0061530,
 111 0153604,0052227,0033147,0105650,
 112 0055121,0055025,0032276,0022015,
 113 };
 114 static unsigned short RQ[32] = {
 115 /*0040200,0000000,0000000,0000000,*/
 116 0042433,0032610,0155604,0033473,
 117 0044572,0173320,0067270,0006616,
 118 0046637,0045246,0162225,0006606,
 119 0050645,0004773,0157577,0053004,
 120 0052612,0033734,0001667,0176501,
 121 0054462,0054121,0173147,0121367,
 122 0056237,0002777,0121451,0176007,
 123 0057623,0136253,0131601,0044710,
 124 };
 125 #endif
 126 #ifdef IBMPC
 127 static unsigned short RP[16] = {
 128 0x4fc5,0x4cda,0xd23c,0xc1ca,
 129 0x8c6b,0x0d43,0x52ba,0x425a,
 130 0xf175,0xe6cc,0x8a92,0xc2d0,
 131 0xc482,0xa697,0x2b42,0x432a,
 132 };
 133 static unsigned short RQ[32] = {
 134 /*0x0000,0x0000,0x0000,0x3ff0,*/
 135 0x86e7,0x1b70,0x66b1,0x4083,
 136 0x01b2,0x0dd7,0x5eda,0x410f,
 137 0xa1b1,0xdc92,0xe954,0x4193,
 138 0xeac1,0x7bef,0xa13f,0x4214,
 139 0xffa8,0x8076,0x46fb,0x4291,
 140 0xf45f,0x3ecc,0x4b0a,0x4306,
 141 0x3f81,0xf465,0xe0bf,0x4373,
 142 0x2939,0x7670,0x7795,0x43d2,
 143 };
 144 #endif
 145 #ifdef MIEEE
 146 static unsigned short RP[16] = {
 147 0xc1ca,0xd23c,0x4cda,0x4fc5,
 148 0x425a,0x52ba,0x0d43,0x8c6b,
 149 0xc2d0,0x8a92,0xe6cc,0xf175,
 150 0x432a,0x2b42,0xa697,0xc482,
 151 };
 152 static unsigned short RQ[32] = {
 153 /*0x3ff0,0x0000,0x0000,0x0000,*/
 154 0x4083,0x66b1,0x1b70,0x86e7,
 155 0x410f,0x5eda,0x0dd7,0x01b2,
 156 0x4193,0xe954,0xdc92,0xa1b1,
 157 0x4214,0xa13f,0x7bef,0xeac1,
 158 0x4291,0x46fb,0x8076,0xffa8,
 159 0x4306,0x4b0a,0x3ecc,0xf45f,
 160 0x4373,0xe0bf,0xf465,0x3f81,
 161 0x43d2,0x7795,0x7670,0x2939,
 162 };
 163 #endif
 164
 165 #ifdef UNK
 166 static double PP[7] = {
 167  7.62125616208173112003E-4,
 168  7.31397056940917570436E-2,
 169  1.12719608129684925192E0,
 170  5.11207951146807644818E0,
 171  8.42404590141772420927E0,
 172  5.21451598682361504063E0,
 173  1.00000000000000000254E0,
 174 };
 175 static double PQ[7] = {
 176  5.71323128072548699714E-4,
 177  6.88455908754495404082E-2,
 178  1.10514232634061696926E0,
 179  5.07386386128601488557E0,
 180  8.39985554327604159757E0,
 181  5.20982848682361821619E0,
 182  9.99999999999999997461E-1,
 183 };
 184 #endif
 185 #ifdef DEC
 186 static unsigned short PP[28] = {
 187 0035507,0144542,0061543,0024326,
 188 0037225,0145105,0017766,0022661,
 189 0040220,0043766,0010254,0133255,
 190 0040643,0113047,0142611,0151521,
 191 0041006,0144344,0055351,0074261,
 192 0040646,0156520,0120574,0006416,
 193 0040200,0000000,0000000,0000000,
 194 };
 195 static unsigned short PQ[28] = {
 196 0035425,0142330,0115041,0165514,
 197 0037214,0177352,0145105,0052026,
 198 0040215,0072515,0141207,0073255,
 199 0040642,0056427,0137222,0106405,
 200 0041006,0062716,0166427,0165450,
 201 0040646,0133352,0035425,0123304,
 202 0040200,0000000,0000000,0000000,
 203 };
 204 #endif
 205 #ifdef IBMPC
 206 static unsigned short PP[28] = {
 207 0x651b,0x4c6c,0xf92c,0x3f48,
 208 0xc4b6,0xa3fe,0xb948,0x3fb2,
 209 0x96d6,0xc215,0x08fe,0x3ff2,
 210 0x3a6a,0xf8b1,0x72c4,0x4014,
 211 0x2f16,0x8b5d,0xd91c,0x4020,
 212 0x81a2,0x142f,0xdbaa,0x4014,
 213 0x0000,0x0000,0x0000,0x3ff0,
 214 };
 215 static unsigned short PQ[28] = {
 216 0x3d69,0x1344,0xb89b,0x3f42,
 217 0xaa83,0x5948,0x9fdd,0x3fb1,
 218 0xeed6,0xb850,0xaea9,0x3ff1,
 219 0x51a1,0xf7d2,0x4ba2,0x4014,
 220 0xfd65,0xdda2,0xccb9,0x4020,
 221 0xb4d9,0x4762,0xd6dd,0x4014,
 222 0x0000,0x0000,0x0000,0x3ff0,
 223 };
 224 #endif
 225 #ifdef MIEEE
 226 static unsigned short PP[28] = {
 227 0x3f48,0xf92c,0x4c6c,0x651b,
 228 0x3fb2,0xb948,0xa3fe,0xc4b6,
 229 0x3ff2,0x08fe,0xc215,0x96d6,
 230 0x4014,0x72c4,0xf8b1,0x3a6a,
 231 0x4020,0xd91c,0x8b5d,0x2f16,
 232 0x4014,0xdbaa,0x142f,0x81a2,
 233 0x3ff0,0x0000,0x0000,0x0000,
 234 };
 235 static unsigned short PQ[28] = {
 236 0x3f42,0xb89b,0x1344,0x3d69,
 237 0x3fb1,0x9fdd,0x5948,0xaa83,
 238 0x3ff1,0xaea9,0xb850,0xeed6,
 239 0x4014,0x4ba2,0xf7d2,0x51a1,
 240 0x4020,0xccb9,0xdda2,0xfd65,
 241 0x4014,0xd6dd,0x4762,0xb4d9,
 242 0x3ff0,0x0000,0x0000,0x0000,
 243 };
 244 #endif
 245
 246 #ifdef UNK
 247 static double QP[8] = {
 248  5.10862594750176621635E-2,
 249  4.98213872951233449420E0,
 250  7.58238284132545283818E1,
 251  3.66779609360150777800E2,
 252  7.10856304998926107277E2,
 253  5.97489612400613639965E2,
 254  2.11688757100572135698E2,
 255  2.52070205858023719784E1,
 256 };
 257 static double QQ[7] = {
 258 /* 1.00000000000000000000E0,*/
 259  7.42373277035675149943E1,
 260  1.05644886038262816351E3,
 261  4.98641058337653607651E3,
 262  9.56231892404756170795E3,
 263  7.99704160447350683650E3,
 264  2.82619278517639096600E3,
 265  3.36093607810698293419E2,
 266 };
 267 #endif
 268 #ifdef DEC
 269 static unsigned short QP[32] = {
 270 0037121,0037723,0055605,0151004,
 271 0040637,0066656,0031554,0077264,
 272 0041627,0122714,0153170,0161466,
 273 0042267,0061712,0036520,0140145,
 274 0042461,0133315,0131573,0071176,
 275 0042425,0057525,0147500,0013201,
 276 0042123,0130122,0061245,0154131,
 277 0041311,0123772,0064254,0172650,
 278 };
 279 static unsigned short QQ[28] = {
 280 /*0040200,0000000,0000000,0000000,*/
 281 0041624,0074603,0002112,0101670,
 282 0042604,0007135,0010162,0175565,
 283 0043233,0151510,0157757,0172010,
 284 0043425,0064506,0112006,0104276,
 285 0043371,0164125,0032271,0164242,
 286 0043060,0121425,0122750,0136013,
 287 0042250,0005773,0053472,0146267,
 288 };
 289 #endif
 290 #ifdef IBMPC
 291 static unsigned short QP[32] = {
 292 0xba40,0x6b70,0x27fa,0x3faa,
 293 0x8fd6,0xc66d,0xedb5,0x4013,
 294 0x1c67,0x9acf,0xf4b9,0x4052,
 295 0x180d,0x47aa,0xec79,0x4076,
 296 0x6e50,0xb66f,0x36d9,0x4086,
 297 0x02d0,0xb9e8,0xabea,0x4082,
 298 0xbb0b,0x4c54,0x760a,0x406a,
 299 0x9eb5,0x4d15,0x34ff,0x4039,
 300 };
 301 static unsigned short QQ[28] = {
 302 /*0x0000,0x0000,0x0000,0x3ff0,*/
 303 0x5077,0x6089,0x8f30,0x4052,
 304 0x5f6f,0xa20e,0x81cb,0x4090,
 305 0xfe81,0x1bfd,0x7a69,0x40b3,
 306 0xd118,0xd280,0xad28,0x40c2,
 307 0x3d14,0xa697,0x3d0a,0x40bf,
 308 0x1781,0xb4bd,0x1462,0x40a6,
 309 0x5997,0x6ae7,0x017f,0x4075,
 310 };
 311 #endif
 312 #ifdef MIEEE
 313 static unsigned short QP[32] = {
 314 0x3faa,0x27fa,0x6b70,0xba40,
 315 0x4013,0xedb5,0xc66d,0x8fd6,
 316 0x4052,0xf4b9,0x9acf,0x1c67,
 317 0x4076,0xec79,0x47aa,0x180d,
 318 0x4086,0x36d9,0xb66f,0x6e50,
 319 0x4082,0xabea,0xb9e8,0x02d0,
 320 0x406a,0x760a,0x4c54,0xbb0b,
 321 0x4039,0x34ff,0x4d15,0x9eb5,
 322 };
 323 static unsigned short QQ[28] = {
 324 /*0x3ff0,0x0000,0x0000,0x0000,*/
 325 0x4052,0x8f30,0x6089,0x5077,
 326 0x4090,0x81cb,0xa20e,0x5f6f,
 327 0x40b3,0x7a69,0x1bfd,0xfe81,
 328 0x40c2,0xad28,0xd280,0xd118,
 329 0x40bf,0x3d0a,0xa697,0x3d14,
 330 0x40a6,0x1462,0xb4bd,0x1781,
 331 0x4075,0x017f,0x6ae7,0x5997,
 332 };
 333 #endif
 334
 335 #ifdef UNK
 336 static double YP[6] = {
 337  1.26320474790178026440E9,
 338 -6.47355876379160291031E11,
 339  1.14509511541823727583E14,
 340 -8.12770255501325109621E15,
 341  2.02439475713594898196E17,
 342 -7.78877196265950026825E17,
 343 };
 344 static double YQ[8] = {
 345 /* 1.00000000000000000000E0,*/
 346  5.94301592346128195359E2,
 347  2.35564092943068577943E5,
 348  7.34811944459721705660E7,
 349  1.87601316108706159478E10,
 350  3.88231277496238566008E12,
 351  6.20557727146953693363E14,
 352  6.87141087355300489866E16,
 353  3.97270608116560655612E18,
 354 };
 355 #endif
 356 #ifdef DEC
 357 static unsigned short YP[24] = {
 358 0047626,0112763,0013715,0133045,
 359 0152026,0134552,0142033,0024411,
 360 0053720,0045245,0102210,0077565,
 361 0155347,0000321,0136415,0102031,
 362 0056463,0146550,0055633,0032605,
 363 0157054,0171012,0167361,0054265,
 364 };
 365 static unsigned short YQ[32] = {
 366 /*0040200,0000000,0000000,0000000,*/
 367 0042424,0111515,0044773,0153014,
 368 0044546,0005405,0171307,0075774,
 369 0046614,0023575,0047105,0063556,
 370 0050613,0143034,0101533,0156026,
 371 0052541,0175367,0166514,0114257,
 372 0054415,0014466,0134350,0171154,
 373 0056164,0017436,0025075,0022101,
 374 0057534,0103614,0103663,0121772,
 375 };
 376 #endif
 377 #ifdef IBMPC
 378 static unsigned short YP[24] = {
 379 0xb6c5,0x62f9,0xd2be,0x41d2,
 380 0x6521,0x5883,0xd72d,0xc262,
 381 0x0fef,0xb091,0x0954,0x42da,
 382 0xb083,0x37a1,0xe01a,0xc33c,
 383 0x66b1,0x0b73,0x79ad,0x4386,
 384 0x2b17,0x5dde,0x9e41,0xc3a5,
 385 };
 386 static unsigned short YQ[32] = {
 387 /*0x0000,0x0000,0x0000,0x3ff0,*/
 388 0x7ac2,0xa93f,0x9269,0x4082,
 389 0xef7f,0xbe58,0xc160,0x410c,
 390 0xacee,0xa9c8,0x84ef,0x4191,
 391 0x7b83,0x906b,0x78c3,0x4211,
 392 0x9316,0xfda9,0x3f5e,0x428c,
 393 0x1e4e,0xd71d,0xa326,0x4301,
 394 0xa488,0xc547,0x83e3,0x436e,
 395 0x747f,0x90f6,0x90f1,0x43cb,
 396 };
 397 #endif
 398 #ifdef MIEEE
 399 static unsigned short YP[24] = {
 400 0x41d2,0xd2be,0x62f9,0xb6c5,
 401 0xc262,0xd72d,0x5883,0x6521,
 402 0x42da,0x0954,0xb091,0x0fef,
 403 0xc33c,0xe01a,0x37a1,0xb083,
 404 0x4386,0x79ad,0x0b73,0x66b1,
 405 0xc3a5,0x9e41,0x5dde,0x2b17,
 406 };
 407 static unsigned short YQ[32] = {
 408 /*0x3ff0,0x0000,0x0000,0x0000,*/
 409 0x4082,0x9269,0xa93f,0x7ac2,
 410 0x410c,0xc160,0xbe58,0xef7f,
 411 0x4191,0x84ef,0xa9c8,0xacee,
 412 0x4211,0x78c3,0x906b,0x7b83,
 413 0x428c,0x3f5e,0xfda9,0x9316,
 414 0x4301,0xa326,0xd71d,0x1e4e,
 415 0x436e,0x83e3,0xc547,0xa488,
 416 0x43cb,0x90f1,0x90f6,0x747f,
 417 };
 418 #endif
 419
 420
 421 #ifdef UNK
 422 static double Z1 = 1.46819706421238932572E1;
 423 static double Z2 = 4.92184563216946036703E1;
 424 #endif
 425
 426 #ifdef DEC
 427 static unsigned short DZ1[] = {0041152,0164532,0006114,0010540};
 428 static unsigned short DZ2[] = {0041504,0157663,0001625,0020621};
 429 #define Z1 (*(double *)DZ1)
 430 #define Z2 (*(double *)DZ2)
 431 #endif
 432
 433 #ifdef IBMPC
 434 static unsigned short DZ1[] = {0x822c,0x4189,0x5d2b,0x402d};
 435 static unsigned short DZ2[] = {0xa432,0x6072,0x9bf6,0x4048};
 436 #define Z1 (*(double *)DZ1)
 437 #define Z2 (*(double *)DZ2)
 438 #endif
 439
 440 #ifdef MIEEE
 441 static unsigned short DZ1[] = {0x402d,0x5d2b,0x4189,0x822c};
 442 static unsigned short DZ2[] = {0x4048,0x9bf6,0x6072,0xa432};
 443 #define Z1 (*(double *)DZ1)
 444 #define Z2 (*(double *)DZ2)
 445 #endif
 446
 447 #ifdef ANSIPROT
 448 extern double polevl ( double, void *, int );
 449 extern double p1evl ( double, void *, int );
 450 extern double log ( double );
 451 extern double sin ( double );
 452 extern double cos ( double );
 453 extern double sqrt ( double );
 454 double j1 ( double );
 455 #else
 456 double polevl(), p1evl(), log(), sin(), cos(), sqrt();
 457 double j1();
 458 #endif
 459 extern double TWOOPI, THPIO4, SQ2OPI;
 460
 461 double j1(x)
 462 double x;
 463 {
 464 double w, z, p, q, xn;
 465
 466 w = x;
 467 if( x < 0 )
 468         w = -x;
 469
 470 if( w <= 5.0 )
 471         {
 472         z = x * x;
 473         w = polevl( z, RP, 3 ) / p1evl( z, RQ, 8 );
 474         w = w * x * (z - Z1) * (z - Z2);
 475         return( w );
 476         }
 477
 478 w = 5.0/x;
 479 z = w * w;
 480 p = polevl( z, PP, 6)/polevl( z, PQ, 6 );
 481 q = polevl( z, QP, 7)/p1evl( z, QQ, 7 );
 482 xn = x - THPIO4;
 483 p = p * cos(xn) - w * q * sin(xn);
 484 return( p * SQ2OPI / sqrt(x) );
 485 }
 486
 487
 488 extern double MAXNUM;
 489
 490 double y1(x)
 491 double x;
 492 {
 493 double w, z, p, q, xn;
 494
 495 if( x <= 5.0 )
 496         {
 497         if( x <= 0.0 )
 498                 {
 499                 mtherr( "y1", DOMAIN );
 500                 return( -MAXNUM );
 501                 }
 502         z = x * x;
 503         w = x * (polevl( z, YP, 5 ) / p1evl( z, YQ, 8 ));
 504         w += TWOOPI * ( j1(x) * log(x)  -  1.0/x );
 505         return( w );
 506         }
 507
 508 w = 5.0/x;
 509 z = w * w;
 510 p = polevl( z, PP, 6)/polevl( z, PQ, 6 );
 511 q = polevl( z, QP, 7)/p1evl( z, QQ, 7 );
 512 xn = x - THPIO4;
 513 p = p * sin(xn) + w * q * cos(xn);
 514 return( p * SQ2OPI / sqrt(x) );
 515 }