User:Peter Luschny/CentralizedBernoulli

How I once found a
Guy Steele sequence

KEYWORDS: Bernoulli polynomial, Clausen numbers, Guy Steele sequence.
Concerned with sequences:  A064538, A135517, A091090, A160014, A141056.

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

$\sum _{j=0}^{k}{\frac {(-1)^{k-j}}{k+1}}{\binom {k+1}{j+1}}(j+1)^{n+1}$

then I think these are very complicated numbers. However, when I see their algorithmic description

T(n, k) = k*T(n-1, k-1) + (k+1)*T(n-1, 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 n-ten 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 n-ten 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 Staudt-Clausen 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.