cp's OEIS Frontend

This is a front-end for the Online Encyclopedia of Integer Sequences, made by Christian Perfect. The idea is to provide OEIS entries in non-ancient HTML, and then to think about how they're presented visually. The source code is on GitHub.

Showing 1-8 of 8 results.

A324370 Product of all primes p not dividing n such that the sum of the base-p digits of n is at least p, or 1 if no such prime exists.

Original entry on oeis.org

1, 1, 2, 1, 6, 1, 6, 3, 10, 1, 6, 1, 210, 15, 2, 3, 30, 5, 210, 21, 110, 15, 30, 5, 546, 21, 14, 1, 30, 1, 462, 231, 1190, 105, 6, 1, 51870, 1365, 70, 21, 2310, 55, 2310, 105, 322, 105, 210, 35, 6630, 663, 286, 33, 330, 55, 798, 57, 290, 15, 30, 1, 930930, 15015, 1430, 2145, 1122, 85, 82110, 2415, 70, 3, 330, 55, 21111090, 285285
Offset: 1

Views

Author

Keywords

Comments

The product is finite, as the sum of the base-p digits of n is n if p > n.
a(198) = 2465 is the only term below 10^6 that is a Carmichael number (A002997).
It appears that a(n)=1 if and only if n is in A094960. - Robert Israel, Mar 30 2020
It turns out that a(n) equals the denominator of the first derivative of the Bernoulli polynomial B(n,x). So a(n)=1 if and only if n is in A094960, also impyling that n+1 is prime. A324370 is also involved in such formulas regarding higher derivatives. See Kellner 2023. - Bernd C. Kellner, Oct 12 2023

Examples

			For p = 2, 3, and 5, the sum of the base p digits of 7 is 1+1+1 = 3 >= 2, 2+1 = 3 >= 3, and 1+2 = 3 < 5, respectively, so a(7) = 2*3 = 6.
		

Crossrefs

Programs

  • Maple
    N:= 100: # for a(1)..a(N)
    V:= Vector(N,1):
    p:= 1:
    for iter from 1 do
       p:= nextprime(p);
       if p >= N then break fi;
       for n from p+1 to N do
         if n mod p <> 0 and convert(convert(n,base,p),`+`)>= p then
           V[n]:= V[n]*p
         fi
    od od:
    convert(V,list); # Robert Israel, Mar 30 2020
    # Alternatively, note that this formula is suggesting offset 0 and a(0) = 1:
    seq(denom(diff(bernoulli(n, x), x)), n = 1..51); # Peter Luschny, Oct 13 2023
  • Mathematica
    SD[n_, p_] := If[n < 1 || p < 2, 0, Plus @@ IntegerDigits[n, p]];
    DD2[n_] := Times @@ Select[Prime[Range[PrimePi[(n+1)/(2+Mod[n+1, 2])]]], !Divisible[n, #] && SD[n, #] >= # &];
    Table[DD2[n], {n, 1, 100}]
    (* From Bernd C. Kellner, Oct 12 2023 (Start) *)
    (* Denominator of first derivative of BP *)
    k = 1; Table[Denominator[Together[D[BernoulliB[n, x], {x, k}]]], {n, 1, 100}]
    (* End *)
  • Python
    from math import prod
    from sympy.ntheory import digits
    from sympy import primefactors, primerange
    def a(n):
        nonpf = set(primerange(1, n+1)) - set(primefactors(n))
        return prod(p for p in nonpf if sum(digits(n, p)[1:]) >= p)
    print([a(n) for n in range(1, 75)]) # Michael S. Branicky, Jul 03 2022

Formula

a(n) * A324369(n) = A195441(n-1) = denominator(Bernoulli_n(x) - Bernoulli_n).
a(n) * A324369(n) * A324371(n) = A144845(n-1) = denominator(Bernoulli_{n-1}(x)).
a(n+1) = A195441(n)/A324369(n+1) = A144845(n)/A007947(n+1) = A318256(n). Essentially the same as A318256. - Peter Luschny, Mar 05 2019
From Bernd C. Kellner, Oct 12 2023: (Start)
a(n) = denominator(Bernoulli_n(x)').
k-th derivative: let (n)_m be the falling factorial.
For n > k, a(n-k+1)/gcd(a(n-k+1), (n)_{k-1}) = denominator(Bernoulli_n(x)^(k)). Otherwise, the denominator equals 1. (End)

A195441 a(n) = denominator(Bernoulli_{n+1}(x) - Bernoulli_{n+1}).

Original entry on oeis.org

1, 1, 2, 1, 6, 2, 6, 3, 10, 2, 6, 2, 210, 30, 6, 3, 30, 10, 210, 42, 330, 30, 30, 30, 546, 42, 14, 2, 30, 2, 462, 231, 3570, 210, 6, 2, 51870, 2730, 210, 42, 2310, 330, 2310, 210, 4830, 210, 210, 210, 6630, 1326, 858, 66, 330, 110, 798, 114, 870, 30, 30, 6
Offset: 0

Views

Author

Peter Luschny, Sep 18 2011

Keywords

Comments

If s(n) is the smallest number such that s(n)*(1^n + 2^n + ... + x^n) is a polynomial in x with integer coefficients then a(n)=s(n)/(n+1) (see A064538).
a(n) is squarefree, by the von Staudt-Clausen theorem on the denominators of Bernoulli numbers. - Kieren MacMillan and Jonathan Sondow, Nov 20 2015
Kellner and Sondow give a detailed analysis of this sequence and provide a simple way to compute the terms without using Bernoulli polynomials and numbers. They prove that a(n) is the product of the primes less than or equal to (n+2)/(2+(n mod 2)) such that the sum of digits of n+1 in base p is at least p. - Peter Luschny, May 14 2017
The equation a(n-1) = denominator(Bernoulli_n(x) - Bernoulli_n) = rad(n+1) has only finitely many solutions, where rad(n) = A007947(n) is the radical of n. It is conjectured that S = {3, 5, 8, 9, 11, 27, 29, 35, 59} is the full set of all such solutions. Note that (S\{8})+1 joined with {1,2} equals A094960. More precisely, the set S implies the finite sequence of A094960. See Kellner 2023. - Bernd C. Kellner, Oct 18 2023
As was observed in the example section of A318256: denominator(B_n(x)) = rad(n+1) if n is in {0, 1, 3, 5, 9, 11, 27, 29, 35, 59} = {A094960(n) - 1: 1 <= n <= 10}. - Peter Luschny, Oct 18 2023

Crossrefs

Programs

  • Julia
    using Nemo, Primes
    function A195441(n::Int)
        n < 4 && return ZZ([1,1,2,1][n+1])
        P = primes(2, div(n+2, 2+n%2))
        prod([ZZ(p) for p in P if p <= sum(digits(n+1, base=p))])
    end
    println([A195441(n) for n in 0:59]) # Peter Luschny, May 14 2017
    
  • Maple
    A195441 := n -> denom(bernoulli(n+1, x)-bernoulli(n+1)):
    seq(A195441(i),i=0..59);
    # Formula of Kellner and Sondow:
    a := proc(n) local s; s := (p,n) -> add(i,i=convert(n,base,p));
    select(isprime,[$2..(n+2)/(2+irem(n,2))]); mul(i,i=select(p->s(p,n+1)>=p,%)) end: seq(a(n), n=0..59); # Peter Luschny, May 14 2017
  • Mathematica
    a[n_] := Denominator[Together[(BernoulliB[n + 1, x] - BernoulliB[n + 1])]]; Table[a[n], {n, 0, 59}] (* Jonathan Sondow, Nov 20 2015 *)
    SD[n_, p_] := If[n < 1 || p < 2, 0, Plus @@ IntegerDigits[n, p]]; DD[n_] := Times @@ Select[Prime[Range[PrimePi[(n+2)/(2+Mod[n, 2])]]], SD[n+1, #] >= # &]; Table[DD[n], {n, 0, 59}] (* Bernd C. Kellner, Oct 18 2023 *)
  • PARI
    a(n) = {my(vp = Vec(bernpol(n+1, x)-bernfrac(n+1))); lcm(vector(#vp, k, denominator(vp[k])));} \\ Michel Marcus, Feb 08 2016
    
  • Python
    from math import prod
    from sympy.ntheory.factor_ import primerange, digits
    def A195441(n): return prod(p for p in primerange((n+2)//(2|n&1)+1) if sum(digits(n+1,p)[1:])>=p) # Chai Wah Wu, Oct 04 2023
  • Sage
    A195441 = lambda n: mul([p for p in (2..(n+2)//(2+n%2)) if is_prime(p) and sum((n+1).digits(base=p))>=p])
    print([A195441(n) for n in (0..59)]) # Peter Luschny, May 14 2017
    

Formula

a(n) = A064538(n)/(n+1). - Jonathan Sondow, Nov 12 2015
A001221(a(n)) = A001222(a(n)). - Kieren MacMillan and Jonathan Sondow, Nov 20 2015
a(2*n)/a(2*n+1) = A286516(n+1). - Bernd C. Kellner and Jonathan Sondow, May 24 2017
a(n) = A007947(A338025(n+1)). - Harald Hofstätter, Oct 10 2020
From Bernd C. Kellner, Oct 18 2023: (Start)
Note that the formulas here are shifted in index by 1 due to the definition of a(n) using index n+1!
a(n) = A324369(n+1) * A324370(n+1).
a(n) = A144845(n) / A324371(n+1).
a(n-1) = lcm(a(n), rad(n+1)), if n >= 3 is odd.
If n+1 is composite, then rad(n+1) divides a(n-1).
If m is a Carmichael number (A002997), then m divides both a(m-1) and a(m-2).
See papers of Kellner and Kellner & Sondow. (End)

Extensions

Definition simplified by Jonathan Sondow, Nov 20 2015

A366168 Denominator of the second derivative of the n-th Bernoulli polynomial B(n,x).

Original entry on oeis.org

1, 1, 1, 1, 1, 1, 1, 3, 1, 1, 1, 1, 1, 15, 1, 1, 3, 5, 5, 21, 1, 5, 15, 5, 1, 21, 7, 1, 1, 1, 1, 231, 7, 35, 3, 1, 1, 1365, 35, 7, 21, 55, 55, 105, 7, 7, 105, 35, 5, 663, 13, 11, 33, 55, 1, 57, 1, 5, 15, 1, 1, 15015, 715, 715, 33, 17, 85, 2415, 35, 1, 3, 55, 55, 285285, 19019, 1001
Offset: 1

Views

Author

Bernd C. Kellner, Oct 02 2023

Keywords

Comments

The sequence consists only of odd numbers. The denominators are connected with A324370, from which an explicit formula follows as given below. See Kellner 2023.

Examples

			B(5,x) = x^5 - (5x^4)/2 + (5 x^3)/3 - x/6 and B''(5,x) = 20x^3 - 30x^2 + 10x, so a(5) = 1.
a(14) = A324370(13)/gcd(A324370(13), 14) = 210/gcd(210, 14) = 15.
		

Crossrefs

Programs

  • Mathematica
    (* k-th derivative of BP *)
    k = 2; Table[Denominator[Together[D[BernoulliB[n, x], {x, k}]]], {n, 1, 100}]
    (* 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; Table[DBP[n, k], {n, 1, 100}]
  • Python
    from math import lcm
    from sympy import Poly, diff, bernoulli
    from sympy.abc import x
    def A366168(n): return lcm(*(c.q for c in Poly(diff(bernoulli(n,x),x,2)).coeffs())) if n>=3 else 1 # Chai Wah Wu, Oct 04 2023

Formula

Let (n)_k be the falling factorial. Let s_p(n) be the sum of the p-adic digits of n.
a(1) = 1, and for n > 1, a(n) = A324370(n-1)/gcd(A324370(n-1), n) = Product_{prime p <= n/(2+(n mod 2)): gcd(p,(n)_2)=1, s_p(n-1) >= p} p.

A366169 Positive integers k such that the second derivative of the k-th Bernoulli polynomial B(k,x) contains only integer coefficients.

Original entry on oeis.org

1, 2, 3, 4, 5, 6, 7, 9, 10, 11, 12, 13, 15, 16, 21, 25, 28, 29, 30, 31, 36, 37, 55, 57, 60, 61, 70, 121, 190
Offset: 1

Views

Author

Bernd C. Kellner, Oct 02 2023

Keywords

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.

Examples

			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.
		

Crossrefs

Programs

  • 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&]
  • 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))
    A366169_list = list(islice(A366169_gen(),20)) # Chai Wah Wu, Oct 03 2023

Formula

k is a term if A366168(k) = 1.

A366186 Positive integers k such that the third derivative of the k-th Bernoulli polynomial B(k, x) contains only integer coefficients.

Original entry on oeis.org

1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 20, 21, 22, 25, 26, 28, 29, 30, 31, 32, 35, 36, 37, 38, 42, 50, 52, 55, 56, 57, 58, 60, 61, 62, 66, 70, 71, 72, 78, 80, 92, 110, 121, 122, 156, 176, 177, 190, 191, 210, 392
Offset: 1

Views

Author

Peter Luschny, Oct 03 2023

Keywords

Comments

From Bernd C. Kellner, Oct 04 2023: (Start)
As a published result on Oct 02 2023 (cf. A366169), all such sequences regarding higher derivatives of the Bernoulli polynomials having only integer coefficients are finite. We have an infinite chain of subsets: A094960 subset of A366169 subset of A366186 subset of A366187 subset of A366188 subset of ... . See Kellner 2023 (Theorem 13, Conjecture 14, and S_3 (this sequence)).
The sequence is finite and is a supersequence of A366169. It remains to show that 392 is the last term. This is very likely, since the terms depend on the estimate 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. (End)

Crossrefs

Cf. A094960 (m=1), A366169 (m=2), this sequence (m=3), A366187 (m=4), A366188 (m=5), A366189.

Programs

  • Maple
    aList := len -> select(n -> denom(diff(bernoulli(n, x), `$`(x, 3))) = 1, [seq(1..len)]): aList(400);
  • Mathematica
    (* From Bernd C. Kellner, Oct 04 2023 (Start) *)
    (* k-th derivative of BP *)
    k = 3; 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 = 3; Select[Range[1000], DBP[#, k] == 1&]
    (* End *)
  • PARI
    isok(k) = #select(x->denominator(x)>1, Vec(deriv(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 A366186_gen(): # generator of terms
        return filter(lambda k:k<=3 or all(c.is_integer for c in Poly(diff(bernoulli(k,x),x,3)).coeffs()),count(1))
    A366186_list = list(islice(A366186_gen(),30)) # Chai Wah Wu, Oct 03 2023

Formula

From Bernd C. Kellner, Oct 04 2023: (Start)
Let (n)_m be the falling factorial. Let s_p(n) be the sum of the p-adic digits of n.
The denominator of the third derivative of the n-th Bernoulli polynomial B(n, x) is given as follows (Kellner 2023, Theorem 12).
D_3(n) = 1 for 1 <= n <= 3. For n > 3, D_3(n) = A324370(n-2)/gcd(A324370(n-2), (n)2) = Product{prime p <= (n-1)/(2+((n-1) mod 2)): gcd(p,(n)_3)=1, s_p(n-2) >= p} p.
Then k is a term if and only if D_3(k) = 1. (End)

A366187 Positive integers k such that the fourth derivative of the k-th Bernoulli polynomial B(k, x) contains only integer coefficients.

Original entry on oeis.org

1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 25, 26, 27, 28, 29, 30, 31, 32, 33, 35, 36, 37, 38, 39, 42, 43, 45, 50, 51, 52, 53, 55, 56, 57, 58, 59, 60, 61, 62, 63, 66, 67, 68, 70, 71, 72, 73, 78, 79, 80, 81, 91, 92, 93, 110, 111, 121, 122, 123, 143, 147, 156, 157, 171, 176, 177, 178, 190, 191, 192, 210, 211, 255, 392, 393
Offset: 1

Views

Author

Peter Luschny, Oct 03 2023

Keywords

Comments

From Bernd C. Kellner, Oct 04 2023: (Start)
As a published result on Oct 02 2023 (cf. A366169), all such sequences regarding higher derivatives of the Bernoulli polynomials having only integer coefficients are finite. We have an infinite chain of subsets: A094960 subset of A366169 subset of A366186 subset of A366187 subset of A366188 subset of ... . See Kellner 2023 (Theorem 13).
The sequence is finite and is a supersequence of A366186. It remains to show that 393 is the last term. This is very likely, since the terms depend on the estimate 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. (End)

Crossrefs

Cf. A094960 (m=1), A366169 (m=2), A366186 (m=3), this sequence (m=4), A366188 (m=5), A366189.

Programs

  • Maple
    aList := len -> select(n -> denom(diff(bernoulli(n, x), `$`(x, 4))) = 1, [seq(1..len)]): aList(400);
  • Mathematica
    (* From Bernd C. Kellner, Oct 04 2023 (Start) *)
    (* k-th derivative of BP *)
    k = 4; 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 = 4; Select[Range[1000], DBP[#, k] == 1&]
    (* End *)
  • PARI
    isok(k) = #select(x->denominator(x)>1, Vec(deriv(deriv(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 A366187_gen(): # generator of terms
        return filter(lambda k:k<=4 or all(c.is_integer for c in Poly(diff(bernoulli(k,x),x,4)).coeffs()),count(1))
    A366187_list = list(islice(A366187_gen(),40)) # Chai Wah Wu, Oct 03 2023

Formula

From Bernd C. Kellner, Oct 04 2023: (Start)
Let (n)_m be the falling factorial. Let s_p(n) be the sum of the p-adic digits of n.
The denominator of the fourth derivative of the n-th Bernoulli polynomial B(n, x) is given as follows (Kellner 2023, Theorem 12).
D_4(n) = 1 for 1 <= n <= 4. For n > 4, D_4(n) = A324370(n-3)/gcd(A324370(n-3), (n)3) = Product{prime p <= (n-2)/(2+((n-2) mod 2)): gcd(p,(n)_4)=1, s_p(n-3) >= p} p.
Then k is a term if and only if D_4(k) = 1. (End)

A366188 Positive integers k such that the fifth derivative of the k-th Bernoulli polynomial B(k, x) contains only integer coefficients.

Original entry on oeis.org

1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 42, 43, 44, 45, 46, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 66, 67, 68, 69, 70, 71
Offset: 1

Views

Author

Peter Luschny, Oct 03 2023

Keywords

Comments

From Bernd C. Kellner, Oct 04 2023: (Start)
As a published result on Oct 02 2023 (cf. A366169), all such sequences regarding higher derivatives of the Bernoulli polynomials having only integer coefficients are finite. We have an infinite chain of subsets: A094960 subset of A366169 subset of A366186 subset of A366187 subset of A366188 subset of ... . See Kellner 2023 (Theorem 13).
The sequence is finite and is a supersequence of A366187. It remains to show that 904 is the last term. This is very likely, since the terms depend on the estimate 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. (End)

Crossrefs

Cf. A094960 (m=1), A366169 (m=2), A366186 (m=3), A366187 (m=4), (this sequence) (m=5), A366189.

Programs

  • Maple
    aList := len -> select(n -> denom(diff(bernoulli(n, x), `$`(x, 5))) = 1, [seq(1..len)]): aList(1000);
  • Mathematica
    (* From Bernd C. Kellner, Oct 04 2023 (Start) *)
    (* k-th derivative of BP *)
    k = 5; 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 = 5; Select[Range[1000], DBP[#, k] == 1&]
    (* End *)
  • PARI
    isok(k) = #select(x->denominator(x)>1, Vec(deriv(deriv(deriv(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 A366188_gen(): # generator of terms
        return filter(lambda k:k<=5 or all(c.is_integer for c in Poly(diff(bernoulli(k,x),x,5)).coeffs()),count(1))
    A366188_list = list(islice(A366188_gen(),30)) # Chai Wah Wu, Oct 03 2023

Formula

From Bernd C. Kellner, Oct 04 2023: (Start)
Let (n)_m be the falling factorial. Let s_p(n) be the sum of the p-adic digits of n.
The denominator of the fifth derivative of the n-th Bernoulli polynomial B(n, x) is given as follows (Kellner 2023, Theorem 12).
D_5(n) = 1 for 1 <= n <= 5. For n > 5, D_5(n) = A324370(n-4)/gcd(A324370(n-4), (n)4) = Product{prime p <= (n-3)/(2+((n-3) mod 2)): gcd(p,(n)_5)=1, s_p(n-4) >= p} p.
Then k is a term if and only if D_5(k) = 1. (End)

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.

Original entry on oeis.org

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, 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, 5, 2, 3, 3, 4, 5, 6, 7, 6, 3, 4, 3, 4, 5, 6, 5, 6, 7
Offset: 1

Views

Author

Peter Luschny, Oct 03 2023

Keywords

Comments

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.
Conjecture: Every integer appears in this sequence only a finite number of times. (This generalizes the conjectures made in A366186-A366188.)

Examples

			B = Bernoulli(8, x).
B = -(1/30) + (2/3)*x^2 - (7/3)*x^4 + (14/3)*x^6 - 4*x^7 + x^8;
B' = (4/3)*x - (28/3)*x^3 + 28*x^5 - 28*x^6 + 8*x^7;
B'' = (4/3) - 28*x^2 + 140*x^4 - 168*x^5 + 56*x^6;
B''' = -56*x + 560*x^3 - 840*x^4 + 336*x^5.
Thus the integral index of B is a(8) = 3.
		

Crossrefs

Bernoulli polynomials: A196838/A196839 (with rising powers).
Cf. A094960 (m=1), A366169 (m=2), A366186 (m=3), A366187 (m=4), A366188 (m=5).

Programs

  • Maple
    aList := proc(len) local n, k, d, A;
        A := Array([seq(0, n = 0..len-1)]);
        for n from 1 to len do
           k := 0: d:= 0;
           while d <> 1 do
              k := k + 1;
              d := denom(diff(bernoulli(n, x), `$`(x, k)));
           od;
           A[n] := k;
        od;
    convert(A, list) end:
    aList(86);
  • Mathematica
    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]]]];
    Table[a[n], {n, 1, 100}] (* Jean-François Alcover, Oct 23 2023 *)
  • Python
    from itertools import count
    from sympy import Poly, bernoulli, diff
    from sympy.abc import x
    def A366189(n):
        p = Poly(bernoulli(n,x))
        for i in count(1):
            p = diff(p)
            if all(c.is_integer for c in p.coeffs()):
                return i # Chai Wah Wu, Oct 03 2023
    
  • SageMath
    def A366189List(len):
        A = [0 for _ in range(len)]
        P. = ZZ[]
        for n in range(len):
            ber = bernoulli_polynomial(x, n + 1)
            k = 0
            while True:
                k = k + 1
                ber = diff(ber, x)
                if ber.denominator() == 1:
                    A[n] = k; break
        return A
    print(A366189List(86))  # Peter Luschny, Oct 04 2023
Showing 1-8 of 8 results.