Skip to content
Permalink

Comparing changes

Choose two branches to see what’s changed or to start a new pull request. If you need to, you can also or learn more about diff comparisons.

Open a pull request

Create a new pull request by comparing changes across two branches. If you need to, you can also . Learn more about diff comparisons here.
base repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected base ref is valid, then try again.
Loading
base: af0daf09eaf6
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: 596aff238456
Choose a head ref

Commits on Aug 8, 2020

  1. Declaration.

    pleroy committed Aug 8, 2020
    Copy the full SHA
    ec1a04a View commit details
  2. Bibliography.

    pleroy committed Aug 8, 2020

    Unverified

    This commit is not signed, but one or more authors requires that any commit attributed to them is signed.
    Copy the full SHA
    5358007 View commit details
  3. Helper templates.

    pleroy committed Aug 8, 2020
    Copy the full SHA
    35361d1 View commit details
  4. Ordering.

    pleroy committed Aug 8, 2020
    Copy the full SHA
    a536b6f View commit details
  5. Tracing for polynomials.

    pleroy committed Aug 8, 2020
    Copy the full SHA
    4cd5d3b View commit details
  6. Nicer formatting.

    pleroy committed Aug 8, 2020
    Copy the full SHA
    779f5d7 View commit details
  7. Poisson series output.

    pleroy committed Aug 8, 2020
    Copy the full SHA
    96ca66e View commit details
  8. Poisson series output.

    pleroy committed Aug 8, 2020
    Copy the full SHA
    c8088f6 View commit details
  9. Code restructuring.

    pleroy committed Aug 8, 2020
    Copy the full SHA
    e22436c View commit details
  10. Remove One.

    pleroy committed Aug 8, 2020
    Copy the full SHA
    9030703 View commit details
  11. Transcribing the algorithm.

    pleroy committed Aug 8, 2020
    Copy the full SHA
    aabe33f View commit details

Commits on Aug 9, 2020

  1. Merge.

    pleroy committed Aug 9, 2020
    Copy the full SHA
    68dd6ed View commit details
  2. Copy the full SHA
    6e8f9f6 View commit details
  3. Use types in the test.

    pleroy committed Aug 9, 2020
    Copy the full SHA
    51c88ab View commit details
  4. Copy the full SHA
    ff558dc View commit details
  5. Return projection.

    pleroy committed Aug 9, 2020
    Copy the full SHA
    c1c0fce View commit details
  6. Copy the full SHA
    339164d View commit details
  7. += and -= for Poisson series.

    pleroy committed Aug 9, 2020
    Copy the full SHA
    0344bbd View commit details
  8. A test that fails.

    pleroy committed Aug 9, 2020
    Copy the full SHA
    fcaae29 View commit details
  9. A test that passes.

    pleroy committed Aug 9, 2020
    Copy the full SHA
    b176a54 View commit details
  10. Unicode FTW.

    pleroy committed Aug 9, 2020
    Copy the full SHA
    9cddad9 View commit details
  11. Comments and code cleanup.

    pleroy committed Aug 9, 2020
    Copy the full SHA
    de3ebc6 View commit details
  12. Test cleanup.

    pleroy committed Aug 9, 2020
    Copy the full SHA
    0130925 View commit details
  13. A TODO.

    pleroy committed Aug 9, 2020
    Copy the full SHA
    c1951c7 View commit details
  14. Another TODO.

    pleroy committed Aug 9, 2020
    Copy the full SHA
    4eb8eaf View commit details
  15. Lint.

    pleroy committed Aug 9, 2020
    Copy the full SHA
    422a431 View commit details
  16. After egg's review.

    pleroy committed Aug 9, 2020
    Copy the full SHA
    a0aaca3 View commit details

Commits on Aug 10, 2020

  1. Merge pull request #2669 from pleroy/Projection

    A first implementation of Кудрявцев projection
    pleroy authored Aug 10, 2020
    Copy the full SHA
    596aff2 View commit details
12 changes: 12 additions & 0 deletions documentation/bibliography.bib
Original file line number Diff line number Diff line change
@@ -498,6 +498,18 @@ @article{IAUWGCC2009
volume = {109},
}

@article{Kudryavtsev2007,
author = {Kudryavtsev, S.M.},
url = {http://dx.doi.org/10.1051/0004-6361:20077568},
date = {2007},
doi = {10.1051/0004-6361:20077568},
journaltitle = {Astronomy \& Astrophysics},
keywords = {celestial mechanics; method: analytical; Moon},
pages = {1069--1075},
title = {Long-term harmonic development of lunar ephemeris},
volume = {471},
}

@article{LaineyDuriezVienne2004,
author = {Lainey, V. and Duriez, L. and Vienne, A.},
url = {http://dx.doi.org/10.1051/0004-6361:20034565},
Binary file modified documentation/bibliography.pdf
Binary file not shown.
42 changes: 32 additions & 10 deletions numerics/frequency_analysis.hpp
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@


#pragma once

#include <type_traits>
@@ -20,26 +20,48 @@ using quantities::Primitive;
using quantities::Product;
using quantities::Time;

// Computes the precise mode of a quasi-periodic function, assuming that the
// mode is over the interval fft_mode (so named because it has presumably been
// obtained using FFT). The function weight is the apodization function. The
// function dot is the dot product between function and Poisson series. See
// [Cha95].
// A function that implements the dot product between a |Function| and a Poisson
// series with an apodization function |weight| function that is itself a
// Poisson series. |Function| must be a functor taking an Instant.
template<typename Function,
typename RValue, int rdegree_, int wdegree_,
template<typename, typename, int> class Evaluator>
using DotProduct =
std::function<Product<std::invoke_result_t<Function, Instant>, RValue>(
Function const& left,
PoissonSeries<RValue, rdegree_, Evaluator> const& right,
PoissonSeries<double, wdegree_, Evaluator> const& weight)>;

// Computes the precise mode of a quasi-periodic |function|, assuming that the
// mode is over the interval |fft_mode| (so named because it has presumably been
// obtained using FFT). See [Cha95].
template<typename Function,
int wdegree_,
template<typename, typename, int> class Evaluator>
AngularFrequency PreciseMode(
Interval<AngularFrequency> const& fft_mode,
Function const& function,
PoissonSeries<double, wdegree_, Evaluator> const& weight,
std::function<Product<std::invoke_result_t<Function, Instant>, RValue>(
Function const& left,
PoissonSeries<RValue, rdegree_, Evaluator> const& right,
PoissonSeries<double, wdegree_, Evaluator> const& weight)> const& dot);
DotProduct<Function, double, 0, wdegree_, Evaluator> const& dot);

// Computes the Кудрявцев projection of |function| on a basis with angular
// frequency ω and maximum degree |degree_|. See [Kud07].
// TODO(phl): We really need multiple angular frequencies.
template<typename Function,
int degree_, int wdegree_,
template<typename, typename, int> class Evaluator>
PoissonSeries<std::invoke_result_t<Function, Instant>, degree_, Evaluator>
Projection(
AngularFrequency const& ω,
Function const& function,
PoissonSeries<double, wdegree_, Evaluator> const& weight,
DotProduct<Function, std::invoke_result_t<Function, Instant>,
degree_, wdegree_, Evaluator> const& dot);

} // namespace internal_frequency_analysis

using internal_frequency_analysis::PreciseMode;
using internal_frequency_analysis::Projection;

} // namespace frequency_analysis
} // namespace numerics
212 changes: 205 additions & 7 deletions numerics/frequency_analysis_body.hpp
Original file line number Diff line number Diff line change
@@ -3,29 +3,134 @@

#include "numerics/frequency_analysis.hpp"

#include <algorithm>
#include <functional>
#include <vector>

#include "numerics/fixed_arrays.hpp"
#include "numerics/root_finders.hpp"
#include "quantities/elementary_functions.hpp"
#include "quantities/si.hpp"

namespace principia {
namespace numerics {
namespace frequency_analysis {
namespace internal_frequency_analysis {

using quantities::Inverse;
using quantities::Sqrt;
using quantities::Square;
using quantities::SquareRoot;
namespace si = quantities::si;

// A helper struct for generating the Poisson series tⁿ sin ω t and tⁿ cos ω t.
template<typename Series, int n>
struct SeriesGenerator {
// The series tⁿ sin ω t.
static Series Sin(AngularFrequency const& ω, Instant const& origin);
// The series tⁿ cos ω t.
static Series Cos(AngularFrequency const& ω, Instant const& origin);

private:
// The polynomial tⁿ.
static typename Series::Polynomial Unit(Instant const& origin);
};

// A helper struct for generating the Кудрявцев basis, i.e., functions of the
// form tⁿ sin ω t and tⁿ cos ω t properly ordered.
template<typename Series,
typename = std::make_index_sequence<Series::degree + 1>>
struct BasisGenerator;

template<typename Series, std::size_t... indices>
struct BasisGenerator<Series, std::index_sequence<indices...>> {
// TODO(phl): Will need to properly handle ω = 0.
static std::array<Series, 2 * Series::degree + 2> Basis(
AngularFrequency const& ω,
Instant const& origin);
};


template<typename Series, int n>
Series SeriesGenerator<Series, n>::Sin(AngularFrequency const& ω,
Instant const& origin) {
typename Series::Polynomial::Coefficients const zeros;
typename Series::Polynomial const zero{zeros, origin};
return Series(zero,
{{ω,
{/*sin=*/Unit(origin),
/*cos=*/zero}}});
}

template<typename Series, int n>
Series SeriesGenerator<Series, n>::Cos(AngularFrequency const& ω,
Instant const& origin) {
typename Series::Polynomial::Coefficients const zeros;
typename Series::Polynomial const zero{zeros, origin};
return Series(zero,
{{ω,
{/*sin=*/zero,
/*cos=*/Unit(origin)}}});
}

template<typename Series, int n>
typename Series::Polynomial SeriesGenerator<Series, n>::Unit(
Instant const& origin) {
typename Series::Polynomial::Coefficients coefficients;
std::get<n>(coefficients) = si::Unit<
std::tuple_element_t<n, typename Series::Polynomial::Coefficients>>;
return Series::Polynomial(coefficients, origin);
}

template<typename Series, std::size_t... indices>
std::array<Series, 2 * Series::degree + 2>
BasisGenerator<Series, std::index_sequence<indices...>>::Basis(
AngularFrequency const& ω,
Instant const& origin) {
// This has the elements {Sin(ωt), t Sin(ωt), t² Sin(ωt), ..., Cos(ωt), ...}
// which is not the order we want (we want lower-degree polynomials first).
std::array<Series, 2 * Series::degree + 2> all_series = {
SeriesGenerator<Series, indices>::Sin(ω, origin)...,
SeriesGenerator<Series, indices>::Cos(ω, origin)...};

// Order all_series by repeatedly swapping its elements.
if (Series::degree >= 2) {
// The index of this array is the current index of a series in all_series.
// The value is the index of the final resting place of that series in
// all_series. The elements at indices 0 and 2 * Series::degree + 1 are
// unused.
std::array<int, 2 * Series::degree + 2> permutation;
for (int i = 1; i < 2 * Series::degree + 1; ++i) {
permutation[i] =
i <= Series::degree ? 2 * i : 2 * (i - Series::degree) - 1;
}
for (int i = 1; i < 2 * Series::degree + 1;) {
// Swap the series currently at index i to its final resting place.
// Iterate until the series at index i is at its final resting place
// (i.e., after we have executed an entire cycle of the permutation).
// Then move to the next series.
if (i == permutation[i]) {
++i;
} else {
int const j = permutation[i];
std::swap(all_series[i], all_series[j]);
std::swap(permutation[i], permutation[j]);
}
}
}
return all_series;
}


template<typename Function,
typename RValue, int rdegree_, int wdegree_,
int wdegree_,
template<typename, typename, int> class Evaluator>
AngularFrequency PreciseMode(
Interval<AngularFrequency> const& fft_mode,
Function const& function,
PoissonSeries<double, wdegree_, Evaluator> const& weight,
std::function<Product<std::invoke_result_t<Function, Instant>, RValue>(
Function const& left,
PoissonSeries<RValue, rdegree_, Evaluator> const& right,
PoissonSeries<double, wdegree_, Evaluator> const& weight)> const& dot) {
using DotResult = Product<std::invoke_result_t<Function, Instant>, RValue>;
DotProduct<Function, double, 0, wdegree_, Evaluator> const& dot) {
using Value = std::invoke_result_t<Function, Instant>;
using Degree0 = PoissonSeries<double, 0, Evaluator>;

auto amplitude = [&dot, &function, &weight](AngularFrequency const& ω) {
@@ -46,7 +151,100 @@ AngularFrequency PreciseMode(
return GoldenSectionSearch(amplitude,
fft_mode.min,
fft_mode.max,
std::greater<Square<DotResult>>());
std::greater<Square<Value>>());
}

template<typename Function,
int degree_, int wdegree_,
template<typename, typename, int> class Evaluator>
PoissonSeries<std::invoke_result_t<Function, Instant>, degree_, Evaluator>
Projection(
AngularFrequency const& ω,
Function const& function,
PoissonSeries<double, wdegree_, Evaluator> const& weight,
DotProduct<Function, std::invoke_result_t<Function, Instant>,
degree_, wdegree_, Evaluator> const& dot) {
using Value = std::invoke_result_t<Function, Instant>;
using Series = PoissonSeries<Value, degree_, Evaluator>;

Instant const& t0 = weight.origin();
auto const basis = BasisGenerator<Series>::Basis(ω, t0);
constexpr int basis_size = std::tuple_size_v<decltype(basis)>;

// This code follows [Kud07], section 2. Our indices start at 0, unlike those
// of Кудрявцев which start at 1.
FixedLowerTriangularMatrix<Inverse<Value>, basis_size> α;
std::vector<Function> f;

// Only indices 0 to m - 1 are used in this array. At the beginning of
// iteration m it contains Aⱼ⁽ᵐ⁻¹⁾.
std::array<double, basis_size> A;

auto const F₀ = dot(function, basis[0], weight);
// TODO(phl): This does not work if basis does not have the same type as
// Function, i.e., if the degrees don't match.
auto const Q₀₀ = dot(basis[0], basis[0], weight);
α[0][0] = 1 / Sqrt(Q₀₀);
A[0] = F₀ / Q₀₀;
f.emplace_back(function - A[0] * basis[0]);
for (int m = 1; m < basis_size; ++m) {
// Contains Fₘ.
auto const F = dot(f[m - 1], basis[m], weight);

// Only indices 0 to m are used in this array. It contains Qₘⱼ.
std::array<Square<Value>, basis_size> Q;
for (int j = 0; j <= m; ++j) {
Q[j] = dot(basis[m], basis[j], weight);
}

// Only indices 0 to m - 1 are used in this array. It contains Bⱼ⁽ᵐ⁾.
std::array<Value, basis_size> B;
for (int j = 0; j < m; ++j) {
Value Σ_αⱼₛ_Qₘₛ{};
for (int s = 0; s <= j; ++s) {
Σ_αⱼₛ_Qₘₛ += α[j][s] * Q[s];
}
B[j] = -Σ_αⱼₛ_Qₘₛ;
}

{
Square<Value> Σ_Bₛ⁽ᵐ⁾²{};
for (int s = 0; s < m; ++s) {
Σ_Bₛ⁽ᵐ⁾² += B[s] * B[s];
}
DCHECK_LE(Σ_Bₛ⁽ᵐ⁾², Q[m]);
α[m][m] = 1 / Sqrt(Q[m] - Σ_Bₛ⁽ᵐ⁾²);
}

for (int j = 0; j < m; ++j) {
double Σ_Bₛ⁽ᵐ⁾_αₛⱼ = 0;
for (int s = j; s < m; ++s) {
Σ_Bₛ⁽ᵐ⁾_αₛⱼ += B[s] * α[s][j];
}
α[m][j] = α[m][m] * Σ_Bₛ⁽ᵐ⁾_αₛⱼ;
}

A[m] = α[m][m] * α[m][m] * F;

for (int j = 0; j < m; ++j) {
A[j] += α[m][m] * α[m][j] * F;
}

{
PoissonSeries<double, degree_, Evaluator> Σ_αₘᵢ_eᵢ =
α[m][0] * basis[0];
for (int i = 1; i <= m; ++i) {
Σ_αₘᵢ_eᵢ += α[m][i] * basis[i];
}
f.emplace_back(f[m - 1] - α[m][m] * F * Σ_αₘᵢ_eᵢ);
}
}

PoissonSeries<Value, degree_, Evaluator> result = A[0] * basis[0];
for (int i = 1; i < basis_size; ++i) {
result += A[i] * basis[i];
}
return result;
}

} // namespace internal_frequency_analysis
86 changes: 70 additions & 16 deletions numerics/frequency_analysis_test.cpp
Original file line number Diff line number Diff line change
@@ -14,6 +14,7 @@
#include "quantities/named_quantities.hpp"
#include "quantities/quantities.hpp"
#include "quantities/si.hpp"
#include "testing_utilities/almost_equals.hpp"
#include "testing_utilities/approximate_quantity.hpp"
#include "testing_utilities/is_near.hpp"
#include "testing_utilities/numerics_matchers.hpp"
@@ -25,14 +26,18 @@ namespace frequency_analysis {
using geometry::Instant;
using quantities::AngularFrequency;
using quantities::Length;
using quantities::Pow;
using quantities::Sin;
using quantities::Square;
using quantities::Time;
using quantities::si::Metre;
using quantities::si::Radian;
using quantities::si::Second;
using testing_utilities::AlmostEquals;
using testing_utilities::IsNear;
using testing_utilities::RelativeErrorFrom;
using testing_utilities::operator""_⑴;
namespace si = quantities::si;

class FrequencyAnalysisTest : public ::testing::Test {};

@@ -45,32 +50,32 @@ TEST_F(FrequencyAnalysisTest, PreciseMode) {
std::uniform_real_distribution<> amplitude_noise(-0.1, 0.1);
std::uniform_real_distribution<> frequency_noise(-100.0, 100.0);

using Degree0 = PoissonSeries<double, 0, HornerEvaluator>;
Degree0::PolynomialsByAngularFrequency polynomials;
using Series = PoissonSeries<Length, 0, HornerEvaluator>;
Series::PolynomialsByAngularFrequency polynomials;

// Main harmonic.
polynomials.emplace(
ω,
Degree0::Polynomials{/*sin=*/Degree0::Polynomial({1}, t0),
/*cos=*/Degree0::Polynomial({0}, t0)});
Series::Polynomials{/*sin=*/Series::Polynomial({1 * Metre}, t0),
/*cos=*/Series::Polynomial({0 * Metre}, t0)});

// Noise with lower amplitude and higher frequency.
for (int i = 0; i < 10; ++i) {
auto const sin_amplitude = amplitude_noise(random);
auto const cos_amplitude = amplitude_noise(random);
auto const sin_amplitude = amplitude_noise(random) * Metre;
auto const cos_amplitude = amplitude_noise(random) * Metre;
polynomials.emplace(
ω * frequency_noise(random),
Degree0::Polynomials{/*sin=*/Degree0::Polynomial({sin_amplitude}, t0),
/*cos=*/Degree0::Polynomial({cos_amplitude}, t0)});
Series::Polynomials{/*sin=*/Series::Polynomial({sin_amplitude}, t0),
/*cos=*/Series::Polynomial({cos_amplitude}, t0)});
}
Degree0 const sin(Degree0::Polynomial({amplitude_noise(random)}, t0),
Series const sin(Series::Polynomial({amplitude_noise(random) * Metre}, t0),
polynomials);

Instant const t_min = t0;
Instant const t_max = t0 + (FFT::size - 1) * Δt;
std::vector<Length> signal;
for (int n = 0; n < FFT::size; ++n) {
signal.push_back(sin(t0 + n * Δt) * Metre);
signal.push_back(sin(t0 + n * Δt));
}

// Won't fit on the stack.
@@ -81,12 +86,14 @@ TEST_F(FrequencyAnalysisTest, PreciseMode) {
EXPECT_THAT(mode.midpoint(),
RelativeErrorFrom(ω, IsNear(8.1e-4_⑴)));

std::function<double(Degree0 const& left,
Degree0 const& right,
Degree0 const& weight)> const dot =
[t_min, t_max](Degree0 const& left,
Degree0 const& right,
Degree0 const& weight) {
using Double = PoissonSeries<double, 0, HornerEvaluator>;

std::function<Length(Series const& left,
Double const& right,
Double const& weight)> const dot =
[t_min, t_max](Series const& left,
Double const& right,
Double const& weight) {
return Dot(left, right, weight, t_min, t_max);
};

@@ -98,6 +105,53 @@ TEST_F(FrequencyAnalysisTest, PreciseMode) {
RelativeErrorFrom(ω, IsNear(6.4e-11_⑴)));
}

TEST_F(FrequencyAnalysisTest, Projection) {
Instant const t0;
AngularFrequency const ω = 666.543 * π * Radian / Second;
std::mt19937_64 random(42);
std::uniform_real_distribution<> noise(-10.0, 10.0);

using Series = PoissonSeries<Length, 4, HornerEvaluator>;

auto random_polynomial = [&noise, &random, &t0]() {
auto const c0 = noise(random) * Metre;
auto const c1 = noise(random) * Metre / Second;
auto const c2 = noise(random) * Metre / Pow<2>(Second);
auto const c3 = noise(random) * Metre / Pow<3>(Second);
auto const c4 = noise(random) * Metre / Pow<4>(Second);

return Series::Polynomial({c0, c1, c2, c3, c4}, t0);
};

auto const sin = random_polynomial();
auto const cos = random_polynomial();
Series const series(
Series::Polynomial(Series::Polynomial::Coefficients{}, t0),
{{ω, Series::Polynomials{sin, cos}}});

Instant const t_min = t0;
Instant const t_max = t0 + 100 * Radian / ω;

using Double = PoissonSeries<double, 0, HornerEvaluator>;

std::function<Square<Length>(Series const& left,
Series const& right,
Double const& weight)> const dot =
[t_min, t_max](Series const& left,
Series const& right,
Double const& weight) {
return Dot(left, right, weight, t_min, t_max);
};

auto const projection = Projection(
ω, series, apodization::Dirichlet<HornerEvaluator>(t_min, t_max), dot);

for (int i = 0; i <= 100; ++i) {
EXPECT_THAT(projection(t0 + i * Radian / ω),
AlmostEquals(series(t0 + i * Radian / ω), 0, 2432));
}
}

} // namespace frequency_analysis
} // namespace numerics
} // namespace principia
4 changes: 4 additions & 0 deletions numerics/poisson_series.hpp
Original file line number Diff line number Diff line change
@@ -32,6 +32,7 @@ template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
class PoissonSeries {
public:
static const int degree = degree_;
using Polynomial =
numerics::PolynomialInMonomialBasis<Value, Instant, degree_, Evaluator>;

@@ -55,6 +56,9 @@ class PoissonSeries {
PoissonSeries<quantities::Primitive<Value, Time>, degree_ + 1, Evaluator>
Primitive() const;

PoissonSeries& operator+=(PoissonSeries const& right);
PoissonSeries& operator-=(PoissonSeries const& right);

private:
Instant origin_; // Common to all polynomials.
Polynomial aperiodic_;
22 changes: 20 additions & 2 deletions numerics/poisson_series_body.hpp
Original file line number Diff line number Diff line change
@@ -131,6 +131,24 @@ PoissonSeries<Value, degree_, Evaluator>::Primitive() const {
return Result{aperiodic, std::move(periodic)};
}

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
PoissonSeries<Value, degree_, Evaluator>&
PoissonSeries<Value, degree_, Evaluator>::operator+=(
PoissonSeries const& right) {
*this = *this + right;
return *this;
}

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
PoissonSeries<Value, degree_, Evaluator>&
PoissonSeries<Value, degree_, Evaluator>::operator-=(
PoissonSeries const& right) {
*this = *this - right;
return *this;
}

template<typename Value, int rdegree_,
template<typename, typename, int> class Evaluator>
PoissonSeries<Value, rdegree_, Evaluator> operator+(
@@ -377,15 +395,15 @@ std::ostream& operator<<(
if (!is_start_of_output) {
out << " + ";
}
out << polynomials.sin << " * Sin(" << quantities::DebugString(ω)
out <<"(" << polynomials.sin << ") * Sin(" << quantities::DebugString(ω)
<< " * (T - " << series.origin_ << "))";
is_start_of_output = false;
}
if (!polynomials.cos.is_zero()) {
if (!is_start_of_output) {
out << " + ";
}
out << polynomials.cos << " * Cos(" << quantities::DebugString(ω)
out << "(" << polynomials.cos << ") * Cos(" << quantities::DebugString(ω)
<< " * (T - " << series.origin_ << "))";
is_start_of_output = false;
}