#include <string>
#include <iostream>
#include <iomanip>
#include <cmath>

#include "TMath.h"
#include "Math/PdfFuncMathCore.h"
#include "Math/ProbFuncMathCore.h"

#include "JMath/JMathSupportkit.hh"

#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
#include "Jeep/JPrint.hh"


/**
 * \file
 * Test application of consistency between ROOT and C++ or Jpp functions.
 *
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  double       precision;
  int          debug;

  try {

    JParser<> zap;

    zap['p'] = make_field(precision)    = 1.0e-10;
    zap['d'] = make_field(debug)        = 3;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  for (int n = 1; n != 12; ++n) {
    for (int m = 1; m <= n; ++m) {

      DEBUG(setw(2) << n << ' ' << setw(2) << m << ' ' << TMath::Binomial(n,m) << ' ' << binomial(n,m) << endl);

      ASSERT(TMath::Binomial(n,m) == binomial(n,m), "Test of binomial function");
    }
  }

  for (const double x : { -1.0e10, -10.0, -1.0, 0.0, +1.0, +10.0, +1.0e10 } ) {

    DEBUG("TMath::Erf  (" << SCIENTIFIC(12,5) << x << ") " << SCIENTIFIC(12,5) << TMath::Erf(x)   << ' ' <<
	  "erf("          << SCIENTIFIC(12,5) << x << ") " << SCIENTIFIC(12,5) << erf(x)          << endl);

    ASSERT(fabs(TMath::Erf(x) - erf(x))      <= precision, "Test of erf function");
  }

  for (const double x : { -1.0e10, -10.0, -1.0, 0.0, +1.0, +10.0, +1.0e10 } ) {

    DEBUG("TMath::Erfc (" << SCIENTIFIC(12,5) << x << ") " << SCIENTIFIC(12,5) << TMath::Erfc (x) << ' ' << 
	  "erfc("         << SCIENTIFIC(12,5) << x << ") " << SCIENTIFIC(12,5) << erfc(x)         << endl);

    ASSERT(fabs(TMath::Erfc(x) - erfc(x))    <= precision, "Test of erfc function");
  }

  for (const double x : { 1.0e-10, +1.0, +10.0, 100.0 } ) {

    DEBUG("TMath::Gamma(" << SCIENTIFIC(12,5) << x << ") " << SCIENTIFIC(12,5) << TMath::Gamma(x) << ' ' <<
	  "tgamma("       << SCIENTIFIC(12,5) << x << ") " << SCIENTIFIC(12,5) << tgamma(x)       << endl);

    ASSERT(fabs(TMath::Gamma(x) - tgamma(x)) <= precision * tgamma(x), "Test of Gamma function");
  }

  for (const double a : { 1.0e-10, +1.0, +10.0, +1.0e10 } ) {
    for (const double x : { 1.0e-10, +1.0, +10.0, +1.0e10 } ) {

      DEBUG("TMath::Gamma(" << SCIENTIFIC(12,5) << a << "," << SCIENTIFIC(12,5) << x << ") " << SCIENTIFIC(12,5) << TMath::Gamma(a, x) << ' ' << 
	    "Gamma("        << SCIENTIFIC(12,5) << a << "," << SCIENTIFIC(12,5) << x << ") " << SCIENTIFIC(12,5) << Gamma(a, x)        << endl);

      ASSERT(fabs(TMath::Gamma(a,x) - Gamma(a,x)) <= precision, "Test of Gamma function");
    }
  }

  for (int n : { 0, 1, 10, 100 }) {
    for (const double mu : { 1.0e-10, +1.0, +10.0, +1.0e10 } ) {

      DEBUG("ROOT::Math::poisson_pdf(" << n << "," << SCIENTIFIC(12,5) << mu << ") " << SCIENTIFIC(12,5) << ROOT::Math::poisson_pdf(n, mu) << ' ' << 
	    "poisson("                 << n << "," << SCIENTIFIC(12,5) << mu << ") " << SCIENTIFIC(12,5) << poisson(n, mu)                 << endl);

      ASSERT(fabs(ROOT::Math::poisson_cdf(n, mu) - Poisson(n,mu)) <= precision, "Test of Poisson cumulative density function");
    }
  }

  for (int n : { 0, 1, 10, 100 }) {
    for (const double mu : { 1.0e-10, +1.0, +10.0, +1.0e10 } ) {

      DEBUG("ROOT::Math::poisson_cdf(" << n << "," << SCIENTIFIC(12,5) << mu << ") " << SCIENTIFIC(12,5) << ROOT::Math::poisson_cdf(n, mu) << ' ' << 
	    "Poisson("                 << n << "," << SCIENTIFIC(12,5) << mu << ") " << SCIENTIFIC(12,5) << Poisson(n, mu)                 << endl);

      ASSERT(fabs(ROOT::Math::poisson_cdf(n, mu) - Poisson(n,mu)) <= precision, "Test of Poisson cumulative density function");
    }
  }

  return 0;
}