Re: computing Bernoulli numbers

On 2008-01-31, Bart Vandewoestyne <MyFirstName.MyLastName@xxxxxxxxxx> wrote:

else

mysum = 0
i = 1
do
term = 1.0_qp/(i**n)
if (term < spacing(mysum)) then
exit
end if
mysum = mysum + term
i = i+1
end do

bn = 2*mysum*real_factorial(n)/((2*pi)**n)

Notice that I have to add

if (modulo(n,4)==0) then
bn = -bn
end if

to get the sign right. But that doesn't change my original question.

Regards,
Bart

--
"Share what you know. Learn what you don't."
.