// Gnurbs - A curve and surface library
// Copyright (C) 2008-2026 Eric Bechet
//
// See the LICENSE file for contributions and license information.
// Please report all bugs and problems to <bechet@cadxfem.org>.
//
// Contributor: Leblanc Christophe.
// cleblancad@gmail.com
//
// Reference:
// "The Nurbs Book" Les Piegl, Wayne Tiller,
// ISBN 3-540-61545-8 2nd ed. Springer-Verlarg Berlin heidelberg New York.
// pp. 241-263

#include "reparametrization.h"

#include <limits>
#include <stdexcept>
#include <algorithm>
#include <assert.h>
#include <cmath>
#include <stdio.h>

/// \brief Compares to floatting-point numbers for near equality.
/// \source: http://en.cppreference.com/w/cpp/types/numeric_limits/epsilon
template<class T>
inline bool nearly_equal(const T &x, const T &y, const int ulp = 2)
{
  // the machine epsilon has to be scaled to the magnitude of the values used
  // and multiplied by the desired precision in ULPs (units in the last place)
  return std::abs(x-y) <= std::numeric_limits<T>::epsilon() * std::abs(x+y) * ulp
    // unless the result is subnormal
           || std::abs(x-y) < std::numeric_limits<T>::min();
}

///\brief Binomial Coefficients (n k) = n! / ((n-k)! k!)
class BinomialCoefficients
{
private:
  int NMAX;
  size_t** table;

  size_t compute(const size_t n, const size_t k)
  {
    if(n < k) return 0;
    if(n == k) return 1;
    if(n == 0) return 1;

    if(table[n][k] > 0) return table[n][k];
    else table[n][k] = compute(n-1, k-1) + compute(n-1, k);
    return table[n][k];
  }

public:
  BinomialCoefficients(const int max): NMAX(max+1)
  {
    table = new size_t*[NMAX];

    for(int i = 0; i < NMAX; ++i)
    {
      table[i] = new size_t[NMAX];

      for(int j = 0; j < NMAX; ++j)
      {
        table[i][j] = 0;
      }
    }
  }

  ~BinomialCoefficients()
  {
    for(int i = 0; i < NMAX; ++i)
      delete[] table[i];
    delete[] table;
  }

  inline size_t apply(const size_t n, const size_t k)
  { if(n < k) return 0; else return compute(n, k); }

  inline size_t coeff(const size_t n, const size_t k)
  { if(n < k) return 0; else return compute(n, k); }
};

/* Using bisection, return the root of a function or functor func known to lie between x1 and x2.
The root will be refined until its accuracy is  ̇xacc. */
template<typename FunctorType>
bool rootbissection(const FunctorType &func, const double x1, const double x2, const double xacc,
  double &rtb)
{
  double dx, xmid;
  double f = func(x1);
  double fmid = func(x2);

  if(nearly_equal(f, 0.0)) 
  {
    rtb = x1;
    return true;
  }

  if(nearly_equal(fmid, 0.0))
  {
    rtb = x2;
    return true;
  }

  if(f*fmid >= 0.0) return false;

  rtb = (f < 0.0) ? (dx = x2-x1, x1) : (dx = x1-x2, x2);

  while(true)
  {
    fmid = func(xmid = rtb + (dx *= 0.5));
    if(fmid <= 0.0) rtb = xmid;
    if(fabs(dx) < xacc || nearly_equal(fmid, 0.0)) return true;
  }

  return false;
}

/// \brief Functor for knot 'inversion'.
class ReparametrizationFunctor
{
public:
  ReparametrizationFunctor(const nbspline &c):
    curve(c), target(std::numeric_limits<double>::quiet_NaN()) {}

  double operator()(const double s) const
  {
    npoint p;
    curve.P(s, p);
    return (p[0]-target); // No perspective division.
  }

  double get_target() const
  { return target; }

  void set_target(const double t)
  { target = t; }

private:
  const nbspline& curve;
  double target;
};


//-----------------------------------------------------------------------------


bool Reparametrization::is_bspline(const nbspline &curve) const
{
  for(int i = 0; i < curve.nb_CP(); ++i)
    if(!nearly_equal(curve.CP(i)[3], 1.0))
      return false;

  return true;
}

void Reparametrization::apply()
{
  if(m_curve.nb_CP() < 2)
    throw std::runtime_error("Reparametrization requires a curve with at least two control points.");

  if(m_curve.degree() < 1)
    throw std::runtime_error("Reparametrization requires a curve of degree at least 1.");

  if(m_reparametrization.degree() < 1)
    throw std::runtime_error("Reparametrization requires a reparametrizatiobn curve of degree at least 1.");

  // Steps 1 and 2: set nodal sequence (Nurbs Book p251).
  nbspline reparametrized_curve(m_curve), saturated_curve;
  set_nodal_sequence(reparametrized_curve, saturated_curve);

  // Step 3: compute the control points (Nurbs Book p251).
  set_control_points(reparametrized_curve, saturated_curve);

  // Step 4: knots removal (Nurbs Book p251).
  reparametrized_curve.clean(std::numeric_limits<double>::epsilon());
  m_curve = reparametrized_curve;

  m_done = true;
}

void Reparametrization::set_nodal_sequence(nbspline &curve,
  nbspline &saturated_curve) const
{
  if(m_use_inverse)
  {
  // Adapted code for the case here where the roles of g and f are reversed from Nurbs Book p251.
  // Step 1.
  {
    // Loop on internal knots of the reparametrization.
    for(int i = 1; i < m_reparametrization.nb_knots()-1; ++i)
      curve.insertknot(m_reparametrization.u(i), 1);

    saturated_curve = curve;
    saturated_curve.saturate();
  }

  // Step 2.
  {
    // Form the S' nodal sequence.
    std::set<double> ui;
    std::vector<double> si;
    get_distinct_knots(curve, ui); // Here saturated_curve and curve have the same distinct knots.

    const int new_degree = curve.degree() * m_reparametrization.degree();
    const int new_number_of_knots = 2*(new_degree+1) // extremities
                                  + new_degree*(ui.size()-2); // internal knots
    const int new_number_of_cps = new_number_of_knots-1-new_degree;

    curve.reset(new_number_of_cps, new_degree);

    // Compute distinct knots si.
    si.reserve(ui.size());
    for(typename std::set<double>::const_iterator it_u = ui.begin();
        it_u != ui.end(); ++it_u)
    {
      npoint p;
      m_reparametrization.P(*it_u, p);
      si.push_back(p[0]);
    }

    // Set p*q+1 beginning knots.
    int count = 0;
    typename std::vector<double>::const_iterator it_si = si.begin();

    for(; count <= new_degree; ++count)
      curve.u(count) = *it_si;

    // Set the internal knots qi with multiplicities p*q.
    ++it_si;
    for(; count < curve.nb_knots()-1-new_degree; ++count)
    {
      curve.u(count) = *it_si;
      if(count % new_degree == 0) ++it_si;
    }

    // Set p*q+1 ending knots.
    for(; count < curve.nb_knots(); ++count)
      curve.u(count) = *it_si;
  }
  }
  else
  {
  //
  // Step 1 in Nurbs Book p251.
  //
  {
    std::set<double> si, ui;
    get_distinct_internal_knots(m_reparametrization, si);
    get_distinct_knots(curve, ui);

    // Insert knots ui = f(si) p times in Cw(u)
    for(typename std::set<double>::const_iterator it_si = si.begin();
        it_si != si.end(); ++it_si)
    {
      npoint ui;
      m_reparametrization.P(*it_si, ui);
      curve.insertknot(ui[0], curve.degree());
    }

    // Set the internal knots of U to multiplicity p.
    for(typename std::set<double>::const_iterator it_ui = ui.begin();
        it_ui != ui.end(); ++it_ui)
    {
      curve.insertknot(*it_ui, curve.degree());
    }

    saturated_curve = curve;
  }

  //
  // Step 2 in Nurbs Book p251.
  //
  {
    // Form the S' nodal sequence of knots si = g(ui) = f^{(-1)}(si).
    std::set<double> ui;
    std::vector<double> si;
    get_distinct_knots(curve, ui);

    const int new_degree = curve.degree() * m_reparametrization.degree();
    const int new_number_of_knots = 2*(new_degree+1) // extremities
                                  + new_degree*(ui.size()-2); // internal knots
    const int new_number_of_cps = new_number_of_knots-1-new_degree;

    curve.reset(new_number_of_cps, new_degree);

    // Compute distinct knots si.
    get_inverse_knots(ui, si);

    // Set p*q+1 beginning knots.
    int count = 0;
    typename std::vector<double>::const_iterator it_si = si.begin();

    for(; count <= new_degree; ++count)
      curve.u(count) = *it_si;

    // Set the internal knots qi with multiplicities p*q.
    ++it_si;
    for(; count < curve.nb_knots()-1-new_degree; ++count)
    {
      curve.u(count) = *it_si;
      if(count % new_degree == 0) ++it_si;
    }

    // Set p*q+1 ending knots.
    for(; count < curve.nb_knots(); ++count)
      curve.u(count) = *it_si;
  }
  }
}

void Reparametrization::get_distinct_internal_knots(const nbspline &curve,
  std::set<double> &knots) const
{
  knots.clear();
  int knots_index = 0;

  // Skip first knot.
  while(nearly_equal(curve.min_u(), curve.u(knots_index))) ++knots_index;

  // Insert the other knots except the last one.
  while(!nearly_equal(curve.max_u(), curve.u(knots_index)))
  {
    knots.insert(curve.u(knots_index));
    ++knots_index;
  }
}

void Reparametrization::get_distinct_knots(const nbspline &curve,
  std::set<double> &knots) const
{
  knots.clear();
  if(curve.nb_knots() < 1) return; // Empty nodal sequence.

  knots.insert(curve.u(0));
  for(int i = 1; i < curve.nb_knots(); ++i)
  {
    if(!nearly_equal(curve.u(i), *knots.rbegin()))
      knots.insert(curve.u(i));
  }
}

void Reparametrization::get_inverse_knots(const std::set<double> &ui,
  std::vector<double> &si) const
{
  si.clear();
  si.reserve(ui.size());

  ReparametrizationFunctor functor(m_reparametrization);
  const double tol = std::numeric_limits<double>::epsilon();

  for(typename std::set<double>::const_iterator it_ui = ui.begin();
      it_ui != ui.end(); ++it_ui)
  {
    double ret;
    npoint p;

    m_reparametrization.P(*it_ui, p);
    functor.set_target(p[0]);

    const bool flag = rootbissection(functor, m_reparametrization.min_u(),
      m_reparametrization.max_u(), tol*fabs(m_reparametrization.max_u()), ret);

    if(!flag)
      throw std::runtime_error("Reparametrization: unable to invert nodal sequence.");

    si.push_back(ret);
  }
}

void Reparametrization::set_control_points(nbspline &reparametrized_curve,
  const nbspline &saturated_curve) const
{
  const npoint nan_point = npoint(0.0, 0.0, 0.0, 0.0)
                         * std::numeric_limits<double>::quiet_NaN();
  const int p = m_curve.degree();
  const int q = m_reparametrization.degree();

  // Initially, set all the control points of the reparametrized curve to NaN.
  for(int i = 0; i < reparametrized_curve.nb_CP(); ++i)
    reparametrized_curve.CP(i) = nan_point;

  // Step 1 in Nurbs Book p247.
  // Compute the interpolant control points for each Bezier segment.
  // (First an last CPs of each segment).
  {
    int i = 0;
    while((i*p) < saturated_curve.nb_CP())
    {
      reparametrized_curve.CP(i*p*q) = saturated_curve.CP(i*p);
      ++i;
    }
  }

  // Done for degree 1 reparametrizations.
  if(m_reparametrization.degree() == 1) return;

  // Steps 2 and 3 in Nurbs Book p247.
  // Compute the needed derivatives of the curve and the reparametrization.
  std::vector< std::vector<npoint> > C; // C[i][j] is the control point of the i-th hodograph at the j-th distinct knot.
  std::vector< std::vector<double> > f; // f[i][j] is the i-th derivative at the j-th distinct knot.

  std::vector<double> Sprime_vector;

  {
    const int n = p*q; // Maximum derivative order to compute.

    std::set<double> Uprime, Sprime; // Distinct knots and their 'inverse'.
    get_distinct_knots(saturated_curve, Uprime);
    get_distinct_knots(reparametrized_curve, Sprime);
    assert(Uprime.size() == Sprime.size());

    C.resize(n+1, std::vector<npoint>(Uprime.size()));
    f.resize(n+1, std::vector<double>(Sprime.size()));

    int i, j;
    typename std::set<double>::const_iterator it_Uprime, it_Sprime;

    // Fill-in Sprime_vector with Sprime.
    Sprime_vector.reserve(Sprime.size());
    for(it_Sprime = Sprime.begin(); it_Sprime != Sprime.end(); ++it_Sprime)
      Sprime_vector.push_back(*it_Sprime);

    for(i = 0; i <= n; ++i)
    {
      for(j = 0, it_Uprime = Uprime.begin(), it_Sprime = Sprime.begin();
          j < static_cast<int>(Uprime.size()); ++j, ++it_Uprime, ++it_Sprime)
      {
        npoint pt;

        if(i > p)
          C[i][j] = npoint(0.0, 0.0, 0.0, 0.0);
        else
        {
          m_curve.dercurvepoint(*it_Uprime, i, pt);
          C[i][j] = pt;
        }

        if(i > q)
          f[i][j] = 0.0;
        else
        {
          if(m_use_inverse)
            throw std::runtime_error("Not implemented (complex problem)");
          else
            m_reparametrization.dercurvepoint(*it_Sprime, i, pt);

          f[i][j] = pt[0];
        }
      }
    }
  }

  // Step 4 in Nurbs Book p247.
  // Compute the interior control points of each Bézier segment
  // using equations 6.40, 6.41 and 6.42.
  {
    const int ml = static_cast<int>(floor((p*q+1.0)/2.0));
    const int mr = p*q-1-static_cast<int>(floor((p*q)/2.0));
    const int nb_bezier_segments = (reparametrized_curve.nb_CP()-1)/(p*q);
    BinomialCoefficients binomial(std::max(ml, mr));

    // Loop over Bézier segments.
    for(int k = 0; k < nb_bezier_segments; ++k)
    {
      const double c = Sprime_vector[k];
      const double d = Sprime_vector[k+1];
      double delta_s = (d-c);
      int sign1 = 1;
      int sign2 = 1;

      // Compute left interior control points (equation 6.41).
      for(int i = 1; i <= ml; ++i)
      {
        npoint &Qi = reparametrized_curve.CP(i+k*p*q);

        // Left term.
        double factor = delta_s;
        for(int j = 0; j < i; ++j)
          factor /= (p*q-j);

        // Compute equation 6.40 for s = c in Nurbs Book, p248.
        Qi = compute_equation640(C, f, k, i) * factor;

        // Right terms (sum over j).
        sign2 = sign1;
        for(int j = 0; j < i; ++j)
        {
          const npoint &Qj = reparametrized_curve.CP(j+k*p*q);
          const int coeff = sign2*static_cast<int>(binomial.coeff(i, j));

          Qi += Qj*coeff;
          sign2 = -sign2;
        }

        sign1 = -sign1;
        delta_s *= (d-c); // Raise by one power.
      }

      // Compute right interior control points (equation 6.42).
      sign1 = sign2 = 1;
      delta_s = (d-c);

      for(int i = 1; i <= mr; ++i)
      {
        npoint &Qi = reparametrized_curve.CP((k+1)*p*q-i);

        // Left term.
        double factor = delta_s;
        for(int j = 0; j < i; ++j)
          factor /= (p*q-j);

        // Compute equation 6.40 for s = d in Nurbs Book, p248.
        Qi = compute_equation640(C, f, k+1, i) * (factor*(-sign1));

        // Right terms (sum over j).
        sign2 = sign1;
        for(int j = 0; j < i; ++j)
        {
          const npoint &Qj = reparametrized_curve.CP((k+1)*p*q-j);
          const int coeff = sign2*static_cast<int>(binomial.coeff(i, j));

          Qi += Qj*coeff;
          sign2 = -sign2;
        }

        sign1 = -sign1;
        delta_s *= (d-c); // Raise by one power.
      }
    }
  }
}

npoint Reparametrization::compute_equation640(const std::vector< std::vector<npoint> > &C,
  const std::vector< std::vector<double> > &f, const size_t knot_index,
  const size_t n) const
{
  npoint ret(0.0, 0.0, 0.0, 0.0);

  // Simple implementation of formula 6.40 by "brute force"
  // (explore and test all the combinations) + some optimization hints.
  for(size_t j = 1; j <= n; ++j) // Optimization: j = 0 never contributes.
  {
    double sum_k = 0.0; // Result from the sum over the k's.
    std::vector<size_t> k(n);
    k.assign(k.size(), 0);

    size_t k_increment_index = 0;
    if(nearly_equal(C[j][knot_index].norm(), 0.0)) continue; // Optimization: if C(j)==(0,0): no contribution.

    while(true)
    {
      size_t sumn = 0;
      for(size_t i = 0; i < k.size(); ++i)
        sumn += (i+1)*k[i];

      {
      // Optimization:
      // For any value of k[k_increment_index] in [0, j] the constrain k1+2*k2+...+n*kn = n
      // can not be satisfied.
      if((sumn > n)  ||
         ((sumn - (k_increment_index+1)*k[k_increment_index] + (k_increment_index+1)*j) < n))
      {
        k[k_increment_index] = j; // Put current k to value j, so
                                  // k_increment_index will be increment here below.
      }
      else // The constraint k1+2*k2+...+n*kn = n can be satisfied.
      {
        while(k[k_increment_index] <= j) // Optimization: each k_i can not be bigger than j.
        {
          bool smaller_than_n_constraint_satisfied = true;

          if(constraints640_satisfied(k, j, smaller_than_n_constraint_satisfied))
          {
            sum_k += compute_equation640_helper(f, knot_index, j, k);
          }

          // Optimization:
          // If the relaxed constraint k_1+2*k_2+...+n*k_n <= n
          // (from constraint k_1+2*k_2+...+n*k_n = n) is
          // no more satisfied. It makes no sense to continue to loop
          // trying to increase the current k[k_increment_index] value
          // for violating the relaxed constraint more and more...
          if(!smaller_than_n_constraint_satisfied)
          {
            k[k_increment_index] = j; // Put current k to value j, so
                                      // k_increment_index will be increment here below.
            break;
          }
          else
          {
            if(sumn < n-1) // Optimization: Try to 'jump' closer to sumn. (k increment index is always zero)
            {
              const size_t diff = n - sumn;
              sumn -= k[k_increment_index];
              k[k_increment_index] += diff;
              sumn += k[k_increment_index];
            }
            else // Simple '+1 jump'.
              ++k[k_increment_index];
          }
        } // while(k[k_increment_index] <= j)
      } // else
      }

      // Increase by one next non-saturated counter (if any).
      ++k_increment_index;

      while((k_increment_index < n) && (k[k_increment_index] >= j))
        ++k_increment_index;

      if(k_increment_index == n) break; // All counters saturated: stop.

      // Optimization:
      // For ceil(n/2) <= k_increment_index <= n-1, the only possible values of
      // k[k_increment_index] are 0 and 1. So, they are saturated at value 1.
      if(k_increment_index >= ceil(n/2.0) && k_increment_index <= n-1)
      {
        while((k_increment_index < n) && (k[k_increment_index] >= 1))
          ++k_increment_index;

        if(k_increment_index == n) break; // All counters saturated: stop.
      }

      ++k[k_increment_index];

      // Zeroes all counters below k_increment_index
      // and set k_increment_index to zero.
      for(size_t i = 0; i < k_increment_index; ++i)
        k[i] = 0;

      k_increment_index = 0;
    } // while true

    ret += C[j][knot_index] * sum_k; // j-th derivative at given knot index.
  } // for j

  return ret;
}

double Reparametrization::compute_equation640_helper(
    const std::vector< std::vector<double> > &f, const int knot_index, const int j,
    const std::vector<size_t> &k) const
{
  const size_t n = k.size();

  //
  // Compute the factor involving factorials in equation 6.40, Nurbs Book p248.
  // Compute it under log form.
  //
  long double factor = 0.0;
  if(n <= 170)
  {
    factor = log(factorial(n));
    for(size_t i = 1; i <= n; ++i)
      factor -= (log(factorial(k[i-1])) + k[i-1]*log(factorial(i)));

    factor = exp(factor);
  }
  else // Case n > 170 (in practice should seldom happen, except if one uses a curve
       //               of degree p and a reparametrization of degree q such that p*q > 170).
  {
    factor = 0.0;
    for(size_t i = 1; i <= n; ++i)
    {
      factor += log(i);
      long double sum = 0.0;

      for(size_t j = 1; j <= k[i-1]; ++j)
        sum += log(j);

      for(size_t j = 1; j <= i; ++j)
        sum += k[i-1]*log(i);

      factor -= sum;
    }

    factor = exp(factor);
  }

  //
  // Compute one term of the sum over k's of equation 6.40.
  //
  for(size_t i = 1; i <= n; ++i)
  {
    if(k[i-1] != 0)
      factor *= pow(f[i][knot_index], k[i-1]);
  }

  return static_cast<double>(factor);
}

bool Reparametrization::constraints640_satisfied(const std::vector<size_t> &k,
  const size_t j, bool &smaller_than_n_constraint_satisfied) const
{
  const size_t n = k.size();
  size_t sumj = 0, sumn = 0;
  smaller_than_n_constraint_satisfied = true;

  size_t count = 1;
  for(typename std::vector<size_t>::const_iterator it_k = k.begin();
      it_k != k.end(); ++it_k, ++count)
  {
    sumj += (*it_k);
    sumn += (*it_k)*count;
  }

  if(sumn > n) smaller_than_n_constraint_satisfied = false;
  return (sumj == j && sumn == n);
}

long double Reparametrization::factorial(const int n) const
{
  // Factorials table for the long double type (from 0! to 170!).
  // (From the boost library).
  static const long double factorials[171] = {
      1L,
      1L,
      2L,
      6L,
      24L,
      120L,
      720L,
      5040L,
      40320L,
      362880.0L,
      3628800.0L,
      39916800.0L,
      479001600.0L,
      6227020800.0L,
      87178291200.0L,
      1307674368000.0L,
      20922789888000.0L,
      355687428096000.0L,
      6402373705728000.0L,
      121645100408832000.0L,
      0.243290200817664e19L,
      0.5109094217170944e20L,
      0.112400072777760768e22L,
      0.2585201673888497664e23L,
      0.62044840173323943936e24L,
      0.15511210043330985984e26L,
      0.403291461126605635584e27L,
      0.10888869450418352160768e29L,
      0.304888344611713860501504e30L,
      0.8841761993739701954543616e31L,
      0.26525285981219105863630848e33L,
      0.822283865417792281772556288e34L,
      0.26313083693369353016721801216e36L,
      0.868331761881188649551819440128e37L,
      0.29523279903960414084761860964352e39L,
      0.103331479663861449296666513375232e41L,
      0.3719933267899012174679994481508352e42L,
      0.137637530912263450463159795815809024e44L,
      0.5230226174666011117600072241000742912e45L,
      0.203978820811974433586402817399028973568e47L,
      0.815915283247897734345611269596115894272e48L,
      0.3345252661316380710817006205344075166515e50L,
      0.1405006117752879898543142606244511569936e52L,
      0.6041526306337383563735513206851399750726e53L,
      0.265827157478844876804362581101461589032e55L,
      0.1196222208654801945619631614956577150644e57L,
      0.5502622159812088949850305428800254892962e58L,
      0.2586232415111681806429643551536119799692e60L,
      0.1241391559253607267086228904737337503852e62L,
      0.6082818640342675608722521633212953768876e63L,
      0.3041409320171337804361260816606476884438e65L,
      0.1551118753287382280224243016469303211063e67L,
      0.8065817517094387857166063685640376697529e68L,
      0.427488328406002556429801375338939964969e70L,
      0.2308436973392413804720927426830275810833e72L,
      0.1269640335365827592596510084756651695958e74L,
      0.7109985878048634518540456474637249497365e75L,
      0.4052691950487721675568060190543232213498e77L,
      0.2350561331282878571829474910515074683829e79L,
      0.1386831185456898357379390197203894063459e81L,
      0.8320987112741390144276341183223364380754e82L,
      0.507580213877224798800856812176625227226e84L,
      0.3146997326038793752565312235495076408801e86L,
      0.1982608315404440064116146708361898137545e88L,
      0.1268869321858841641034333893351614808029e90L,
      0.8247650592082470666723170306785496252186e91L,
      0.5443449390774430640037292402478427526443e93L,
      0.3647111091818868528824985909660546442717e95L,
      0.2480035542436830599600990418569171581047e97L,
      0.1711224524281413113724683388812728390923e99L,
      0.1197857166996989179607278372168909873646e101L,
      0.8504785885678623175211676442399260102886e102L,
      0.6123445837688608686152407038527467274078e104L,
      0.4470115461512684340891257138125051110077e106L,
      0.3307885441519386412259530282212537821457e108L,
      0.2480914081139539809194647711659403366093e110L,
      0.188549470166605025498793226086114655823e112L,
      0.1451830920282858696340707840863082849837e114L,
      0.1132428117820629783145752115873204622873e116L,
      0.8946182130782975286851441715398316520698e117L,
      0.7156945704626380229481153372318653216558e119L,
      0.5797126020747367985879734231578109105412e121L,
      0.4753643337012841748421382069894049466438e123L,
      0.3945523969720658651189747118012061057144e125L,
      0.3314240134565353266999387579130131288001e127L,
      0.2817104114380550276949479442260611594801e129L,
      0.2422709538367273238176552320344125971528e131L,
      0.210775729837952771721360051869938959523e133L,
      0.1854826422573984391147968456455462843802e135L,
      0.1650795516090846108121691926245361930984e137L,
      0.1485715964481761497309522733620825737886e139L,
      0.1352001527678402962551665687594951421476e141L,
      0.1243841405464130725547532432587355307758e143L,
      0.1156772507081641574759205162306240436215e145L,
      0.1087366156656743080273652852567866010042e147L,
      0.103299784882390592625997020993947270954e149L,
      0.9916779348709496892095714015418938011582e150L,
      0.9619275968248211985332842594956369871234e152L,
      0.942689044888324774562618574305724247381e154L,
      0.9332621544394415268169923885626670049072e156L,
      0.9332621544394415268169923885626670049072e158L,
      0.9425947759838359420851623124482936749562e160L,
      0.9614466715035126609268655586972595484554e162L,
      0.990290071648618040754671525458177334909e164L,
      0.1029901674514562762384858386476504428305e167L,
      0.1081396758240290900504101305800329649721e169L,
      0.1146280563734708354534347384148349428704e171L,
      0.1226520203196137939351751701038733888713e173L,
      0.132464181945182897449989183712183259981e175L,
      0.1443859583202493582204882102462797533793e177L,
      0.1588245541522742940425370312709077287172e179L,
      0.1762952551090244663872161047107075788761e181L,
      0.1974506857221074023536820372759924883413e183L,
      0.2231192748659813646596607021218715118256e185L,
      0.2543559733472187557120132004189335234812e187L,
      0.2925093693493015690688151804817735520034e189L,
      0.339310868445189820119825609358857320324e191L,
      0.396993716080872089540195962949863064779e193L,
      0.4684525849754290656574312362808384164393e195L,
      0.5574585761207605881323431711741977155627e197L,
      0.6689502913449127057588118054090372586753e199L,
      0.8094298525273443739681622845449350829971e201L,
      0.9875044200833601362411579871448208012564e203L,
      0.1214630436702532967576624324188129585545e206L,
      0.1506141741511140879795014161993280686076e208L,
      0.1882677176888926099743767702491600857595e210L,
      0.237217324288004688567714730513941708057e212L,
      0.3012660018457659544809977077527059692324e214L,
      0.3856204823625804217356770659234636406175e216L,
      0.4974504222477287440390234150412680963966e218L,
      0.6466855489220473672507304395536485253155e220L,
      0.8471580690878820510984568758152795681634e222L,
      0.1118248651196004307449963076076169029976e225L,
      0.1487270706090685728908450891181304809868e227L,
      0.1992942746161518876737324194182948445223e229L,
      0.269047270731805048359538766214698040105e231L,
      0.3659042881952548657689727220519893345429e233L,
      0.5012888748274991661034926292112253883237e235L,
      0.6917786472619488492228198283114910358867e237L,
      0.9615723196941089004197195613529725398826e239L,
      0.1346201247571752460587607385894161555836e242L,
      0.1898143759076170969428526414110767793728e244L,
      0.2695364137888162776588507508037290267094e246L,
      0.3854370717180072770521565736493325081944e248L,
      0.5550293832739304789551054660550388118e250L,
      0.80479260574719919448490292577980627711e252L,
      0.1174997204390910823947958271638517164581e255L,
      0.1727245890454638911203498659308620231933e257L,
      0.2556323917872865588581178015776757943262e259L,
      0.380892263763056972698595524350736933546e261L,
      0.571338395644585459047893286526105400319e263L,
      0.8627209774233240431623188626544191544816e265L,
      0.1311335885683452545606724671234717114812e268L,
      0.2006343905095682394778288746989117185662e270L,
      0.308976961384735088795856467036324046592e272L,
      0.4789142901463393876335775239063022722176e274L,
      0.7471062926282894447083809372938315446595e276L,
      0.1172956879426414428192158071551315525115e279L,
      0.1853271869493734796543609753051078529682e281L,
      0.2946702272495038326504339507351214862195e283L,
      0.4714723635992061322406943211761943779512e285L,
      0.7590705053947218729075178570936729485014e287L,
      0.1229694218739449434110178928491750176572e290L,
      0.2004401576545302577599591653441552787813e292L,
      0.3287218585534296227263330311644146572013e294L,
      0.5423910666131588774984495014212841843822e296L,
      0.9003691705778437366474261723593317460744e298L,
      0.1503616514864999040201201707840084015944e301L,
      0.2526075744973198387538018869171341146786e303L,
      0.4269068009004705274939251888899566538069e305L,
      0.7257415615307998967396728211129263114717e307L,
   };

  return factorials[n];
}