Re: computing Bernoulli numbers
- From: "robin" <robin_v@xxxxxxxxxxx>
- Date: Fri, 08 Feb 2008 23:49:33 GMT
"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.
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
I think that you have compiled these as separate units.
You should use CONTAINS, so that the kind of real_factorial
in the subroutine is visible to the subroutine.
Alternatively, explicitly type the function.
As well, you need to use IMPLICIT NONE in your procedures.
Had you done so, you would have received an error message
for subroutine bernoulli_series
telling you that real_factorial had not been given a type.
... 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:
The fact that they are accurate to single precision should give a clue ...
.
- Prev by Date: Re: What is Standards Conformance? (Was: Is it possible to invoke a Unix command within Fortran 90 source)
- Next by Date: Re: OT Launched - pictures from Andoya, Norway - Andenes, Bleik, and the rocket range
- Previous by thread: Re: computing Bernoulli numbers
- Next by thread: Re: passing a variable number of args/vectors to a routine
- Index(es):
Relevant Pages
|
|
