'*********************************************************************
'* 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   1                                                          *
'*    1  -.5                                                         *
'*    2   .1666666666666667                                          *
'*    4  -3.333333333333335D-02                                      *
'*    6   2.380952380952382D-02                                      *
'*    8  -3.333333333333325D-02                                      *
'*   10   7.575757575757569D-02                                      *
'*  ----------------------------                                     *
'*                                                                   *
'* ----------------------------------------------------------------- *
'* REFERENCE: "Fortran Routines for Computation of Special Functions,*
'*             jin.ece.uiuc.edu/routines/routines.html".             *
'*                                                                   *
'*                               Basic Release By J-P Moreau, Paris. *
'*                                       (www.jpmoreau.fr)           *
'*********************************************************************
' PROGRAM MBERNOA;
DEFDBL A-H, O-Z
DEFINT H-I

OPTION BASE 0

    CLS
    PRINT
    INPUT "  Please enter Nmax: ", N

    DIM BN(N)

    GOSUB 1000 'call BERNOA(N,BN)

    F$ = "####"

    PRINT
    PRINT "   n            Bn"
    PRINT " ----------------------------"
    PRINT USING F$; 0; : PRINT "  "; BN(0)
    PRINT USING F$; 1; : PRINT "  "; BN(1)
FOR K = 2 TO N STEP 2
      PRINT USING F$; K; : PRINT "  "; BN(K)
    NEXT K
    PRINT " ----------------------------"
    PRINT

END 'of main program


1000 'Subroutine BERNOA(N, BN)
' ======================================
'  Purpose: Compute Bernoulli number Bn
'  Input :  n --- Serial number
'  Output:  BN(n) --- Bn
' ======================================
    BN(0) = 1!
    BN(1) = -.5
    FOR M = 2 TO N
      S = -(1! / (M + 1!) - .5)
      FOR K = 2 TO M - 1
        R = 1!
        FOR J = 2 TO K
          R = R * (J + M - K) / J
        NEXT J
        S = S - R * BN(K)
      NEXT K
      BN(M) = S
    NEXT M
    FOR M = 3 TO N STEP 2
      BN(M) = 0!
    NEXT M
RETURN

'end of file mbernoa.bas