32 * K.C. Ng, August, 1989.
33 *
34 * We choose [1.5,2.5] to be the primary interval. Our algorithms
35 * are mainly derived from
36 *
37 *
38 * zeta(2)-1 2 zeta(3)-1 3
39 * lgamma(2+s) = s*(1-euler) + --------- * s - --------- * s + ...
40 * 2 3
41 *
42 *
43 * Note 1. Since gamma(1+s)=s*gamma(s), hence
44 * lgamma(1+s) = log(s) + lgamma(s), or
45 * lgamma(s) = lgamma(1+s) - log(s).
46 * When s is really tiny (like roundoff), lgamma(1+s) ~ s(1-enler)
47 * Hence lgamma(s) ~ -log(s) for tiny s
48 *
49 */
50
51 #include "libm.h"
52 #include "libm_synonyms.h"
53 #include "longdouble.h"
54
55 static long double neg(long double, int *);
56 static long double poly(long double, const long double *, int);
57 static long double polytail(long double);
58 static long double primary(long double);
59
60 static const long double
61 c0 = 0.0L,
62 ch = 0.5L,
63 c1 = 1.0L,
64 c2 = 2.0L,
65 c3 = 3.0L,
66 c4 = 4.0L,
67 c5 = 5.0L,
68 c6 = 6.0L,
69 pi = 3.1415926535897932384626433832795028841971L,
70 tiny = 1.0e-40L;
71
72 long double
|
32 * K.C. Ng, August, 1989.
33 *
34 * We choose [1.5,2.5] to be the primary interval. Our algorithms
35 * are mainly derived from
36 *
37 *
38 * zeta(2)-1 2 zeta(3)-1 3
39 * lgamma(2+s) = s*(1-euler) + --------- * s - --------- * s + ...
40 * 2 3
41 *
42 *
43 * Note 1. Since gamma(1+s)=s*gamma(s), hence
44 * lgamma(1+s) = log(s) + lgamma(s), or
45 * lgamma(s) = lgamma(1+s) - log(s).
46 * When s is really tiny (like roundoff), lgamma(1+s) ~ s(1-enler)
47 * Hence lgamma(s) ~ -log(s) for tiny s
48 *
49 */
50
51 #include "libm.h"
52 #include "longdouble.h"
53
54 static long double neg(long double, int *);
55 static long double poly(long double, const long double *, int);
56 static long double polytail(long double);
57 static long double primary(long double);
58
59 static const long double
60 c0 = 0.0L,
61 ch = 0.5L,
62 c1 = 1.0L,
63 c2 = 2.0L,
64 c3 = 3.0L,
65 c4 = 4.0L,
66 c5 = 5.0L,
67 c6 = 6.0L,
68 pi = 3.1415926535897932384626433832795028841971L,
69 tiny = 1.0e-40L;
70
71 long double
|