cp's OEIS Frontend

A366189 a(n) is the positive integer k such that the k-th derivative of the n-th Bernoulli polynomial B(n, x) contains only integer coefficients but no lower derivative of B(n, x) has this property.

%I A366189 #28 Oct 23 2023 07:17:06
%S A366189 1,1,2,1,2,1,2,3,2,1,2,1,2,3,2,2,3,3,4,3,2,3,4,5,2,3,4,1,2,1,2,3,4,5,
%T A366189 3,1,2,3,4,5,6,3,4,5,4,5,6,7,5,3,4,3,4,5,2,3,2,3,4,1,2,3,4,5,6,3,4,4,
%U A366189 5,2,3,3,4,5,6,7,6,3,4,3,4,5,6,5,6,7
%N A366189 a(n) is the positive integer k such that the k-th derivative of the n-th Bernoulli polynomial B(n, x) contains only integer coefficients but no lower derivative of B(n, x) has this property.
%C A366189 The 'integral index' k of a rational polynomial p(x) is the smallest integer k such that p^[k](x) is an integer polynomial, where p^[k](x) denotes the k-th derivative of p. (Integer polynomials have integral index 0.) Using this way of speaking, the a(n) are the integral indices of the Bernoulli polynomials.
%C A366189 Conjecture: Every integer appears in this sequence only a finite number of times. (This generalizes the conjectures made in A366186-A366188.)
%H A366189 Peter Luschny, <a href="/A366189/b366189.txt">Table of n, a(n) for n = 1..4000</a>
%e A366189 B = Bernoulli(8, x).
%e A366189 B = -(1/30) + (2/3)*x^2 - (7/3)*x^4 + (14/3)*x^6 - 4*x^7 + x^8;
%e A366189 B' = (4/3)*x - (28/3)*x^3 + 28*x^5 - 28*x^6 + 8*x^7;
%e A366189 B'' = (4/3) - 28*x^2 + 140*x^4 - 168*x^5 + 56*x^6;
%e A366189 B''' = -56*x + 560*x^3 - 840*x^4 + 336*x^5.
%e A366189 Thus the integral index of B is a(8) = 3.
%p A366189 aList := proc(len) local n, k, d, A;
%p A366189     A := Array([seq(0, n = 0..len-1)]);
%p A366189     for n from 1 to len do
%p A366189        k := 0: d:= 0;
%p A366189        while d <> 1 do
%p A366189           k := k + 1;
%p A366189           d := denom(diff(bernoulli(n, x), `$`(x, k)));
%p A366189        od;
%p A366189        A[n] := k;
%p A366189     od;
%p A366189 convert(A, list) end:
%p A366189 aList(86);
%t A366189 a[n_] := Module[{b, k, x}, b = BernoulliB[n, x]; For[k = 1, True, k++, b = D[b, x]; If[AllTrue[CoefficientList[b, x], IntegerQ], Return[k]]]];
%t A366189 Table[a[n], {n, 1, 100}] (* _Jean-François Alcover_, Oct 23 2023 *)
%o A366189 (Python)
%o A366189 from itertools import count
%o A366189 from sympy import Poly, bernoulli, diff
%o A366189 from sympy.abc import x
%o A366189 def A366189(n):
%o A366189     p = Poly(bernoulli(n,x))
%o A366189     for i in count(1):
%o A366189         p = diff(p)
%o A366189         if all(c.is_integer for c in p.coeffs()):
%o A366189             return i # _Chai Wah Wu_, Oct 03 2023
%o A366189 (SageMath)
%o A366189 def A366189List(len):
%o A366189     A = [0 for _ in range(len)]
%o A366189     P.<x> = ZZ[]
%o A366189     for n in range(len):
%o A366189         ber = bernoulli_polynomial(x, n + 1)
%o A366189         k = 0
%o A366189         while True:
%o A366189             k = k + 1
%o A366189             ber = diff(ber, x)
%o A366189             if ber.denominator() == 1:
%o A366189                 A[n] = k; break
%o A366189     return A
%o A366189 print(A366189List(86))  # _Peter Luschny_, Oct 04 2023
%Y A366189 Bernoulli polynomials: A196838/A196839 (with rising powers).
%Y A366189 Cf. A094960 (m=1), A366169 (m=2), A366186 (m=3), A366187 (m=4), A366188 (m=5).
%K A366189 nonn
%O A366189 1,3
%A A366189 _Peter Luschny_, Oct 03 2023