KEYWORDS: Bernoulli, Euler, Tangent, Harmonic, Binomial,
SwissKnife, Worpitzky, Akiyama–Tanigawa.
Concerned with sequences: A027641, A027642, A131689, A019538, A028246, A141056,
A027760, A176276, A176277.
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:
f_{1} 
xf_{1}+f_{1}−f_{2} 
x^{2}f_{1}+(2f_{1}−2f_{2}) x+2f_{3}+f_{1}−3f_{2} 
x^{3}f_{1}+(3f_{1}−3f_{2})x^{2}+(6f_{3}+3f_{1} −9f_{2})x−6f_{4}+12f_{3}+f_{1}−7f_{2} 
x^{4}f_{1}+(4f_{1}−4f_{2})x^{3}+(12f_{3} +6f_{1}−18f_{2})x^{2} +(−24f_{4}+48f_{3}+4f_{1}−28f_{2})x−60f_{4} +50f_{3}+f_{1}15f_{2}+24f_{5} 
The sequence of coefficients of these polynomials, sorted in descending order,
is

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(m1,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 mathparlance.
Let B_{n}(x) be the Bernoulli polynomials and H_{n} 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...) .
I am amazed that this formula seems to be unknown. Much more exotic relations between the harmonic numbers and the Bernoulli numbers can be easily found (look for example at the section Assorted identities in the Wikipedia article). Indeed the formula looks like an exercise to the chapter Harmonic Summation in GKP's Concrete Mathematics (which is followed by the chapter Bernoulli numbers).
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 B_{n} = B_{n}(1) is to be preferred over the definition B_{n} = B_{n}(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 B_{1} = 1/2 as you can see from the facsimile on Wikipedia, and Jakob was a very concrete mathematician.
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 B_{1} = −1/2 and in the second table B_{1} = 1/2. So the first case refers to a definition B_{n} = B_{n}(0) and the second case to a definition B_{n} = B_{n}(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+1k+1}. Here {nk} denotes the Stirling set numbers (aka. S_{n,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 B_{n} = 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);
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 now 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.

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 B_{n} = 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):

VB is up to ten times faster! This is fun. ;) However, Maple's Bernoulli function is even faster. So I remembered the Knuth−Buckholz method, which is based on the computation of the tangent numbers (see CM, 6.5) and only needs a single division.
SeidelTangentBernoulli := proc(n) local i, j: if n = 0 then 1 elif n = 1 then 1/2 elif n mod 2 = 1 then 0 else i := 2^n; i := i*(i1); (1)^iquo(n1,2)*n*SeidelEulerTangent(n+1)/i fi end:
To compute the tangent numbers I used Seidel's boustrophedon transform (which is described on the Wikipedia page of the Bernoulli numbers, and of course in the paper by Millar, Sloane and Young, see example 3). The timings of the resulting algorithm are nearly the same as those of Maple's function (MV R5).
SeidelEulerTangent := proc(N) local A,Am,n,i,h,k,e,j; n := 2*N  5; if n < 0 then RETURN(1) fi; h := iquo(n+1,4); A := array[h..h]; A[1] := 0; A[0] := 1; k := 0; e := 1; for i from 0 to iquo(n+1,2) do Am := 0; A[k + e] := 0; e := e; for j from 0 to i do Am := Am + A[k]; A[k] := Am; k := k + e; od od; # print([seq(A[z], z=h..h)]); if N mod 2 = 0 then A[h] else A[h] fi end:
Consider the signless triangle Worpitzky(n, k)Harmonic(k) (which is also Abs(Stirling1(k+1,2)Stirling2(n+1,k+1))) for n ≥ 0 and k ≥ 0. The fact that the alternating sum of the rows is 0,1,0,0,0,... was a crucial step in some part of the proof above. This 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 (unsigned Stirling numbers of first kind).
Based on the above considerations I submitted the following sequences to OEIS.










Copyright: Peter Luschny, 2010. This article was written online on my user page of the OEISwiki. License: Creative Commons AttributionShareAlike 3.0 Unported License (the same license which Wikipedia uses).