A while ago I found this algorithm. Today I read in wikipedia that Euler zig zag numbers can be used for computing the Bernoulli numbers. This Mathematica program computes the Euler zig zag numbers and there after the Bernoulli numbers.
Clear[nn, A1, A2, A3, A4, A5, A6, A7, A8, A9]
nn = 42;
Clear[t, n, k];
t[n_, 1] = 1;
t[n_, k_] :=
t[n, k] =
If[n >= k, Sum[t[n - i, k - 1] + t[n - i, k], {i, 1, 1}], 0];
A1 = Table[Table[t[n, k], {k, 1, nn}], {n, 1, nn}];
MatrixForm[A1];
Clear[t, n, k];
t[n_, 1] = If[Or[Mod[n, 4] == 1, Mod[n, 4] == 0], 1, -1];
t[n_, k_] := t[n, k] = If[n >= k, t[n - 1, k - 1], 0];
A2 = Table[Table[t[n, k], {k, 1, nn}], {n, 1, nn}];
MatrixForm[A2];
Clear[t, n, k];
t[n_, k_] :=
t[n, k] = If[n >= k, If[n == k, 1, If[Mod[k, 2] == 1, 0, 1]], 0];
A3 = Table[Table[t[n, k], {k, 1, nn}], {n, 1, nn}];
MatrixForm[A3];
Clear[t, n, k];
t[n_, k_] :=
t[n, k] = If[n >= k, If[n == k, 1, If[Mod[k, 2] == 1, 1, 0]], 0];
A4 = Table[Table[t[n, k], {k, 1, nn}], {n, 1, nn}];
MatrixForm[A4];
A5 = A1*A2*A3;
MatrixForm[A5];
A6 = A1*A2*A4;
MatrixForm[A6];
A7 = Inverse[A5];
MatrixForm[A7];
A8 = Inverse[A6];
MatrixForm[A8];
A9 = A7 + A8 - IdentityMatrix[nn];
MatrixForm[A9];
Total[Transpose[A9]];
aa = A9[[All, 1]];
Table[(-1)^Floor[n/2]*n/(2^n - 4^n)*aa[[n]], {n, 2, Length[aa], 2}];
Denominator[%];
Table[(-1)^(2*n)*(n*2)/(2^(n*2) - 4^(n*2))*aa[[2*n]], {n, 1,
Length[aa]/2}]
Denominator[%]
Output: {-(1/6), -(1/30), -(1/42), -(1/30), -(5/66), -(691/2730), -(7/6), -( 3617/510), -(43867/798), -(174611/330), -(854513/138), -(236364091/ 2730), -(8553103/6), -(23749461029/870), -(8615841276005/14322), -( 7709321041217/510), -(2577687858367/6), -(26315271553053477373/ 1919190), -(2929993913841559/6), -(261082718496449122051/13530), -( 1520097643918070802691/1806)}
Denominators: {6, 30, 42, 30, 66, 2730, 6, 510, 798, 330, 138, 2730, 6, 870, 14322, \ 510, 6, 1919190, 6, 13530, 1806}
Since it is basically matrix inversion of the Pascal triangle, can we find a closed form for the Bernoulli numbers?
The closed-form formula \begin{equation}\label{Higgins-Gould-B}\tag{1} B_n=\sum_{k=0}^n\frac1{k+1}\sum_{j=0}^k(-1)^j\binom{k}{j}j^n,\quad n\ge0, \end{equation} has a long history, it appeared in the paper
It is a special case of the formula (2.5) in the paper
Its equivalent form is \begin{equation}\label{Bernoulli-Stirling-eq}\tag{2} B_n=\sum_{k=0}^n(-1)^k\frac{k!}{k+1}S(n,k), \quad n\ge0, \end{equation} where $S(n,k)$ is the Stirling numbers of the second kind. See page 29, Remark 2 in the paper
There existed at least seven alternative proofs of the formulas \eqref{Higgins-Gould-B} and \eqref{Bernoulli-Stirling-eq} in the paper [1, 2] above and in the following monographs and papers:
To the best of my knowledge, except the formulas \eqref{Higgins-Gould-B} and \eqref{Bernoulli-Stirling-eq} above, there are also the following explicit formulas for the Bernoulli numbers $B_n$: \begin{align}\label{Higgins-Gould-B(11)}\tag{3} B_n&=\sum_{j=0}^n(-1)^j\binom{n+1}{j+1}\frac{n!}{(n+j)!}\sum_{k=0}^j(-1)^{j-k}\binom{j}{k}k^{n+j}, \quad n\ge0;\\ B_n&=\sum_{i=0}^n(-1)^{i}\frac{\binom{n+1}{i+1}}{\binom{n+i}{i}}S(n+i,i), \quad n\ge0;\label{Bernoulli-Stirling-formula}\tag{4}\\ B_{2k}&=1+\sum_{m=1}^{2k-1}\frac{S(2k+1,m+1) S(2k,2k-m)}{\binom{2k}{m}}\\ &\quad-\frac{2k}{2k+1}\sum_{m=1}^{2k}\frac{S(2k,m)S(2k+1,2k-m+1)}{\binom{2k}{m-1}}, \quad k\in\mathbb{N};\tag{5}\\ B_{2k}&=\frac{(-1)^{k-1}k}{2^{2(k-1)}(2^{2k}-1)}\sum_{i=0}^{k-1}\sum_{\ell=0}^{k-i-1} (-1)^{i+\ell}\binom{2k}{\ell}(k-i-\ell)^{2k-1}, \quad k\in\mathbb{N};\tag{6}\\ B_{2m}&=(-1)^{m-1}\frac{m} {2^{2m-1}\bigl(2^{2m}-1\bigr)}\Biggl[\sum_{k=0}^{m-1} (-1)^k\binom{2m}{k}(m-k)^{2m-1}\\ &\quad+2\sum_{k=1}^{m-1}(-1)^k\sum_{\ell=0}^{m-k-1} (-1)^{\ell}\binom{2m}{\ell}(m-k-\ell)^{2m-1}\Biggr],\quad m\in\mathbb{N};\\ B_{2m}&=\frac{m} {2^{2m-1}\bigl(2^{2m}-1\bigr)}\sum_{\ell=1}^{2m}\frac{(-1)^{\ell-1}}{2^\ell} \biggl(\frac1\ell-\frac1{m+1}\biggr) \binom{2m+1}{\ell} \sum_{q=0}^\ell\binom{\ell}{q}(2q-\ell)^{2m}, \quad m\in\mathbb{N};\\ B_{2k}&= \frac12 - \frac1{2k+1} - 2k \sum_{i=1}^{k-1} \frac{A_{2(k-i)}}{2(k - i) + 1},\quad k\in\mathbb{N};\tag{7} \end{align} where $A_m$ is defined by \begin{equation*} \sum_{m=1}^nm^k=\sum_{m=0}^{k+1}A_mn^{m}. \end{equation*} The formulas \eqref{Higgins-Gould-B(11)} and \eqref{Bernoulli-Stirling-formula} are also equivalent to eah other.
By the way, I would like to mention two intereting double inequalities related to the Bernoulli numbers $B_{2n}$ as follows.
More related references