computing Bernoulli numbers

I would like to implement a routine to compute floating point
Bernoulli numbers.

After reading through the thread

http://groups.google.be/group/sci.math.num-analysis/browse_thread/thread/88cab7bfaf115f6d/671edcdf8ebc5523?hl=nl&lnk=st&q=#671edcdf8ebc5523

it seemed to me that equation (23.2.16) that uses the
Riemann-Zeta function is the way to go.

I tried implementing it as follows:

subroutine bernoulli_series(n, bn)
integer(kind=i4b), intent(in) :: n
real(kind=qp), intent(out) :: bn

integer(kind=i4b) :: i
real(kind=qp) :: mysum, term


if (n<0) then
write(unit=*, fmt="(A)") &
"ERROR: cannot compute Bernoulli number for negative n!"
stop
end if

if (n==0) then

bn = 1.0_qp

else if (n==1) then

bn = -0.5_qp

else if (modulo(n,2)==1) then

bn = 0.0_qp

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)

end if

end subroutine bernoulli_series


.... where ...


function real_factorial(n) result(res)
integer(kind=i4b), intent(in) :: n
real(kind=qp) :: res

integer(kind=i4b) :: i

res = 1.0_qp
do i = 1,n
res = res*i
end do

end function real_factorial


.... and qp is the usual double precision.


When I test this subroutine and compare the values against what one
should expect theoretically, i already notice quite some accuracy loss
from the beginning:

=> Testing bernoulli(n)...
bernoulli_series(0)
bernoulli_series(1)
bernoulli_series(2)

ERROR:
found : 0.1666644802167240
expecting: 0.1666666666666667

bernoulli_series(3)
bernoulli_series(4)

ERROR:
found : 3.3333332307553823E-02
expecting: -3.3333333333333333E-02

bernoulli_series(5)
bernoulli_series(6)

ERROR:
found : 2.3809523726588188E-02
expecting: 2.3809523809523808E-02

bernoulli_series(7)
bernoulli_series(8)

ERROR:
found : 3.3333333298527612E-02
expecting: -3.3333333333333333E-02

bernoulli_series(9)
bernoulli_series(10)

ERROR:
found : 7.5757575723059675E-02
expecting: 7.5757575757575760E-02

bernoulli_series(11)


As the calculated value is always smaller than the expected value,
I assume I am not adding enough terms. The condition that checks
whether an extra term should be added or not is

if (term < spacing(mysum)) then

Is this causing the loss in accuracy? Or does anybody else see another
reason?

How do I change my code so that my Bernoulli numbers are computed more accurate?

Thanks,
Bart

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

Relevant Pages

  • Re: computing Bernoulli numbers
    ... subroutine bernoulli_series ... real:: mysum, term ... cannot compute Bernoulli number for negative n!" ... expecting: -3.3333333333333333E-02 ...
    (comp.lang.fortran)
  • Re: computing Bernoulli numbers
    ... subroutine bernoulli_series ... real:: mysum, term ... cannot compute Bernoulli number for negative n!" ... expecting: -3.3333333333333333E-02 ...
    (comp.lang.fortran)
  • Re: computing Bernoulli numbers
    ... subroutine bernoulli_series ... real:: mysum, term ... cannot compute Bernoulli number for negative n!" ... i already notice quite some accuracy loss ...
    (comp.lang.fortran)
  • Re: Get date between two dates
    ... Instead of posting all the code in your subroutine and expecting the ... people who may be able to help you to sort it out, ...
    (microsoft.public.office.developer.outlook.vba)