1 /*
2 Copyright (C) 2016 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 "acb.h"
13
14 void
acb_bernoulli_poly_ui(acb_t res,ulong n,const acb_t x,slong prec)15 acb_bernoulli_poly_ui(acb_t res, ulong n, const acb_t x, slong prec)
16 {
17 acb_t s, x2;
18 arb_t t, c;
19 ulong k;
20
21 if (n == 0)
22 {
23 acb_one(res);
24 return;
25 }
26
27 if (n == 1)
28 {
29 acb_mul_2exp_si(res, x, 1);
30 acb_sub_ui(res, res, 1, prec);
31 acb_mul_2exp_si(res, res, -1);
32 return;
33 }
34
35 if (acb_is_real(x))
36 {
37 arb_bernoulli_poly_ui(acb_realref(res), n, acb_realref(x), prec);
38 arb_zero(acb_imagref(res));
39 return;
40 }
41
42 /* assuming small n simplifies the code that follows */
43 if (n >> (FLINT_BITS / 2) || !acb_is_finite(x))
44 {
45 acb_indeterminate(res);
46 return;
47 }
48
49 acb_init(s);
50 acb_init(x2);
51 arb_init(t);
52 arb_init(c);
53
54 acb_mul(x2, x, x, prec);
55
56 /* s = x^2 - x n / 2 */
57 acb_mul_ui(s, x, n, prec);
58 acb_mul_2exp_si(s, s, -1);
59 acb_sub(s, x2, s, prec);
60
61 /* c = n (n-1) / 2; s = s + c / 6 */
62 arb_set_ui(c, n * (n - 1));
63 arb_mul_2exp_si(c, c, -1);
64 arb_div_ui(t, c, 6, prec);
65 acb_add_arb(s, s, t, prec);
66
67 for (k = 4; k <= n; k += 2)
68 {
69 /* c = binomial(n,k) */
70 arb_mul_ui(c, c, (n + 1 - k) * (n + 2 - k), prec);
71 arb_div_ui(c, c, k * (k - 1), prec);
72
73 /* s = s x^2 + b_k c */
74 acb_mul(s, s, x2, prec);
75 arb_bernoulli_ui(t, k, prec);
76 arb_mul(t, t, c, prec);
77 acb_add_arb(s, s, t, prec);
78 }
79
80 if (n >= 3 && n % 2)
81 acb_mul(s, s, x, prec);
82
83 acb_swap(res, s);
84
85 acb_clear(s);
86 acb_clear(x2);
87 arb_clear(t);
88 arb_clear(c);
89 }
90
91