1 /*
2 Copyright (C) 2015 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
14 static void
bsplit(arb_t P,arb_t Q,const fmpz_t n,const fmpz_t a,const fmpz_t b,slong prec)15 bsplit(arb_t P, arb_t Q, const fmpz_t n, const fmpz_t a, const fmpz_t b, slong prec)
16 {
17 fmpz_t t;
18 fmpz_init(t);
19 fmpz_sub(t, b, a);
20
21 if (fmpz_sgn(t) <= 0)
22 {
23 arb_zero(P);
24 arb_one(Q);
25 }
26 else if (fmpz_cmp_ui(t, 20) < 0)
27 {
28 slong steps, k;
29 arb_t u;
30 arb_init(u);
31
32 arb_zero(P);
33 arb_one(Q);
34
35 steps = fmpz_get_si(t);
36
37 for (k = steps - 1; k >= 0; k--)
38 {
39 fmpz_add_ui(t, a, k);
40
41 arb_set_round_fmpz(u, t, prec);
42 arb_pow_fmpz(u, u, n, prec);
43 arb_addmul(P, Q, u, prec);
44
45 if (!fmpz_is_zero(t))
46 arb_mul_fmpz(Q, Q, t, prec);
47 }
48
49 arb_clear(u);
50 }
51 else
52 {
53 fmpz_t m;
54 arb_t P1, Q2;
55
56 fmpz_init(m);
57 arb_init(P1);
58 arb_init(Q2);
59
60 fmpz_add(m, a, b);
61 fmpz_tdiv_q_2exp(m, m, 1);
62
63 bsplit(P1, Q, n, a, m, prec);
64 bsplit(P, Q2, n, m, b, prec);
65
66 arb_mul(Q, Q, Q2, prec);
67 arb_addmul(P, P1, Q2, prec);
68
69 fmpz_clear(m);
70 arb_clear(P1);
71 arb_clear(Q2);
72 }
73
74 fmpz_clear(t);
75 }
76
77 void
arb_bell_sum_bsplit(arb_t res,const fmpz_t n,const fmpz_t a,const fmpz_t b,const fmpz_t mmag,slong prec)78 arb_bell_sum_bsplit(arb_t res, const fmpz_t n,
79 const fmpz_t a, const fmpz_t b, const fmpz_t mmag, slong prec)
80 {
81 if (fmpz_cmp(a, b) >= 0)
82 {
83 arb_zero(res);
84 }
85 else
86 {
87 slong wp;
88 arb_t P, Q;
89
90 wp = _fmpz_sub_small(b, a);
91 wp = FLINT_BIT_COUNT(FLINT_ABS(wp));
92 wp = prec + fmpz_bits(n) + fmpz_bits(a) + wp;
93
94 arb_init(P);
95 arb_init(Q);
96
97 bsplit(P, Q, n, a, b, wp);
98 arb_div(res, P, Q, wp);
99
100 if (!fmpz_is_zero(a))
101 {
102 arb_gamma_fmpz(P, a, wp);
103 arb_div(res, res, P, wp);
104 }
105
106 arb_set_round(res, res, prec);
107
108 arb_clear(P);
109 arb_clear(Q);
110 }
111 }
112
113