/********************************************************************
* 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.000000000000                                        *
*    1       -0.500000000000                                        *
*    2        0.166666666667                                        *
*    4       -0.033333333333                                        *
*    6        0.023809523810                                        *
*    8       -0.033333333333                                        *
*   10        0.075757575758                                        *
*  --------------------------                                       *
*                                                                   *
* ----------------------------------------------------------------- *
* REFERENCE: "Fortran Routines for Computation of Special Functions,*
*             jin.ece.uiuc.edu/routines/routines.html".             *
*                                                                   *
*                                 C++ Release By J-P Moreau, Paris. *
*                                        (www.jpmoreau.fr)          *
********************************************************************/
#include <stdio.h>
#include <math.h>

#define  NMAX  100

  double B[NMAX];
  int    K,N;


  void BERNOA(int N, double *BN)  {
/*======================================
   Purpose: Compute Bernoulli number Bn
   Input :  n --- Serial number
   Output:  BN(n) --- Bn
  ======================================*/
    int J,K,M;
    double R,S;

    BN[0]=1.0;
    BN[1]=-0.5;
    for (M=2; M<=N; M++) {
      S=-(1.0/(M+1.0)-0.5);
      for (K=2; K<M; K++) {
        R=1.0;
        for (J=2; J<=K; J++) R *= 1.0*(J+M-K)/J;
        S -= R*BN[K];
      }
      BN[M]=S;
    } 
    for (M=3; M<=N; M+=2)  BN[M]=0.0;
  }


  void main()  {

    printf("\n  Please enter Nmax: ");
    scanf("%d", &N);

    BERNOA(N,B);

    printf("\n");
    printf("   n            Bn\n");
    printf(" --------------------------\n");
    printf("   0  %20.12f\n", B[0]);
    printf("   1  %20.12f\n", B[1]);
    for (K=2; K<=N; K+=2)
      printf("%4d  %20.12f\n", K, B[K]);
    printf(" --------------------------\n\n");
  }

// end of file mbernoa.cpp