Re: computing Bernoulli numbers

Fortran code for evaluating a sequence of Bernoulli numbers and other special functions:
http://pagesperso-orange.fr/jean-pierre.moreau/f_function.html
and
http://jin.ece.uiuc.edu/routines/routines.html

Skip Knoble

On Thu, 31 Jan 2008 14:14:51 GMT, Bart Vandewoestyne <MyFirstName.MyLastName@xxxxxxxxxx>
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:
-|
-| 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

.

Relevant Pages

  • Re: Fortran Compling problem
    ... >> I got a Fortran code called APEX1310.f and two codes called ... or write your own GETDAT and GETTIM wrapper functions. ... subroutine GetDat ...
    (comp.lang.fortran)
  • Re: matlab mex file (fortran)
    ... > I send you my hypothetical fortran code too. ... Program Celsius Table: Prints simple Fahrenheit-Celsius table ... >>Subroutine mexFunction(nlhs, plhs, nrhs, prhs) ...
    (comp.soft-sys.matlab)
  • Re: Recompilation of mex-files fails
    ... I used 'function' as the fortran code to be called is a ... changing 'function' to 'subroutine' avoids this warning. ... C The gateway routine ... subroutine mexFunction(nlhs, plhs, nrhs, prhs) ...
    (comp.soft-sys.matlab)
  • Re: computing Bernoulli numbers
    ... subroutine bernoulli_series ... real:: mysum, term ... cannot compute Bernoulli number for negative n!" ... expecting: -3.3333333333333333E-02 ...
    (comp.lang.fortran)
  • Re: computing Bernoulli numbers
    ... subroutine bernoulli_series ... real:: mysum, term ... cannot compute Bernoulli number for negative n!" ... i already notice quite some accuracy loss ...
    (comp.lang.fortran)
Loading