Bernoulli numbers are used in some series expansions of several functions   (trigonometric, hyperbolic, gamma, etc.),   and are extremely important in number theory and analysis.

Note that there are two definitions of Bernoulli numbers;   this task will be using the modern usage   (as per   The National Institute of Standards and Technology convention).

The   nth   Bernoulli number is expressed as   Bn.



Task
   show the Bernoulli numbers   B0   through   B60.
   suppress the output of values which are equal to zero.   (Other than   B1 , all   odd   Bernoulli numbers have a value of zero.)
   express the Bernoulli numbers as fractions  (most are improper fractions).
   the fractions should be reduced.
   index each number in some way so that it can be discerned which Bernoulli number is being displayed.
   align the solidi   (/)   if used  (extra credit).


An algorithm
The Akiyama–Tanigawa algorithm for the "second Bernoulli numbers" as taken from wikipedia is as follows:

 for m from 0 by 1 to n do
    A[m] ← 1/(m+1)
    for j from m by -1 to 1 do
      A[j-1] ← j×(A[j-1] - A[j])
  return A[0] (which is Bn)

See also
 Sequence A027641 Numerator of Bernoulli number B_n on The On-Line Encyclopedia of Integer Sequences.
 Sequence A027642 Denominator of Bernoulli number B_n on The On-Line Encyclopedia of Integer Sequences.
 Entry Bernoulli number on The Eric Weisstein's World of Mathematics (TM).
 Luschny's The Bernoulli Manifesto for a discussion on    B1   =   -½   versus   +½. 


 
(* Taken from the 'Ada 99' project, https://marquisdegeek.com/code_ada99 *)
 
program BernoulliForAda99;
 
uses BigDecimalMath; {library for arbitary high precision BCD numbers}
 
type
  Fraction = object
  private
    numerator, denominator: BigDecimal;
 
  public
    procedure assign(n, d: Int64);
    procedure subtract(rhs: Fraction);
    procedure multiply(value: Int64);
    procedure reduce();
    procedure writeOutput();
end;
 
 
function gcd(a, b: BigDecimal):BigDecimal;
begin
  if (b = 0) then begin
    gcd := a;
    end
  else begin
    gcd := gcd(b, a mod b);
 end;
end;
 
 
procedure Fraction.writeOutput();
var sign : char;
begin
  sign := ' ';
  if (numerator<0) then sign := '-';
  if (denominator<0) then sign := '-';
  write(sign + BigDecimalToStr(abs(numerator)):45);
  write(' / ');
  write(BigDecimalToStr(abs(denominator)));
end;
 
 
procedure Fraction.assign(n, d: Int64);
begin
 
  numerator := n;
  denominator := d;
end;
 
 
procedure Fraction.subtract(rhs: Fraction);
begin
  numerator := numerator * rhs.denominator;
  numerator := numerator - (rhs.numerator * denominator);
  denominator := denominator * rhs.denominator;
end;
 
 
procedure Fraction.multiply(value: Int64);
var
  temp :BigDecimal;
begin
  temp := value;
  numerator := numerator * temp;
end;
 
 
procedure Fraction.reduce();
var gcdResult: BigDecimal;
begin
  gcdResult := gcd(numerator, denominator);
  begin
    numerator := numerator div gcdResult;     (* div is Int64 division *)
    denominator := denominator div gcdResult; (* could also use round(d/r) *)
  end;  
end;
 
 
function calculateBernoulli(n: Int64) : Fraction;
var
  m, j: Int64;
  results: array of Fraction;
 
  begin
    setlength(results, 60) ; {largest value 60}
    for m:= 0 to n do
    begin
      results[m].assign(1, m+1);
 
      for j:= m downto 1 do
        begin
          results[j-1].subtract(results[j]);
          results[j-1].multiply(j);
          results[j-1].reduce();
        end;  
    end;
 
    calculateBernoulli := results[0];
end;
 
 
(* Main program starts here *)
 
var
  b: Int64;
  result: Fraction;
 
begin
  writeln('Calculating Bernoulli numbers...');
  writeln('B( 0) :                                             1 / 1');
  for b:= 1 to 60  do
    begin 
	if (b<3) or ((b mod 2) = 0) then begin
          result := calculateBernoulli(b);
          write('B(',b:2,')');
          write(' : ');
          result.writeOutput();
        writeln;
      end;
  end;
end.
 