1389 return (xx);
1390 }
1391
1392 /* INDENT OFF */
1393 static const double
1394 /* 0.134861805732790769689793935774652917006 */
1395 t0z1 = 0.1348618057327907737708,
1396 t0z1_l = -4.0810077708578299022531e-18,
1397 /* 0.461632144968362341262659542325721328468 */
1398 t0z2 = 0.4616321449683623567850,
1399 t0z2_l = -1.5522348162858676890521e-17,
1400 /* 0.819773101100500601787868704921606996312 */
1401 t0z3 = 0.8197731011005006118708,
1402 t0z3_l = -1.0082945122487103498325e-17;
1403 /* 1.134861805732790769689793935774652917006 */
1404 /* INDENT ON */
1405
1406 /* gamma(x+i) for 0 <= x < 1 */
1407 static struct Double
1408 gam_n(int i, double x) {
1409 struct Double rr, yy;
1410 double r1, r2, t2, z, xh, xl, yh, yl, zh, z1, z2, zl, x5, wh, wl;
1411
1412 /* compute yy = gamma(x+1) */
1413 if (x > 0.2845) {
1414 if (x > 0.6374) {
1415 r1 = x - t0z3;
1416 r2 = (double) ((float) (r1 - t0z3_l));
1417 t2 = r1 - r2;
1418 yy = GT3(r2, t2 - t0z3_l);
1419 } else {
1420 r1 = x - t0z2;
1421 r2 = (double) ((float) (r1 - t0z2_l));
1422 t2 = r1 - r2;
1423 yy = GT2(r2, t2 - t0z2_l);
1424 }
1425 } else {
1426 r1 = x - t0z1;
1427 r2 = (double) ((float) (r1 - t0z1_l));
1428 t2 = r1 - r2;
1429 yy = GT1(r2, t2 - t0z1_l);
1541 lx = __LO(x);
1542 ix = hx & 0x7fffffff;
1543 y = x;
1544
1545 if (ix < 0x3ca00000)
1546 return (one / x); /* |x| < 2**-53 */
1547 if (ix >= 0x7ff00000)
1548 /* +Inf -> +Inf, -Inf or NaN -> NaN */
1549 return (x * ((hx < 0)? 0.0 : x));
1550 if (hx > 0x406573fa || /* x > 171.62... overflow to +inf */
1551 (hx == 0x406573fa && lx > 0xE561F647)) {
1552 z = x / tiny;
1553 return (z * z);
1554 }
1555 if (hx >= 0x40200000) { /* x >= 8 */
1556 ww = large_gam(x, &m);
1557 w = ww.h + ww.l;
1558 __HI(w) += m << 20;
1559 return (w);
1560 }
1561 if (hx > 0) { /* x from 0 to 8 */
1562 i = (int) x;
1563 ww = gam_n(i, x - (double) i);
1564 return (ww.h + ww.l);
1565 }
1566
1567 /* negative x */
1568 /* INDENT OFF */
1569 /*
1570 * compute: xk =
1571 * -2 ... x is an even int (-inf is even)
1572 * -1 ... x is an odd int
1573 * +0 ... x is not an int but chopped to an even int
1574 * +1 ... x is not an int but chopped to an odd int
1575 */
1576 /* INDENT ON */
1577 xk = 0;
1578 if (ix >= 0x43300000) {
1579 if (ix >= 0x43400000)
1580 xk = -2;
1581 else
|
1389 return (xx);
1390 }
1391
1392 /* INDENT OFF */
1393 static const double
1394 /* 0.134861805732790769689793935774652917006 */
1395 t0z1 = 0.1348618057327907737708,
1396 t0z1_l = -4.0810077708578299022531e-18,
1397 /* 0.461632144968362341262659542325721328468 */
1398 t0z2 = 0.4616321449683623567850,
1399 t0z2_l = -1.5522348162858676890521e-17,
1400 /* 0.819773101100500601787868704921606996312 */
1401 t0z3 = 0.8197731011005006118708,
1402 t0z3_l = -1.0082945122487103498325e-17;
1403 /* 1.134861805732790769689793935774652917006 */
1404 /* INDENT ON */
1405
1406 /* gamma(x+i) for 0 <= x < 1 */
1407 static struct Double
1408 gam_n(int i, double x) {
1409 struct Double rr = {0.0L, 0.0L}, yy;
1410 double r1, r2, t2, z, xh, xl, yh, yl, zh, z1, z2, zl, x5, wh, wl;
1411
1412 /* compute yy = gamma(x+1) */
1413 if (x > 0.2845) {
1414 if (x > 0.6374) {
1415 r1 = x - t0z3;
1416 r2 = (double) ((float) (r1 - t0z3_l));
1417 t2 = r1 - r2;
1418 yy = GT3(r2, t2 - t0z3_l);
1419 } else {
1420 r1 = x - t0z2;
1421 r2 = (double) ((float) (r1 - t0z2_l));
1422 t2 = r1 - r2;
1423 yy = GT2(r2, t2 - t0z2_l);
1424 }
1425 } else {
1426 r1 = x - t0z1;
1427 r2 = (double) ((float) (r1 - t0z1_l));
1428 t2 = r1 - r2;
1429 yy = GT1(r2, t2 - t0z1_l);
1541 lx = __LO(x);
1542 ix = hx & 0x7fffffff;
1543 y = x;
1544
1545 if (ix < 0x3ca00000)
1546 return (one / x); /* |x| < 2**-53 */
1547 if (ix >= 0x7ff00000)
1548 /* +Inf -> +Inf, -Inf or NaN -> NaN */
1549 return (x * ((hx < 0)? 0.0 : x));
1550 if (hx > 0x406573fa || /* x > 171.62... overflow to +inf */
1551 (hx == 0x406573fa && lx > 0xE561F647)) {
1552 z = x / tiny;
1553 return (z * z);
1554 }
1555 if (hx >= 0x40200000) { /* x >= 8 */
1556 ww = large_gam(x, &m);
1557 w = ww.h + ww.l;
1558 __HI(w) += m << 20;
1559 return (w);
1560 }
1561 if (hx > 0) { /* 0 < x < 8 */
1562 i = (int) x;
1563 ww = gam_n(i, x - (double) i);
1564 return (ww.h + ww.l);
1565 }
1566
1567 /* negative x */
1568 /* INDENT OFF */
1569 /*
1570 * compute: xk =
1571 * -2 ... x is an even int (-inf is even)
1572 * -1 ... x is an odd int
1573 * +0 ... x is not an int but chopped to an even int
1574 * +1 ... x is not an int but chopped to an odd int
1575 */
1576 /* INDENT ON */
1577 xk = 0;
1578 if (ix >= 0x43300000) {
1579 if (ix >= 0x43400000)
1580 xk = -2;
1581 else
|