bernob Subroutine

subroutine bernob ( n, bn )

  !*****************************************************************************80
  !
  !! BERNOB computes the Bernoulli number Bn.
  !
  !  Licensing:
  !
  !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However, 
  !    they give permission to incorporate this routine into a user program 
  !    provided that the copyright is acknowledged.
  !
  !  Modified:
  !
  !    11 July 2012
  !
  !  Author:
  !
  !    Shanjie Zhang, Jianming Jin
  !
  !  Reference:
  !
  !    Shanjie Zhang, Jianming Jin,
  !    Computation of Special Functions,
  !    Wiley, 1996,
  !    ISBN: 0-471-11963-6,
  !    LC: QA351.C45.
  !
  !  Parameters:
  !
  !    Input, integer ( kind = 4 ) N, the index.
  !
  !    Output, real ( kind = 8 ) BN, the value of the N-th Bernoulli number.
  !
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) bn(0:n)
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  real ( kind = 8 ) r1
  real ( kind = 8 ) r2
  real ( kind = 8 ) s
  real ( kind = 8 ) tpi

  tpi = 6.283185307179586D+00
  bn(0) = 1.0D+00
  bn(1) = -0.5D+00
  bn(2) = 1.0D+00 / 6.0D+00
  r1 = ( 2.0D+00 / tpi )**2

  do m = 4, n, 2

     r1 = - r1 * ( m - 1 ) * m / ( tpi * tpi )
     r2 = 1.0D+00

     do k = 2, 10000
        s = ( 1.0D+00 / k ) ** m
        r2 = r2 + s
        if ( s < 1.0D-15 ) then
           exit
        end if
     end do

     bn(m) = r1 * r2

  end do

  return
end subroutine bernob