login
a(n) = denominator(b_n(x)), where b_n(x) are the polynomials defined in A335947.
2

%I #13 Jul 02 2020 13:05:11

%S 1,1,12,4,240,48,1344,192,3840,1280,33792,3072,5591040,430080,245760,

%T 49152,16711680,983040,522977280,27525120,1211105280,173015040,

%U 1447034880,62914560,22900899840,4580179968,1409286144,469762048,116769423360,4026531840,7689065201664

%N a(n) = denominator(b_n(x)), where b_n(x) are the polynomials defined in A335947.

%C The sequence can also be computed without reference to the Bernoulli polynomials (ultimately thanks to the von Staudt-Clausen theorem) by the method of Kellner and Sondow (2019). Compare the SageMath program.

%H Bernd C. Kellner and Jonathan Sondow, <a href="https://arxiv.org/abs/1902.10672">On Carmichael and polygonal numbers, Bernoulli polynomials, and sums of base-p digits</a>, arXiv:1902.10672 [math.NT], 2019.

%F a(n) = min {m | m*([x^k] b(n, x)) is integer for all k = 0..n}.

%F The odd part of a(n) is squarefree (A000265).

%F a(n) and A144845(n) have the same odd prime factors.

%F a(n)/A144845(n) = 4^floor(n/2)/2 for n >= 1.

%F a(n)/rad(a(n)) = A158302(n+1), (rad=A007947).

%o (SageMath)

%o def A335949(n):

%o a = set(prime_divisors(n + 1)) - set([2])

%o b = (

%o p

%o for p in prime_range(3, (n + 2) // (2 + n % 2))

%o if not p.divides(n + 1) and sum((n + 1).digits(base=p)) >= p

%o )

%o p = list(a.union(set(b)))

%o return 4 ^ (n // 2) * mul(p)

%o print([A335949(n) for n in range(31)])

%Y Cf. A335947/A335948, A144845, A158302.

%K nonn

%O 0,3

%A _Peter Luschny_, Jul 01 2020