1 /*
2 Copyright (C) 2013 Fredrik Johansson
3
4 This file is part of Arb.
5
6 Arb is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include "bernoulli.h"
13 #include "acb.h"
14 #include "acb_hypgeom.h"
15
16 void
17 acb_gamma_stirling_choose_param(int * reflect, slong * r, slong * n,
18 const acb_t z, int use_reflect, int digamma, slong prec);
19
20 void acb_gamma_stirling_bound(mag_ptr err, const acb_t z, slong k0, slong knum, slong n);
21
22 void arb_gamma_stirling_coeff(arb_t b, ulong k, int digamma, slong prec);
23
24
25 void
acb_gamma_stirling_eval(acb_t s,const acb_t z,slong nterms,int digamma,slong prec)26 acb_gamma_stirling_eval(acb_t s, const acb_t z, slong nterms, int digamma, slong prec)
27 {
28 acb_t t, logz, zinv, zinv2;
29 arb_t b;
30 mag_t err;
31
32 slong k, term_prec;
33 double z_mag, term_mag;
34
35 acb_init(t);
36 acb_init(logz);
37 acb_init(zinv);
38 acb_init(zinv2);
39 arb_init(b);
40
41 acb_log(logz, z, prec);
42 acb_inv(zinv, z, prec);
43
44 nterms = FLINT_MAX(nterms, 1);
45
46 acb_zero(s);
47 if (nterms > 1)
48 {
49 acb_mul(zinv2, zinv, zinv, prec);
50
51 z_mag = arf_get_d(arb_midref(acb_realref(logz)), ARF_RND_UP) * 1.44269504088896;
52
53 for (k = nterms - 1; k >= 1; k--)
54 {
55 term_mag = bernoulli_bound_2exp_si(2 * k);
56 term_mag -= (2 * k - 1) * z_mag;
57 term_prec = prec + term_mag;
58 term_prec = FLINT_MIN(term_prec, prec);
59 term_prec = FLINT_MAX(term_prec, 10);
60
61 arb_gamma_stirling_coeff(b, k, digamma, term_prec);
62
63 if (prec > 2000)
64 {
65 acb_set_round(t, zinv2, term_prec);
66 acb_mul(s, s, t, term_prec);
67 }
68 else
69 acb_mul(s, s, zinv2, term_prec);
70
71 arb_add(acb_realref(s), acb_realref(s), b, term_prec);
72 }
73
74 if (digamma)
75 acb_mul(s, s, zinv2, prec);
76 else
77 acb_mul(s, s, zinv, prec);
78 }
79
80 /* remainder bound */
81 mag_init(err);
82 acb_gamma_stirling_bound(err, z, digamma ? 1 : 0, 1, nterms);
83 mag_add(arb_radref(acb_realref(s)), arb_radref(acb_realref(s)), err);
84 mag_add(arb_radref(acb_imagref(s)), arb_radref(acb_imagref(s)), err);
85 mag_clear(err);
86
87 if (digamma)
88 {
89 acb_neg(s, s);
90 acb_mul_2exp_si(zinv, zinv, -1);
91 acb_sub(s, s, zinv, prec);
92 acb_add(s, s, logz, prec);
93 }
94 else
95 {
96 /* (z-0.5)*log(z) - z + log(2*pi)/2 */
97 arb_one(b);
98 arb_mul_2exp_si(b, b, -1);
99 arb_set(acb_imagref(t), acb_imagref(z));
100 arb_sub(acb_realref(t), acb_realref(z), b, prec);
101 acb_mul(t, logz, t, prec);
102 acb_add(s, s, t, prec);
103 acb_sub(s, s, z, prec);
104 arb_const_log_sqrt2pi(b, prec);
105 arb_add(acb_realref(s), acb_realref(s), b, prec);
106 }
107
108 acb_clear(t);
109 acb_clear(logz);
110 acb_clear(zinv);
111 acb_clear(zinv2);
112 arb_clear(b);
113 }
114
115 void
acb_gamma(acb_t y,const acb_t x,slong prec)116 acb_gamma(acb_t y, const acb_t x, slong prec)
117 {
118 acb_hypgeom_gamma(y, x, prec);
119 }
120
121 void
acb_rgamma(acb_t y,const acb_t x,slong prec)122 acb_rgamma(acb_t y, const acb_t x, slong prec)
123 {
124 acb_hypgeom_rgamma(y, x, prec);
125 }
126
127 void
acb_lgamma(acb_t y,const acb_t x,slong prec)128 acb_lgamma(acb_t y, const acb_t x, slong prec)
129 {
130 acb_hypgeom_lgamma(y, x, prec);
131 }
132