// Harmonic.java
// Andrew Davison, ad@coe.psu.ac.th, April 2023

/* Calculate a sum of Harmonic numbers up to n, and plot
   a curve of the sums from 1 to n.

   https://en.wikipedia.org/wiki/Harmonic_number

   In addition, use the sum to calculate an approximation for
   Euler's constant, although n has to be quite large to 
   get a reasonable value.

   https://en.wikipedia.org/wiki/Euler%27s_constant
    gamma = 0.57721566490153286060651209008240243104215933593992... 

   h(n) ~= ln(n) + gamma + 1/(2n)
   https://brilliant.org/wiki/euler-mascheroni-constant/

   Usage:
     > compile Harmonic.java
     > run Harmonic <n>
   e.g
     > run Harmonic 200
     > run Harmonic 2000
*/

import java.util.*;
import java.math.*;

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

import de.tilman_neumann.jml.transcendental.Ln;
import de.tilman_neumann.jml.precision.Precision;

import com.github.plot.Plot;
import com.github.plot.Plot.*;


public class Harmonic 
{

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

    Plot.Data plotData = Plot.data();

    BigFraction sum = BigFraction.ZERO;
    for (int i = 1; i <= n; i++) {
      sum = sum.add( BigFraction.of(1, i));
      // System.out.println(i + ". " + sum + " = " + sum.doubleValue());
      plotData.xy(i, sum.doubleValue());
    }
    System.out.println("sum = " + sum.doubleValue());

    // plot configuration
    Plot plot = Plot.plot( Plot.plotOpts().
                      title("Harmonic sum up to " + n));
    plot.xAxis("n", Plot.axisOpts().
                      format(Plot.AxisFormat.NUMBER_INT));

    // add data and draw
    plot.series("", plotData, Plot.seriesOpts().
                                 marker(Plot.Marker.NONE));
    PlotUtils.viewPlot("Harmonic Sum", plot);


    // calculate an approximation to Euler's constant
    // gamma ~=  h(n) - ln(n) - 1/(2n)  

    Precision pre = Precision.valueOf(50);
    BigDecimal hn = sum.bigDecimalValue(50, RoundingMode.HALF_UP);
    BigDecimal lnn = Ln.ln(BigDecimal.valueOf(n), pre);
    BigDecimal frac2n = BigDecimal.ONE.divide( 
                    (BigDecimal.valueOf(2)).multiply( BigDecimal.valueOf(n)));

    BigDecimal gamma = hn.subtract( lnn.add(frac2n));
    System.out.println("\nEuler's constant approx (" + n + "):\n       " + 
                                                         gamma);
    System.out.println("\nGamma: 0.57721566490153286060651209008240243104215933593992...");

  }  // end of main()

}  // end of Harmonic class