eulera Subroutine

subroutine eulera ( n, en )

  !*****************************************************************************80
  !
  !! EULERA computes the Euler number En.
  !
  !  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:
  !
  !    10 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 of the highest value to compute.
  !
  !    Output, real ( kind = 8 ) EN(0:N), the Euler numbers up to the N-th value.
  !
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) en(0:n)
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  real ( kind = 8 ) r
  real ( kind = 8 ) s

  en(0) = 1.0D+00

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

  return
end subroutine eulera