!*********************************************************************
!* This program computes Bernoulli number Bn using subroutine BERNOA *
!* ----------------------------------------------------------------- *
!* SAMPLE RUN:                                                       *
!*                                                                   *
!*  Compute Bernoulli number Bn for n = 0,1,...,10.                  *
!*                                                                   *
!*  Please enter Nmax: 10                                            *
!*                                                                   *
!*    n            Bn                                                *
!*  --------------------------                                       *
!*    0    0.100000000000D+01                                        *
!*    1   -0.500000000000D+00                                        *
!*    2    0.166666666667D+00                                        *
!*    4   -0.333333333333D-01                                        *
!*    6    0.238095238095D-01                                        *
!*    8   -0.333333333333D-01                                        *
!*   10    0.757575757576D-01                                        *
!*  --------------------------                                       *
!*                                                                   *
!* ----------------------------------------------------------------- *
!* REFERENCE: "Fortran Routines for Computation of Special Functions,*
!*             jin.ece.uiuc.edu/routines/routines.html".             *
!*                                                                   *
!*                                 F90 Release By J-P Moreau, Paris. *
!*                                         (www.jpmoreau.fr)         *
!*********************************************************************
  
    PROGRAM MBERNOA
    
    DOUBLE PRECISION B
    
    DIMENSION B(0:200)
    PRINT *,' '
    WRITE(*,30,advance='no')
    READ *, N

    CALL BERNOA(N,B)

    PRINT *,' '
    WRITE(*,*)'   n            Bn'
    
    WRITE(*,*)' --------------------------'
    
    WRITE(*,20)0,B(0)
    
    WRITE(*,20)1,B(1)
    
    DO K=2,N,2
      
      WRITE(*,20)K,B(K)
    END DO
    WRITE(*,*)' --------------------------'
    PRINT *,' '

20  FORMAT(2X,I3,D22.12)
30  FORMAT('  Please enter Nmax: ')
    END
 

    SUBROUTINE BERNOA(N,BN)

! ======================================
!  Purpose: Compute Bernoulli number Bn
!  Input :  n --- Serial number
!  Output:  BN(n) --- Bn
! ======================================
    
    IMPLICIT DOUBLE PRECISION (A-H,O-Z)
    
    DIMENSION BN(0:N)
    
    BN(0)=1.0D0
    
    BN(1)=-0.5D0
    
    DO M=2,N
      
      S=-(1.0D0/(M+1.0D0)-0.5D0)
      
      DO K=2,M-1
        
      R=1.0D0
        
        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.0D0
    END DO
    
    RETURN
    
    END

!end of file mbernoa.f90