383 kpcos(double x) {
384 double z;
385
386 z = x * x;
387 return (kc[0] + z * (kc[1] + z * kc[2] + (z * z) * kc[3]));
388 }
389
390 /* INDENT OFF */
391 static const double
392 t0z1 = 0.134861805732790769689793935774652917006,
393 t0z2 = 0.461632144968362341262659542325721328468,
394 t0z3 = 0.819773101100500601787868704921606996312;
395 /* 1.134861805732790769689793935774652917006 */
396 /* INDENT ON */
397
398 /*
399 * gamma(x+i) for 0 <= x < 1
400 */
401 static double
402 gam_n(int i, double x) {
403 double rr, yy;
404 double z1, z2;
405
406 /* compute yy = gamma(x+1) */
407 if (x > 0.2845) {
408 if (x > 0.6374)
409 yy = GT3(x - t0z3);
410 else
411 yy = GT2(x - t0z2);
412 } else
413 yy = GT1(x - t0z1);
414
415 /* compute gamma(x+i) = (x+i-1)*...*(x+1)*yy, 0<i<8 */
416 switch (i) {
417 case 0: /* yy/x */
418 rr = yy / x;
419 break;
420 case 1: /* yy */
421 rr = yy;
422 break;
423 case 2: /* (x+1)*yy */
455 double ss, ww;
456 double x, y, z;
457 int i, j, k, ix, hx, xk;
458
459 hx = *(int *) &xf;
460 ix = hx & 0x7fffffff;
461
462 x = (double) xf;
463 if (ix < 0x33800000)
464 return (1.0F / xf); /* |x| < 2**-24 */
465
466 if (ix >= 0x7f800000)
467 return (xf * ((hx < 0)? 0.0F : xf)); /* +-Inf or NaN */
468
469 if (hx > 0x420C290F) /* x > 35.040096283... overflow */
470 return (float)(x / tiny);
471
472 if (hx >= 0x41000000) /* x >= 8 */
473 return ((float) large_gam(x));
474
475 if (hx > 0) { /* x from 0 to 8 */
476 i = (int) xf;
477 return ((float) gam_n(i, x - (double) i));
478 }
479
480 /* negative x */
481 /* INDENT OFF */
482 /*
483 * compute xk =
484 * -2 ... x is an even int (-inf is considered even)
485 * -1 ... x is an odd int
486 * +0 ... x is not an int but chopped to an even int
487 * +1 ... x is not an int but chopped to an odd int
488 */
489 /* INDENT ON */
490 xk = 0;
491 if (ix >= 0x4b000000) {
492 if (ix > 0x4b000000)
493 xk = -2;
494 else
495 xk = -2 + (ix & 1);
|
383 kpcos(double x) {
384 double z;
385
386 z = x * x;
387 return (kc[0] + z * (kc[1] + z * kc[2] + (z * z) * kc[3]));
388 }
389
390 /* INDENT OFF */
391 static const double
392 t0z1 = 0.134861805732790769689793935774652917006,
393 t0z2 = 0.461632144968362341262659542325721328468,
394 t0z3 = 0.819773101100500601787868704921606996312;
395 /* 1.134861805732790769689793935774652917006 */
396 /* INDENT ON */
397
398 /*
399 * gamma(x+i) for 0 <= x < 1
400 */
401 static double
402 gam_n(int i, double x) {
403 double rr = 0.0L, yy;
404 double z1, z2;
405
406 /* compute yy = gamma(x+1) */
407 if (x > 0.2845) {
408 if (x > 0.6374)
409 yy = GT3(x - t0z3);
410 else
411 yy = GT2(x - t0z2);
412 } else
413 yy = GT1(x - t0z1);
414
415 /* compute gamma(x+i) = (x+i-1)*...*(x+1)*yy, 0<i<8 */
416 switch (i) {
417 case 0: /* yy/x */
418 rr = yy / x;
419 break;
420 case 1: /* yy */
421 rr = yy;
422 break;
423 case 2: /* (x+1)*yy */
455 double ss, ww;
456 double x, y, z;
457 int i, j, k, ix, hx, xk;
458
459 hx = *(int *) &xf;
460 ix = hx & 0x7fffffff;
461
462 x = (double) xf;
463 if (ix < 0x33800000)
464 return (1.0F / xf); /* |x| < 2**-24 */
465
466 if (ix >= 0x7f800000)
467 return (xf * ((hx < 0)? 0.0F : xf)); /* +-Inf or NaN */
468
469 if (hx > 0x420C290F) /* x > 35.040096283... overflow */
470 return (float)(x / tiny);
471
472 if (hx >= 0x41000000) /* x >= 8 */
473 return ((float) large_gam(x));
474
475 if (hx > 0) { /* 0 < x < 8 */
476 i = (int) xf;
477 return ((float) gam_n(i, x - (double) i));
478 }
479
480 /* negative x */
481 /* INDENT OFF */
482 /*
483 * compute xk =
484 * -2 ... x is an even int (-inf is considered even)
485 * -1 ... x is an odd int
486 * +0 ... x is not an int but chopped to an even int
487 * +1 ... x is not an int but chopped to an odd int
488 */
489 /* INDENT ON */
490 xk = 0;
491 if (ix >= 0x4b000000) {
492 if (ix > 0x4b000000)
493 xk = -2;
494 else
495 xk = -2 + (ix & 1);
|