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   +½. 


/*REXX program calculates  N  number of Bernoulli numbers expressed as vulgar fractions.*/
parse arg N .;     if N=='' | N==","  then N= 60 /*Not specified?  Then use the default.*/
numeric digits max(9, n*2)                       /*increase the decimal digits if needed*/
w= max(length(N), 4);      Nw= N  + w  +  N % 4  /*used for aligning (output) fractions.*/
say 'B(n)'   center("Bernoulli numbers expressed as vulgar fractions",  max(78-w, Nw) )
say copies('─',w)  copies("─", max(78-w,Nw+2*w)) /*display 2nd line of title, separators*/
!.= .;    do #=0  to  N                          /*process the numbers from  0  ──►  N. */
          b= bern(#);      if b==0  then iterate /*calculate Bernoulli number, skip if 0*/
          indent= max(0, nW - pos('/', b) )      /*calculate the alignment (indentation)*/
          say right(#, w)  left('', indent)  b   /*display the indented Bernoulli number*/
          end   /*#*/                            /* [↑]  align the Bernoulli fractions. */
exit 0                                           /*stick a fork in it,  we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
bern: parse arg x; if x==0  then return  '1/1'   /*handle the special case of  zero.    */
                   if x==1  then return '-1/2'   /*   "    "     "      "   "  one.     */
                   if x//2  then return   0      /*   "    "     "      "   "  odds > 1.*/
        do j=2  to x  by 2;      jp= j+1         /*process the positive integers up to X*/
        sn= 1 - j                                /*define the  numerator.               */
        sd= 2                                    /*   "    "   denominator.             */
                   do k=2  to j-1  by 2          /*calculate a  SN/SD  sequence.        */
                   parse var  @.k    bn  '/'  ad /*get a previously calculated fraction.*/
                   an= comb(jp, k) * bn          /*use  COMBination  for the next term. */
                   $lcm= LCM(sd, ad)             /*use Least Common Denominator function*/
                   sn= $lcm % sd * sn;  sd= $lcm /*calculate the   current  numerator.  */
                   an= $lcm % ad * an            /*    "      "      next      "        */
                   sn= sn + an                   /*    "      "    current     "        */
                   end   /*k*/                   /* [↑]  calculate the  SN/SD  sequence.*/
        sn= -sn                                  /*flip the sign for the numerator.     */
        sd= sd * jp                              /*calculate         the denominator.   */
        if sn\==1  then do;  _= GCD(sn, sd)      /*get the  Greatest Common Denominator.*/
                            sn= sn%_;   sd= sd%_ /*reduce the numerator and denominator.*/
                        end                      /* [↑]   done with the reduction(s).   */
        @.j= sn'/'sd                             /*save the result for the next round.  */
        end              /*j*/                   /* [↑]  done calculating Bernoulli #'s.*/
      return sn'/'sd
/*──────────────────────────────────────────────────────────────────────────────────────*/
comb: procedure expose !.; parse arg x,y;   if x==y  then return 1
      if !.C.x.y\==.  then return !.C.x.y                 /*combination computed before?*/
      if   x-y  <  y  then y= x-y                         /*x-y < y?   Then use a new Y.*/
      z= perm(x, y);           do j=2  for y-1;  z= z % j
                               end   /*j*/
      !.C.x.y= z;     return z                            /*assign memoization & return.*/
/*──────────────────────────────────────────────────────────────────────────────────────*/
GCD:  procedure;               parse arg x,y;                     x= abs(x)
           do  until y==0;     parse value  x//y  y    with    y  x;  end;        return x
/*──────────────────────────────────────────────────────────────────────────────────────*/
LCM:  procedure; parse arg x,y      /*X=ABS(X);  Y=ABS(Y)   not needed for Bernoulli #s.*/
                                    /*IF Y==0 THEN RETURN 0  "    "    "      "       " */
      $= x * y                                        /*calculate part of the LCM here. */
                   do  until y==0;   parse  value   x//y  y     with     y  x
                   end   /*until*/                    /* [↑]  this is a short & fast GCD*/
      return $ % x                                    /*divide the pre─calculated value.*/
/*──────────────────────────────────────────────────────────────────────────────────────*/
perm: procedure expose !.;  parse arg x,y;          if !.P.x.y\==.  then return !.P.x.y
      z= 1;       do j=x-y+1  to x;     z= z*j;     end;        !.P.x.y= z;       return z