%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