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!
> I think the importance of computation speed depends on what the user needs.

I agree. I was posting more to explore other options.

--
J. J. Weimer
Chemistry / Chemical & Materials Engineering, UAH

Forum

Support

Gallery

Igor Pro 8

Learn More

Igor XOP Toolkit

Learn More

Igor NIDAQ Tools MX

Learn More