Re: computing Bernoulli numbers
- From: Michel Olagnon <molagnon@xxxxxxxxxxxxxxxxx>
- Date: Thu, 31 Jan 2008 15:24:43 +0100
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.
.
- Follow-Ups:
- Re: computing Bernoulli numbers
- From: Bart Vandewoestyne
- Re: computing Bernoulli numbers
- References:
- computing Bernoulli numbers
- From: Bart Vandewoestyne
- computing Bernoulli numbers
- Prev by Date: Re: computing Bernoulli numbers
- Next by Date: Re: How to quote the " character in a format statement?
- Previous by thread: Re: computing Bernoulli numbers
- Next by thread: Re: computing Bernoulli numbers
- Index(es):
Relevant Pages
|
|
