<commit_before>﻿
#include "numerics/чебышёв_series.hpp"

#include <vector>

#include "geometry/grassmann.hpp"
#include "geometry/r3_element.hpp"
#include "geometry/serialization.hpp"
#include "glog/logging.h"
#include "numerics/fixed_arrays.hpp"
#include "numerics/newhall.mathematica.cpp"

namespace principia {

using geometry::DoubleOrQuantityOrMultivectorSerializer;
using geometry::Multivector;
using geometry::R3Element;

namespace numerics {
namespace internal {

// The compiler does a much better job on an |R3Element<double>| than on a
// |Vector<Quantity>| so we specialize this case.
template<typename Scalar, typename Frame, int rank>
class EvaluationHelper<Multivector<Scalar, Frame, rank>> {
 public:
  EvaluationHelper(
      std::vector<Multivector<Scalar, Frame, rank>> const& coefficients,
      int const degree);
  EvaluationHelper(EvaluationHelper&& other) = default;
  EvaluationHelper& operator=(EvaluationHelper&& other) = default;

  Multivector<Scalar, Frame, rank> EvaluateImplementation(
      double const scaled_t) const;

  Multivector<Scalar, Frame, rank> coefficients(int const index) const;
  int degree() const;

 private:
  std::vector<R3Element<double>> coefficients_;
  int degree_;
};

template<typename Vector>
EvaluationHelper<Vector>::EvaluationHelper(
    std::vector<Vector> const& coefficients,
    int const degree) : coefficients_(coefficients), degree_(degree) {}

template<typename Vector>
Vector EvaluationHelper<Vector>::EvaluateImplementation(
    double const scaled_t) const {
  double const two_scaled_t = scaled_t + scaled_t;
  Vector const c_0 = coefficients_[0];
  switch (degree_) {
    case 0:
      return c_0;
    case 1:
      return c_0 + scaled_t * coefficients_[1];
    default:
      // b_degree   = c_degree.
      Vector b_i = coefficients_[degree_];
      // b_degree-1 = c_degree-1 + 2 t b_degree.
      Vector b_j = coefficients_[degree_ - 1] + two_scaled_t * b_i;
      int k = degree_ - 3;
      for (; k >= 1; k -= 2) {
        // b_k+1 = c_k+1 + 2 t b_k+2 - b_k+3.
        b_i = coefficients_[k + 1] + two_scaled_t * b_j - b_i;
        // b_k   = c_k   + 2 t b_k+1 - b_k+2.
        b_j = coefficients_[k] + two_scaled_t * b_i - b_j;
      }
      if (k == 0) {
        // b_1 = c_1 + 2 t b_2 - b_3.
        b_i = coefficients_[1] + two_scaled_t * b_j - b_i;
        // c_0 + t b_1 - b_2.
        return c_0 + scaled_t * b_i - b_j;
      } else {
        // c_0 + t b_1 - b_2.
        return c_0 + scaled_t * b_j - b_i;
      }
  }
}

template<typename Vector>
Vector EvaluationHelper<Vector>::coefficients(int const index) const {
  return coefficients_[index];
}

template<typename Vector>
int EvaluationHelper<Vector>::degree() const {
  return degree_;
}

template<typename Scalar, typename Frame, int rank>
EvaluationHelper<Multivector<Scalar, Frame, rank>>::EvaluationHelper(
    std::vector<Multivector<Scalar, Frame, rank>> const& coefficients,
    int const degree) : degree_(degree) {
  for (auto const& coefficient : coefficients) {
    coefficients_.push_back(coefficient.coordinates() / SIUnit<Scalar>());
  }
}

template<typename Scalar, typename Frame, int rank>
Multivector<Scalar, Frame, rank>
EvaluationHelper<Multivector<Scalar, Frame, rank>>::EvaluateImplementation(
    double const scaled_t) const {
  double const two_scaled_t = scaled_t + scaled_t;
  R3Element<double> const c0 = coefficients_[0];
  switch (degree_) {
    case 0: {
      return Multivector<double, Frame, rank>(c0) * SIUnit<Scalar>();
    }
    case 1: {
      R3Element<double> const c1 = coefficients_[1];
      return Multivector<double, Frame, rank>(
                 {c0.x + scaled_t * c1.x,
                  c0.y + scaled_t * c1.y,
                  c0.z + scaled_t * c1.z}) * SIUnit<Scalar>();
    }
    default: {
      R3Element<double> const cd = coefficients_[degree_];
      double b_kplus2x = cd.x;
      double b_kplus2y = cd.y;
      double b_kplus2z = cd.z;
      R3Element<double> const cdm1 = coefficients_[degree_ - 1];
      double b_kplus1x = cdm1.x + two_scaled_t * b_kplus2x;
      double b_kplus1y = cdm1.y + two_scaled_t * b_kplus2y;
      double b_kplus1z = cdm1.z + two_scaled_t * b_kplus2z;
      double b_k;
      for (int k = degree_ - 2; k >= 1; --k) {
        R3Element<double> const ck = coefficients_[k];
        b_k = ck.x + two_scaled_t * b_kplus1x - b_kplus2x;
        b_kplus2x = b_kplus1x;
        b_kplus1x = b_k;
        b_k = ck.y + two_scaled_t * b_kplus1y - b_kplus2y;
        b_kplus2y = b_kplus1y;
        b_kplus1y = b_k;
        b_k = ck.z + two_scaled_t * b_kplus1z - b_kplus2z;
        b_kplus2z = b_kplus1z;
        b_kplus1z = b_k;
      }
      return Multivector<double, Frame, rank>(
                 {c0.x + scaled_t * b_kplus1x - b_kplus2x,
                  c0.y + scaled_t * b_kplus1y - b_kplus2y,
                  c0.z + scaled_t * b_kplus1z - b_kplus2z}) * SIUnit<Scalar>();
    }
  }
}

template<typename Scalar, typename Frame, int rank>
Multivector<Scalar, Frame, rank>
EvaluationHelper<Multivector<Scalar, Frame, rank>>::coefficients(
    int const index) const {
  return Multivector<double, Frame, rank>(
             coefficients_[index]) * SIUnit<Scalar>();
}

template<typename Scalar, typename Frame, int rank>
int EvaluationHelper<Multivector<Scalar, Frame, rank>>::degree() const {
  return degree_;
}

}  // namespace internal

template<typename Vector>
ЧебышёвSeries<Vector>::ЧебышёвSeries(std::vector<Vector> const& coefficients,
                                     Instant const& t_min,
                                     Instant const& t_max)<commit_msg>About as fast and cleaner.<commit_after>﻿
#include "numerics/чебышёв_series.hpp"

#include <vector>

#include "geometry/grassmann.hpp"
#include "geometry/r3_element.hpp"
#include "geometry/serialization.hpp"
#include "glog/logging.h"
#include "numerics/fixed_arrays.hpp"
#include "numerics/newhall.mathematica.cpp"

namespace principia {

using geometry::DoubleOrQuantityOrMultivectorSerializer;
using geometry::Multivector;
using geometry::R3Element;

namespace numerics {
namespace internal {

// The compiler does a much better job on an |R3Element<double>| than on a
// |Vector<Quantity>| so we specialize this case.
template<typename Scalar, typename Frame, int rank>
class EvaluationHelper<Multivector<Scalar, Frame, rank>> {
 public:
  EvaluationHelper(
      std::vector<Multivector<Scalar, Frame, rank>> const& coefficients,
      int const degree);
  EvaluationHelper(EvaluationHelper&& other) = default;
  EvaluationHelper& operator=(EvaluationHelper&& other) = default;

  Multivector<Scalar, Frame, rank> EvaluateImplementation(
      double const scaled_t) const;

  Multivector<Scalar, Frame, rank> coefficients(int const index) const;
  int degree() const;

 private:
  std::vector<R3Element<double>> coefficients_;
  int degree_;
};

template<typename Vector>
EvaluationHelper<Vector>::EvaluationHelper(
    std::vector<Vector> const& coefficients,
    int const degree) : coefficients_(coefficients), degree_(degree) {}

template<typename Vector>
Vector EvaluationHelper<Vector>::EvaluateImplementation(
    double const scaled_t) const {
  double const two_scaled_t = scaled_t + scaled_t;
  Vector const c_0 = coefficients_[0];
  switch (degree_) {
    case 0:
      return c_0;
    case 1:
      return c_0 + scaled_t * coefficients_[1];
    default:
      // b_degree   = c_degree.
      Vector b_i = coefficients_[degree_];
      // b_degree-1 = c_degree-1 + 2 t b_degree.
      Vector b_j = coefficients_[degree_ - 1] + two_scaled_t * b_i;
      int k = degree_ - 3;
      for (; k >= 1; k -= 2) {
        // b_k+1 = c_k+1 + 2 t b_k+2 - b_k+3.
        b_i = coefficients_[k + 1] + two_scaled_t * b_j - b_i;
        // b_k   = c_k   + 2 t b_k+1 - b_k+2.
        b_j = coefficients_[k] + two_scaled_t * b_i - b_j;
      }
      if (k == 0) {
        // b_1 = c_1 + 2 t b_2 - b_3.
        b_i = coefficients_[1] + two_scaled_t * b_j - b_i;
        // c_0 + t b_1 - b_2.
        return c_0 + scaled_t * b_i - b_j;
      } else {
        // c_0 + t b_1 - b_2.
        return c_0 + scaled_t * b_j - b_i;
      }
  }
}

template<typename Vector>
Vector EvaluationHelper<Vector>::coefficients(int const index) const {
  return coefficients_[index];
}

template<typename Vector>
int EvaluationHelper<Vector>::degree() const {
  return degree_;
}

template<typename Scalar, typename Frame, int rank>
EvaluationHelper<Multivector<Scalar, Frame, rank>>::EvaluationHelper(
    std::vector<Multivector<Scalar, Frame, rank>> const& coefficients,
    int const degree) : degree_(degree) {
  for (auto const& coefficient : coefficients) {
    coefficients_.push_back(coefficient.coordinates() / SIUnit<Scalar>());
  }
}

template<typename Scalar, typename Frame, int rank>
Multivector<Scalar, Frame, rank>
EvaluationHelper<Multivector<Scalar, Frame, rank>>::EvaluateImplementation(
    double const scaled_t) const {
  double const two_scaled_t = scaled_t + scaled_t;
  R3Element<double> const c_0 = coefficients_[0];
  switch (degree_) {
    case 0:
      return Multivector<double, Frame, rank>(c_0) * SIUnit<Scalar>();
    case 1:
      return Multivector<double, Frame, rank>(
                 c_0 + scaled_t * coefficients_[1]) * SIUnit<Scalar>();
    default:
      // b_degree   = c_degree.
      R3Element<double> b_i = coefficients_[degree_];
      // b_degree-1 = c_degree-1 + 2 t b_degree.
      R3Element<double> b_j = coefficients_[degree_ - 1] + two_scaled_t * b_i;
      int k = degree_ - 3;
      for (; k >= 1; k -= 2) {
        // b_k+1 = c_k+1 + 2 t b_k+2 - b_k+3.
        R3Element<double> const c_kplus1 = coefficients_[k + 1];
        b_i.x = c_kplus1.x + two_scaled_t * b_j.x - b_i.x;
        b_i.y = c_kplus1.y + two_scaled_t * b_j.y - b_i.y;
        b_i.z = c_kplus1.z + two_scaled_t * b_j.z - b_i.z;
        // b_k   = c_k   + 2 t b_k+1 - b_k+2.
        R3Element<double> const c_k = coefficients_[k];
        b_j.x = c_k.x + two_scaled_t * b_i.x - b_j.x;
        b_j.y = c_k.y + two_scaled_t * b_i.y - b_j.y;
        b_j.z = c_k.z + two_scaled_t * b_i.z - b_j.z;
      }
      if (k == 0) {
        // b_1 = c_1 + 2 t b_2 - b_3.
        b_i = coefficients_[1] + two_scaled_t * b_j - b_i;
        // c_0 + t b_1 - b_2.
        return Multivector<double, Frame, rank>(
                   c_0 + scaled_t * b_i - b_j) * SIUnit<Scalar>();
      } else {
        // c_0 + t b_1 - b_2.
        return Multivector<double, Frame, rank>(
                   c_0 + scaled_t * b_j - b_i) * SIUnit<Scalar>();
      }
    }
}

template<typename Scalar, typename Frame, int rank>
Multivector<Scalar, Frame, rank>
EvaluationHelper<Multivector<Scalar, Frame, rank>>::coefficients(
    int const index) const {
  return Multivector<double, Frame, rank>(
             coefficients_[index]) * SIUnit<Scalar>();
}

template<typename Scalar, typename Frame, int rank>
int EvaluationHelper<Multivector<Scalar, Frame, rank>>::degree() const {
  return degree_;
}

}  // namespace internal

template<typename Vector>
ЧебышёвSeries<Vector>::ЧебышёвSeries(std::vector<Vector> const& coefficients,
                                     Instant const& t_min,
                                     Instant const& t_max)
{'fraction_non_alphanumeric': 0.1179615110477548, 'fraction_numerical': 0.016838916607270136, 'mean_word_length': 2.980496453900709, 'pattern_counts': {'":': 0, '<': 113, '<?xml version=': 0, '>': 117, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4925003', 'n_tokens_mistral': 4023, 'n_tokens_neox': 3881, 'n_words': 929}