OFFSET
0,10
COMMENTS
REFERENCES
Johann Faulhaber, Academia Algebra. Darinnen die miraculosische Inventiones zu den höchsten Cossen weiters continuirt und profitiert werden. Johann Ulrich Schönigs, Augsburg, 1631.
LINKS
C. G. J. Jacobi, De usu legitimo formulae summatoriae Maclaurinianae, J. Reine Angew. Math., 12 (1834), 263-272.
Donald E. Knuth, Johann Faulhaber and sums of powers, arXiv:math/9207222 [math.CA], 1992; Math. Comp. 61 (1993), no. 203, 277-294.
Peter Luschny, Illustrating the Faulhaber polynomials for n = 1..7.
FORMULA
Let F_n(x) be the polynomial after substituting (sqrt(8*x + 1) - 1)/2 for x in b_n(x), where b_n(x) = (Bernoulli_{2*n)(x+1) - Bernoulli_{2*n)(1))/(2*n).
F_n(1) = 1 for all n >= 0.
T(n, k) = numerator([x^k] F_n(x)).
EXAMPLE
The first few polynomials are:
[0] 1;
[1] x;
[2] x^2;
[3] (4*x - 1)*x^2*(1/3);
[4] (6*x^2 - 4*x + 1)*x^2*(1/3);
[5] (16*x^3 - 20*x^2 + 12*x - 3)*x^2*(1/5);
[6] (16*x^4 - 32*x^3 + 34*x^2 - 20*x + 5)*x^2*(1/3);
[7] (960*x^5 - 2800*x^4 + 4592*x^3 - 4720*x^2 + 2764*x - 691)*x^2*(1/105);
[8] (48*x^6 - 192*x^5 + 448*x^4 - 704*x^3 + 718*x^2 - 420*x + 105)*x^2*(1/3);
[9] (1280*x^7-6720*x^6+21120*x^5-46880*x^4+72912*x^3-74220*x^2+43404*x-10851)*x^2*(1/45);
.
Triangle starts:
[0] 1;
[1] 0, 1;
[2] 0, 0, 1;
[3] 0, 0, -1, 4;
[4] 0, 0, 1, -4, 6;
[5] 0, 0, -3, 12, -20, 16;
[6] 0, 0, 5, -20, 34, -32, 16;
[7] 0, 0, -691, 2764, -4720, 4592, -2800, 960;
[8] 0, 0, 105, -420, 718, -704, 448, -192, 48;
[9] 0, 0, -10851, 43404, -74220, 72912, -46880, 21120, -6720, 1280;
MAPLE
FaulhaberPolynomial := proc(n) if n = 0 then return 1 fi;
expand((bernoulli(2*n, x+1) - bernoulli(2*n, 1))/(2*n));
sort(simplify(expand(subs(x = (sqrt(8*x+1)-1)/2, %))), [x], ascending) end:
Trow := n -> seq(coeff(numer(FaulhaberPolynomial(n)), x, k), k=0..n):
seq(print(Trow(n)), n=0..9);
PROG
(Python)
from math import lcm
from itertools import count, islice
from sympy import simplify, sqrt, bernoulli
from sympy.abc import x
def A335951_T(n, k):
z = simplify((bernoulli(2*n, (sqrt(8*x+1)+1)/2)-bernoulli(2*n, 1))/(2*n)).as_poly().all_coeffs()
return z[n-k]*lcm(*(d.q for d in z))
def A335951_gen(): # generator of terms
yield from (A335951_T(n, k) for n in count(0) for k in range(n+1))
(SageMath)
def A335951Row(n):
R.<x> = PolynomialRing(QQ)
if n == 0: return [1]
b = expand((bernoulli_polynomial(x + 1, 2*n) -
bernoulli_polynomial(1, 2*n))/(2*n))
s = expand(b.subs(x = (sqrt(8*x+1)-1)/2))
return numerator(s).list()
for n in range(10): print(A335951Row(n)) # Peter Luschny, May 17 2022
CROSSREFS
KEYWORD
AUTHOR
Peter Luschny, Jul 16 2020
STATUS
approved