#include "qc/boys.hpp"

#include <cmath>

namespace qc {

namespace {
constexpr double kPi = 3.14159265358979323846;
}

double boys_fn(int n, double t) {
  if (t < 0.0) {
    t = 0.0;
  }
  if (t < 1e-14) {
    return 1.0 / static_cast<double>(2 * n + 1);
  }
  if (t > 36.0) {
    double f = 0.5 * std::sqrt(kPi / t);
    for (int i = 1; i <= n; ++i) {
      f *= (static_cast<double>(2 * i - 1)) / (2.0 * t);
    }
    return f;
  }
  const double sqrt_t = std::sqrt(t);
  double f0 = 0.5 * std::sqrt(kPi / t) * std::erf(sqrt_t);
  if (n == 0) {
    return f0;
  }
  double fm = f0;
  for (int m = 1; m <= n; ++m) {
    fm = (2.0 * static_cast<double>(m - 1) * fm - std::exp(-t)) / (2.0 * t);
  }
  return fm;
}

}  // namespace qc