864 #else
865 t0z2 = 0.46163214496836234126265954232572132343318L,
866 t0z2_l = 5.03501162329616380465302666480916271611101e-36L,
867 #endif
868 /* 0.81977310110050060178786870492160699631174407846245179119586 */
869 #if defined(__x86)
870 t0z3 = 0.81977310110050060178773362329351925836817L,
871 t0z3_l = 1.350816280877379435658077052534574556256230e-22L
872 #else
873 t0z3 = 0.8197731011005006017878687049216069516957449L,
874 t0z3_l = 4.461599916947014419045492615933551648857380e-35L
875 #endif
876 ;
877 /* INDENT ON */
878
879 /*
880 * gamma(x+i) for 0 <= x < 1
881 */
882 static struct LDouble
883 gam_n(int i, long double x) {
884 struct LDouble rr, yy;
885 long double r1, r2, t2, z, xh, xl, yh, yl, zh, z1, z2, zl, x5, wh, wl;
886
887 /* compute yy = gamma(x+1) */
888 if (x > 0.2845L) {
889 if (x > 0.6374L) {
890 r1 = x - t0z3;
891 r2 = CHOPPED((r1 - t0z3_l));
892 t2 = r1 - r2;
893 yy = GT3(r2, t2 - t0z3_l);
894 } else {
895 r1 = x - t0z2;
896 r2 = CHOPPED((r1 - t0z2_l));
897 t2 = r1 - r2;
898 yy = GT2(r2, t2 - t0z2_l);
899 }
900 } else {
901 r1 = x - t0z1;
902 r2 = CHOPPED((r1 - t0z1_l));
903 t2 = r1 - r2;
904 yy = GT1(r2, t2 - t0z1_l);
1015 int i, j, m, ix, hx, xk;
1016 unsigned lx;
1017
1018 hx = H0_WORD(x);
1019 lx = H3_WORD(x);
1020 ix = hx & 0x7fffffff;
1021 y = x;
1022 if (ix < 0x3f8e0000) { /* x < 2**-113 */
1023 return (one / x);
1024 }
1025 if (ix >= 0x7fff0000)
1026 return (x * ((hx < 0)? zero : x)); /* Inf or NaN */
1027 if (x > overflow) /* overflow threshold */
1028 return (x * 1.0e4932L);
1029 if (hx >= 0x40020000) { /* x >= 8 */
1030 ww = large_gam(x, &m);
1031 w = ww.h + ww.l;
1032 return (scalbnl(w, m));
1033 }
1034
1035 if (hx > 0) { /* x from 0 to 8 */
1036 i = (int) x;
1037 ww = gam_n(i, x - (long double) i);
1038 return (ww.h + ww.l);
1039 }
1040 /* INDENT OFF */
1041 /* negative x */
1042 /*
1043 * compute xk =
1044 * -2 ... x is an even int (-inf is considered an even #)
1045 * -1 ... x is an odd int
1046 * +0 ... x is not an int but chopped to an even int
1047 * +1 ... x is not an int but chopped to an odd int
1048 */
1049 /* INDENT ON */
1050 xk = 0;
1051 #if defined(__x86)
1052 if (ix >= 0x403e0000) { /* x >= 2**63 } */
1053 if (ix >= 0x403f0000)
1054 xk = -2;
1055 else
|
864 #else
865 t0z2 = 0.46163214496836234126265954232572132343318L,
866 t0z2_l = 5.03501162329616380465302666480916271611101e-36L,
867 #endif
868 /* 0.81977310110050060178786870492160699631174407846245179119586 */
869 #if defined(__x86)
870 t0z3 = 0.81977310110050060178773362329351925836817L,
871 t0z3_l = 1.350816280877379435658077052534574556256230e-22L
872 #else
873 t0z3 = 0.8197731011005006017878687049216069516957449L,
874 t0z3_l = 4.461599916947014419045492615933551648857380e-35L
875 #endif
876 ;
877 /* INDENT ON */
878
879 /*
880 * gamma(x+i) for 0 <= x < 1
881 */
882 static struct LDouble
883 gam_n(int i, long double x) {
884 struct LDouble rr = {0.0L, 0.0L}, yy;
885 long double r1, r2, t2, z, xh, xl, yh, yl, zh, z1, z2, zl, x5, wh, wl;
886
887 /* compute yy = gamma(x+1) */
888 if (x > 0.2845L) {
889 if (x > 0.6374L) {
890 r1 = x - t0z3;
891 r2 = CHOPPED((r1 - t0z3_l));
892 t2 = r1 - r2;
893 yy = GT3(r2, t2 - t0z3_l);
894 } else {
895 r1 = x - t0z2;
896 r2 = CHOPPED((r1 - t0z2_l));
897 t2 = r1 - r2;
898 yy = GT2(r2, t2 - t0z2_l);
899 }
900 } else {
901 r1 = x - t0z1;
902 r2 = CHOPPED((r1 - t0z1_l));
903 t2 = r1 - r2;
904 yy = GT1(r2, t2 - t0z1_l);
1015 int i, j, m, ix, hx, xk;
1016 unsigned lx;
1017
1018 hx = H0_WORD(x);
1019 lx = H3_WORD(x);
1020 ix = hx & 0x7fffffff;
1021 y = x;
1022 if (ix < 0x3f8e0000) { /* x < 2**-113 */
1023 return (one / x);
1024 }
1025 if (ix >= 0x7fff0000)
1026 return (x * ((hx < 0)? zero : x)); /* Inf or NaN */
1027 if (x > overflow) /* overflow threshold */
1028 return (x * 1.0e4932L);
1029 if (hx >= 0x40020000) { /* x >= 8 */
1030 ww = large_gam(x, &m);
1031 w = ww.h + ww.l;
1032 return (scalbnl(w, m));
1033 }
1034
1035 if (hx > 0) { /* 0 < x < 8 */
1036 i = (int) x;
1037 ww = gam_n(i, x - (long double) i);
1038 return (ww.h + ww.l);
1039 }
1040 /* INDENT OFF */
1041 /* negative x */
1042 /*
1043 * compute xk =
1044 * -2 ... x is an even int (-inf is considered an even #)
1045 * -1 ... x is an odd int
1046 * +0 ... x is not an int but chopped to an even int
1047 * +1 ... x is not an int but chopped to an odd int
1048 */
1049 /* INDENT ON */
1050 xk = 0;
1051 #if defined(__x86)
1052 if (ix >= 0x403e0000) { /* x >= 2**63 } */
1053 if (ix >= 0x403f0000)
1054 xk = -2;
1055 else
|