# Bernoulli Number generator

Here is a function for computing Bernoulli numbers (https://en.wikipedia.org/wiki/Bernoulli_number). The function is recursive, so it allows for succinct coding. However it is not the fastest method, particularly if you need many sequential Bernoulli numbers. Developed with IP7.05 x64 (Windows), and tested up to N=20. The output is numeric, and not in rational-fraction format.
function fB(nB)     // recursive function for Bernoulli Number
variable nB  // number  index argument

variable  v
if(nB==0)
return 1
else
make/O/D/FREE/N=(nB) wB
wB = binomial(nB+1,p)*fB(p)
v = -sum(wB)/(nB+1)
return v  // Bernoulli number value
endif
end
After further research, I have found a new function in IP7 that simplifies the generation of Bernoulli numbers:
DisplayHelpTopic "zeta"

function fBernoulli(n)
variable n

if(n==0)
return 1
elseif(n==1)
return -0.5
elseif(mod(n,2)!=0) //  n>1 is odd here
return 0
else  // n>1 is even here
return   (-1)^(n/2-1)*factorial(n)*2*zeta(n,1,400)/( ((2*pi)^(n)) )
endif
end

The second argument, 1, in `zeta(n,1,400)` makes it the Riemann zeta function; the third argument improves the precision for small 'n'.
A more compact implementation (perhaps faster)?

function fBernoulli(n)
variable n

variable rtn

make/N=2/FREE rtw = {1,-0.5}

if (n<=1)
return rtw[n]
endif

rtn = mod(n,2) != 0 ? 0 : (-1)^(n/2-1)*factorial(n)*2*zeta(n,1,400)/( ((2*pi)^(n)) )

return rtn
end

--
J. J. Weimer
Chemistry / Chemical & Materials Engineering, UAH
Jeffrey,
I think the importance of computation speed depends on what the user needs. For computing many Bn at the same time a recursive approach is widely used. As Brent notes (https://arxiv.org/abs/1108.0286v3): "However, this is unsatisfactory if floating-point numbers are used, because the recurrence in numerically unstable." Harvey (http://web.maths.unsw.edu.au/~davidharvey/talks/bernoulli.pdf) gives an interesting approach wherein ONE division operation provides all necessary coefficients for the Bn up to a given order. The hitch is that huge numerators and denominators are required, and the division is best done in base 2. I tested the concept using `APmath` in base 10, and it works! Forum Support Gallery