// Bernoulli.java
// https://introcs.cs.princeton.edu/java/96optimization/Bernoulli.java.html
// Andrew Davison, ad@coe.psu.ac.th, April 2023

/* Bernoulli numbers appear in the Taylor expansion 
   of the tangent function, the Euler-MacLaurin 
   summation formula, and the Riemann zeta function. 

   https://en.wikipedia.org/wiki/Bernoulli_number
and
    https://www.bernoulli.org/
and
   https://mathworld.wolfram.com/BernoulliNumber.html

   They can be defined recursively by defining
   B0 = 1, B1 = -1/2
   and using the fact that the sum from k = 0 to n-1 of 
   binomial(n, k)*Bk = 0 (for n >= 2)
  
   Usage:
    > compile Bernoulli.java
    > run bernoulli <n>
   e.g.
    > run Bernoulli 15
    Bernoulli numbers (0 - 14)
      0. = 1
      1. = -1 / 2
      2. = 1 / 6
      3. = 0
      4. = -1 / 30
      5. = 0
      6. = 1 / 42
      7. = 0
      8. = -1 / 30
      9. = 0
      10. = 5 / 66
      11. = 0
      12. = -691 / 2730
      13. = 0
      14. = 7 / 6
*/

import java.math.*;

import org.apache.commons.numbers.fraction.*;

import org.apache.commons.numbers.combinatorics.BinomialCoefficient;
   // or use JML Binomial class to return BigIntegers


public class Bernoulli
{

  public static void main(String[] args)
  {
    if (args.length != 1) {
      System.out.println("Usage: java Bernoulli <n>");
      return;
    }
    int n = Nut.parseInt(args[0]);

    BigFraction[] bern = new BigFraction[n];

    bern[0] = BigFraction.of(1, 1);    // B0 == 1
    bern[1] = BigFraction.of(-1, 2);   // B1 == -1/2

    for (int k = 2; k < n; k++) {
      bern[k] = BigFraction.ZERO;  
      for (int i = 0; i < k; i++) {  // calculate sum to k
        BigFraction coef = BigFraction.of(
                   BinomialCoefficient.value(k+1, k+1-i));
        bern[k] = bern[k].subtract(coef.multiply(bern[i]));
      }
      bern[k] = bern[k].divide(BigFraction.of(k + 1));  // /(k+1)
    }

    System.out.println("Bernoulli numbers (0 - " + (n-1) + ")");
    for (int i = 0; i < n; i++)
      System.out.println("  " + i + ". = " + bern[i]);
  } // end of main()

}  // end of Bernoulli class