OFFSET
1,2
COMMENTS
The sequence is finite and is a supersequence of A094960. The terms are those numbers k where the denominator A366168(k) = 1. It remains to show that 190 is the last term. This is very likely, since the terms depend on the estimation of a product of primes satisfying certain p-adic conditions that is connected with A324370. A proven asymptotic formula related to that product implies that this sequence is finite. See Kellner 2017, 2023, and BLMS 2018.
LINKS
Olivier Bordellès, Florian Luca, Pieter Moree, and Igor E. Shparlinski, Denominators of Bernoulli polynomials, Mathematika 64 (2018), 519-541.
Bernd C. Kellner, On a product of certain primes, J. Number Theory, 179 (2017), 126-141; arXiv:1705.04303 [math.NT], 2017.
Bernd C. Kellner, On the finiteness of Bernoulli polynomials whose derivative has only integral coefficients, J. Integer Seq. 27 (2024), Article 24.2.8, 11 pp.; arXiv:2310.01325 [math.NT], 2023.
Bernd C. Kellner and Jonathan Sondow, Power-Sum Denominators, Amer. Math. Monthly, 124 (2017), 695-709; arXiv:1705.03857 [math.NT], 2017.
Bernd C. Kellner and Jonathan Sondow, The denominators of power sums of arithmetic progressions, Integers 18 (2018), #A95, 17 pp.; arXiv:1705.05331 [math.NT], 2017.
Bernd C. Kellner and Jonathan Sondow, On Carmichael and polygonal numbers, Bernoulli polynomials, and sums of base-p digits, Integers 21 (2021), #A52, 21 pp.; arXiv:1902.10672 [math.NT], 2019.
FORMULA
k is a term if A366168(k) = 1.
EXAMPLE
B(5,x) = x^5 - (5x^4)/2 + (5 x^3)/3 - x/6 and B''(5,x) = 20x^3 - 30x^2 + 10x, so 5 is a term.
MAPLE
aList := len -> select(n -> denom(diff(diff(bernoulli(n, x), x), x)) = 1, [seq(1..len)]): aList(200); # Peter Luschny, Oct 03 2023
MATHEMATICA
(* k-th derivative of BP *)
k = 2; Select[Range[1000], Denominator[Together[D[BernoulliB[#, x], {x, k}]]] == 1&]
(* exact denominator formula *)
SD[n_, p_] := If[n < 1 || p < 2, 0, Plus@@IntegerDigits[n, p]];
DBP[n_, k_] := Module[{m = n-k+1, fac = FactorialPower[n, k]}, If[n < 1 || k < 1 || n <= k, Return[1]]; Times@@Select[Prime[Range[PrimePi[(m+1)/(2 + Mod[m+1, 2])]]], !Divisible[fac, #] && SD[m, #] >= #&]];
k = 2; Select[Range[1000], DBP[#, k] == 1&]
PROG
(PARI) isok(k) = #select(x->denominator(x)>1, Vec(deriv(deriv(bernpol(k))))) == 0; \\ Michel Marcus, Oct 03 2023
(Python)
from itertools import count, islice
from sympy import Poly, diff, bernoulli
from sympy.abc import x
def A366169_gen(): # generator of terms
return filter(lambda k:k<=2 or all(c.is_integer for c in Poly(diff(bernoulli(k, x), x, 2)).coeffs()), count(1))
CROSSREFS
KEYWORD
nonn,fini,hard
AUTHOR
Bernd C. Kellner, Oct 02 2023
STATUS
approved