This site is supported by donations to The OEIS Foundation.
User:Peter Luschny/CentralizedBernoulli
From OeisWiki
How I once found a
Guy Steele sequence
KEYWORDS: Bernoulli polynomial, Clausen numbers, Guy Steele sequence.
Concerned with sequences: A064538, A135517, A091090, A160014, A141056.
Contents 
A Maple code review
When I try to understand a sequence the computational description of the sequence is for me as important as the mathematical formula. For instance when I see the formula

then I think these are very complicated numbers. However, when I see their algorithmic description
T(n, k) = k*T(n1, k1) + (k+1)*T(n1, k), base T(0, 0) = 1, boundary T(n, k) = 0 if k > n,
then I think that mathematical notation sometimes obfuscates inherent formal simplicity. Having said that you might understand why I am sometimes upset reading code which obfuscates simple math. Now my story begins.
Recently I read the Maple code of sequence A064538.
A064538 := proc(n) local t1; t1 := eval((bernoulli(n+1,m+1) − bernoulli(n+1))/(n+1)); denominator(factor(t1)) end;
I can see no reason for the local variable 't1', for 'eval' nor for 'factor'. 'denominator' is unknown to my version of Maple, perhaps it is a typo. And there is also a more subtle point. The function computes (implicit) bernoulli(n+1, m+1) − bernoulli(n+1, 0) whereas for me bernoulli(n+1, m+1) − bernoulli(n+1, 1) is clearly the right thing here. So the code reduces to:
A064538 := n > denom((bernoulli(n+1,x+1) − bernoulli(n+1,1))/(n+1)):
Now observe that the quotient (n+1) can be pulled out of the denominator becoming a factor. [You need to justify this, because it is not true in general that denom(a/(b*c)) = c*denom(a/b) for positive integers a, b, c. For example, 1 = denom(1/1) = denom(2/2) != 2*denom(2/1) = 2. Thanks.  Jonathan Sondow 22:09, 20 November 2015 (UTC)]
A064538 := n > (n+1)*denom(bernoulli(n+1,x+1)  bernoulli(n+1,1)):
This is a further simplification. But what becomes even more obvious is that this is not a natural enumeration n → (n+1). Indeed it would be more natural to start with '0' and write in a version to be published on OEIS
B064538 := n > n*denom(bernoulli(n,x+1) − bernoulli(n,1)): B064538(n) (n ≥ 0) = 0, 1, 2, 6, 4, 30, 12, 42, 24, 90, 20
Note that this is only possible due to the new form if we want n = 0 to be included in the range of arguments. The left shift would have caused a division by zero in the original formula. Clearly the wording of the original interpretation of the formula has to be adjusted.
Still I am not satisfied. The factor n can obfuscate the real meaning of the sequence in some other contexts. So I decided to submit the sequence in its more general and more flexible form.
with(numtheory): a := n > denom(B(n,x+1) − B(n,1)): a(n) (n ≥ 0) = 1, 1, 1, 2, 1, 6, 2, 6, 3, 10, 2, 6, 2, 210
Is this the end of the story?
Guy Steele's sequences
Well, I thought what makes sense for the Bernoulli polynomials should also make sense for the Euler polynomials. So let's see.
a := n > denom(euler(n,x+1) − euler(n,1)): a(n) (n ≥ 0) = 1, 1, 1, 2, 1, 2, 1, 4, 1, 2, 1, 4, 1, 2, 1
Does OEIS know this sequence? Eureka, it does! Here is what I found:
a(n) = A135517(n) = 2^(A091090(n) − 1). Author: N. J. A. Sloane based on a message from Guy Steele and D. E. Knuth, Mar 01 2008.
Our formula shows that the right offset is 0 and a(0) = 1. However, I saw no indication that the connection with the Euler polynomials is known. Certainly this is worth a comment.
So the denominator of E(n,x+1) − E(n,1) is Guy Steele's sequence GS(2,5) (conjecturally).
GS := proc(p, q, len) local a, n; a := [1]; for n from 1 to len do a := a,[0,1,a[n],a[n]+1,2*a[n],2*a[n]+1][p]; a := a,[0,1,a[n],a[n]+1,2*a[n],2*a[n]+1][q]; od end:
GS(2,5,...) > [1],1,[2],1,2,1,[4],1,2,1,4,1,2,1,[8], 1,2,1,4,1,2,1,8,1,2,1,4,1,2,1,[16],
I inserted the squared brackets to make the pattern more visible. And indeed, does this pattern not bring the Matryoshka sequences back to our mind? Perhaps a sequence of this type might be called a fractal Matryoshka sequence?
There is a lot more to find out about fractal Matryoshka sequences and GS sequences. For example fix the parameters p = 2 and q = 5 in the GS sequence but introduce new parameters x and y like this:
xGS := proc(x, y, len) local a, n; a := 1; for n from 1 to len do a := a,[0,x,a[n]+x,a[n],2*a[n]+y,2*a[n]+x][2]; a := a,[0,x,a[n]+x,a[n],2*a[n]+y,2*a[n]+x][5]; od end:
Observe:
xGS(x=0, y=1) → GS(p=1, q=6); xGS(x=1, y=0) → GS(p=2, q=5); xGS(x=1, y=1) → GS(p=2, q=6).
Does anyone know a reference to Steele's sequences? (In this case please leave an answer on the discussion page.)
Centralized Bernoulli polynomials
Let's come back to our sequence a(n) = denom(B(n,x+1) − B(n,1)). We have seen that this sequence starts:
1, 1, 1, 2, 1, 6, 2, 6, 3, 10, 2, 6, 2, 210, 30, 6, 3
However, my personal taste would prefer this sequence, which is more homogeneous.
1, 1, 1, 1, 1, 3, 1, 3, 3, 5, 1, 3, 1, 105, 15, 3, 3
Divide by 2 will not work because of the subsequence a(2^n), which is interesting by itself:
Let n>0, then a(n) is odd <=> n = 2^k for some k≥0.
a(2^n) (n ≥ 0) 1 1 1 3 3 3 * 7 * 11 3 * 5 * 11 * 13 3 * 7 * 13 * 19 * 43 3 * 7 * 13 * 19 * 29 * 37 * 43
So I will try to find my preferred sequence. First I introduce a variant of the Bernoulli polynomials, the centralized Bernoulli polynomials CB_{n}(x).
with(numtheory): CentralizedBernoulli := (n, x) > 2^n*B(n, x/2 + 1/2); alias(CB = CentralizedBernoulli);
This amounts to shifting the polynomials left into the origin (which makes the symmetry of the polynomials B(n, 1 − x) = (−1)n B(n, x) look much more natural).
Next we look at the denominator of CB_{n} (which is of course CB_{n}(1)).
seq(denom(CB(n,1)),n=0..16) = 1, 1, 3, 1, 15, 1, 21, 1
It is now crucial to get things right with Dr. Clausen's famous formula, otherwise no one will buy the new polynomials.
Der Bruch der nten Bernoullischen Zahl wird so gefunden: Man addire zu den Theilern von 2n ... 1, 2, a, a’, a”, ... , 2n die Einheit, wodurch man die Reihe Zahlen 2, 3, a + 1, a’ + 1, ... , 2n+1 bekommt. Aus dieser nimmt man bloß die Primzahlen 2, 3, p, p’ etc. und bildet den Bruch der nten Bernoullischen Zahl ... 
You can find a translation on Wikipedia's Bernoulli number page. Here is the variant for the centralized Bernoulli polynomials: Aus dieser nimmt man bloß die ungeraden Primzahlen 3, p, p’ etc... Which means that we restrict the set of factors to the odd primes.
ClausenOdd := proc(n) local S, m; S := divisors(n); S := map(i > i + 1, S); S := select(isprime, S) minus {2}; mul(m, m = S) end:
A first test looks very promising:
seq(ClausenOdd(n),n=0..32) = 1, 1, 3, 1, 15, 1, 21, 1
And indeed denom(CB(n, x+1) − CB(n, 1)) is a sequence of odd integers, no exceptions anymore.
seq(denom(CB(n, x+1)  CB(n, 1)), n = 0..32) = 1, 1, 1, 1, 1, 3, 1, 3, 3, 5, 1, 3, 1, 105, 15, 3, 3
This is of course no surprise, by the von StaudtClausen theorem the denominator of B_{n} can have at most one 2 as factor and our scaling cancelled this 2.
Anyway, this puzzle is solved, but it leaves me with another very intriguing one. What is the real meaning of this sequence?
1 1,1 1,1,1,1 3,1,1,1,1,1,3,1 3,1,1,1,3,1,3,1,15,1,3,1,1,1,1,1 3,1,3,1,1,1,3,1,3,1,3,1,15,1,3,1,105,1,3,1,3,1,1,1,3,1,3,1,3,1,3,1 15,1,3,1,15,1,3,1,15,1,3,1,3,1,3,1,3,1,1,1,5,1,3,1,15,1,7,1,15,1,3,1,105,.. ..1,3,1,3,1,3,1,15,1,3,1,5,1,3,1,15,1,3,1,15,1,3,1,165,1,3,1,15,1,3,1
New sequences for the OEIS
Based on the above considerations I might submit the following sequences to OEIS.
 Denominator of B(n,x+1) − B(n,1); B(n, x) Bernoulli polynomials.
 Denominator of CB(n,x+1) − CB(n,1); CB(n, x) = 2^n B(n, x/2+1/2);
B(n, x) Bernoulli polynomials.  Denominator of the polynomials 2^n B(n, x/2+1/2),
B(n, x) Bernoulli polynomials.  Denom(CB(n,x+1)−CB(n,1))Denom(CB(n,1))/Denom(CB(n,x+1)),
CB(n,x) = 2^n B(n,x/2+1/2), B(n,x) Bernoulli polynomials.