Erst mal Danke für die schnelle, ausführliche und verständliche Antwort!
An Hand von deinen Ausführungen hat sich für mich aber jetzt noch eine Frage ergeben. So wie du die Formel oben ausgeschrieben hast, ist mir beim nochmaligen Stöbern in Wikipedia gleich die Formel hier ins Auge gesprungen:
$$B_m^{+} = \sum \limits_{k=0}^{m}\sum \limits_{v=0}^{k}(-1)^{v}\begin{pmatrix} k\\v \end{pmatrix}\frac{(v+1)^{m}}{k+1}$$
https://en.wikipedia.org/wiki/Bernoulli_number#Explicit_definition
Statt es rekursiv zu programmieren, habe ich es also einfach mal genau so umgesetzt, wie die explizite Formel es vorgesehen hat. Im Kern entspricht das auch genau deiner Erläuterung mit der Doppelschleife. Ich habe das auch darum gemacht, weil das System, in dem ich die Funktion verwenden will, leider eine maximale Rekursionstiefe hat, die bei großen Bernoulli-Zahlen leider überschritten wird.
Nun rechnet mein Programm mit der expliziten Formel aber leider falsche Zahlen aus.
$$B_0^{+} = \frac{-1}{1}$$
$$B_1^{+} = \frac{-5}{2}$$
$$B_2^{+} = \frac{-19}{2}$$
Ich habe das mal händisch nachgeprüft und komme auf die gleichen Ergebnisse und wundere mich darum, wo der Fehler liegt...
$$B_0^{+} = \sum \limits_{k=0}^{0}\sum \limits_{v=0}^{k}(-1)^{v}\begin{pmatrix} k\\v \end{pmatrix}\frac{(v+1)^{m}}{k+1}\\ = (-1)^{0}\begin{pmatrix} 0\\0 \end{pmatrix}\frac{(0+1)^{0}}{0+1}\\ = (-1)^{0} \cdot \frac{0!}{0!\cdot(0 - 0)!}\cdot\frac{(0+1)^{0}}{0+1}\\ = -1 \cdot \frac{1}{1} \cdot \frac{1}{1}\\ = -1 \cdot 1 \cdot 1\\ B_0^{+} = -1$$
Und für
$$B_1^{+} = \sum \limits_{k=0}^{1}\sum \limits_{v=0}^{k}(-1)^{v}\begin{pmatrix} k\\v \end{pmatrix}\frac{(v+1)^{m}}{k+1}\\ \text{ Ist letztlich die Summe aus drei Elementen:}\\ 0.0 \rightarrow(-1)^{0} \cdot \begin{pmatrix} 0\\0 \end{pmatrix} \cdot \frac{(0+1)^{0}}{0+1} = -1 \cdot \frac{1}{1} \cdot \frac{1}{1} = -1\\ 1.0 \rightarrow(-1)^{0} \cdot \begin{pmatrix} 1\\0 \end{pmatrix} \cdot \frac{(0+1)^{1}}{1+1} = -1 \cdot \frac{1}{1} \cdot \frac{1}{2} = -\frac{1}{2}\\ 1.1 \rightarrow(-1)^{1} \cdot \begin{pmatrix} 1\\1 \end{pmatrix} \cdot \frac{(1+1)^{1}}{1+1} = -1 \cdot \frac{1}{1} \cdot \frac{2}{2} = -1\\ \text{ Aus den zwei "inneren" Schleifen ergeben sich:}\\ \sum -1 = -1\\ \sum -\frac{1}{2}, -1 = -\frac{3}{2}\\ \text{ Beide Summen werden in der "äußeren" Schleife wiederum summiert:}\\ \sum -1, -\frac{3}{2} = -\frac{5}{2}\\ B_1^{+} = -\frac{5}{2}$$
Genau das gibt mein Programm auch für die ersten beiden Zahlen aus und wenn ich mich nicht zweimal gleichmäßig händisch verrechnet habe, scheint das ja auch zu stimmen.
Aber wie gesagt sind das natürlich nicht die ersten beiden Bernoulli Zahlen und ich frage mich was da falsch läuft.
Im Gegensatz zu den rekursiven Definitionen, bei denen $$B_0^{+}$$ und $$B_1^{+}$$ ja vordefiniert sind, müsste sich über die explizite Formel ja jede Bernoulli Zahl entsprechend errechnen lassen.