<commit_before>﻿#pragma once

#include "numerics/newhall.hpp"

#include <vector>

#include "geometry/barycentre_calculator.hpp"
#include "glog/logging.h"
#include "numerics/fixed_arrays.hpp"
#include "quantities/elementary_functions.hpp"
#include "quantities/quantities.hpp"

namespace principia {
namespace numerics {
namespace internal_newhall {

using base::make_not_null_unique;
using geometry::Barycentre;
using quantities::Exponentiation;
using quantities::Frequency;
using quantities::Pow;
using quantities::Time;

// Only supports 8 divisions for now.
constexpr int divisions = 8;

template<typename Vector, int degree,
         template<typename, typename, int> class Evaluator>
PolynomialInMonomialBasis<Vector, Instant, degree, Evaluator> Dehomogeneize(
    FixedVector<Vector, degree + 1> const& homogeneous_coefficients,
    Frequency const& scale,
    Instant const& origin);

template<typename Vector,
         typename DehomogeneizedCoefficients, int degree, int k>
struct Dehomogeneizer {
  static void Convert(
      FixedVector<Vector, degree + 1> const& homogeneous_coefficients,
      Frequency const& scale,
      DehomogeneizedCoefficients& dehomogeneized_coefficients);
};

template<typename Vector,
         typename DehomogeneizedCoefficients, int degree>
struct Dehomogeneizer<Vector, DehomogeneizedCoefficients, degree, degree> {
  static void Convert(
      FixedVector<Vector, degree + 1> const& homogeneous_coefficients,
      Frequency const& scale,
      DehomogeneizedCoefficients& dehomogeneized_coefficients);
};

template<typename Vector, int degree,
         template<typename, typename, int> class Evaluator>
PolynomialInMonomialBasis<Vector, Instant, degree, Evaluator> Dehomogeneize(
    FixedVector<Vector, degree + 1> const& homogeneous_coefficients,
    Frequency const& scale,
    Instant const& origin) {
  using P = PolynomialInMonomialBasis<Vector, Instant, degree, Evaluator>;
  typename P::Coefficients dehomogeneized_coefficients;
  Dehomogeneizer<Vector, typename P::Coefficients, degree, /*k=*/0>::Convert(
      homogeneous_coefficients,
      scale,
      dehomogeneized_coefficients);
  return P(dehomogeneized_coefficients, origin);
}

template<typename Vector,
         typename DehomogeneizedCoefficients, int degree, int k>
void Dehomogeneizer<Vector, DehomogeneizedCoefficients, degree, k>::Convert(
    FixedVector<Vector, degree + 1> const& homogeneous_coefficients,
    Frequency const& scale,
    DehomogeneizedCoefficients& dehomogeneized_coefficients) {
  std::get<k>(dehomogeneized_coefficients) =
      homogeneous_coefficients[k] * Pow<k>(scale);
  Dehomogeneizer<Vector, DehomogeneizedCoefficients, degree, k + 1>::Convert(
      homogeneous_coefficients,
      scale,
      dehomogeneized_coefficients);
}

template<typename Vector,
         typename DehomogeneizedCoefficients, int degree>
void Dehomogeneizer<Vector, DehomogeneizedCoefficients, degree, degree>::
Convert(FixedVector<Vector, degree + 1> const& homogeneous_coefficients,
        Frequency const& scale,
        DehomogeneizedCoefficients& dehomogeneized_coefficients) {
  std::get<degree>(dehomogeneized_coefficients) =
      homogeneous_coefficients[degree] * Pow<degree>(scale);
}

template<typename Vector, int degree,
         template<typename, typename, int> class Evaluator>
struct NewhallAppromixator {
  static FixedVector<Vector, degree + 1> HomogeneousCoefficients(
      FixedVector<Vector, 2 * divisions + 2> const& qv,
      Vector& error_estimate);
};

// TODO(phl): Use macros everywhere in this file.  It will be less error-prone.

template<typename Vector,
         template<typename, typename, int> class Evaluator>
struct NewhallAppromixator<Vector, 3, Evaluator> {
  static FixedVector<Vector, 4> HomogeneousCoefficients(
      FixedVector<Vector, 2 * divisions + 2> const& qv,
      Vector& error_estimate) {
    error_estimate =
        newhall_c_matrix_чебышёв_degree_3_divisions_8_w04.row<3>() * qv;
    return newhall_c_matrix_monomial_degree_3_divisions_8_w04 * qv;
  }
};

template<typename Vector,
         template<typename, typename, int> class Evaluator>
struct NewhallAppromixator<Vector, 4, Evaluator> {
  static FixedVector<Vector, 5> HomogeneousCoefficients(
      FixedVector<Vector, 2 * divisions + 2> const& qv,
      Vector& error_estimate) {
    error_estimate =
        newhall_c_matrix_чебышёв_degree_4_divisions_8_w04.row<4>() * qv;
    return newhall_c_matrix_monomial_degree_4_divisions_8_w04 * qv;
  }
};

template<typename Vector,
         template<typename, typename, int> class Evaluator>
struct NewhallAppromixator<Vector, 5, Evaluator> {
  static FixedVector<Vector, 6> HomogeneousCoefficients(
      FixedVector<Vector, 2 * divisions + 2> const& qv,
      Vector& error_estimate) {
    error_estimate =
        newhall_c_matrix_чебышёв_degree_5_divisions_8_w04.row<5>() * qv;
    return newhall_c_matrix_monomial_degree_5_divisions_8_w04 * qv;
  }
};

template<typename Vector,
         template<typename, typename, int> class Evaluator>
struct NewhallAppromixator<Vector, 6, Evaluator> {
  static FixedVector<Vector, 7> HomogeneousCoefficients(
      FixedVector<Vector, 2 * divisions + 2> const& qv,
      Vector& error_estimate) {
    error_estimate =
        newhall_c_matrix_чебышёв_degree_6_divisions_8_w04.row<6>() * qv;
    return newhall_c_matrix_monomial_degree_6_divisions_8_w04 * qv;
  }
};

template<typename Vector,
         template<typename, typename, int> class Evaluator>
struct NewhallAppromixator<Vector, 7, Evaluator> {
  static FixedVector<Vector, 8> HomogeneousCoefficients(
      FixedVector<Vector, 2 * divisions + 2> const& qv,
      Vector& error_estimate) {
    error_estimate =
        newhall_c_matrix_чебышёв_degree_7_divisions_8_w04.row<7>() * qv;
    return newhall_c_matrix_monomial_degree_7_divisions_8_w04 * qv;
  }
};

template<typename Vector,
         template<typename, typename, int> class Evaluator>
struct NewhallAppromixator<Vector, 8, Evaluator> {
  static FixedVector<Vector, 9> HomogeneousCoefficients(
      FixedVector<Vector, 2 * divisions + 2> const& qv,
      Vector& error_estimate) {
    error_estimate =
        newhall_c_matrix_чебышёв_degree_8_divisions_8_w04.row<8>() * qv;
    return newhall_c_matrix_monomial_degree_8_divisions_8_w04 * qv;
  }
};

template<typename Vector,
         template<typename, typename, int> class Evaluator>
struct NewhallAppromixator<Vector, 9, Evaluator> {
  static FixedVector<Vector, 10> HomogeneousCoefficients(
      FixedVector<Vector, 2 * divisions + 2> const& qv,
      Vector& error_estimate) {
    error_estimate =
        newhall_c_matrix_чебышёв_degree_9_divisions_8_w04.row<9>() * qv;
    return newhall_c_matrix_monomial_degree_9_divisions_8_w04 * qv;
  }
};

template<typename Vector,
         template<typename, typename, int> class Evaluator>
struct NewhallAppromixator<Vector, 10, Evaluator> {
  static FixedVector<Vector, 11> HomogeneousCoefficients(
      FixedVector<Vector, 2 * divisions + 2> const& qv,
      Vector& error_estimate) {
    error_estimate =
        newhall_c_matrix_чебышёв_degree_10_divisions_8_w04.row<10>() * qv;
    return newhall_c_matrix_monomial_degree_10_divisions_8_w04 * qv;
  }
};

template<typename Vector,
         template<typename, typename, int> class Evaluator>
struct NewhallAppromixator<Vector, 11, Evaluator> {
  static FixedVector<Vector, 12> HomogeneousCoefficients(
      FixedVector<Vector, 2 * divisions + 2> const& qv,
      Vector& error_estimate) {
    error_estimate =
        newhall_c_matrix_чебышёв_degree_11_divisions_8_w04.row<11>() * qv;
    return newhall_c_matrix_monomial_degree_11_divisions_8_w04 * qv;
  }
};

template<typename Vector,
         template<typename, typename, int> class Evaluator>
struct NewhallAppromixator<Vector, 12, Evaluator> {
  static FixedVector<Vector, 13> HomogeneousCoefficients(
      FixedVector<Vector, 2 * divisions + 2> const& qv,
      Vector& error_estimate) {
    error_estimate =
        newhall_c_matrix_чебышёв_degree_12_divisions_8_w04.row<12>() * qv;
    return newhall_c_matrix_monomial_degree_12_divisions_8_w04 * qv;
  }
};

template<typename Vector,
         template<typename, typename, int> class Evaluator>
struct NewhallAppromixator<Vector, 13, Evaluator> {
  static FixedVector<Vector, 14> HomogeneousCoefficients(
      FixedVector<Vector, 2 * divisions + 2> const& qv,
      Vector& error_estimate) {
    error_estimate =
        newhall_c_matrix_чебышёв_degree_13_divisions_8_w04.row<13>() * qv;
    return newhall_c_matrix_monomial_degree_13_divisions_8_w04 * qv;
  }
};

template<typename Vector,
         template<typename, typename, int> class Evaluator>
struct NewhallAppromixator<Vector, 14, Evaluator> {
  static FixedVector<Vector, 15> HomogeneousCoefficients(
      FixedVector<Vector, 2 * divisions + 2> const& qv,
      Vector& error_estimate) {
    error_estimate =
        newhall_c_matrix_чебышёв_degree_14_divisions_8_w04.row<14>() * qv;
    return newhall_c_matrix_monomial_degree_14_divisions_8_w04 * qv;
  }
};

template<typename Vector,
         template<typename, typename, int> class Evaluator>
struct NewhallAppromixator<Vector, 15, Evaluator> {
  static FixedVector<Vector, 16> HomogeneousCoefficients(
      FixedVector<Vector, 2 * divisions + 2> const& qv,
      Vector& error_estimate) {
    error_estimate =
        newhall_c_matrix_чебышёв_degree_15_divisions_8_w04.row<15>() * qv;
    return newhall_c_matrix_monomial_degree_15_divisions_8_w04 * qv;
  }
};

template<typename Vector,
         template<typename, typename, int> class Evaluator>
struct NewhallAppromixator<Vector, 16, Evaluator> {
  static FixedVector<Vector, 17> HomogeneousCoefficients(
      FixedVector<Vector, 2 * divisions + 2> const& qv,
      Vector& error_estimate) {
    error_estimate =
        newhall_c_matrix_чебышёв_degree_16_divisions_8_w04.row<16>() * qv;
    return newhall_c_matrix_monomial_degree_16_divisions_8_w04 * qv;
  }
};

template<typename Vector,
         template<typename, typename, int> class Evaluator>
struct NewhallAppromixator<Vector, 17, Evaluator> {
  static FixedVector<Vector, 18> HomogeneousCoefficients(
      FixedVector<Vector, 2 * divisions + 2> const& qv,
      Vector& error_estimate) {
    error_estimate =
        newhall_c_matrix_чебышёв_degree_17_divisions_8_w04.row<17>() * qv;
    return newhall_c_matrix_monomial_degree_17_divisions_8_w04 * qv;
  }
};

template<typename Vector>
ЧебышёвSeries<Vector>
NewhallApproximationInЧебышёвBasis(int degree,
                                   std::vector<Vector> const& q,
                                   std::vector<Variation<Vector>> const& v,
                                   Instant const& t_min,
                                   Instant const& t_max,
                                   Vector& error_estimate) {
  CHECK_EQ(divisions + 1, q.size());
  CHECK_EQ(divisions + 1, v.size());

  Time const duration_over_two = 0.5 * (t_max - t_min);

  // Tricky.  The order in Newhall's matrices is such that the entries for the
  // largest time occur first.
  FixedVector<Vector, 2 * divisions + 2> qv;
  for (int i = 0, j = 2 * divisions;
       i < divisions + 1 && j >= 0;
       ++i, j -= 2) {
    qv[j] = q[i];
    qv[j + 1] = v[i] * duration_over_two;
  }

  std::vector<Vector> coefficients;
  coefficients.reserve(degree);
  switch (degree) {
    case 3:
      coefficients = newhall_c_matrix_чебышёв_degree_3_divisions_8_w04 * qv;
      break;
    case 4:
      coefficients = newhall_c_matrix_чебышёв_degree_4_divisions_8_w04 * qv;
      break;
    case 5:
      coefficients = newhall_c_matrix_чебышёв_degree_5_divisions_8_w04 * qv;
      break;
    case 6:
      coefficients = newhall_c_matrix_чебышёв_degree_6_divisions_8_w04 * qv;
      break;
    case 7:
      coefficients = newhall_c_matrix_чебышёв_degree_7_divisions_8_w04 * qv;
      break;
    case 8:
      coefficients = newhall_c_matrix_чебышёв_degree_8_divisions_8_w04 * qv;
      break;
    case 9:
      coefficients = newhall_c_matrix_чебышёв_degree_9_divisions_8_w04 * qv;
      break;
    case 10:
      coefficients = newhall_c_matrix_чебышёв_degree_10_divisions_8_w04 * qv;
      break;
    case 11:
      coefficients = newhall_c_matrix_чебышёв_degree_11_divisions_8_w04 * qv;
      break;
    case 12:
      coefficients = newhall_c_matrix_чебышёв_degree_12_divisions_8_w04 * qv;
      break;
    case 13:
      coefficients = newhall_c_matrix_чебышёв_degree_13_divisions_8_w04 * qv;
      break;
    case 14:
      coefficients = newhall_c_matrix_чебышёв_degree_14_divisions_8_w04 * qv;
      break;
    case 15:
      coefficients = newhall_c_matrix_чебышёв_degree_15_divisions_8_w04 * qv;
      break;
    case 16:
      coefficients = newhall_c_matrix_чебышёв_degree_16_divisions_8_w04 * qv;
      break;
    case 17:
      coefficients = newhall_c_matrix_чебышёв_degree_17_divisions_8_w04 * qv;
      break;
    default:
      LOG(FATAL) << "Unexpected degree " << degree;
      break;
  }
  CHECK_EQ(degree + 1, coefficients.size());
  error_estimate = coefficients[degree];
  return ЧебышёвSeries<Vector>(coefficients, t_min, t_max);
}

template<typename Vector, int degree,
         template<typename, typename, int> class Evaluator>
PolynomialInMonomialBasis<Vector, Instant, degree, Evaluator>
NewhallApproximationInMonomialBasis(std::vector<Vector> const& q,
                                    std::vector<Variation<Vector>> const& v,
                                    Instant const& t_min,
                                    Instant const& t_max,
                                    Vector& error_estimate) {
  CHECK_EQ(divisions + 1, q.size());
  CHECK_EQ(divisions + 1, v.size());

  Time const duration_over_two = 0.5 * (t_max - t_min);

  // Tricky.  The order in Newhall's matrices is such that the entries for the
  // largest time occur first.
  FixedVector<Vector, 2 * divisions + 2> qv;
  for (int i = 0, j = 2 * divisions;
       i < divisions + 1 && j >= 0;
       ++i, j -= 2) {
    qv[j] = q[i];
    qv[j + 1] = v[i] * duration_over_two;
  }

  Instant const t_mid = Barycentre<Instant, double>({t_min, t_max}, {1, 1});
  return Dehomogeneize<Vector, degree, Evaluator>(
             NewhallAppromixator<Vector, degree, Evaluator>::
                 HomogeneousCoefficients(qv, error_estimate),
             /*scale=*/1.0 / duration_over_two,
             t_mid);
}

template<typename Vector,
         template<typename, typename, int> class Evaluator>
not_null<std::unique_ptr<Polynomial<Vector, Instant>>>
NewhallApproximationInMonomialBasis(int degree,
                                    std::vector<Vector> const& q,
                                    std::vector<Variation<Vector>> const& v,
                                    Instant const& t_min,
                                    Instant const& t_max,
                                    Vector& error_estimate) {
  switch (degree) {
    case 3:
      return make_not_null_unique<
                 PolynomialInMonomialBasis<Vector, Instant, 3, Evaluator>>(
                     NewhallApproximationInMonomialBasis<Vector, 3, Evaluator>(
                         q, v,
                         t_min, t_max,
                         error_estimate));
    case 4:
      return make_not_null_unique<
                 PolynomialInMonomialBasis<Vector, Instant, 4, Evaluator>>(
                     NewhallApproximationInMonomialBasis<Vector, 4, Evaluator>(
                         q, v,
                         t_min, t_max,
                         error_estimate));
    case 5:
      return make_not_null_unique<
                 PolynomialInMonomialBasis<Vector, Instant, 5, Evaluator>>(
                     NewhallApproximationInMonomialBasis<Vector, 5, Evaluator>(
                         q, v,
                         t_min, t_max,
                         error_estimate));
    case 6:
      return make_not_null_unique<
                 PolynomialInMonomialBasis<Vector, Instant, 6, Evaluator>>(
                     NewhallApproximationInMonomialBasis<Vector, 6, Evaluator>(
                         q, v,
                         t_min, t_max,
                         error_estimate));
    case 7:
      return make_not_null_unique<
                 PolynomialInMonomialBasis<Vector, Instant, 7, Evaluator>>(
                     NewhallApproximationInMonomialBasis<Vector, 7, Evaluator>(
                         q, v,
                         t_min, t_max,
                         error_estimate));
    case 8:
      return make_not_null_unique<
                 PolynomialInMonomialBasis<Vector, Instant, 8, Evaluator>>(
                     NewhallApproximationInMonomialBasis<Vector, 8, Evaluator>(
                         q, v,
                         t_min, t_max,
                         error_estimate));
    case 9:
      return make_not_null_unique<
                 PolynomialInMonomialBasis<Vector, Instant, 9, Evaluator>>(
                     NewhallApproximationInMonomialBasis<Vector, 9, Evaluator>(
                         q, v,
                         t_min, t_max,
                         error_estimate));
    case 10:
      return make_not_null_unique<
                 PolynomialInMonomialBasis<Vector, Instant, 10, Evaluator>>(
                     NewhallApproximationInMonomialBasis<Vector, 10, Evaluator>(
                         q, v,
                         t_min, t_max,
                         error_estimate));
    case 11:
      return make_not_null_unique<
                 PolynomialInMonomialBasis<Vector, Instant, 11, Evaluator>>(
                     NewhallApproximationInMonomialBasis<Vector, 11, Evaluator>(
                         q, v,
                         t_min, t_max,
                         error_estimate));
    case 12:
      return make_not_null_unique<
                 PolynomialInMonomialBasis<Vector, Instant, 12, Evaluator>>(
                     NewhallApproximationInMonomialBasis<Vector, 12, Evaluator>(
                         q, v,
                         t_min, t_max,
                         error_estimate));
    case 13:
      return make_not_null_unique<
                 PolynomialInMonomialBasis<Vector, Instant, 13, Evaluator>>(
                     NewhallApproximationInMonomialBasis<Vector, 13, Evaluator>(
                         q, v,
                         t_min, t_max,
                         error_estimate));
    case 14:
      return make_not_null_unique<
                 PolynomialInMonomialBasis<Vector, Instant, 14, Evaluator>>(
                     NewhallApproximationInMonomialBasis<Vector, 14, Evaluator>(
                         q, v,
                         t_min, t_max,
                         error_estimate));
    case 15:
      return make_not_null_unique<
                 PolynomialInMonomialBasis<Vector, Instant, 15, Evaluator>>(
                     NewhallApproximationInMonomialBasis<Vector, 15, Evaluator>(
                         q, v,
                         t_min, t_max,
                         error_estimate));
    case 16:
      return make_not_null_unique<
                 PolynomialInMonomialBasis<Vector, Instant, 16, Evaluator>>(
                     NewhallApproximationInMonomialBasis<Vector, 16, Evaluator>(
                         q, v,
                         t_min, t_max,
                         error_estimate));
    case 17:
      return make_not_null_unique<
                 PolynomialInMonomialBasis<Vector, Instant, 17, Evaluator>>(
                     NewhallApproximationInMonomialBasis<Vector, 17, Evaluator>(
                         q, v,
                         t_min, t_max,
                         error_estimate));
    default:
      LOG(FATAL) << "Unexpected degree " << degree;
      break;
  }
}

}  // namespace internal_newhall
}  // namespace numerics
}  // namespace principia
<commit_msg>Use macros for the repetitive code in Newhall.<commit_after>﻿#pragma once

#include "numerics/newhall.hpp"

#include <vector>

#include "geometry/barycentre_calculator.hpp"
#include "glog/logging.h"
#include "numerics/fixed_arrays.hpp"
#include "quantities/elementary_functions.hpp"
#include "quantities/quantities.hpp"

namespace principia {
namespace numerics {
namespace internal_newhall {

using base::make_not_null_unique;
using geometry::Barycentre;
using quantities::Exponentiation;
using quantities::Frequency;
using quantities::Pow;
using quantities::Time;

// Only supports 8 divisions for now.
constexpr int divisions = 8;

template<typename Vector, int degree,
         template<typename, typename, int> class Evaluator>
PolynomialInMonomialBasis<Vector, Instant, degree, Evaluator> Dehomogeneize(
    FixedVector<Vector, degree + 1> const& homogeneous_coefficients,
    Frequency const& scale,
    Instant const& origin);

template<typename Vector,
         typename DehomogeneizedCoefficients, int degree, int k>
struct Dehomogeneizer {
  static void Convert(
      FixedVector<Vector, degree + 1> const& homogeneous_coefficients,
      Frequency const& scale,
      DehomogeneizedCoefficients& dehomogeneized_coefficients);
};

template<typename Vector,
         typename DehomogeneizedCoefficients, int degree>
struct Dehomogeneizer<Vector, DehomogeneizedCoefficients, degree, degree> {
  static void Convert(
      FixedVector<Vector, degree + 1> const& homogeneous_coefficients,
      Frequency const& scale,
      DehomogeneizedCoefficients& dehomogeneized_coefficients);
};

template<typename Vector, int degree,
         template<typename, typename, int> class Evaluator>
PolynomialInMonomialBasis<Vector, Instant, degree, Evaluator> Dehomogeneize(
    FixedVector<Vector, degree + 1> const& homogeneous_coefficients,
    Frequency const& scale,
    Instant const& origin) {
  using P = PolynomialInMonomialBasis<Vector, Instant, degree, Evaluator>;
  typename P::Coefficients dehomogeneized_coefficients;
  Dehomogeneizer<Vector, typename P::Coefficients, degree, /*k=*/0>::Convert(
      homogeneous_coefficients,
      scale,
      dehomogeneized_coefficients);
  return P(dehomogeneized_coefficients, origin);
}

template<typename Vector,
         typename DehomogeneizedCoefficients, int degree, int k>
void Dehomogeneizer<Vector, DehomogeneizedCoefficients, degree, k>::Convert(
    FixedVector<Vector, degree + 1> const& homogeneous_coefficients,
    Frequency const& scale,
    DehomogeneizedCoefficients& dehomogeneized_coefficients) {
  std::get<k>(dehomogeneized_coefficients) =
      homogeneous_coefficients[k] * Pow<k>(scale);
  Dehomogeneizer<Vector, DehomogeneizedCoefficients, degree, k + 1>::Convert(
      homogeneous_coefficients,
      scale,
      dehomogeneized_coefficients);
}

template<typename Vector,
         typename DehomogeneizedCoefficients, int degree>
void Dehomogeneizer<Vector, DehomogeneizedCoefficients, degree, degree>::
Convert(FixedVector<Vector, degree + 1> const& homogeneous_coefficients,
        Frequency const& scale,
        DehomogeneizedCoefficients& dehomogeneized_coefficients) {
  std::get<degree>(dehomogeneized_coefficients) =
      homogeneous_coefficients[degree] * Pow<degree>(scale);
}

template<typename Vector, int degree,
         template<typename, typename, int> class Evaluator>
struct NewhallAppromixator {
  static FixedVector<Vector, degree + 1> HomogeneousCoefficients(
      FixedVector<Vector, 2 * divisions + 2> const& qv,
      Vector& error_estimate);
};

// TODO(phl): Use macros everywhere in this file.  It will be less error-prone.

#define PRINCIPIA_NEWHALL_APPROXIMATOR_SPECIALIZATION(degree)                  \
  template<typename Vector,                                                    \
           template<typename, typename, int> class Evaluator>                  \
  struct NewhallAppromixator<Vector, (degree), Evaluator> {                    \
    static FixedVector<Vector, ((degree) + 1)> HomogeneousCoefficients(        \
        FixedVector<Vector, 2 * divisions + 2> const& qv,                      \
        Vector& error_estimate) {                                              \
      error_estimate =                                                         \
          newhall_c_matrix_чебышёв_degree_##degree##_divisions_8_w04.row<      \
              (degree)>() * qv;                                                \
      return newhall_c_matrix_monomial_degree_##degree##_divisions_8_w04 * qv; \
    }                                                                          \
  }

PRINCIPIA_NEWHALL_APPROXIMATOR_SPECIALIZATION(3);
PRINCIPIA_NEWHALL_APPROXIMATOR_SPECIALIZATION(4);
PRINCIPIA_NEWHALL_APPROXIMATOR_SPECIALIZATION(5);
PRINCIPIA_NEWHALL_APPROXIMATOR_SPECIALIZATION(6);
PRINCIPIA_NEWHALL_APPROXIMATOR_SPECIALIZATION(7);
PRINCIPIA_NEWHALL_APPROXIMATOR_SPECIALIZATION(8);
PRINCIPIA_NEWHALL_APPROXIMATOR_SPECIALIZATION(9);
PRINCIPIA_NEWHALL_APPROXIMATOR_SPECIALIZATION(10);
PRINCIPIA_NEWHALL_APPROXIMATOR_SPECIALIZATION(11);
PRINCIPIA_NEWHALL_APPROXIMATOR_SPECIALIZATION(12);
PRINCIPIA_NEWHALL_APPROXIMATOR_SPECIALIZATION(13);
PRINCIPIA_NEWHALL_APPROXIMATOR_SPECIALIZATION(14);
PRINCIPIA_NEWHALL_APPROXIMATOR_SPECIALIZATION(15);
PRINCIPIA_NEWHALL_APPROXIMATOR_SPECIALIZATION(16);
PRINCIPIA_NEWHALL_APPROXIMATOR_SPECIALIZATION(17);

#undef PRINCIPIA_NEWHALL_APPROXIMATOR_SPECIALIZATION

#define PRINCIPIA_NEWHALL_APPROXIMATION_IN_ЧЕБЫШЁВ_BASIS_CASE(degree)    \
  case (degree):                                                         \
    coefficients =                                                       \
        newhall_c_matrix_чебышёв_degree_##degree##_divisions_8_w04 * qv; \
    break

template<typename Vector>
ЧебышёвSeries<Vector>
NewhallApproximationInЧебышёвBasis(int degree,
                                   std::vector<Vector> const& q,
                                   std::vector<Variation<Vector>> const& v,
                                   Instant const& t_min,
                                   Instant const& t_max,
                                   Vector& error_estimate) {
  CHECK_EQ(divisions + 1, q.size());
  CHECK_EQ(divisions + 1, v.size());

  Time const duration_over_two = 0.5 * (t_max - t_min);

  // Tricky.  The order in Newhall's matrices is such that the entries for the
  // largest time occur first.
  FixedVector<Vector, 2 * divisions + 2> qv;
  for (int i = 0, j = 2 * divisions;
       i < divisions + 1 && j >= 0;
       ++i, j -= 2) {
    qv[j] = q[i];
    qv[j + 1] = v[i] * duration_over_two;
  }

  std::vector<Vector> coefficients;
  coefficients.reserve(degree);
  switch (degree) {
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_ЧЕБЫШЁВ_BASIS_CASE(3);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_ЧЕБЫШЁВ_BASIS_CASE(4);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_ЧЕБЫШЁВ_BASIS_CASE(5);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_ЧЕБЫШЁВ_BASIS_CASE(6);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_ЧЕБЫШЁВ_BASIS_CASE(7);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_ЧЕБЫШЁВ_BASIS_CASE(8);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_ЧЕБЫШЁВ_BASIS_CASE(9);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_ЧЕБЫШЁВ_BASIS_CASE(10);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_ЧЕБЫШЁВ_BASIS_CASE(11);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_ЧЕБЫШЁВ_BASIS_CASE(12);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_ЧЕБЫШЁВ_BASIS_CASE(13);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_ЧЕБЫШЁВ_BASIS_CASE(14);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_ЧЕБЫШЁВ_BASIS_CASE(15);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_ЧЕБЫШЁВ_BASIS_CASE(16);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_ЧЕБЫШЁВ_BASIS_CASE(17);
    default:
      LOG(FATAL) << "Unexpected degree " << degree;
      break;
  }
  CHECK_EQ(degree + 1, coefficients.size());
  error_estimate = coefficients[degree];
  return ЧебышёвSeries<Vector>(coefficients, t_min, t_max);
}

#undef PRINCIPIA_NEWHALL_APPROXIMATION_IN_ЧЕБЫШЁВ_BASIS_CASE

template<typename Vector, int degree,
         template<typename, typename, int> class Evaluator>
PolynomialInMonomialBasis<Vector, Instant, degree, Evaluator>
NewhallApproximationInMonomialBasis(std::vector<Vector> const& q,
                                    std::vector<Variation<Vector>> const& v,
                                    Instant const& t_min,
                                    Instant const& t_max,
                                    Vector& error_estimate) {
  CHECK_EQ(divisions + 1, q.size());
  CHECK_EQ(divisions + 1, v.size());

  Time const duration_over_two = 0.5 * (t_max - t_min);

  // Tricky.  The order in Newhall's matrices is such that the entries for the
  // largest time occur first.
  FixedVector<Vector, 2 * divisions + 2> qv;
  for (int i = 0, j = 2 * divisions;
       i < divisions + 1 && j >= 0;
       ++i, j -= 2) {
    qv[j] = q[i];
    qv[j + 1] = v[i] * duration_over_two;
  }

  Instant const t_mid = Barycentre<Instant, double>({t_min, t_max}, {1, 1});
  return Dehomogeneize<Vector, degree, Evaluator>(
             NewhallAppromixator<Vector, degree, Evaluator>::
                 HomogeneousCoefficients(qv, error_estimate),
             /*scale=*/1.0 / duration_over_two,
             t_mid);
}

#define PRINCIPIA_NEWHALL_APPROXIMATION_IN_MONOMIAL_BASIS_CASE(degree)    \
  case (degree):                                                          \
    return make_not_null_unique<                                          \
        PolynomialInMonomialBasis<Vector, Instant, (degree), Evaluator>>( \
        NewhallApproximationInMonomialBasis<Vector, (degree), Evaluator>( \
            q, v,                                                         \
            t_min, t_max,                                                 \
            error_estimate))

template<typename Vector, template<typename, typename, int> class Evaluator>
not_null<std::unique_ptr<Polynomial<Vector, Instant>>>
NewhallApproximationInMonomialBasis(int degree,
                                    std::vector<Vector> const& q,
                                    std::vector<Variation<Vector>> const& v,
                                    Instant const& t_min,
                                    Instant const& t_max,
                                    Vector& error_estimate) {
  switch (degree) {
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_MONOMIAL_BASIS_CASE(3);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_MONOMIAL_BASIS_CASE(4);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_MONOMIAL_BASIS_CASE(5);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_MONOMIAL_BASIS_CASE(6);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_MONOMIAL_BASIS_CASE(7);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_MONOMIAL_BASIS_CASE(8);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_MONOMIAL_BASIS_CASE(9);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_MONOMIAL_BASIS_CASE(10);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_MONOMIAL_BASIS_CASE(11);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_MONOMIAL_BASIS_CASE(12);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_MONOMIAL_BASIS_CASE(13);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_MONOMIAL_BASIS_CASE(14);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_MONOMIAL_BASIS_CASE(15);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_MONOMIAL_BASIS_CASE(16);
    PRINCIPIA_NEWHALL_APPROXIMATION_IN_MONOMIAL_BASIS_CASE(17);
    default:
      LOG(FATAL) << "Unexpected degree " << degree;
      break;
  }
}

#undef PRINCIPIA_NEWHALL_APPROXIMATION_IN_MONOMIAL_BASIS_CASE

}  // namespace internal_newhall
}  // namespace numerics
}  // namespace principia

{'fraction_non_alphanumeric': 0.08235858101629914, 'fraction_numerical': 0.017929050814956854, 'mean_word_length': 2.9040548970679976, 'pattern_counts': {'":': 0, '<': 281, '<?xml version=': 0, '>': 265, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 51, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5906292', 'n_tokens_mistral': 10852, 'n_tokens_neox': 10144, 'n_words': 1862}