Print this page
11210 libm should be cstyle(1ONBLD) clean
*** 20,37 ****
*/
/*
* Copyright 2011 Nexenta Systems, Inc. All rights reserved.
*/
/*
* Copyright 2006 Sun Microsystems, Inc. All rights reserved.
* Use is subject to license terms.
*/
#pragma weak __tgamma = tgamma
! /* INDENT OFF */
/*
* True gamma function
* double tgamma(double x)
*
* Error:
--- 20,38 ----
*/
/*
* Copyright 2011 Nexenta Systems, Inc. All rights reserved.
*/
+
/*
* Copyright 2006 Sun Microsystems, Inc. All rights reserved.
* Use is subject to license terms.
*/
#pragma weak __tgamma = tgamma
! /* BEGIN CSTYLED */
/*
* True gamma function
* double tgamma(double x)
*
* Error:
*** 745,759 ****
* expm1(r) = r + Et1*r^2 + Et2*r^3 + ... + Et5*r^6, and
* 2**(j/32) is obtained by table look-up S[j]+S_trail[j].
* Remez error bound:
* |exp(r) - (1+r+Et1*r^2+...+Et5*r^6)| <= 2^(-63).
*/
#include "libm.h"
! #define __HI(x) ((int *) &x)[HIWORD]
! #define __LO(x) ((unsigned *) &x)[LOWORD]
struct Double {
double h;
double l;
};
--- 746,761 ----
* expm1(r) = r + Et1*r^2 + Et2*r^3 + ... + Et5*r^6, and
* 2**(j/32) is obtained by table look-up S[j]+S_trail[j].
* Remez error bound:
* |exp(r) - (1+r+Et1*r^2+...+Et5*r^6)| <= 2^(-63).
*/
+ /* END CSTYLED */
#include "libm.h"
! #define __HI(x) ((int *)&x)[HIWORD]
! #define __LO(x) ((unsigned *)&x)[LOWORD]
struct Double {
double h;
double l;
};
*** 1131,1223 ****
#define Q32 cr[29]
#define Q33 cr[30]
#define Q34 cr[31]
#define Q35 cr[32]
! static const double
! GZ1_h = +0.938204627909682398190,
GZ1_l = +5.121952600248205157935e-17,
GZ2_h = +0.885603194410888749921,
GZ2_l = -4.964236872556339810692e-17,
GZ3_h = +0.936781411463652347038,
GZ3_l = -2.541923110834479415023e-17,
TZ1 = -0.3517214357852935791015625,
TZ3 = +0.280530631542205810546875;
- /* INDENT ON */
! /* compute gamma(y=yh+yl) for y in GT1 = [1.0000, 1.2845] */
! /* assume yh got 20 significant bits */
static struct Double
! GT1(double yh, double yl) {
double t3, t4, y, z;
struct Double r;
y = yh + yl;
z = y * y;
! t3 = (z * (P10 + y * ((P11 + y * P12) + z * (P13 + y * P14)))) /
! (Q10 + y * ((Q11 + y * Q12) + z * ((Q13 + Q14 * y) + z * Q15)));
t3 += (TZ1 * yl + GZ1_l);
t4 = TZ1 * yh;
! r.h = (double) ((float) (t4 + GZ1_h + t3));
t3 += (t4 - (r.h - GZ1_h));
r.l = t3;
return (r);
}
! /* compute gamma(y=yh+yl) for y in GT2 = [1.2844, 1.6374] */
! /* assume yh got 20 significant bits */
static struct Double
! GT2(double yh, double yl) {
double t3, y, z;
struct Double r;
y = yh + yl;
z = y * y;
! t3 = (z * (P20 + y * P21 + z * (P22 + y * P23))) /
! (Q20 + (y * ((Q21 + Q22 * y) + z * Q23) +
! (z * z) * ((Q24 + Q25 * y) + z * Q26))) + GZ2_l;
! r.h = (double) ((float) (GZ2_h + t3));
r.l = t3 - (r.h - GZ2_h);
return (r);
}
! /* compute gamma(y=yh+yl) for y in GT3 = [1.6373, 2.0000] */
! /* assume yh got 20 significant bits */
static struct Double
! GT3(double yh, double yl) {
double t3, t4, y, z;
struct Double r;
y = yh + yl;
z = y * y;
! t3 = (z * (P30 + y * ((P31 + y * P32) + z * (P33 + y * P34)))) /
! (Q30 + y * ((Q31 + y * Q32) + z * ((Q33 + Q34 * y) + z * Q35)));
t3 += (TZ3 * yl + GZ3_l);
t4 = TZ3 * yh;
! r.h = (double) ((float) (t4 + GZ3_h + t3));
t3 += (t4 - (r.h - GZ3_h));
r.l = t3;
return (r);
}
! /* INDENT OFF */
/*
* return tgamma(x) scaled by 2**-m for 8<x<=171.62... using Stirling's formula
* log(G(x)) ~= (x-.5)*(log(x)-1) + .5(log(2*pi)-1) + (1/x)*P(1/(x*x))
* = L1 + L2 + L3,
*/
- /* INDENT ON */
static struct Double
! large_gam(double x, int *m) {
! double z, t1, t2, t3, z2, t5, w, y, u, r, z4, v, t24 = 16777216.0,
! p24 = 1.0 / 16777216.0;
int n2, j2, k, ix, j;
unsigned lx;
struct Double zz;
double u2, ss_h, ss_l, r_h, w_h, w_l, t4;
! /* INDENT OFF */
/*
* compute ss = ss.h+ss.l = log(x)-1 (see tgamma_log.h for details)
*
* log(x) - 1 = T1(n) + T2(j) + T3(s), where x = 2**n * y, 1<=y<2,
* j=[64*y], z[j]=1+j/64+1/128, s = (y-z[j])/(y+z[j]), and
--- 1133,1232 ----
#define Q32 cr[29]
#define Q33 cr[30]
#define Q34 cr[31]
#define Q35 cr[32]
! static const double GZ1_h = +0.938204627909682398190,
GZ1_l = +5.121952600248205157935e-17,
GZ2_h = +0.885603194410888749921,
GZ2_l = -4.964236872556339810692e-17,
GZ3_h = +0.936781411463652347038,
GZ3_l = -2.541923110834479415023e-17,
TZ1 = -0.3517214357852935791015625,
TZ3 = +0.280530631542205810546875;
! /*
! * compute gamma(y=yh+yl) for y in GT1 = [1.0000, 1.2845]
! * assume yh got 20 significant bits
! */
static struct Double
! GT1(double yh, double yl)
! {
double t3, t4, y, z;
struct Double r;
y = yh + yl;
z = y * y;
! t3 = (z * (P10 + y * ((P11 + y * P12) + z * (P13 + y * P14)))) / (Q10 +
! y * ((Q11 + y * Q12) + z * ((Q13 + Q14 * y) + z * Q15)));
t3 += (TZ1 * yl + GZ1_l);
t4 = TZ1 * yh;
! r.h = (double)((float)(t4 + GZ1_h + t3));
t3 += (t4 - (r.h - GZ1_h));
r.l = t3;
return (r);
}
! /*
! * compute gamma(y=yh+yl) for y in GT2 = [1.2844, 1.6374]
! * assume yh got 20 significant bits
! */
static struct Double
! GT2(double yh, double yl)
! {
double t3, y, z;
struct Double r;
y = yh + yl;
z = y * y;
! t3 = (z * (P20 + y * P21 + z * (P22 + y * P23))) / (Q20 + (y * ((Q21 +
! Q22 * y) + z * Q23) + (z * z) * ((Q24 + Q25 * y) + z * Q26))) +
! GZ2_l;
! r.h = (double)((float)(GZ2_h + t3));
r.l = t3 - (r.h - GZ2_h);
return (r);
}
! /*
! * compute gamma(y=yh+yl) for y in GT3 = [1.6373, 2.0000]
! * assume yh got 20 significant bits
! */
static struct Double
! GT3(double yh, double yl)
! {
double t3, t4, y, z;
struct Double r;
y = yh + yl;
z = y * y;
! t3 = (z * (P30 + y * ((P31 + y * P32) + z * (P33 + y * P34)))) / (Q30 +
! y * ((Q31 + y * Q32) + z * ((Q33 + Q34 * y) + z * Q35)));
t3 += (TZ3 * yl + GZ3_l);
t4 = TZ3 * yh;
! r.h = (double)((float)(t4 + GZ3_h + t3));
t3 += (t4 - (r.h - GZ3_h));
r.l = t3;
return (r);
}
!
/*
* return tgamma(x) scaled by 2**-m for 8<x<=171.62... using Stirling's formula
* log(G(x)) ~= (x-.5)*(log(x)-1) + .5(log(2*pi)-1) + (1/x)*P(1/(x*x))
* = L1 + L2 + L3,
*/
static struct Double
! large_gam(double x, int *m)
! {
! double z, t1, t2, t3, z2, t5, w, y, u, r, z4, v, t24 = 16777216.0, p24 =
! 1.0 / 16777216.0;
int n2, j2, k, ix, j;
unsigned lx;
struct Double zz;
double u2, ss_h, ss_l, r_h, w_h, w_l, t4;
!
/*
* compute ss = ss.h+ss.l = log(x)-1 (see tgamma_log.h for details)
*
* log(x) - 1 = T1(n) + T2(j) + T3(s), where x = 2**n * y, 1<=y<2,
* j=[64*y], z[j]=1+j/64+1/128, s = (y-z[j])/(y+z[j]), and
*** 1237,1247 ****
* __________________________
* + T3(s)-2s: |__________________________|
* -------------------------------------------
* [leading] + [Trailing]
*/
- /* INDENT ON */
ix = __HI(x);
lx = __LO(x);
n2 = (ix >> 20) - 0x3ff; /* exponent of x, range:3-7 */
n2 += n2; /* 2n */
ix = (ix & 0x000fffff) | 0x3ff00000; /* y = scale x to [1,2] */
--- 1246,1255 ----
*** 1251,1300 ****
__LO(z) = 0;
j2 = (ix >> 13) & 0x7e; /* 2j */
t1 = y + z;
t2 = y - z;
r = one / t1;
! t1 = (double) ((float) t1);
u = r * t2; /* u = (y-z)/(y+z) */
t4 = T2[j2 + 1] + T1[n2 + 1];
z2 = u * u;
k = __HI(u) & 0x7fffffff;
t3 = T2[j2] + T1[n2];
if ((k >> 20) < 0x3ec) { /* |u|<2**-19 */
t2 = t4 + u * ((two + z2 * A1) + (z2 * z2) * (A2 + z2 * A3));
} else {
t5 = t4 + u * (z2 * A1 + (z2 * z2) * (A2 + z2 * A3));
u2 = u + u;
! v = (double) ((int) (u2 * t24)) * p24;
t2 = t5 + r * ((two * t2 - v * t1) - v * (y - (t1 - z)));
t3 += v;
}
! ss_h = (double) ((float) (t2 + t3));
ss_l = t2 - (ss_h - t3);
/*
* compute ww = (x-.5)*(log(x)-1) + .5*(log(2pi)-1) + 1/x*(P(1/x^2)))
* where ss = log(x) - 1 in already in extra precision
*/
z = one / x;
r = x - half;
! r_h = (double) ((float) r);
w_h = r_h * ss_h + hln2pi_h;
z2 = z * z;
w = (r - r_h) * ss_h + r * ss_l;
z4 = z2 * z2;
t1 = z2 * (GP1 + z4 * (GP3 + z4 * (GP5 + z4 * GP7)));
t2 = z4 * (GP2 + z4 * (GP4 + z4 * GP6));
t1 += t2;
w += hln2pi_l;
w_l = z * (GP0 + t1) + w;
! k = (int) ((w_h + w_l) * invln2_32 + half);
/* compute the exponential of w_h+w_l */
j = k & 0x1f;
*m = (k >> 5);
! t3 = (double) k;
/* perform w - k*ln2_32 (represent as w_h - w_l) */
t1 = w_h - t3 * ln2_32hi;
t2 = t3 * ln2_32lo;
w = w_l - t2;
--- 1259,1310 ----
__LO(z) = 0;
j2 = (ix >> 13) & 0x7e; /* 2j */
t1 = y + z;
t2 = y - z;
r = one / t1;
! t1 = (double)((float)t1);
u = r * t2; /* u = (y-z)/(y+z) */
t4 = T2[j2 + 1] + T1[n2 + 1];
z2 = u * u;
k = __HI(u) & 0x7fffffff;
t3 = T2[j2] + T1[n2];
+
if ((k >> 20) < 0x3ec) { /* |u|<2**-19 */
t2 = t4 + u * ((two + z2 * A1) + (z2 * z2) * (A2 + z2 * A3));
} else {
t5 = t4 + u * (z2 * A1 + (z2 * z2) * (A2 + z2 * A3));
u2 = u + u;
! v = (double)((int)(u2 * t24)) * p24;
t2 = t5 + r * ((two * t2 - v * t1) - v * (y - (t1 - z)));
t3 += v;
}
!
! ss_h = (double)((float)(t2 + t3));
ss_l = t2 - (ss_h - t3);
/*
* compute ww = (x-.5)*(log(x)-1) + .5*(log(2pi)-1) + 1/x*(P(1/x^2)))
* where ss = log(x) - 1 in already in extra precision
*/
z = one / x;
r = x - half;
! r_h = (double)((float)r);
w_h = r_h * ss_h + hln2pi_h;
z2 = z * z;
w = (r - r_h) * ss_h + r * ss_l;
z4 = z2 * z2;
t1 = z2 * (GP1 + z4 * (GP3 + z4 * (GP5 + z4 * GP7)));
t2 = z4 * (GP2 + z4 * (GP4 + z4 * GP6));
t1 += t2;
w += hln2pi_l;
w_l = z * (GP0 + t1) + w;
! k = (int)((w_h + w_l) * invln2_32 + half);
/* compute the exponential of w_h+w_l */
j = k & 0x1f;
*m = (k >> 5);
! t3 = (double)k;
/* perform w - k*ln2_32 (represent as w_h - w_l) */
t1 = w_h - t3 * ln2_32hi;
t2 = t3 * ln2_32lo;
w = w_l - t2;
*** 1310,1320 ****
zz.l = S_trail[j] * (one + t3) + S[j] * t3;
zz.h = S[j];
return (zz);
}
! /* INDENT OFF */
/*
* kpsin(x)= sin(pi*x)/pi
* 3 5 7 9 11 13 15
* = x+ks[0]*x +ks[1]*x +ks[2]*x +ks[3]*x +ks[4]*x +ks[5]*x +ks[6]*x
*/
--- 1320,1330 ----
zz.l = S_trail[j] * (one + t3) + S[j] * t3;
zz.h = S[j];
return (zz);
}
!
/*
* kpsin(x)= sin(pi*x)/pi
* 3 5 7 9 11 13 15
* = x+ks[0]*x +ks[1]*x +ks[2]*x +ks[3]*x +ks[4]*x +ks[5]*x +ks[6]*x
*/
*** 1325,1354 ****
+2.61478477632554278317289628332654539353521911570e-0002,
-2.34607978510202710377617190278735525354347705866e-0003,
+1.48413292290051695897242899977121846763824221705e-0004,
-6.87730769637543488108688726777687262485357072242e-0006,
};
- /* INDENT ON */
/* assume x is not tiny and positive */
static struct Double
! kpsin(double x) {
double z, t1, t2, t3, t4;
struct Double xx;
z = x * x;
xx.h = x;
t1 = z * x;
t2 = z * z;
t4 = t1 * ks[0];
! t3 = (t1 * z) * ((ks[1] + z * ks[2] + t2 * ks[3]) + (z * t2) *
! (ks[4] + z * ks[5] + t2 * ks[6]));
xx.l = t4 + t3;
return (xx);
}
! /* INDENT OFF */
/*
* kpcos(x)= cos(pi*x)/pi
* 2 4 6 8 10 12
* = 1/pi +kc[0]*x +kc[1]*x +kc[2]*x +kc[3]*x +kc[4]*x +kc[5]*x
*/
--- 1335,1364 ----
+2.61478477632554278317289628332654539353521911570e-0002,
-2.34607978510202710377617190278735525354347705866e-0003,
+1.48413292290051695897242899977121846763824221705e-0004,
-6.87730769637543488108688726777687262485357072242e-0006,
};
/* assume x is not tiny and positive */
static struct Double
! kpsin(double x)
! {
double z, t1, t2, t3, t4;
struct Double xx;
z = x * x;
xx.h = x;
t1 = z * x;
t2 = z * z;
t4 = t1 * ks[0];
! t3 = (t1 * z) * ((ks[1] + z * ks[2] + t2 * ks[3]) + (z * t2) * (ks[4] +
! z * ks[5] + t2 * ks[6]));
xx.l = t4 + t3;
return (xx);
}
!
/*
* kpcos(x)= cos(pi*x)/pi
* 2 4 6 8 10 12
* = 1/pi +kc[0]*x +kc[1]*x +kc[2]*x +kc[3]*x +kc[4]*x +kc[5]*x
*/
*** 1363,1459 ****
-4.25027339940149518500158850753393173519732149213e-0001,
+7.49080625187015312373925142219429422375556727752e-0002,
-8.21442040906099210866977352284054849051348692715e-0003,
+6.10411356829515414575566564733632532333904115968e-0004,
};
- /* INDENT ON */
/* assume x is not tiny and positive */
static struct Double
! kpcos(double x) {
double z, t1, t2, t3, t4, x4, x8;
struct Double xx;
z = x * x;
xx.h = one_pi_h;
! t1 = (double) ((float) x);
x4 = z * z;
t2 = npi_2_l * z + npi_2_h * (x + t1) * (x - t1);
! t3 = one_pi_l + x4 * ((kc[1] + z * kc[2]) + x4 * (kc[3] + z *
! kc[4] + x4 * kc[5]));
t4 = t1 * t1; /* 48 bits mantissa */
x8 = t2 + t3;
t4 *= npi_2_h; /* npi_2_h is 5 bits const. The product is exact */
xx.l = x8 + t4; /* that will minimized the rounding error in xx.l */
return (xx);
}
- /* INDENT OFF */
static const double
! /* 0.134861805732790769689793935774652917006 */
t0z1 = 0.1348618057327907737708,
t0z1_l = -4.0810077708578299022531e-18,
! /* 0.461632144968362341262659542325721328468 */
t0z2 = 0.4616321449683623567850,
t0z2_l = -1.5522348162858676890521e-17,
! /* 0.819773101100500601787868704921606996312 */
t0z3 = 0.8197731011005006118708,
t0z3_l = -1.0082945122487103498325e-17;
! /* 1.134861805732790769689793935774652917006 */
! /* INDENT ON */
/* gamma(x+i) for 0 <= x < 1 */
static struct Double
! gam_n(int i, double x) {
! struct Double rr = {0.0L, 0.0L}, yy;
double r1, r2, t2, z, xh, xl, yh, yl, zh, z1, z2, zl, x5, wh, wl;
/* compute yy = gamma(x+1) */
if (x > 0.2845) {
if (x > 0.6374) {
r1 = x - t0z3;
! r2 = (double) ((float) (r1 - t0z3_l));
t2 = r1 - r2;
yy = GT3(r2, t2 - t0z3_l);
} else {
r1 = x - t0z2;
! r2 = (double) ((float) (r1 - t0z2_l));
t2 = r1 - r2;
yy = GT2(r2, t2 - t0z2_l);
}
} else {
r1 = x - t0z1;
! r2 = (double) ((float) (r1 - t0z1_l));
t2 = r1 - r2;
yy = GT1(r2, t2 - t0z1_l);
}
/* compute gamma(x+i) = (x+i-1)*...*(x+1)*yy, 0<i<8 */
switch (i) {
case 0: /* yy/x */
r1 = one / x;
! xh = (double) ((float) x); /* x is not tiny */
! rr.h = (double) ((float) ((yy.h + yy.l) * r1));
! rr.l = r1 * (yy.h - rr.h * xh) -
! ((r1 * rr.h) * (x - xh) - r1 * yy.l);
break;
case 1: /* yy */
rr.h = yy.h;
rr.l = yy.l;
break;
case 2: /* (x+1)*yy */
z = x + one; /* may not be exact */
! zh = (double) ((float) z);
rr.h = zh * yy.h;
rr.l = z * yy.l + (x - (zh - one)) * yy.h;
break;
case 3: /* (x+2)*(x+1)*yy */
z1 = x + one;
z2 = x + 2.0;
z = z1 * z2;
! xh = (double) ((float) z);
! zh = (double) ((float) z1);
xl = (x - (zh - one)) * (z2 + zh) - (xh - zh * (zh + one));
rr.h = xh * yy.h;
rr.l = z * yy.l + xl * yy.h;
break;
--- 1373,1471 ----
-4.25027339940149518500158850753393173519732149213e-0001,
+7.49080625187015312373925142219429422375556727752e-0002,
-8.21442040906099210866977352284054849051348692715e-0003,
+6.10411356829515414575566564733632532333904115968e-0004,
};
/* assume x is not tiny and positive */
static struct Double
! kpcos(double x)
! {
double z, t1, t2, t3, t4, x4, x8;
struct Double xx;
z = x * x;
xx.h = one_pi_h;
! t1 = (double)((float)x);
x4 = z * z;
t2 = npi_2_l * z + npi_2_h * (x + t1) * (x - t1);
! t3 = one_pi_l + x4 * ((kc[1] + z * kc[2]) + x4 * (kc[3] + z * kc[4] +
! x4 * kc[5]));
t4 = t1 * t1; /* 48 bits mantissa */
x8 = t2 + t3;
t4 *= npi_2_h; /* npi_2_h is 5 bits const. The product is exact */
xx.l = x8 + t4; /* that will minimized the rounding error in xx.l */
return (xx);
}
static const double
! /* 0.134861805732790769689793935774652917006 */
t0z1 = 0.1348618057327907737708,
t0z1_l = -4.0810077708578299022531e-18,
! /* 0.461632144968362341262659542325721328468 */
t0z2 = 0.4616321449683623567850,
t0z2_l = -1.5522348162858676890521e-17,
! /* 0.819773101100500601787868704921606996312 */
t0z3 = 0.8197731011005006118708,
t0z3_l = -1.0082945122487103498325e-17;
!
! /*
! * 1.134861805732790769689793935774652917006
! */
/* gamma(x+i) for 0 <= x < 1 */
static struct Double
! gam_n(int i, double x)
! {
! struct Double rr = { 0.0L, 0.0L }, yy;
double r1, r2, t2, z, xh, xl, yh, yl, zh, z1, z2, zl, x5, wh, wl;
/* compute yy = gamma(x+1) */
if (x > 0.2845) {
if (x > 0.6374) {
r1 = x - t0z3;
! r2 = (double)((float)(r1 - t0z3_l));
t2 = r1 - r2;
yy = GT3(r2, t2 - t0z3_l);
} else {
r1 = x - t0z2;
! r2 = (double)((float)(r1 - t0z2_l));
t2 = r1 - r2;
yy = GT2(r2, t2 - t0z2_l);
}
} else {
r1 = x - t0z1;
! r2 = (double)((float)(r1 - t0z1_l));
t2 = r1 - r2;
yy = GT1(r2, t2 - t0z1_l);
}
/* compute gamma(x+i) = (x+i-1)*...*(x+1)*yy, 0<i<8 */
switch (i) {
case 0: /* yy/x */
r1 = one / x;
! xh = (double)((float)x); /* x is not tiny */
! rr.h = (double)((float)((yy.h + yy.l) * r1));
! rr.l = r1 * (yy.h - rr.h * xh) - ((r1 * rr.h) * (x - xh) - r1 *
! yy.l);
break;
case 1: /* yy */
rr.h = yy.h;
rr.l = yy.l;
break;
case 2: /* (x+1)*yy */
z = x + one; /* may not be exact */
! zh = (double)((float)z);
rr.h = zh * yy.h;
rr.l = z * yy.l + (x - (zh - one)) * yy.h;
break;
case 3: /* (x+2)*(x+1)*yy */
z1 = x + one;
z2 = x + 2.0;
z = z1 * z2;
! xh = (double)((float)z);
! zh = (double)((float)z1);
xl = (x - (zh - one)) * (z2 + zh) - (xh - zh * (zh + one));
rr.h = xh * yy.h;
rr.l = z * yy.l + xl * yy.h;
break;
*** 1463,1537 ****
zh = z1;
__LO(zh) = 0;
__HI(zh) &= 0xfffffff8; /* zh 18 bits mantissa */
zl = x - (zh - 2.0);
z = z1 * z2;
! xh = (double) ((float) z);
xl = zl * (z2 + zh * (z1 + zh)) - (xh - zh * (zh * zh - one));
rr.h = xh * yy.h;
rr.l = z * yy.l + xl * yy.h;
break;
case 5: /* ((x+1)*(x+4)*(x+2)*(x+3))*yy */
z1 = x + 2.0;
z2 = x + 3.0;
z = z1 * z2;
! zh = (double) ((float) z1);
! yh = (double) ((float) z);
yl = (x - (zh - 2.0)) * (z2 + zh) - (yh - zh * (zh + one));
z2 = z - 2.0;
z *= z2;
! xh = (double) ((float) z);
xl = yl * (z2 + yh) - (xh - yh * (yh - 2.0));
rr.h = xh * yy.h;
rr.l = z * yy.l + xl * yy.h;
break;
case 6: /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5))*yy */
z1 = x + 2.0;
z2 = x + 3.0;
z = z1 * z2;
! zh = (double) ((float) z1);
! yh = (double) ((float) z);
z1 = x - (zh - 2.0);
yl = z1 * (z2 + zh) - (yh - zh * (zh + one));
z2 = z - 2.0;
x5 = x + 5.0;
z *= z2;
! xh = (double) ((float) z);
zh += 3.0;
xl = yl * (z2 + yh) - (xh - yh * (yh - 2.0));
! /* xh+xl=(x+1)*...*(x+4) */
! /* wh+wl=(x+5)*yy */
! wh = (double) ((float) (x5 * (yy.h + yy.l)));
wl = (z1 * yy.h + x5 * yy.l) - (wh - zh * yy.h);
rr.h = wh * xh;
rr.l = z * wl + xl * wh;
break;
case 7: /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5)*(x+6))*yy */
z1 = x + 3.0;
z2 = x + 4.0;
z = z2 * z1;
! zh = (double) ((float) z1);
! yh = (double) ((float) z); /* yh+yl = (x+3)(x+4) */
yl = (x - (zh - 3.0)) * (z2 + zh) - (yh - (zh * (zh + one)));
z1 = x + 6.0;
z2 = z - 2.0; /* z2 = (x+2)*(x+5) */
z *= z2;
! xh = (double) ((float) z);
xl = yl * (z2 + yh) - (xh - yh * (yh - 2.0));
! /* xh+xl=(x+2)*...*(x+5) */
! /* wh+wl=(x+1)(x+6)*yy */
z2 -= 4.0; /* z2 = (x+1)(x+6) */
! wh = (double) ((float) (z2 * (yy.h + yy.l)));
wl = (z2 * yy.l + yl * yy.h) - (wh - (yh - 6.0) * yy.h);
rr.h = wh * xh;
rr.l = z * wl + xl * wh;
}
return (rr);
}
double
! tgamma(double x) {
struct Double ss, ww;
double t, t1, t2, t3, t4, t5, w, y, z, z1, z2, z3, z5;
int i, j, k, m, ix, hx, xk;
unsigned lx;
--- 1475,1557 ----
zh = z1;
__LO(zh) = 0;
__HI(zh) &= 0xfffffff8; /* zh 18 bits mantissa */
zl = x - (zh - 2.0);
z = z1 * z2;
! xh = (double)((float)z);
xl = zl * (z2 + zh * (z1 + zh)) - (xh - zh * (zh * zh - one));
rr.h = xh * yy.h;
rr.l = z * yy.l + xl * yy.h;
break;
case 5: /* ((x+1)*(x+4)*(x+2)*(x+3))*yy */
z1 = x + 2.0;
z2 = x + 3.0;
z = z1 * z2;
! zh = (double)((float)z1);
! yh = (double)((float)z);
yl = (x - (zh - 2.0)) * (z2 + zh) - (yh - zh * (zh + one));
z2 = z - 2.0;
z *= z2;
! xh = (double)((float)z);
xl = yl * (z2 + yh) - (xh - yh * (yh - 2.0));
rr.h = xh * yy.h;
rr.l = z * yy.l + xl * yy.h;
break;
case 6: /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5))*yy */
z1 = x + 2.0;
z2 = x + 3.0;
z = z1 * z2;
! zh = (double)((float)z1);
! yh = (double)((float)z);
z1 = x - (zh - 2.0);
yl = z1 * (z2 + zh) - (yh - zh * (zh + one));
z2 = z - 2.0;
x5 = x + 5.0;
z *= z2;
! xh = (double)((float)z);
zh += 3.0;
xl = yl * (z2 + yh) - (xh - yh * (yh - 2.0));
!
! /*
! * xh+xl=(x+1)*...*(x+4)
! * wh+wl=(x+5)*yy
! */
! wh = (double)((float)(x5 * (yy.h + yy.l)));
wl = (z1 * yy.h + x5 * yy.l) - (wh - zh * yy.h);
rr.h = wh * xh;
rr.l = z * wl + xl * wh;
break;
case 7: /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5)*(x+6))*yy */
z1 = x + 3.0;
z2 = x + 4.0;
z = z2 * z1;
! zh = (double)((float)z1);
! yh = (double)((float)z); /* yh+yl = (x+3)(x+4) */
yl = (x - (zh - 3.0)) * (z2 + zh) - (yh - (zh * (zh + one)));
z1 = x + 6.0;
z2 = z - 2.0; /* z2 = (x+2)*(x+5) */
z *= z2;
! xh = (double)((float)z);
xl = yl * (z2 + yh) - (xh - yh * (yh - 2.0));
!
! /*
! * xh+xl=(x+2)*...*(x+5)
! * wh+wl=(x+1)(x+6)*yy
! */
z2 -= 4.0; /* z2 = (x+1)(x+6) */
! wh = (double)((float)(z2 * (yy.h + yy.l)));
wl = (z2 * yy.l + yl * yy.h) - (wh - (yh - 6.0) * yy.h);
rr.h = wh * xh;
rr.l = z * wl + xl * wh;
}
+
return (rr);
}
double
! tgamma(double x)
! {
struct Double ss, ww;
double t, t1, t2, t3, t4, t5, w, y, z, z1, z2, z3, z5;
int i, j, k, m, ix, hx, xk;
unsigned lx;
*** 1540,1686 ****
ix = hx & 0x7fffffff;
y = x;
if (ix < 0x3ca00000)
return (one / x); /* |x| < 2**-53 */
if (ix >= 0x7ff00000)
/* +Inf -> +Inf, -Inf or NaN -> NaN */
! return (x * ((hx < 0)? 0.0 : x));
if (hx > 0x406573fa || /* x > 171.62... overflow to +inf */
(hx == 0x406573fa && lx > 0xE561F647)) {
z = x / tiny;
return (z * z);
}
if (hx >= 0x40200000) { /* x >= 8 */
ww = large_gam(x, &m);
w = ww.h + ww.l;
__HI(w) += m << 20;
return (w);
}
if (hx > 0) { /* 0 < x < 8 */
! i = (int) x;
! ww = gam_n(i, x - (double) i);
return (ww.h + ww.l);
}
! /* negative x */
! /* INDENT OFF */
/*
* compute: xk =
* -2 ... x is an even int (-inf is even)
* -1 ... x is an odd int
* +0 ... x is not an int but chopped to an even int
* +1 ... x is not an int but chopped to an odd int
*/
- /* INDENT ON */
xk = 0;
if (ix >= 0x43300000) {
if (ix >= 0x43400000)
xk = -2;
else
xk = -2 + (lx & 1);
} else if (ix >= 0x3ff00000) {
k = (ix >> 20) - 0x3ff;
if (k > 20) {
j = lx >> (52 - k);
if ((j << (52 - k)) == lx)
xk = -2 + (j & 1);
else
xk = j & 1;
} else {
j = ix >> (20 - k);
if ((j << (20 - k)) == ix && lx == 0)
xk = -2 + (j & 1);
else
xk = j & 1;
}
}
if (xk < 0)
/* ideally gamma(-n)= (-1)**(n+1) * inf, but c99 expect NaN */
return ((x - x) / (x - x)); /* 0/0 = NaN */
-
/* negative underflow thresold */
if (ix > 0x4066e000 || (ix == 0x4066e000 && lx > 11)) {
/* x < -183.0 - 11ulp */
z = tiny / x;
if (xk == 1)
z = -z;
return (z * tiny);
}
/* now compute gamma(x) by -1/((sin(pi*y)/pi)*gamma(1+y)), y = -x */
/*
* First compute ss = -sin(pi*y)/pi , so that
* gamma(x) = 1/(ss*gamma(1+y))
*/
y = -x;
! j = (int) y;
! z = y - (double) j;
! if (z > 0.3183098861837906715377675)
if (z > 0.6816901138162093284622325)
ss = kpsin(one - z);
else
ss = kpcos(0.5 - z);
! else
ss = kpsin(z);
if (xk == 0) {
ss.h = -ss.h;
ss.l = -ss.l;
}
/* Then compute ww = gamma(1+y), note that result scale to 2**m */
m = 0;
if (j < 7) {
ww = gam_n(j + 1, z);
} else {
w = y + one;
if ((lx & 1) == 0) { /* y+1 exact (note that y<184) */
ww = large_gam(w, &m);
} else {
t = w - one;
if (t == y) { /* y+one exact */
ww = large_gam(w, &m);
} else { /* use y*gamma(y) */
if (j == 7)
ww = gam_n(j, z);
else
ww = large_gam(y, &m);
t4 = ww.h + ww.l;
! t1 = (double) ((float) y);
! t2 = (double) ((float) t4);
/* t4 will not be too large */
ww.l = y * (ww.l - (t2 - ww.h)) + (y - t1) * t2;
ww.h = t1 * t2;
}
}
}
/* compute 1/(ss*ww) */
t3 = ss.h + ss.l;
t4 = ww.h + ww.l;
! t1 = (double) ((float) t3);
! t2 = (double) ((float) t4);
z1 = ss.l - (t1 - ss.h); /* (t1,z1) = ss */
z2 = ww.l - (t2 - ww.h); /* (t2,z2) = ww */
t3 = t3 * t4; /* t3 = ss*ww */
z3 = one / t3; /* z3 = 1/(ss*ww) */
t5 = t1 * t2;
z5 = z1 * t4 + t1 * z2; /* (t5,z5) = ss*ww */
! t1 = (double) ((float) t3); /* (t1,z1) = ss*ww */
z1 = z5 - (t1 - t5);
! t2 = (double) ((float) z3); /* leading 1/(ss*ww) */
z2 = z3 * (t2 * z1 - (one - t2 * t1));
z = t2 - z2;
/* check whether z*2**-m underflow */
if (m != 0) {
hx = __HI(z);
i = hx & 0x80000000;
ix = hx ^ i;
j = ix >> 20;
if (j > m) {
ix -= m << 20;
__HI(z) = ix ^ i;
} else if ((m - j) > 52) {
/* underflow */
--- 1560,1725 ----
ix = hx & 0x7fffffff;
y = x;
if (ix < 0x3ca00000)
return (one / x); /* |x| < 2**-53 */
+
if (ix >= 0x7ff00000)
/* +Inf -> +Inf, -Inf or NaN -> NaN */
! return (x * ((hx < 0) ? 0.0 : x));
!
if (hx > 0x406573fa || /* x > 171.62... overflow to +inf */
(hx == 0x406573fa && lx > 0xE561F647)) {
z = x / tiny;
return (z * z);
}
+
if (hx >= 0x40200000) { /* x >= 8 */
ww = large_gam(x, &m);
w = ww.h + ww.l;
__HI(w) += m << 20;
return (w);
}
+
if (hx > 0) { /* 0 < x < 8 */
! i = (int)x;
! ww = gam_n(i, x - (double)i);
return (ww.h + ww.l);
}
! /*
! * negative x
! */
!
/*
* compute: xk =
* -2 ... x is an even int (-inf is even)
* -1 ... x is an odd int
* +0 ... x is not an int but chopped to an even int
* +1 ... x is not an int but chopped to an odd int
*/
xk = 0;
+
if (ix >= 0x43300000) {
if (ix >= 0x43400000)
xk = -2;
else
xk = -2 + (lx & 1);
} else if (ix >= 0x3ff00000) {
k = (ix >> 20) - 0x3ff;
+
if (k > 20) {
j = lx >> (52 - k);
+
if ((j << (52 - k)) == lx)
xk = -2 + (j & 1);
else
xk = j & 1;
} else {
j = ix >> (20 - k);
+
if ((j << (20 - k)) == ix && lx == 0)
xk = -2 + (j & 1);
else
xk = j & 1;
}
}
+
if (xk < 0)
/* ideally gamma(-n)= (-1)**(n+1) * inf, but c99 expect NaN */
return ((x - x) / (x - x)); /* 0/0 = NaN */
/* negative underflow thresold */
if (ix > 0x4066e000 || (ix == 0x4066e000 && lx > 11)) {
/* x < -183.0 - 11ulp */
z = tiny / x;
+
if (xk == 1)
z = -z;
+
return (z * tiny);
}
/* now compute gamma(x) by -1/((sin(pi*y)/pi)*gamma(1+y)), y = -x */
/*
* First compute ss = -sin(pi*y)/pi , so that
* gamma(x) = 1/(ss*gamma(1+y))
*/
y = -x;
! j = (int)y;
! z = y - (double)j;
!
! if (z > 0.3183098861837906715377675) {
if (z > 0.6816901138162093284622325)
ss = kpsin(one - z);
else
ss = kpcos(0.5 - z);
! } else {
ss = kpsin(z);
+ }
+
if (xk == 0) {
ss.h = -ss.h;
ss.l = -ss.l;
}
/* Then compute ww = gamma(1+y), note that result scale to 2**m */
m = 0;
+
if (j < 7) {
ww = gam_n(j + 1, z);
} else {
w = y + one;
+
if ((lx & 1) == 0) { /* y+1 exact (note that y<184) */
ww = large_gam(w, &m);
} else {
t = w - one;
+
if (t == y) { /* y+one exact */
ww = large_gam(w, &m);
} else { /* use y*gamma(y) */
if (j == 7)
ww = gam_n(j, z);
else
ww = large_gam(y, &m);
+
t4 = ww.h + ww.l;
! t1 = (double)((float)y);
! t2 = (double)((float)t4);
/* t4 will not be too large */
ww.l = y * (ww.l - (t2 - ww.h)) + (y - t1) * t2;
ww.h = t1 * t2;
}
}
}
/* compute 1/(ss*ww) */
t3 = ss.h + ss.l;
t4 = ww.h + ww.l;
! t1 = (double)((float)t3);
! t2 = (double)((float)t4);
z1 = ss.l - (t1 - ss.h); /* (t1,z1) = ss */
z2 = ww.l - (t2 - ww.h); /* (t2,z2) = ww */
t3 = t3 * t4; /* t3 = ss*ww */
z3 = one / t3; /* z3 = 1/(ss*ww) */
t5 = t1 * t2;
z5 = z1 * t4 + t1 * z2; /* (t5,z5) = ss*ww */
! t1 = (double)((float)t3); /* (t1,z1) = ss*ww */
z1 = z5 - (t1 - t5);
! t2 = (double)((float)z3); /* leading 1/(ss*ww) */
z2 = z3 * (t2 * z1 - (one - t2 * t1));
z = t2 - z2;
/* check whether z*2**-m underflow */
if (m != 0) {
hx = __HI(z);
i = hx & 0x80000000;
ix = hx ^ i;
j = ix >> 20;
+
if (j > m) {
ix -= m << 20;
__HI(z) = ix ^ i;
} else if ((m - j) > 52) {
/* underflow */
*** 1696,1702 ****
--- 1735,1742 ----
ix -= m << 20;
__HI(z) = ix ^ i;
z *= t;
}
}
+
return (z);
}