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: 83dd8c3d7798
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: 9fc4ea64c0b1
Choose a head ref
  • 6 commits
  • 3 files changed
  • 2 contributors

Commits on Jun 24, 2020

  1. Implement operator*.

    pleroy committed Jun 24, 2020
    Copy the full SHA
    97906c0 View commit details

Commits on Jun 27, 2020

  1. Copy the full SHA
    c67a481 View commit details
  2. Refactor the test.

    pleroy committed Jun 27, 2020
    Copy the full SHA
    6df6684 View commit details
  3. A test for the algebra.

    pleroy committed Jun 27, 2020
    Copy the full SHA
    0ec2d27 View commit details
  4. Lint.

    pleroy committed Jun 27, 2020
    Copy the full SHA
    5d13f7d View commit details
  5. Merge pull request #2616 from pleroy/PoissonAlgebra

    Algebra of Poisson series
    pleroy authored Jun 27, 2020
    Copy the full SHA
    9fc4ea6 View commit details
Showing with 161 additions and 57 deletions.
  1. +32 −16 numerics/poisson_series.hpp
  2. +68 −0 numerics/poisson_series_body.hpp
  3. +61 −41 numerics/poisson_series_test.cpp
48 changes: 32 additions & 16 deletions numerics/poisson_series.hpp
Original file line number Diff line number Diff line change
@@ -49,33 +49,40 @@ class PoissonSeries {
PolynomialsByAngularFrequency periodic_;

template<typename V, int r, template<typename, typename, int> class E>
PoissonSeries<V, r, E> friend operator-(PoissonSeries<V, r, E> const& right);
PoissonSeries<V, r, E>
friend operator-(PoissonSeries<V, r, E> const& right);
template<typename V, int l, int r, template<typename, typename, int> class E>
PoissonSeries<V, std::max(l, r), E> friend operator+(
PoissonSeries<V, l, E> const& left,
PoissonSeries<V, r, E> const& right);
PoissonSeries<V, std::max(l, r), E>
friend operator+(PoissonSeries<V, l, E> const& left,
PoissonSeries<V, r, E> const& right);
template<typename V, int l, int r, template<typename, typename, int> class E>
PoissonSeries<V, std::max(l, r), E> friend operator-(
PoissonSeries<V, l, E> const& left,
PoissonSeries<V, r, E> const& right);
PoissonSeries<V, std::max(l, r), E>
friend operator-(PoissonSeries<V, l, E> const& left,
PoissonSeries<V, r, E> const& right);
template<typename Scalar,
typename V, int d,
template<typename, typename, int> class E>
PoissonSeries<Product<Scalar, V>, d, E> friend operator*(
Scalar const& left,
PoissonSeries<V, d, E> const& right);
PoissonSeries<Product<Scalar, V>, d, E>
friend operator*(Scalar const& left,
PoissonSeries<V, d, E> const& right);
template<typename Scalar,
typename V, int d,
template<typename, typename, int> class E>
PoissonSeries<Product<V, Scalar>, d, E> friend operator*(
PoissonSeries<V, d, E> const& left,
Scalar const& right);
PoissonSeries<Product<V, Scalar>, d, E>
friend operator*(PoissonSeries<V, d, E> const& left,
Scalar const& right);
template<typename Scalar,
typename V, int d,
template<typename, typename, int> class E>
PoissonSeries<Quotient<V, Scalar>, d, E> friend operator/(
PoissonSeries<V, d, E> const& left,
Scalar const& right);
PoissonSeries<Quotient<V, Scalar>, d, E>
friend operator/(PoissonSeries<V, d, E> const& left,
Scalar const& right);
template<typename L, typename R,
int l, int r,
template<typename, typename, int> class E>
PoissonSeries<Product<L, R>, l + r, E>
friend operator*(PoissonSeries<L, l, E> const& left,
PoissonSeries<R, r, E> const& right);
};

// Vector space of Poisson series.
@@ -123,6 +130,15 @@ PoissonSeries<Quotient<Value, Scalar>, degree_, Evaluator>
operator/(PoissonSeries<Value, degree_, Evaluator> const& left,
Scalar const& right);

// Algebra of Poisson series.

template<typename LValue, typename RValue,
int ldegree_, int rdegree_,
template<typename, typename, int> class Evaluator>
PoissonSeries<Product<LValue, RValue>, ldegree_ + rdegree_, Evaluator>
operator*(PoissonSeries<LValue, ldegree_, Evaluator> const& left,
PoissonSeries<RValue, rdegree_, Evaluator> const& right);

} // namespace internal_poisson_series

using internal_poisson_series::PoissonSeries;
68 changes: 68 additions & 0 deletions numerics/poisson_series_body.hpp
Original file line number Diff line number Diff line change
@@ -5,6 +5,8 @@

#include <algorithm>
#include <functional>
#include <map>
#include <optional>

#include "quantities/elementary_functions.hpp"

@@ -221,6 +223,72 @@ operator/(PoissonSeries<Value, degree_, Evaluator> const& left,
return {aperiodic, std::move(periodic)};
}

template<typename LValue, typename RValue,
int ldegree_, int rdegree_,
template<typename, typename, int> class Evaluator>
PoissonSeries<Product<LValue, RValue>, ldegree_ + rdegree_, Evaluator>
operator*(PoissonSeries<LValue, ldegree_, Evaluator> const& left,
PoissonSeries<RValue, rdegree_, Evaluator> const& right) {
using Result =
PoissonSeries<Product<LValue, RValue>, ldegree_ + rdegree_, Evaluator>;

// Compute all the individual terms using elementary trigonometric identities
// and put them in a multimap, because the same frequency may appear multiple
// times.
std::multimap<AngularFrequency, Result::Polynomials> terms;
auto const aperiodic = left.aperiodic_ * right.aperiodic_;
for (auto const& [ω, polynomials] : left.periodic_) {
terms.emplace(
ω,
Result::Polynomials{/*sin=*/polynomials.sin * right.aperiodic_,
/*cos=*/polynomials.cos * right.aperiodic_});
}
for (auto const& [ω, polynomials] : right.periodic_) {
terms.emplace(
ω,
Result::Polynomials{/*sin=*/left.aperiodic_ * polynomials.sin,
/*cos=*/left.aperiodic_ * polynomials.cos});
}
for (auto const& [ωl, polynomials_left] : left.periodic_) {
for (auto const& [ωr, polynomials_right] : right.periodic_) {
terms.emplace(ωl - ωr,
Result::Polynomials{
/*sin=*/(-polynomials_left.cos * polynomials_right.sin +
polynomials_left.sin * polynomials_right.cos) /
2,
/*cos=*/(polynomials_left.sin * polynomials_right.sin +
polynomials_left.cos * polynomials_right.cos) /
2});
terms.emplace(ωl + ωr,
Result::Polynomials{
/*sin=*/(polynomials_left.cos * polynomials_right.sin +
polynomials_left.sin * polynomials_right.cos) /
2,
/*cos=*/(-polynomials_left.sin * polynomials_right.sin +
polynomials_left.cos * polynomials_right.cos) /
2});
}
}

// Now group the terms together by frequency.
Result::PolynomialsByAngularFrequency periodic;
std::optional<AngularFrequency> previous_ω;
for (auto it = terms.cbegin(); it != terms.cend(); ++it) {
auto const& ω = it->first;
auto const& polynomials = it->second;
if (previous_ω.has_value() && previous_ω.value() == ω) {
auto& previous_polynomials = periodic.rbegin()->second;
previous_polynomials.sin += polynomials.sin;
previous_polynomials.cos += polynomials.cos;
} else {
periodic.insert(*it);
}
previous_ω = ω;
}

return {aperiodic, std::move(periodic)};
}

} // namespace internal_poisson_series
} // namespace numerics
} // namespace principia
102 changes: 61 additions & 41 deletions numerics/poisson_series_test.cpp
Original file line number Diff line number Diff line change
@@ -1,6 +1,8 @@

#include "numerics/poisson_series.hpp"

#include <memory>

#include "gtest/gtest.h"
#include "numerics/polynomial_evaluators.hpp"
#include "quantities/elementary_functions.hpp"
@@ -21,89 +23,107 @@ using testing_utilities::AlmostEquals;
class PoissonSeriesTest : public ::testing::Test {
protected:
using Degree1 = PoissonSeries<double, 1, HornerEvaluator>;
};

TEST_F(PoissonSeriesTest, VectorSpace) {
Degree1::Polynomial pa0({0, 0 / Second});
Degree1::Polynomial psa0({100, 200 / Second});
Degree1::Polynomial pca0({1, 2 / Second});
Degree1::Polynomial pb0({3, 4 / Second});
PoissonSeriesTest() {
Degree1::Polynomial pa0({0, 0 / Second});
Degree1::Polynomial psa0({100, 200 / Second});
Degree1::Polynomial pca0({1, 2 / Second});
Degree1::Polynomial pb0({3, 4 / Second});

Degree1::Polynomial psa1({5, 6 / Second});
Degree1::Polynomial pca1({7, 8 / Second});
Degree1::Polynomial psb1({9, 10 / Second});
Degree1::Polynomial pcb1({11, 12 / Second});
Degree1::Polynomial psa1({5, 6 / Second});
Degree1::Polynomial pca1({7, 8 / Second});
Degree1::Polynomial psb1({9, 10 / Second});
Degree1::Polynomial pcb1({11, 12 / Second});

Degree1::Polynomial psa2({13, 14 / Second});
Degree1::Polynomial pca2({15, 16 / Second});
Degree1::Polynomial psa2({13, 14 / Second});
Degree1::Polynomial pca2({15, 16 / Second});

Degree1::Polynomial psb3({-17, -18 / Second});
Degree1::Polynomial pcb3({19, 20 / Second});
Degree1::Polynomial psb3({-17, -18 / Second});
Degree1::Polynomial pcb3({19, 20 / Second});

Degree1::Polynomials psca0{/*sin=*/psa0, /*cos=*/pca0};
Degree1::Polynomials psca0{/*sin=*/psa0, /*cos=*/pca0};

Degree1::Polynomials psca1{/*sin=*/psa1, /*cos=*/pca1};
Degree1::Polynomials pscb1{/*sin=*/psb1, /*cos=*/pcb1};
Degree1::Polynomials psca1{/*sin=*/psa1, /*cos=*/pca1};
Degree1::Polynomials pscb1{/*sin=*/psb1, /*cos=*/pcb1};

Degree1::Polynomials psca2{/*sin=*/psa2, /*cos=*/pca2};
Degree1::Polynomials psca2{/*sin=*/psa2, /*cos=*/pca2};

Degree1::Polynomials pscb3{/*sin=*/psb3, /*cos=*/pcb3};
Degree1::Polynomials pscb3{/*sin=*/psb3, /*cos=*/pcb3};

AngularFrequency ω0 = 0 * Radian / Second;
AngularFrequency ω1 = 1 * Radian / Second;
AngularFrequency ω2 = 2 * Radian / Second;
AngularFrequency ω3 = -3 * Radian / Second;
AngularFrequency ω0 = 0 * Radian / Second;
AngularFrequency ω1 = 1 * Radian / Second;
AngularFrequency ω2 = 2 * Radian / Second;
AngularFrequency ω3 = -3 * Radian / Second;

pa_ = std::make_unique<Degree1>(pa0,
Degree1::PolynomialsByAngularFrequency{
0, psca0}, {ω1, psca1}, {ω2, psca2}});
pb_ = std::make_unique<Degree1>(pb0,
Degree1::PolynomialsByAngularFrequency{
1, pscb1}, {ω3, pscb3}});
}

Degree1 pa(pa0, {{ω0, psca0}, {ω1, psca1}, {ω2, psca2}});
Degree1 pb(pb0, {{ω1, pscb1}, {ω3, pscb3}});
std::unique_ptr<Degree1> pa_;
std::unique_ptr<Degree1> pb_;
};

EXPECT_THAT(pa.Evaluate(1 * Second),
TEST_F(PoissonSeriesTest, Evaluate) {
EXPECT_THAT(pa_->Evaluate(1 * Second),
AlmostEquals(3 + 11 * Sin(1 * Radian) + 15 * Cos(1 * Radian) +
27 * Sin(2 * Radian) + 31 * Cos(2 * Radian),
0, 1));
EXPECT_THAT(pb.Evaluate(1 * Second),
EXPECT_THAT(pb_->Evaluate(1 * Second),
AlmostEquals(7 + 19 * Sin(1 * Radian) + 23 * Cos(1 * Radian) +
35 * Sin(3 * Radian) + 39 * Cos(3 * Radian),
32));
}

TEST_F(PoissonSeriesTest, VectorSpace) {
{
auto const identity = +pa;
auto const identity = +*pa_;
EXPECT_THAT(identity.Evaluate(1 * Second),
AlmostEquals(pa.Evaluate(1 * Second), 0));
AlmostEquals(pa_->Evaluate(1 * Second), 0));
}
{
auto const negated = -pb;
auto const negated = -*pb_;
EXPECT_THAT(negated.Evaluate(1 * Second),
AlmostEquals(-pb.Evaluate(1 * Second), 0));
AlmostEquals(-pb_->Evaluate(1 * Second), 0));
}
{
auto const sum = pa + pb;
auto const sum = *pa_ + *pb_;
EXPECT_THAT(
sum.Evaluate(1 * Second),
AlmostEquals(pa.Evaluate(1 * Second) + pb.Evaluate(1 * Second), 1));
AlmostEquals(pa_->Evaluate(1 * Second) + pb_->Evaluate(1 * Second), 1));
}
{
auto const difference = pa - pb;
auto const difference = *pa_ - *pb_;
EXPECT_THAT(
difference.Evaluate(1 * Second),
AlmostEquals(pa.Evaluate(1 * Second) - pb.Evaluate(1 * Second), 0));
AlmostEquals(pa_->Evaluate(1 * Second) - pb_->Evaluate(1 * Second), 0));
}
{
auto const left_product = 3 * pa;
auto const left_product = 3 * *pa_;
EXPECT_THAT(left_product.Evaluate(1 * Second),
AlmostEquals(3 * pa.Evaluate(1 * Second), 1));
AlmostEquals(3 * pa_->Evaluate(1 * Second), 1));
}
{
auto const right_product = pb * 4;
auto const right_product = *pb_ * 4;
EXPECT_THAT(right_product.Evaluate(1 * Second),
AlmostEquals(pb.Evaluate(1 * Second) * 4, 0));
AlmostEquals(pb_->Evaluate(1 * Second) * 4, 0));
}
{
auto const quotient = pb / 1.5;
auto const quotient = *pb_ / 1.5;
EXPECT_THAT(quotient.Evaluate(1 * Second),
AlmostEquals(pb.Evaluate(1 * Second) / 1.5, 0, 32));
AlmostEquals(pb_->Evaluate(1 * Second) / 1.5, 0, 32));
}
}

TEST_F(PoissonSeriesTest, Algebra) {
auto const product = *pa_ * *pb_;
EXPECT_THAT(
product.Evaluate(1 * Second),
AlmostEquals(pa_->Evaluate(1 * Second) * pb_->Evaluate(1 * Second), 6));
}

} // namespace numerics
} // namespace principia