This site is supported by donations to The OEIS Foundation.
User:Peter Luschny/SeqTransformation
Contents
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+f1−f2 |
| 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 | f1−f2 | |||
| 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.
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]);
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)));
- Definition
- 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)));
- Definition
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):
| 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.
