eulerb Subroutine

subroutine eulerb ( n, en )

  !*****************************************************************************80
  !
  !! EULERB 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:
  !
  !    09 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)
  real ( kind = 8 ) hpi
  real ( kind = 8 ) isgn
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  real ( kind = 8 ) r1
  real ( kind = 8 ) r2
  real ( kind = 8 ) s

  hpi = 2.0D+00 / 3.141592653589793D+00
  en(0) = 1.0D+00
  en(2) = -1.0D+00
  r1 = -4.0D+00 * hpi ** 3

  do m = 4, n, 2
     r1 = - r1 * ( m - 1 ) * m * hpi * hpi
     r2 = 1.0D+00
     isgn = 1.0D+00
     do k = 3, 1000, 2
        isgn = - isgn
        s = ( 1.0D+00 / k ) ** ( m + 1 )
        r2 = r2 + isgn * s
        if ( s < 1.0D-15 ) then
           exit
        end if
     end do

     en(m) = r1 * r2

  end do

  return
end subroutine eulerb