hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-05 17:31:38 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

add arb_power_sum_vec

Browse source
This commit is contained in:
Fredrik Johansson 2015-08-04 12:22:38 +02:00
parent a38bc984f8
commit 38a116157f
4 changed files with 178 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
arb.h
View file

@ -616,6 +616,8 @@ void arb_zeta_ui(arb_t z, ulong n, long prec);
void arb_hurwitz_zeta(arb_t z, const arb_t s, const arb_t a, long prec);
void arb_bernoulli_ui(arb_t z, ulong n, long prec);
void arb_power_sum_vec(arb_ptr res, const arb_t a, const arb_t b, long len, long prec);
void arb_rising_ui_bs(arb_t y, const arb_t x, ulong n, long prec);
void arb_rising_ui_rs(arb_t y, const arb_t x, ulong n, ulong m, long prec);
void arb_rising_ui_rec(arb_t y, const arb_t x, ulong n, long prec);

68
arb/power_sum_vec.c Normal file
View file

@ -0,0 +1,68 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2015 Fredrik Johansson
******************************************************************************/
#include "arb.h"
#include "bernoulli.h"
#include "arb_poly.h"
/* todo: don't use exact bernoulli numbers for large len */
/* todo: output exact integers when precise enough */
void
arb_power_sum_vec(arb_ptr res, const arb_t a, const arb_t b, long len, long prec)
{
arb_ptr t, u, v;
long k;
if (len < 1)
return;
t = _arb_vec_init(len + 1);
u = _arb_vec_init(len + 1);
v = _arb_vec_init(len + 1);
/* exp(ax), exp(bx) */
arb_set(t + 1, a);
arb_set(u + 1, b);
_arb_poly_exp_series(t, t, 2, len + 1, prec);
_arb_poly_exp_series(u, u, 2, len + 1, prec);
_arb_vec_sub(t, u, t, len + 1, prec);
/* x/(exp(x)-1) */
BERNOULLI_ENSURE_CACHED(len + 1);
for (k = 0; k <= len; k++)
arb_set_fmpq(u + k, bernoulli_cache + k, prec);
_arb_poly_borel_transform(u, u, len + 1, prec);
_arb_poly_mullow(v, t, len + 1, u, len + 1, len + 1, prec);
_arb_poly_inv_borel_transform(v, v, len + 1, prec);
for (k = 0; k < len; k++)
arb_div_ui(res + k, v + k + 1, k + 1, prec);
_arb_vec_clear(t, len + 1);
_arb_vec_clear(u, len + 1);
_arb_vec_clear(v, len + 1);
}

91
arb/test/t-power_sum_vec.c Normal file
View file

@ -0,0 +1,91 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2015 Fredrik Johansson
******************************************************************************/
#include "arb.h"
int main()
{
long iter;
flint_rand_t state;
printf("power_sum_vec....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 1000; iter++)
{
arb_t a, b, s, t;
arb_ptr res;
long aa, bb, k, n, len;
long prec;
len = n_randint(state, 30);
prec = 2 + n_randint(state, 500);
aa = n_randint(state, 50) - 50;
bb = aa + n_randint(state, 50);
arb_init(a);
arb_init(b);
arb_init(s);
arb_init(t);
res = _arb_vec_init(len);
arb_set_si(a, aa);
arb_set_si(b, bb);
arb_power_sum_vec(res, a, b, len, prec);
for (n = 0; n < len; n++)
{
arb_zero(s);
for (k = aa; k < bb; k++)
{
arb_set_si(t, k);
arb_pow_ui(t, t, n, prec);
arb_add(s, s, t, prec);
}
if (!arb_overlaps(res + n, s))
{
printf("FAIL: overlap\n\n");
printf("a = %ld, b = %ld, n = %ld\n\n", aa, bb, n);
printf("res = "); arb_printd(res + n, 30); printf("\n\n");
printf("s = "); arb_printd(s, 30); printf("\n\n");
abort();
}
}
arb_clear(a);
arb_clear(b);
arb_clear(s);
arb_clear(t);
_arb_vec_clear(res, len);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

17
doc/source/arb.rst
View file

@ -1162,6 +1162,23 @@ Bernoulli numbers
enough and `n` large enough for the Euler product to converge
rapidly (otherwise this function will effectively hang).
.. function:: void arb_power_sum_vec(arb_ptr res, const arb_t a, const arb_t b, long len, long prec)
For *n* from 0 to *len* - 1, sets entry *n* in the output vector *res* to
.. math ::
S_n(a,b) = \frac{1}{n+1}\left(B_{n+1}(b) - B_{n+1}(a)\right)
where `B_n(x)` is a Bernoulli polynomial. If *a* and *b* are integers
and `b \ge a`, this is equivalent to
.. math ::
S_n(a,b) = \sum_{k=a}^{b-1} k^n.
The computation uses the generating function for Bernoulli polynomials.
Polylogarithms
-------------------------------------------------------------------------------