Re: computing Bernoulli numbers
- From: Reinhold Bader <Bader@xxxxxx>
- Date: Thu, 31 Jan 2008 16:58:32 +0100
-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1
Bart Vandewoestyne schrieb:
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)
integer overflow in i**n before converting to floating point? Will impact
slowest converging series (small n) first ...
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
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.2 (GNU/Linux)
Comment: Using GnuPG with SUSE - http://enigmail.mozdev.org
iD8DBQFHofAoFVLhKuD7VgsRAtpLAKDGz6D828KpbtBG/xMxQ7PSmRqLnwCcC5R0
yMy5NdpIQINc+X7leEOT13o=
=ic1K
-----END PGP SIGNATURE-----
.
- References:
- computing Bernoulli numbers
- From: Bart Vandewoestyne
- computing Bernoulli numbers
- Prev by Date: Re: computing Bernoulli numbers
- Next by Date: Re: computing Bernoulli numbers
- Previous by thread: Re: computing Bernoulli numbers
- Next by thread: Re: computing Bernoulli numbers
- Index(es):
Relevant Pages
|
|
