// Visit http://www.johndcook.com/stand_alone_code.html for the source of this code and more like it.

#include "S12_Gamma.hpp"
#include <cmath>
#include <sstream>
#include <iostream>
#include <stdexcept>
#include <cfloat>

namespace SNL {
namespace S12 {

// Can delete these functions and the #include for <iostream> if not testing.

void TestGamma()
{
  struct TestCase
  {
    double input;
    double expected;
  };

  TestCase test[] =
  {
    // Test near branches in code for (0, 0.001), [0.001, 12), (12, infinity)
    {1e-20, 1e+20},
    {2.19824158876e-16, 4.5490905327e+15},  // 0.99*DBL_EPSILON
    {2.24265050974e-16, 4.45900953205e+15}, // 1.01*DBL_EPSILON
    {0.00099, 1009.52477271},
    {0.00100, 999.423772485},
    {0.00101, 989.522792258},
    {6.1, 142.451944066},
    {11.999, 39819417.4793},
    {12, 39916800.0},
    {12.001, 40014424.1571},
    {15.2, 149037380723.0}
  };

  size_t numTests = sizeof(test) / sizeof(TestCase);

  double worst_absolute_error = 0.0;
  double worst_relative_error = 0.0;
  size_t worst_absolute_error_case = 0;
  size_t worst_relative_error_case = 0;

  for (size_t t = 0; t < numTests; t++)
  {
    double computed = Gamma(test[t].input);
    double absolute_error = fabs(computed - test[t].expected);
    double relative_error = absolute_error / test[t].expected;

    if (absolute_error > worst_absolute_error)
    {
      worst_absolute_error = absolute_error;
      worst_absolute_error_case = t;
    }

    if (relative_error > worst_relative_error)
    {
      worst_relative_error = absolute_error;
      worst_relative_error_case = t;
    }
  }

  size_t t = worst_absolute_error_case;
  double x = test[t].input;
  double y = test[t].expected;
  std::cout << "Worst absolute error: "
            << fabs(Gamma(x) - y)
            << "\nGamma( "
            <<  x
            << ") computed as "
            << Gamma(x)
            << " but exact value is "
            << y
            << "\n";

  t = worst_relative_error_case;
  x = test[t].input;
  y = test[t].expected;
  std::cout << "Worst relative error: "
            << (Gamma(x) - y) / y
            << "\nGamma( "
            <<  x
            << ") computed as "
            << Gamma(x)
            << " but exact value is "
            << y
            << "\n";
}

void TestLogGamma()
{
  struct TestCase
  {
    double input;
    double expected;
  };

  TestCase test[] =
  {
    {1e-12, 27.6310211159},
    {0.9999, 5.77297915613e-05},
    {1.0001, -5.77133422205e-05},
    {3.1, 0.787375083274},
    {6.3, 5.30734288962},
    {11.9999, 17.5020635801},
    {12, 17.5023078459},
    {12.0001, 17.5025521125},
    {27.4, 62.5755868211}
  };

  size_t numTests = sizeof(test) / sizeof(TestCase);

  double worst_absolute_error = 0.0;
  double worst_relative_error = 0.0;
  size_t worst_absolute_error_case = 0;
  size_t worst_relative_error_case = 0;

  for (size_t t = 0; t < numTests; t++)
  {
    double computed = LogGamma(test[t].input);
    double absolute_error = fabs(computed - test[t].expected);
    double relative_error = absolute_error / test[t].expected;

    if (absolute_error > worst_absolute_error)
    {
      worst_absolute_error = absolute_error;
      worst_absolute_error_case = t;
    }

    if (relative_error > worst_relative_error)
    {
      worst_relative_error = absolute_error;
      worst_relative_error_case = t;
    }
  }

  size_t t = worst_absolute_error_case;
  double x = test[t].input;
  double y = test[t].expected;
  std::cout << "Worst absolute error: "
            << fabs(LogGamma(x) - y)
            << "\nGamma( "
            <<  x
            << ") computed as "
            << LogGamma(x)
            << " but exact value is "
            << y
            << "\n";

  t = worst_relative_error_case;
  x = test[t].input;
  y = test[t].expected;
  std::cout << "Worst relative error: "
            << (LogGamma(x) - y) / y
            << "\nGamma( "
            <<  x
            << ") computed as "
            << LogGamma(x)
            << " but exact value is "
            << y
            << "\n";
}

} // namespace S12
} // namespace SNL
