Re: computing Bernoulli numbers
- From: "Gerry Ford" <invalid@xxxxxxxxxxx>
- Date: Thu, 31 Jan 2008 18:34:11 -0700
"Bart Vandewoestyne" <MyFirstName.MyLastName@xxxxxxxxxx> wrote in message
news:vJkoj.683$6n5.478@xxxxxxxxxxxxxxxx
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?
You speak of Bernoulli numbers as if they were as familiar as binomial
coefficients. I identify the Bernoulli's with fluid dynamics and doubt they
would have known anything of Riemann's zeta function. A google search turns
up a lot of information, but my baud rate has sunk to less than a tenth of a
dial-up, so tja.
What are Bernoulli numbers?
--
Gerry Ford
http://mathworld.wolfram.com/BernoulliNumber.html
.
- Follow-Ups:
- Re: computing Bernoulli numbers
- From: Gerry Ford
- Re: computing Bernoulli numbers
- Prev by Date: Re: installing fruit
- Next by Date: Re: Automatically transform or expand do loop in a subroutine
- Previous by thread: Re: installing fruit
- Next by thread: Re: computing Bernoulli numbers
- Index(es):
Relevant Pages
|
|
