OFFSET
0,1
COMMENTS
Note also that Sum_{i=1..m} i^n is a polynomial in m of degree n+1.
One can prove that all terms of the sequence are integers.
The sequence appears to possess an astonishing property: every odd prime p is the maximal prime divisor of a((p-1)/2).
EXAMPLE
Let n=1. Since lcmd(Sum_{i=1..m} i ) = 2, lcmd(Sum_{i=1..m} i^2) = 6, lcmd(Sum_{i=1..m} i^3) = 4, then lcmd(Sum_{i=1..m} i*Sum_{j=1..i} j ) = lcmd(Sum_{i=1..m} (i^2*(i+1)/2) = 24, therefore, a(1) = 24/4 = 6.
Let p=53. Then a(26) = 5830 = 2*5*11*53 has maximal prime divisor 53.
MATHEMATICA
LCMD[P_, m_] := LCM @@ Denominator[CoefficientList[P // FunctionExpand, m] ]; a[n_] := LCMD[Sum[i^n*HarmonicNumber[i, -n], {i, 1, m}], m]/ LCMD[ HarmonicNumber[m, -n]^2, m]; Table[an = a[n]; Print["a(", n, ") = ", an]; an, {n, 0, 60}] (* Jean-François Alcover, Feb 18 2016 *)
PROG
(PARI) sp(p) = x * Polrev(vector(p+1, k, (-1)^(k==p)*binomial(p+1, k)*bernfrac(p+1-k))/(p+1));
lcmd(pol) = lcm(apply(x->denominator(x), Vec(pol)));
a(n) = {pol = x^n*sp(n); pnum = sum(i=0, poldegree(pol), polcoeff(pol, i)*sp(i)); lcmd(pnum)/lcmd(sp(n)^2); } \\ Michel Marcus, Feb 17 2016
CROSSREFS
KEYWORD
nonn
AUTHOR
Vladimir Shevelev and Peter J. C. Moses, Dec 20 2011
EXTENSIONS
More terms from Michel Marcus, Feb 17 2016
STATUS
approved