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 "arb.h"
13
14 void
arb_bernoulli_poly_ui(arb_t res,ulong n,const arb_t x,slong prec)15 arb_bernoulli_poly_ui(arb_t res, ulong n, const arb_t x, slong prec)
16 {
17 arb_t s, t, c, x2;
18 ulong m, k;
19 int negate;
20
21 if (n == 0)
22 {
23 arb_one(res);
24 return;
25 }
26
27 if (n == 1)
28 {
29 arb_mul_2exp_si(res, x, 1);
30 arb_sub_ui(res, res, 1, prec);
31 arb_mul_2exp_si(res, res, -1);
32 return;
33 }
34
35 /* small integer x */
36 if (arb_is_int(x) && arf_cmpabs_ui(arb_midref(x), n) < 0 && n < WORD_MAX)
37 {
38 if (arf_sgn(arb_midref(x)) >= 0)
39 {
40 m = arf_get_si(arb_midref(x), ARF_RND_DOWN);
41 negate = 0;
42 }
43 else
44 {
45 m = UWORD(1) - arf_get_si(arb_midref(x), ARF_RND_DOWN);
46 negate = n % 2;
47 }
48
49 arb_init(t);
50 arb_zero(res);
51
52 /* todo: use a dedicated power sum function */
53 for (k = 1; k < m; k++)
54 {
55 arb_ui_pow_ui(t, k, n - 1, prec);
56 arb_add(res, res, t, prec);
57 }
58
59 arb_mul_ui(res, res, n, prec);
60 arb_bernoulli_ui(t, n, prec);
61 arb_add(res, res, t, prec);
62 if (negate)
63 arb_neg(res, res);
64
65 arb_clear(t);
66 return;
67 }
68
69 /* assuming small n simplifies the code that follows */
70 if (n >> (FLINT_BITS / 2) || !arb_is_finite(x))
71 {
72 arb_indeterminate(res);
73 return;
74 }
75
76 arb_init(s);
77 arb_init(t);
78 arb_init(c);
79 arb_init(x2);
80
81 arb_mul(x2, x, x, prec);
82
83 /* s = x^2 - x n / 2 */
84 arb_mul_ui(s, x, n, prec);
85 arb_mul_2exp_si(s, s, -1);
86 arb_sub(s, x2, s, prec);
87
88 /* c = n (n-1) / 2; s = s + c / 6 */
89 arb_set_ui(c, n * (n - 1));
90 arb_mul_2exp_si(c, c, -1);
91 arb_div_ui(t, c, 6, prec);
92 arb_add(s, s, t, prec);
93
94 for (k = 4; k <= n; k += 2)
95 {
96 /* c = binomial(n,k) */
97 arb_mul_ui(c, c, (n + 1 - k) * (n + 2 - k), prec);
98 arb_div_ui(c, c, k * (k - 1), prec);
99
100 /* s = s x^2 + b_k c */
101 arb_mul(s, s, x2, prec);
102 arb_bernoulli_ui(t, k, prec);
103 arb_addmul(s, t, c, prec);
104 }
105
106 if (n >= 3 && n % 2)
107 arb_mul(s, s, x, prec);
108
109 arb_swap(res, s);
110
111 arb_clear(s);
112 arb_clear(t);
113 arb_clear(c);
114 arb_clear(x2);
115 }
116
117