bernoa Subroutine

subroutine bernoa ( n, bn )

  !*****************************************************************************80
  !
  !! BERNOA 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 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  real ( kind = 8 ) r
  real ( kind = 8 ) s

  bn(0) = 1.0D+00
  bn(1) = -0.5D+00

  do m = 2, n
     s = - ( 1.0D+00 / ( m + 1.0D+00 ) - 0.5D+00 )
     do k = 2, m - 1
        r = 1.0D+00
        do j = 2, k
           r = r * ( j + m - k ) / j
        end do
        s = s - r * bn(k)
     end do
     bn(m) = s
  end do

  do m = 3, n, 2
     bn(m) = 0.0D+00
  end do

  return
end subroutine bernoa