This site is supported by donations to The OEIS Foundation.

User:Peter Luschny/SeqTransformation

From OeisWiki
Jump to: navigation, search

A sequence transformation
and the Bernoulli numbers

KEYWORDS:Bernoulli, Euler, Tangent, Harmonic, Binomial, Swiss-Knife, Worpitzky, Akiyama–Tanigawa.

Concerned with sequences: A027641, A027642, A131689, A028246, A141056, A027760, A176276, A176277.

A sequence to sequence transformation

Given an sequence f, an integer n ≥ 0 and a formal symbol x let me define a sequence of polynomials

T := proc(f,n,x) local k,v; 
add(add((-1)^v*binomial(k,v)*f(k+1)*(x+v+1)^n, 
v=0..k), k=0..n) end:
f1
xf1+f1f2
x2f1+(2f1−2f2) x+2f3+f1−3f2
x3f1+(3f1−3f2)x2+(6f3+3f1−9f2)x−6f4+12f3+f1−7f2
x4f1+(4f1−4f2)x3+(12f3+6f1−18f2)x2 +(−24f4+48f3+4f1−28f2)x−60f4+50f3+f1−15f2+24f5

The sequence of coefficients of these polynomials, sorted in descending order, is

f1        
f1 f1f2      
f1 2f1−2f2 f1−3f2+2f3    
f1 3f1−3f2 3f1−9f2+6f3 f1−7f2+12f3−6f4  
f1 4f1−4f2 6f1−18f2+12f3 4f1−28f2+48f3−24f4 f1−15f2+50f3−60f4+24f5

It is this sequence (the triangle read by rows) which I call the transform of f. f will usually be a sequence of integers or a sequence of rationals numbers. As a first check let us see on which sequence the identity function f: x -> x is mapped.

1
1, -1
1, -2,  1
1, -3,  3,  -1
1, -4,  6,  -4,  1
1, -5, 10, -10,  5, -1

A signed Pascal triangle, no big surprise, but a good start. Now let us take some more interesting input sequences. What about f(n) the harmonic numbers?

H := proc(n) local i; add(1/i, i=1..n) end;
T(H, n, x) starts:
1
x - 1/2
x^2 - x + 1/6
x^3 - 3/2 x^2 + 1/2 x

We get the Bernoulli polynomials!

And feeding the transformation with twice the sum of the inverse powers of 2?

G := proc(n) local i; 2*add(1/2^i, i=1..n) end:
T(G, n, x) starts:
1
x - 1/2
x^2 - x
x^3 - 3/2 x^2 + 1/4

We get the Euler polynomials!

I give one more example, perhaps the most interesting one. Let

C := m -> if irem(m,4) = 0 then 0 
          else 1/((-1)^iquo(m,4)*2^iquo(m-1,2)) fi:

Interesting because the polynomials generated have integer coefficients:

               1
               x
           x^2  - 1
          x^3  - 3 x
       x^4  - 6 x^2  + 5
      x^5  - 10 x^3  + 25 x
    x^6  - 15 x^4  + 75 x^2  - 61

The Swiss knife polynomials!

If we evaluate these polynomials at 0 and 1 we get, respectively,

1, 0, -1,  0, 5, 0, -61, 0, 1385, 0, -50521, 0, 2702765, ...
1, 1, 0, -2, 0, 16, 0, -272, 0, 7936, 0, -353792, 0, ...

I am sure you know these sequences. For the purpose of reference I restyle the most important special case in math-parlance.

Let Bn(x) be the Bernoulli polynomials and Hn the harmonic numbers.

I searched in the literature but could not find this formula. So I showed the formula in the newsgroup de.sci.mathematik; however, no one could provide a reference. If you know a reference for this formula, please let me know (the proof is easy once you see the formula on the blackboard...).

The case x = 1 can also be written with the Riemann zeta function as

This formula is valid for all n ≥ 0 provided for n = 0,1 the left hand side is understood as a limit value. By the way, this formula for the zeta values shows that the definition Bn = Bn(1) is to be preferred over the definition Bn = Bn(0) if you intend to go into the number theory business.

Bernoulli numbers.

On the other hand, some concrete mathematicians refrain from this definition being fearful to introduce thereby confusion to the world. Clearly this is a very strange argument since Jakob Bernoulli himself used B1 = 1/2 as you can see from the facsimile on Wikipedia, and Jakob was a very concrete mathematician.

Bernoulli and Worpitzky numbers

But this is not the end of the story. Looking at the table given above we can equate the constant coefficients of the polynomials with the Bernoulli numbers.

B_0 = 1*H(1);
B_1 = 1*H(1)-  1*H(2);                              [Table 1]
B_2 = 1*H(1)-  3*H(2)   +2*H(3);
B_3 = 1*H(1)-  7*H(2)  +12*H(3)-   6*H(4);
B_4 = 1*H(1)- 15*H(2)  +50*H(3)-  60*H(4)+  24*H(5);
B_5 = 1*H(1)- 31*H(2) +180*H(3)- 390*H(4)+ 360*H(5)- 120*H(6);

Now this is shocking: the coefficients on the right hand side of this table are the Worpitzky numbers; however, it is not Worpitzky's representation! Worpitzky gave the following representation with J(n) = 1/n.

B_0 = 1*J(1);
B_1 = 1*J(1)-  1*J(2);                               [Table 2] 
B_2 = 1*J(1)-  3*J(2)   +2*J(3);
B_3 = 1*J(1)-  7*J(2)  +12*J(3)-   6*J(4);
B_4 = 1*J(1)- 15*J(2)  +50*J(3)-  60*J(4)+  24*J(5);
B_5 = 1*J(1)- 31*J(2) +180*J(3)- 390*J(4)+ 360*J(5)- 120*J(6);

So what happens if we expand table 1? We get yet another sum representation of the Bernoulli numbers, with even simpler coefficients!

B_0 = 1*J(1)
B_1 = 0*J(1)- 1*J(2)                                 [Table 3]
B_2 = 0*J(1)- 1*J(2)+  2*J(3)
B_3 = 0*J(1)- 1*J(2)+  6*J(3)-   6*J(4)
B_4 = 0*J(1)- 1*J(2)+ 14*J(3)-  36*J(4)+   24*J(5)
B_5 = 0*J(1)- 1*J(2)+ 30*J(3)- 150*J(4)+  240*J(5)-  120*J(6)

But wait a moment. Does table 1 and table 2 really amount to the same thing? Well, almost. In the first table B1 = −1/2 and in the second table B1 = 1/2. So the first case refers to a definition Bn = Bn(0) and the second case to a definition Bn = Bn(1). Apart from this the values in the tables are identical (as the two definitions differ only in this particular case).

Let us denote the coefficients in table 3 by V(n, k).

             V(n, k)
                1
              0, -1
             0, -1, 2
           0, -1, 6, -6
        0, -1, 14, -36, 24
    0, -1, 30, -150, 240, -120
0, -1, 62, -540, 1560, -1800, 720

What we have to find is a formula for the V(n, k). This is not difficult if we look at the relationship with the signed Worpitzky numbers W(n, k) = (-1)^k k! {n+1|k+1}. Here {n|k} denotes the Stirling set numbers (aka. Sn,k of the second kind). Our formal definition is V(n, k) = Sum(j=k..n) W(n,j). Now we can restate our findings:

For n > 1
Bn = Sum(k=0..n) W(n, k) J(k+1)
    = Sum(k=0..n) W(n, k) H(k+1)
    = Sum(k=0..n) V(n, k)  J(k+1).

The first identity is Wropitzky's, which we take as granted. J(k+1) = H(k+1) − H(k).  So we can rewrite Wropitzky's sum as Sum(k=0..n) W(n, k)(H(k+1) − H(k)). Observe that Sum(k=0..n)  W(n, k) H(k) = 0,-1,0,0,0,0,0,.. (starting at n = 0). Thus if n <> 1 we can simplify to Sum(k=0..n) W(n, k) H(k+1), which is the second identity. The third identity follows from W(n, k) = V(n, k) − V(n, k+1).

For those who would like to check the formulas with Maple:

W := proc(n,k) (-1)^k*combinat[stirling2](n+1,k+1)*k! end:
V := proc(n,k) add(W(n,j),j=k..n) end:
H := proc(n) add(1/k,k=1..n) end:

seq(add(W(n,k)*(H(k+1)-H(k)),k=0..n),n=0..23);
seq(add(W(n,k)/(k+1),k=0..n),n=0..23);
seq(bernoulli(n, 1),n=0..23);

seq(add(W(n,k)*H(k+1),k=0..n),n=0..23);
seq(add(V(n,k)/(k+1),k=0..n),n=0..23);
seq(bernoulli(n, 0),n=0..23);

seq(add(W(n,k)*H(k),k=0..n),n=0..23);

Inverse polynomial harmony for Jakob

Putting the transformation into a plotting bag ...

MyPlot := proc(f,R) local T,i;
T := proc(f,n,x) local k,v; add(add((-1)^v*binomial(k,v)
*f(k+1)*(x+v+1)^n,v=0..k),k=0..n) end:
plot([seq(T(f,i,x),i=2..7)],x=R,thickness=2,axes=boxed) end:

... we can visualize the effect of the harmonic numbers versus the effect of the inverse numbers:

A := MyPlot(x->add(1/i, i=1..x),0..1):
B := MyPlot(x->1/x, -1..0):
plots[display]([A,B]);
Bernoulli polynomials.

Be careful when interpreting this plot: these are not the Bernoulli polynomials. This are two plots in one. The left hand side (the interval [-1,0]) is generated by feeding J(n) into the transformation, the right hand side (the interval [0,1]) by feeding H(n) into the transformation. Only on [0,1] this coincides with the Bernoulli polynomials. 

A primer on Worpitzky numbers

  • Worpitzky (n ≥ 0 and k ≥ 0)
     
    • Definition
      W := proc(n, k) stirling2(n+1,k+1)*k! end:
    • Recursion
      Wrec := proc(n, k) option remember; 
      if k > n then 0 elif n = 0 then 1 
      else k*Wrec(n-1, k-1) + (k+1)*Wrec(n-1, k) fi end:
    • Egf
      w := (x, y) -> exp(x)/(1+y*(1-exp(x)));
  • V-Numbers (n ≥ 0 and k ≥ 0)
     
    • Definition
      V := proc(n, k) add(W(n, j), j=k..n) end:
    • Recursion
      Vrec := proc(n, k) option remember;
      if k > n then 0 elif n = 0 then 1
      else k*(Vrec(n-1, k-1) + Vrec(n-1, k)) fi end:
      
    • Egf
      v := (x, y) -> 1/(1+y*(1-exp(x)));

A challenger for the Akiyama–Tanigawa algorithm?

The Akiyama–Tanigawa algorithm is a cute way to compute the Bernoulli numbers.

AT := proc(n) local m, j, A;
for m from 0 by 1 to n do
   A[m] := 1/(m + 1);
   for j from m by -1 to 1 do
      A[j - 1] := j * (A[j - 1] - A[j])
od od; A[0] end:

The sum representation Bn = Sum(k=0..n) V(n, k)/(k+1) and the recursion for V(n,k) suggest a similar algorithm.

VB := proc(n) local m, j, A;
if n = 0 then 1 else A[0] := 0; A[1] := 1;
for m from 2 by 1 to n do A[m] := 0;
   for j from m by -1 to 1 do
      A[j] := j * (A[j - 1] + A[j])
od od; 
add((-1)^j*A[j]/(j+1),j=1..n) fi end:

An advantage of the VB algorithm over the AT algorithm is that it postpones rational arithmetic to the computation of the sum (the A[j]'s in AT are rational numbers, the A[j]'s in VB are integers). Does this make a noticeable difference? Let us benchmark! I used the two procedures as given above; however, I replaced 'end:' by 'NULL; end:' to suppress the output. On my personal computer with Maple V Release 5 I got the following results (times in seconds):

Bernoulli Benchmark
n AT VB
100 0.047 0.016
200 0.594 0.187
400 7.656 1.203
800 92.766 9.376

VB is up to ten times faster! This is fun. ;-)

New sequences for the OEIS

Based on the above considerations I will submit the following sequences to the OEIS.

                                  0,
                                 0, 1,
                               0, 3, -3,
                            0, 7, -18, 11,
                         0, 15, -75, 110, -50,
                     0, 31, -270, 715, -750, 274,
                0, 63, -903, 3850, -7000, 5754, -1764,
          0, 127, -2898, 18711, -52500, 72884, -49392, 13068,

This triangle is (up to sign) Worpitzky(n,k)Harmonic(k) as well as Stirling1(k+1,2)Stirling2(n+1,k+1) for n ≥ 0 and k ≥ 0. The fact that the sum of the rows is 0,1,0,0,0,... was a crucial step in some part of the proof above. The triangle should be supplemented by the sequence 0,1,3,18,125,1020,9667,104790,.. which is the sum over the odd k's in a row. (Who comes up with a combinatorial interpretation?!) Note that the right hand side of the triangle is A000254 resp. A081048.
The sequences now are : A176276 and A176277.