<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 {

template<typename Vector>
ЧебышёвSeries<Vector>::ЧебышёвSeries(std::vector<Vector> const& coefficients,
                                     Instant const& t_min,
                                     Instant const& t_max)
    : wtf_(coefficients),
      coefficients_(coefficients),
      degree_(static_cast<int>(coefficients_.size()) - 1),
      t_min_(t_min),
      t_max_(t_max) {
  CHECK_LE(0, degree_) << "Degree must be at least 0";
  CHECK_LT(t_min_, t_max_) << "Time interval must not be empty";
  // Precomputed to save operations at the expense of some accuracy loss.
  Time const duration = t_max_ - t_min_;
  t_mean_ = t_min_ + 0.5 * duration;
  two_over_duration_ = 2 / duration;
}


template<typename Vector>
ЧебышёвSeries<Vector>::ЧебышёвSeries(ЧебышёвSeries&& other)
    : wtf_(other.coefficients_),
      coefficients_(std::move(other.coefficients_)),
      degree_(std::move(other.degree_)),
      t_min_(std::move(other.t_min_)),
      t_max_(std::move(other.t_max_)),
      t_mean_(std::move(other.t_mean_)),
      two_over_duration_(std::move(other.two_over_duration_)) {}

template<typename Vector>
ЧебышёвSeries<Vector>& ЧебышёвSeries<Vector>::operator=(
    ЧебышёвSeries&& other) {
  coefficients_ = std::move(other.coefficients_);
  degree_ = other.degree_;
  t_min_ = std::move(other.t_min_);
  t_max_ = std::move(other.t_max_);
  t_mean_ = std::move(other.t_mean_);
  two_over_duration_ = std::move(other.two_over_duration_);
  return *this;
}

template<typename Vector>
bool ЧебышёвSeries<Vector>::operator==(ЧебышёвSeries const& right) const {
  return coefficients_ == right.coefficients_ &&
         t_min_ == right.t_min_ &&
         t_max_ == right.t_max_;
}

template<typename Vector>
bool ЧебышёвSeries<Vector>::operator!=(ЧебышёвSeries const& right) const {
  return !ЧебышёвSeries<Vector>::operator==(right);
}

template<typename Vector>
Instant const& ЧебышёвSeries<Vector>::t_min() const {
  return t_min_;
}

template<typename Vector>
Instant const& ЧебышёвSeries<Vector>::t_max() const {
  return t_max_;
}

template<typename Vector>
Vector const& ЧебышёвSeries<Vector>::last_coefficient() const {
  return coefficients_[degree_];
}

template<typename Vector>
WTF<Vector>::WTF(std::vector<Vector> const& coefficients)
    : coefficients_(coefficients) {}

template<typename Vector>
Vector WTF<Vector>::EvaluateImplementation(
    int const degree,
    double const scaled_t) const {
  double const two_scaled_t = scaled_t + scaled_t;
  // We have to allow |scaled_t| to go slightly out of [-1, 1] because of
  // computation errors.  But if it goes too far, something is broken.
  // TODO(phl): This should use DCHECK but these macros don't work because the
  // Principia projects don't define NDEBUG.
  #ifdef _DEBUG
  CHECK_LE(scaled_t, 1.1);
  CHECK_GE(scaled_t, -1.1);
  #endif

  Vector const c0 = coefficients_[0];
  switch (degree) {
    case 0:
      return c0;
    case 1:
      return c0 + scaled_t * coefficients_[1];
    default:
      Vector b_kplus2 = coefficients_[degree];
      Vector b_kplus1 = coefficients_[degree - 1] + two_scaled_t * b_kplus2;
      Vector b_k;
      for (int k = degree - 2; k >= 1; --k) {
        b_k = coefficients_[k] + two_scaled_t * b_kplus1 - b_kplus2;
        b_kplus2 = b_kplus1;
        b_kplus1 = b_k;
      }
      return c0 + scaled_t * b_kplus1 - b_kplus2;
  }
}

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

template<typename Scalar, typename Frame, int rank>
Multivector<Scalar, Frame, rank> WTF<Multivector<Scalar, Frame, rank>>::
EvaluateImplementation(
    int const degree,
    double const scaled_t) const {
  double const two_scaled_t = scaled_t + scaled_t;
  // We have to allow |scaled_t| to go slightly out of [-1, 1] because of
  // computation errors.  But if it goes too far, something is broken.
  // TODO(phl): This should use DCHECK but these macros don't work because the
  // Principia projects don't define NDEBUG.
#ifdef _DEBUG
  CHECK_LE(scaled_t, 1.1);
  CHECK_GE(scaled_t, -1.1);
#endif

  R3Element<double> const c0 = coefficients_[0];
  double const c0x = c0.x;
  double const c0y = c0.y;
  double const c0z = c0.z;
  switch (degree) {
    case 0: {
      return Multivector<double, Frame, rank>({c0x, c0y, c0z}) * SIUnit<Scalar>();
    }
    case 1: {
      R3Element<double> const c1 = coefficients_[1];
      double const c1x = c1.x;
      double const c1y = c1.y;
      double const c1z = c1.z;
      return Multivector<double, Frame, rank>({c0x + scaled_t * c1x,
                                               c0y + scaled_t * c1y,
                                               c0z + scaled_t * c1z}) * 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>(
            {c0x + scaled_t * b_kplus1x - b_kplus2x,
              c0y + scaled_t * b_kplus1y - b_kplus2y,
              c0z + scaled_t * b_kplus1z - b_kplus2z}) * SIUnit<Scalar>();
    }
  }<commit_msg>What changed?<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 {

template<typename Vector>
ЧебышёвSeries<Vector>::ЧебышёвSeries(std::vector<Vector> const& coefficients,
                                     Instant const& t_min,
                                     Instant const& t_max)
    : wtf_(coefficients),
      coefficients_(coefficients),
      degree_(static_cast<int>(coefficients_.size()) - 1),
      t_min_(t_min),
      t_max_(t_max) {
  CHECK_LE(0, degree_) << "Degree must be at least 0";
  CHECK_LT(t_min_, t_max_) << "Time interval must not be empty";
  // Precomputed to save operations at the expense of some accuracy loss.
  Time const duration = t_max_ - t_min_;
  t_mean_ = t_min_ + 0.5 * duration;
  two_over_duration_ = 2 / duration;
}


template<typename Vector>
ЧебышёвSeries<Vector>::ЧебышёвSeries(ЧебышёвSeries&& other)
    : wtf_(other.coefficients_),
      coefficients_(std::move(other.coefficients_)),
      degree_(std::move(other.degree_)),
      t_min_(std::move(other.t_min_)),
      t_max_(std::move(other.t_max_)),
      t_mean_(std::move(other.t_mean_)),
      two_over_duration_(std::move(other.two_over_duration_)) {}

template<typename Vector>
ЧебышёвSeries<Vector>& ЧебышёвSeries<Vector>::operator=(
    ЧебышёвSeries&& other) {
  coefficients_ = std::move(other.coefficients_);
  degree_ = other.degree_;
  t_min_ = std::move(other.t_min_);
  t_max_ = std::move(other.t_max_);
  t_mean_ = std::move(other.t_mean_);
  two_over_duration_ = std::move(other.two_over_duration_);
  return *this;
}

template<typename Vector>
bool ЧебышёвSeries<Vector>::operator==(ЧебышёвSeries const& right) const {
  return coefficients_ == right.coefficients_ &&
         t_min_ == right.t_min_ &&
         t_max_ == right.t_max_;
}

template<typename Vector>
bool ЧебышёвSeries<Vector>::operator!=(ЧебышёвSeries const& right) const {
  return !ЧебышёвSeries<Vector>::operator==(right);
}

template<typename Vector>
Instant const& ЧебышёвSeries<Vector>::t_min() const {
  return t_min_;
}

template<typename Vector>
Instant const& ЧебышёвSeries<Vector>::t_max() const {
  return t_max_;
}

template<typename Vector>
Vector const& ЧебышёвSeries<Vector>::last_coefficient() const {
  return coefficients_[degree_];
}

template<typename Vector>
WTF<Vector>::WTF(std::vector<Vector> const& coefficients)
    : coefficients_(coefficients) {}

template<typename Vector>
Vector WTF<Vector>::EvaluateImplementation(
    int const degree,
    double const scaled_t) const {
  double const two_scaled_t = scaled_t + scaled_t;
  // We have to allow |scaled_t| to go slightly out of [-1, 1] because of
  // computation errors.  But if it goes too far, something is broken.
  // TODO(phl): This should use DCHECK but these macros don't work because the
  // Principia projects don't define NDEBUG.
  #ifdef _DEBUG
  CHECK_LE(scaled_t, 1.1);
  CHECK_GE(scaled_t, -1.1);
  #endif

  Vector const c0 = coefficients_[0];
  switch (degree) {
    case 0:
      return c0;
    case 1:
      return c0 + scaled_t * coefficients_[1];
    default:
      Vector b_kplus2 = coefficients_[degree];
      Vector b_kplus1 = coefficients_[degree - 1] + two_scaled_t * b_kplus2;
      Vector b_k;
      for (int k = degree - 2; k >= 1; --k) {
        b_k = coefficients_[k] + two_scaled_t * b_kplus1 - b_kplus2;
        b_kplus2 = b_kplus1;
        b_kplus1 = b_k;
      }
      return c0 + scaled_t * b_kplus1 - b_kplus2;
  }
}

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

template<typename Scalar, typename Frame, int rank>
Multivector<Scalar, Frame, rank> WTF<Multivector<Scalar, Frame, rank>>::
EvaluateImplementation(
    int const degree,
    double const scaled_t) const {
  double const two_scaled_t = scaled_t + scaled_t;
  // We have to allow |scaled_t| to go slightly out of [-1, 1] because of
  // computation errors.  But if it goes too far, something is broken.
  // TODO(phl): This should use DCHECK but these macros don't work because the
  // Principia projects don't define NDEBUG.
#ifdef _DEBUG
  CHECK_LE(scaled_t, 1.1);
  CHECK_GE(scaled_t, -1.1);
#endif

  R3Element<double> const c0 = coefficients_[0];
  double const c0x = c0.x;
  double const c0y = c0.y;
  double const c0z = c0.z;
  switch (degree) {
    case 0: {
      return Multivector<double, Frame, rank>({c0x, c0y, c0z}) * SIUnit<Scalar>();
    }
    case 1: {
      R3Element<double> const c1 = coefficients_[1];
      double const c1x = c1.x;
      double const c1y = c1.y;
      double const c1z = c1.z;
      return Multivector<double, Frame, rank>({c0x + scaled_t * c1x,
                                               c0y + scaled_t * c1y,
                                               c0z + scaled_t * c1z}) * 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>(
            {c0x + scaled_t * b_kplus1x - b_kplus2x,
              c0y + scaled_t * b_kplus1y - b_kplus2y,
              c0z + scaled_t * b_kplus1z - b_kplus2z}) * SIUnit<Scalar>();
    }
  }
{'fraction_non_alphanumeric': 0.1156980319803198, 'fraction_numerical': 0.017373923739237393, 'mean_word_length': 3.2347005208333335, 'pattern_counts': {'":': 0, '<': 105, '<?xml version=': 0, '>': 101, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6705589', 'n_tokens_mistral': 4763, 'n_tokens_neox': 4703, 'n_words': 1055}