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 "arb.h"
13 #include "arb_hypgeom.h"
14 #include "bernoulli.h"
15
16 /* todo: move/cleanup helper functions */
17
18 void
acb_gamma_bound_phase(mag_t bound,const acb_t z)19 acb_gamma_bound_phase(mag_t bound, const acb_t z)
20 {
21 arf_t x, y, t, u;
22 int xsign;
23 slong prec;
24
25 arf_init(x);
26 arf_init(y);
27 arf_init(t);
28 arf_init(u);
29
30 prec = MAG_BITS;
31
32 /* first compute x, y such that |arg(z)| <= arg(x+yi) */
33
34 /* argument increases with smaller real parts */
35 arf_set_mag(x, arb_radref(acb_realref(z)));
36 arf_sub(x, arb_midref(acb_realref(z)), x, prec, ARF_RND_FLOOR);
37
38 xsign = arf_sgn(x);
39
40 if (xsign >= 0) /* argument increases away from the real axis */
41 arb_get_abs_ubound_arf(y, acb_imagref(z), prec);
42 else /* argument increases closer to the real axis */
43 arb_get_abs_lbound_arf(y, acb_imagref(z), prec);
44
45 if (arf_is_zero(y))
46 {
47 if (xsign > 0)
48 mag_one(bound);
49 else
50 mag_inf(bound);
51 }
52 else
53 {
54 if (xsign >= 0)
55 {
56 /* compute upper bound for t = y / (sqrt(x^2 + y^2) + x) */
57 arf_mul(t, x, x, prec, ARF_RND_DOWN);
58 arf_mul(u, y, y, prec, ARF_RND_DOWN);
59 arf_add(t, t, u, prec, ARF_RND_DOWN);
60 arf_sqrt(t, t, prec, ARF_RND_DOWN);
61 arf_add(t, t, x, prec, ARF_RND_DOWN);
62 arf_div(t, y, t, prec, ARF_RND_UP);
63 }
64 else
65 {
66 /* compute upper bound for t = (sqrt(x^2 + y^2) - x) / y */
67 arf_mul(t, x, x, prec, ARF_RND_UP);
68 arf_mul(u, y, y, prec, ARF_RND_UP);
69 arf_add(t, t, u, prec, ARF_RND_UP);
70 arf_sqrt(t, t, prec, ARF_RND_UP);
71 arf_sub(t, t, x, prec, ARF_RND_UP);
72 arf_div(t, t, y, prec, ARF_RND_UP);
73 }
74
75 /* compute upper bound for sqrt(1 + t^2) */
76 arf_mul(t, t, t, prec, ARF_RND_UP);
77 arf_add_ui(t, t, 1, prec, ARF_RND_UP);
78 arf_sqrt(t, t, prec, ARF_RND_UP);
79
80 arf_get_mag(bound, t);
81 }
82
83 arf_clear(x);
84 arf_clear(y);
85 arf_clear(t);
86 arf_clear(u);
87 }
88
89 /*
90 2 |B_{2n}| G(2n+k-1) / (G(k+1) G(2n+1)) |z| (T |z|^{-1})^(2n+k)
91 TODO: CHECK n >= 1 ?
92 */
93 void
acb_gamma_stirling_bound(mag_ptr err,const acb_t z,slong k0,slong knum,slong n)94 acb_gamma_stirling_bound(mag_ptr err, const acb_t z, slong k0, slong knum, slong n)
95 {
96 mag_t c, t, u, v;
97 slong i, k;
98
99 if (arb_contains_zero(acb_imagref(z)) &&
100 arb_contains_nonpositive(acb_realref(z)))
101 {
102 for (i = 0; i < knum; i++)
103 mag_inf(err + i);
104 return;
105 }
106
107 mag_init(c);
108 mag_init(t);
109 mag_init(u);
110 mag_init(v);
111
112 /* t = lower bound for |z| */
113 acb_get_mag_lower(t, z);
114 /* v = upper bound for |z| */
115 acb_get_mag(v, z);
116
117 /* c = upper bound for 1/(cos(arg(z)/2) |z|) */
118 acb_gamma_bound_phase(c, z);
119 mag_div(c, c, t);
120
121 /* numerator: 2 B_{2n} gamma(2n+k-1) |z| */
122 mag_bernoulli_div_fac_ui(err, 2 * n);
123 mag_mul_2exp_si(err, err, 1);
124 mag_fac_ui(u, 2 * n + k0 - 2);
125 mag_mul(err, err, u);
126 mag_mul(err, err, v);
127
128 /* denominator gamma(k+1) gamma(2n+1) */
129 mag_rfac_ui(t, k0);
130 mag_mul(err, err, t);
131
132 /* multiply by c^(2n+k) */
133 mag_pow_ui(t, c, 2 * n + k0);
134 mag_mul(err, err, t);
135
136 for (i = 1; i < knum; i++)
137 {
138 /* recurrence factor: c * (2n+k-2) / k */
139 k = k0 + i;
140 mag_mul(err + i, err + i - 1, c);
141 mag_mul_ui(err + i, err + i, 2 * n + k - 2);
142 mag_div_ui(err + i, err + i, k);
143 }
144
145 mag_clear(c);
146 mag_clear(t);
147 mag_clear(u);
148 mag_clear(v);
149 }
150
151 void
arb_gamma_stirling_bound(mag_ptr err,const arb_t x,slong k0,slong knum,slong n)152 arb_gamma_stirling_bound(mag_ptr err, const arb_t x, slong k0, slong knum, slong n)
153 {
154 acb_t z;
155 acb_init(z);
156 acb_set_arb(z, x);
157 acb_gamma_stirling_bound(err, z, k0, knum, n);
158 acb_clear(z);
159 }
160
161 void
arb_gamma_stirling_coeff(arb_t b,ulong k,int digamma,slong prec)162 arb_gamma_stirling_coeff(arb_t b, ulong k, int digamma, slong prec)
163 {
164 fmpz_t d;
165 fmpz_init(d);
166
167 BERNOULLI_ENSURE_CACHED(2 * k);
168
169 arb_set_round_fmpz(b, fmpq_numref(bernoulli_cache + 2 * k), prec);
170
171 if (digamma)
172 fmpz_mul_ui(d, fmpq_denref(bernoulli_cache + 2 * k), 2 * k);
173 else
174 fmpz_mul2_uiui(d, fmpq_denref(bernoulli_cache + 2 * k), 2 * k, 2 * k - 1);
175
176 arb_div_fmpz(b, b, d, prec);
177 fmpz_clear(d);
178 }
179
180
181
182 void
arb_gamma_stirling_eval(arb_t s,const arb_t z,slong nterms,int digamma,slong prec)183 arb_gamma_stirling_eval(arb_t s, const arb_t z, slong nterms, int digamma, slong prec)
184 {
185 arb_t b, t, logz, zinv, zinv2;
186 mag_t err;
187
188 slong k, term_prec;
189 double z_mag, term_mag;
190
191 arb_init(b);
192 arb_init(t);
193 arb_init(logz);
194 arb_init(zinv);
195 arb_init(zinv2);
196
197 arb_log(logz, z, prec);
198 arb_inv(zinv, z, prec);
199
200 nterms = FLINT_MAX(nterms, 1);
201
202 arb_zero(s);
203
204 if (nterms > 1)
205 {
206 arb_mul(zinv2, zinv, zinv, prec);
207
208 z_mag = arf_get_d(arb_midref(logz), ARF_RND_UP) * 1.44269504088896;
209
210 for (k = nterms - 1; k >= 1; k--)
211 {
212 term_mag = bernoulli_bound_2exp_si(2 * k);
213 term_mag -= (2 * k - 1) * z_mag;
214 term_prec = prec + term_mag;
215 term_prec = FLINT_MIN(term_prec, prec);
216 term_prec = FLINT_MAX(term_prec, 10);
217
218 if (prec > 2000)
219 {
220 arb_set_round(t, zinv2, term_prec);
221 arb_mul(s, s, t, term_prec);
222 }
223 else
224 arb_mul(s, s, zinv2, term_prec);
225
226 arb_gamma_stirling_coeff(b, k, digamma, term_prec);
227 arb_add(s, s, b, term_prec);
228 }
229
230 if (digamma)
231 arb_mul(s, s, zinv2, prec);
232 else
233 arb_mul(s, s, zinv, prec);
234 }
235
236 /* remainder bound */
237 mag_init(err);
238 arb_gamma_stirling_bound(err, z, digamma ? 1 : 0, 1, nterms);
239 mag_add(arb_radref(s), arb_radref(s), err);
240 mag_clear(err);
241
242 if (digamma)
243 {
244 arb_neg(s, s);
245 arb_mul_2exp_si(zinv, zinv, -1);
246 arb_sub(s, s, zinv, prec);
247 arb_add(s, s, logz, prec);
248 }
249 else
250 {
251 /* (z-0.5)*log(z) - z + log(2*pi)/2 */
252 arb_one(t);
253 arb_mul_2exp_si(t, t, -1);
254 arb_sub(t, z, t, prec);
255 arb_mul(t, logz, t, prec);
256 arb_add(s, s, t, prec);
257 arb_sub(s, s, z, prec);
258 arb_const_log_sqrt2pi(t, prec);
259 arb_add(s, s, t, prec);
260 }
261
262 arb_clear(t);
263 arb_clear(b);
264 arb_clear(zinv);
265 arb_clear(zinv2);
266 arb_clear(logz);
267 }
268
269 void
arb_gamma_fmpq(arb_t y,const fmpq_t x,slong prec)270 arb_gamma_fmpq(arb_t y, const fmpq_t x, slong prec)
271 {
272 arb_hypgeom_gamma_fmpq(y, x, prec);
273 }
274
275 void
arb_gamma_fmpz(arb_t y,const fmpz_t x,slong prec)276 arb_gamma_fmpz(arb_t y, const fmpz_t x, slong prec)
277 {
278 arb_hypgeom_gamma_fmpz(y, x, prec);
279 }
280
281 void
arb_gamma(arb_t y,const arb_t x,slong prec)282 arb_gamma(arb_t y, const arb_t x, slong prec)
283 {
284 arb_hypgeom_gamma(y, x, prec);
285 }
286
287 void
arb_rgamma(arb_t y,const arb_t x,slong prec)288 arb_rgamma(arb_t y, const arb_t x, slong prec)
289 {
290 arb_hypgeom_rgamma(y, x, prec);
291 }
292
293 void
arb_lgamma(arb_t y,const arb_t x,slong prec)294 arb_lgamma(arb_t y, const arb_t x, slong prec)
295 {
296 arb_hypgeom_lgamma(y, x, prec);
297 }
298
299