computing Bernoulli numbers [Fortran]

Prev: Fortran 2003: abstract interfaces
Next: Type specification and initialization expressions

From: Bart Vandewoestyne on 1 Feb 2008 09:14

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."

From: Bart Vandewoestyne on 1 Feb 2008 09:17

On 2008-01-31, Bart Vandewoestyne <MyFirstName.MyLastName(a)telenet.be> 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."

From: Michel Olagnon on 1 Feb 2008 09:24

Bart Vandewoestyne wrote:
> 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:

Can't see the interface declaring real_factorial() as (kind=qp) in
the calling routine.

>
> 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.
>

From: Herman D. Knoble on 1 Feb 2008 09:41

From: Bart Vandewoestyne on 1 Feb 2008 10:06

On 2008-01-31, Michel Olagnon <molagnon(a)ifremer-a-oter.fr> wrote:
>
> Can't see the interface declaring real_factorial() as (kind=qp) in
> the calling routine.

Michel,

I'm afraid I don't understand what you mean with your response.
I did show how real_factorial(n) looks like:

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 I use it like

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

which looks OK to me.

What would you like to know/see more?

Regards,
Bart

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

| Next | Last
Pages: 1 2 3
Prev: Fortran 2003: abstract interfaces
Next: Type specification and initialization expressions