1 /*
2     Copyright (C) 2021 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_hypgeom.h"
13 #include "acb_hypgeom.h"
14 
15 static void
evaluate_rect(acb_t res,const short * term_prec,slong len,const acb_t x,slong prec)16 evaluate_rect(acb_t res, const short * term_prec, slong len, const acb_t x, slong prec)
17 {
18     slong i, j, m, r, n1, n2;
19     acb_ptr xs;
20     acb_t s, t;
21     arb_struct c[17];
22 
23     m = n_sqrt(len) + 1;
24     m = FLINT_MIN(m, 16);
25     r = (len + m - 1) / m;
26 
27     xs = _acb_vec_init(m + 1);
28     acb_init(s);
29     acb_init(t);
30     _acb_vec_set_powers(xs, x, m + 1, prec);
31 
32     acb_zero(res);
33 
34     for (i = r - 1; i >= 0; i--)
35     {
36         n1 = m * i;
37         n2 = FLINT_MIN(len, n1 + m);
38 
39         for (j = n1; j < n2; j++)
40         {
41             if (j == 0)
42             {
43                 arb_init(c);
44                 arb_one(c);
45             }
46             else
47             {
48                 if (!_arb_hypgeom_gamma_coeff_shallow(arb_midref(c + j - n1), arb_radref(c + j - n1), j, term_prec[j]))
49                     flint_abort();
50             }
51         }
52 
53         arb_dot(acb_realref(s), NULL, 0, acb_realref(xs), 2, c, 1, n2 - n1, prec);
54         arb_dot(acb_imagref(s), NULL, 0, acb_imagref(xs), 2, c, 1, n2 - n1, prec);
55 
56 #if 0
57         acb_set_round(t, xs + m, term_prec[n1]);
58         acb_mul(res, res, t, term_prec[n1]);
59         acb_add(res, res, s, term_prec[n1]);
60 #else
61         acb_mul(res, res, xs + m, term_prec[n1]);
62         acb_add(res, res, s, term_prec[n1]);
63 #endif
64     }
65 
66     _acb_vec_clear(xs, m + 1);
67     acb_clear(s);
68     acb_clear(t);
69 }
70 
71 /* Bound requires: |u| <= 20, N <= 10000, N != (1443, 2005, 9891). */
72 static void
error_bound(mag_t err,const acb_t u,slong N)73 error_bound(mag_t err, const acb_t u, slong N)
74 {
75     mag_t t;
76     mag_init(t);
77 
78     acb_get_mag(t, u);
79 
80     if (N >= 1443 || mag_cmp_2exp_si(t, 4) > 0)
81     {
82         mag_inf(err);
83     }
84     else
85     {
86         mag_pow_ui(err, t, N);
87         mag_mul_2exp_si(err, err, arb_hypgeom_gamma_coeffs[N].exp);
88 
89         if (mag_cmp_2exp_si(t, -1) > 0)
90             mag_mul(err, err, t);
91         else
92             mag_mul_2exp_si(err, err, -1);
93 
94         mag_mul_2exp_si(err, err, 3);
95 
96         if (mag_cmp_2exp_si(err, -8) > 0)
97             mag_inf(err);
98     }
99 
100     mag_clear(t);
101 }
102 
103 static double
want_taylor(double x,double y,slong prec)104 want_taylor(double x, double y, slong prec)
105 {
106     if (y < 0.0) y = -y;
107     if (x < 0.0) x = -x;
108 
109     if ((prec < 128 && y > 4.0) || (prec < 256 && y > 5.0) ||
110         (prec < 512 && y > 8.0) || (prec < 1024 && y > 9.0) || y > 10.0)
111     {
112         return 0;
113     }
114 
115     if (x * (1.0 + 0.75 * y) > 8 + 0.15 * prec)
116     {
117         return 0;
118     }
119 
120     return 1;
121 }
122 
123 int
acb_hypgeom_gamma_taylor(acb_t res,const acb_t z,int reciprocal,slong prec)124 acb_hypgeom_gamma_taylor(acb_t res, const acb_t z, int reciprocal, slong prec)
125 {
126     acb_t s, u;
127     int success;
128     double dua, dub, du2, log2u;
129     slong i, r, n, wp, tail_bound, goal;
130     short term_prec[ARB_HYPGEOM_GAMMA_TAB_NUM];
131     mag_t err;
132 
133     if (!acb_is_finite(z) ||
134         arf_cmp_2exp_si(arb_midref(acb_imagref(z)), 4) >= 0 ||
135         arf_cmp_2exp_si(arb_midref(acb_realref(z)), 10) >= 0)
136     {
137         return 0;
138     }
139 
140     dua = arf_get_d(arb_midref(acb_realref(z)), ARF_RND_UP);
141     dub = arf_get_d(arb_midref(acb_imagref(z)), ARF_RND_UP);
142     dub = fabs(dub);
143 
144     if (!want_taylor(dua, dub, prec))
145         return 0;
146 
147     if (dua >= 0.0)
148         r = (slong) (dua + 0.5);
149     else
150         r = -(slong) (-dua + 0.5);
151 
152     acb_init(s);
153     acb_init(u);
154     mag_init(err);
155 
156     success = 0;
157 
158     /* Argument reduction: u = z - r */
159     acb_sub_si(u, z, r, 2 * prec + 10);
160     dua -= r;
161 
162     goal = acb_rel_accuracy_bits(u);
163 
164     /* not designed for wide intervals (yet) */
165     if (goal < 8)
166     {
167         success = 0;
168         goto cleanup;
169     }
170 
171     goal = FLINT_MIN(goal, prec - MAG_BITS) + MAG_BITS;
172     goal = FLINT_MAX(goal, 5);
173     goal = goal + 5;
174     wp = goal + 4 + FLINT_BIT_COUNT(FLINT_ABS(r));
175 
176     if (wp > ARB_HYPGEOM_GAMMA_TAB_PREC)
177     {
178         success = 0;
179         goto cleanup;
180     }
181 
182     if (!want_taylor(r, dub, goal))
183     {
184         success = 0;
185         goto cleanup;
186     }
187 
188     du2 = dua * dua + dub * dub;
189 
190     if (du2 > 1e-8)
191     {
192         log2u = 0.5 * mag_d_log_upper_bound(du2) * 1.4426950408889634074 * (1 + 1e-14);
193     }
194     else
195     {
196         slong aexp, bexp;
197 
198         aexp = arf_cmpabs_2exp_si(arb_midref(acb_realref(u)), -wp) >= 0 ? ARF_EXP(arb_midref(acb_realref(u))) : -wp;
199         bexp = arf_cmpabs_2exp_si(arb_midref(acb_imagref(u)), -wp) >= 0 ? ARF_EXP(arb_midref(acb_imagref(u))) : -wp;
200         log2u = FLINT_MAX(aexp, bexp) + 1;
201     }
202 
203     term_prec[0] = wp;
204     n = 0;
205 
206     for (i = 1; i < ARB_HYPGEOM_GAMMA_TAB_NUM; i++)
207     {
208         tail_bound = arb_hypgeom_gamma_coeffs[i].exp + i * log2u + 5;
209 
210         if (tail_bound <= -goal)
211         {
212             n = i;
213             break;
214         }
215 
216         term_prec[i] = FLINT_MIN(FLINT_MAX(wp + tail_bound, 2), wp);
217 
218         if (term_prec[i] > arb_hypgeom_gamma_coeffs[i].nlimbs * FLINT_BITS)
219         {
220             success = 0;
221             goto cleanup;
222         }
223     }
224 
225     if (n != 0)
226         error_bound(err, u, n);
227 
228     if (n == 0 || mag_is_inf(err))
229     {
230         success = 0;
231         goto cleanup;
232     }
233 
234     evaluate_rect(s, term_prec, n, u, wp);
235     acb_add_error_mag(s, err);
236 
237     if (r == 0 || r == 1)
238     {
239         if (r == 0)
240             acb_mul(s, s, u, wp);
241 
242         if (reciprocal)
243         {
244             acb_set_round(res, s, prec);
245         }
246         else
247         {
248             acb_one(u);
249             acb_div(res, u, s, prec);
250         }
251     }
252     else if (r >= 2)
253     {
254         acb_add_ui(u, u, 1, wp);
255         acb_hypgeom_rising_ui_rec(u, u, r - 1, wp);
256 
257         if (reciprocal)
258             acb_div(res, s, u, prec);
259         else
260             acb_div(res, u, s, prec);
261     }
262     else
263     {
264         /* gamma(x) = (-1)^r / (rgamma(1+x-r)*rf(1+r-x,-r)*(x-r)) */
265         /* 1/gamma(x) = (-1)^r * rgamma(1+x-r) * rf(1+r-x,-r) * (x-r) */
266 
267         acb_neg(res, z);
268         acb_add_si(res, res, 1 + r, wp);
269         acb_hypgeom_rising_ui_rec(res, res, -r, wp);
270         acb_mul(u, res, u, wp);
271 
272         if (reciprocal)
273         {
274             acb_mul(res, s, u, prec);
275         }
276         else
277         {
278             acb_mul(u, s, u, wp);
279             acb_inv(res, u, prec);
280         }
281 
282         if (r % 2)
283             acb_neg(res, res);
284     }
285 
286     success = 1;
287 
288 cleanup:
289     acb_clear(s);
290     acb_clear(u);
291     mag_clear(err);
292 
293     return success;
294 }
295 
296