cp's OEIS Frontend

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

%I A335949 #17 Jan 11 2025 13:47:37
%S A335949 1,1,12,4,240,48,1344,192,3840,1280,33792,3072,5591040,430080,245760,
%T A335949 49152,16711680,983040,522977280,27525120,1211105280,173015040,
%U A335949 1447034880,62914560,22900899840,4580179968,1409286144,469762048,116769423360,4026531840,7689065201664
%N A335949 a(n) = denominator(b_n(x)), where b_n(x) are the polynomials defined in A335947.
%C A335949 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 A335949 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 A335949 a(n) = min {m | m*([x^k] b(n, x)) is an integer for all k = 0..n}.
%F A335949 The odd part of a(n) is squarefree (A000265).
%F A335949 a(n) and A144845(n) have the same odd prime factors.
%F A335949 a(n)/A144845(n) = 4^floor(n/2)/2 for n >= 1.
%F A335949 a(n)/rad(a(n)) = A158302(n+1), (rad=A007947).
%o A335949 (SageMath)
%o A335949 def A335949(n):
%o A335949     a = set(prime_divisors(n + 1)) - set([2])
%o A335949     b = (
%o A335949         p
%o A335949         for p in prime_range(3, (n + 2) // (2 + n % 2))
%o A335949         if not p.divides(n + 1) and sum((n + 1).digits(base=p)) >= p
%o A335949     )
%o A335949     p = list(a.union(set(b)))
%o A335949     return 4 ^ (n // 2) * mul(p)
%o A335949 print([A335949(n) for n in range(31)])
%Y A335949 Cf. A335947/A335948, A144845, A158302.
%K A335949 nonn
%O A335949 0,3
%A A335949 _Peter Luschny_, Jul 01 2020