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 void
acb_hypgeom_rising_ui_rs(acb_t res,const acb_t x,ulong n,ulong m,slong prec)16 acb_hypgeom_rising_ui_rs(acb_t res, const acb_t x, ulong n, ulong m, slong prec)
17 {
18 slong i, k, l, m0, climbs, climbs_max, wp;
19 acb_ptr xpow;
20 acb_t t, u;
21 mp_ptr c;
22 TMP_INIT;
23
24 if (n <= 1)
25 {
26 if (n == 0)
27 acb_one(res);
28 else
29 acb_set_round(res, x, prec);
30 return;
31 }
32
33 TMP_START;
34
35 if (m == 0 || m == -1)
36 {
37 if (n <= 6)
38 m = 2;
39 else if (n <= 16)
40 m = 4;
41 else if (n <= 40)
42 m = 6;
43 else
44 {
45 m0 = n_sqrt(n);
46 m = 8 + 0.27 * pow(1.5 * FLINT_MAX(0, prec - 1024), 0.4);
47 m = FLINT_MIN(m, m0);
48 m = FLINT_MIN(m, 64);
49 }
50 }
51
52 wp = ARF_PREC_ADD(prec, FLINT_BIT_COUNT(n));
53
54 climbs_max = FLINT_BIT_COUNT(n - 1) * m;
55 c = TMP_ALLOC(sizeof(mp_limb_t) * climbs_max * m);
56
57 xpow = _acb_vec_init(m + 1);
58 _acb_vec_set_powers(xpow, x, m + 1, wp);
59 acb_init(t);
60 acb_init(u);
61
62 for (k = 0; k < n; k += m)
63 {
64 l = FLINT_MIN(m, n - k);
65 climbs = FLINT_BIT_COUNT(k + l - 1) * l;
66 climbs = (climbs + FLINT_BITS - 1) / FLINT_BITS;
67
68 /* assumes l >= 2 */
69 if (l == 1)
70 {
71 acb_add_ui(u, x, k, wp);
72 }
73 else
74 {
75 if (climbs == 1)
76 {
77 _arb_hypgeom_rising_coeffs_1(c, k, l);
78 acb_dot_ui(u, xpow + l, 0, xpow, 1, c, 1, l, wp);
79 }
80 else if (climbs == 2)
81 {
82 _arb_hypgeom_rising_coeffs_2(c, k, l);
83 acb_dot_uiui(u, xpow + l, 0, xpow, 1, c, 1, l, wp);
84 }
85 else
86 {
87 fmpz * f = (fmpz *) c;
88
89 for (i = 0; i < l; i++)
90 fmpz_init(f + i);
91
92 _arb_hypgeom_rising_coeffs_fmpz(f, k, l);
93
94 acb_dot_fmpz(u, xpow + l, 0, xpow, 1, f, 1, l, wp);
95
96 for (i = 0; i < l; i++)
97 fmpz_clear(f + i);
98 }
99 }
100
101 if (k == 0)
102 acb_swap(t, u);
103 else
104 acb_mul(t, t, u, wp);
105 }
106
107 acb_set_round(res, t, prec);
108
109 acb_clear(t);
110 acb_clear(u);
111 _acb_vec_clear(xpow, m + 1);
112 TMP_END;
113 }
114
115