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.

Previous Showing 11-20 of 26 results. Next

A324405 Squarefree integers m > 1 such that if prime p divides m, then s_p(m) >= p and s_p(m) == 3 (mod p-1), where s_p(m) is the sum of the base p digits of m.

Original entry on oeis.org

3003, 3315, 5187, 7395, 8463, 14763, 19803, 26733, 31755, 47523, 50963, 58035, 62403, 88023, 105339, 106113, 123123, 139971, 152643, 157899, 166611, 178923, 183183, 191919
Offset: 1

Views

Author

Keywords

Comments

For d >= 1 define S_d = (terms m in A324315 such that s_p(m) == d (mod p-1) if prime p divides m). Then S_1 is precisely the Carmichael numbers (A002997), S_2 is A324404, S_3 is A324405, and the union of all S_d for d >= 1 is A324315.
Subsequence of the 3-Knödel numbers (A033553). Generally, for d > 1 the terms of S_d that are greater than d form a subsequence of the d-Knödel numbers.
See Kellner and Sondow 2019.

Examples

			3003 = 3*7*11*13 is squarefree and equals 11010020_3, 11520_7, 2290_11, and 14a0_13 in base p = 3, 7, 11, and 13. Then s_3(3003) = 1+1+1+2 = 5 >= 3, s_7(3003) = 1+1+5+2 = 9 >= 7, s_11(3003) = 2+2+9 = 13 >= 11, and s_13(3003) = 1+4+a = 1+4+10 = 15 >= 13. Also, s_3(3003) = 5 == 3 (mod 2), s_7(3003) = 9 == 3 (mod 6), s_11(3003) = 13 == 3 (mod 10), and s_13(3003) = 15 == 3 (mod 12), so 3003 is a member.
		

Crossrefs

Programs

  • Mathematica
    SD[n_, p_] := If[n < 1 || p < 2, 0, Plus @@ IntegerDigits[n, p]];
    LP[n_] := Transpose[FactorInteger[n]][[1]];
    TestSd[n_, d_] := (n > 1) && (d > 0) && SquareFreeQ[n] && VectorQ[LP[n], SD[n, #] >= # && Mod[SD[n, #] - d, # - 1] == 0 &];
    Select[Range[200000], TestSd[#, 3] &]

A324317 Number of primary Carmichael numbers (A324316) less than 10^n.

Original entry on oeis.org

0, 0, 0, 2, 4, 9, 19, 51, 107, 219, 417, 757, 1470, 2666, 5040, 9280, 17210, 32039, 59762, 111811, 210627, 397968
Offset: 1

Views

Author

Keywords

Comments

The number of Carmichael numbers (A002997) less than 10^n is 0, 0, 1, 7, 16, 43, 105, 255, 646, 1547, 3605, 8241, 19279, 44706, 105212, 246683, 585355, 1401644, ... (see A055553).
The terms up to a(10) are given in Table 1 of Kellner and Sondow 2019. The terms up to a(18) and related results are given in Table 1.5 of Kellner 2019.
All computations depend on Pinch's database.

Examples

			There are two primary Carmichael numbers less than 10^4, namely, 1729 and 2821, so a(4) = 2.
		

Crossrefs

Extensions

a(11)-a(18) from Amiram Eldar, Mar 01 2019
a(19) from Amiram Eldar, Dec 05 2020
a(20)-a(22) calculated using data from Claude Goutier and added by Amiram Eldar, Apr 22 2024

A324318 Number of terms in A324315 (squarefree integers m > 1 such that if prime p divides m, then the sum of the base p digits of m is at least p) less than 10^n.

Original entry on oeis.org

0, 0, 2, 57, 636, 7048, 75150, 801931, 8350039, 86361487
Offset: 1

Views

Author

Keywords

Comments

The number of squarefree integers less than 10^n is 0, 6, 61, 608, 6083, 60794, 607926, 6079291, 60792694, 607927124, ... (see A053462).

Examples

			There are two terms of A324315 less than 10^3, namely, 231 and 561, so a(3) = 2.
		

Crossrefs

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

Original entry on oeis.org

1, 2, 4, 6, 10, 12, 28, 30, 36, 60
Offset: 1

Views

Author

Benoit Cloitre, Jun 19 2004

Keywords

Comments

From Max Alekseyev, Dec 08 2011: (Start)
There are no other terms below 10^9.
k belongs to this sequence if k*binomial(k-1,m)*Bernoulli(m) is an integer for each m in 0..k-1. (End)
From Max Alekseyev, Jun 04 2012: (Start)
If for a prime p >= 3, k ends with base-p digits a,b with a+b >= p, then for m = (a+1)*(p-1), the number k*binomial(k-1,m)*Bernoulli(m) is not an integer (it contains p in the denominator). For p=3, this implies that k == 5, 7, or 8 (mod 9) are not in this sequence; for p=5, this implies that k == 9, 13, 14, 17, 18, 19, 21, 22, 23, or 24 (mod 25) are not in this sequence; and so on.
Conjecture: for every integer k > 78, there exists a prime p >= 3 such that the sum of last two base-p digits of k is at least p. This conjecture would imply that this sequence is finite and 60 is the last term. (End)
The conjecture is true for all k such that k+1 is not a prime, a power of 2, or a Giuga number (A007850). In this case, there exists a prime p >= 3 such that the base-p representation of k ends in a,p-1 with a > 0. - Max Alekseyev, Feb 16 2021
The sequence is finite and is a subsequence of A366169. The terms are those numbers k where A324370(k) = 1. It remains to show that 60 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. - Bernd C. Kellner, Oct 02 2023

Examples

			B(6,x) = x^6 - 3*x^5 + (5/2)*x^4 - (1/2)*x^2 + 1/42 so B'(6,x) contains only integer coefficients and 6 is in the sequence.
		

Crossrefs

Programs

  • Maple
    p := n -> if denom(diff(bernoulli(n, x), x)) = 1 then n else fi:
    seq(p(n), n=1..100); # Emeric Deutsch
  • Mathematica
    (* From Bernd C. Kellner, Oct 02 2023. (Start) *)
    (* k-th derivative of BP: *)
    k = 1; 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 = 1; Select[Range[1000], DBP[#, k] == 1&]
    (* End *)
  • PARI
    is_A094960(k) = !#select(x->(denominator(x)!=1), Vec(deriv(bernpol(k)))); \\ Michel Marcus, Feb 15 2021
    
  • Python
    from itertools import count, islice
    from sympy import Poly, diff, bernoulli
    from sympy.abc import x
    def A094960_gen(): # generator of terms
        return filter(lambda k:k<=1 or all(c.is_integer for c in Poly(diff(bernoulli(k,x),x)).coeffs()),count(1))
    A094960_list = list(islice(A094960_gen(),10)) # Chai Wah Wu, Oct 03 2023

Formula

k is a term if A324370(k) = 1. - Bernd C. Kellner, Oct 02 2023
k is a term <=> 0 = Sum_{j=0..k-1} k*binomial(k - 1, j) mod Clausen(j), where Clausen(n) = A160014(n, 1). - Peter Luschny, Oct 04 2023

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.

A318256 a(n) = (denominator of B(n,x)) / (the squarefree kernel of n+1), where B(n,x) is the n-th Bernoulli polynomial.

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
Offset: 0

Views

Author

Peter Luschny, Sep 12 2018

Keywords

Examples

			a(59) = 1 because there exist no number which satisfies the definition (and the product of an empty set is 1).
a(60) = 930930 because {2, 3, 5, 7, 11, 13, 31} are the only primes which satisfy the definition.
The denominator of the Bernoulli polynomial B_n(x) equals the squarefree kernel of n+1 if n is in {0, 1, 3, 5, 9, 11, 27, 29, 35, 59}. These might be the only numbers with this property.
		

Crossrefs

a(n) = A144845(n) / A007947(n+1).
Cf. A324370 (same sequence with offset 1).

Programs

  • Maple
    a := n -> denom(bernoulli(n, x)) / mul(p, p in numtheory:-factorset(n+1)):
    seq(a(n), n=0..61);
  • Mathematica
    sfk[n_] := Times @@ FactorInteger[n][[All, 1]];
    a[n_] := (BernoulliB[n, x] // Together // Denominator)/sfk[n+1];
    Table[a[n], {n, 0, 61}] (* Jean-François Alcover, Feb 14 2019 *)
  • Sage
    def A318256(n): return mul([p for p in (2..(n+2)//(2+n%2))
                    if is_prime(p)
                    and not p.divides(n+1)
                    and sum((n+1).digits(base=p)) >= p])
    print([A318256(n) for n in (0..61)])

Formula

Let Q(n) = {p <= floor((n + 2)/(2 + n mod 2)) and p is prime and p does not divide n + 1 and the sum of the digits in base p of n+1 is at least p} then a(n) = Product_{p in Q(n)} p. (See the Kellner & Sondow links.)
a(n) = denominator(Bernoulli'(n+1, x)), where ' denotes d/dx. - Peter Luschny, Oct 15 2023

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)
Previous Showing 11-20 of 26 results. Next