!********************************************************************* !* 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