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

Commits on Aug 20, 2020

  1. Try to implement AtOrigin.

    pleroy committed Aug 20, 2020
    Copy the full SHA
    a654012 View commit details
  2. Copy the full SHA
    c8ea95f View commit details
  3. Something that compiles.

    pleroy committed Aug 20, 2020
    Copy the full SHA
    801481d View commit details

Commits on Aug 21, 2020

  1. A test that passes.

    pleroy committed Aug 21, 2020
    Copy the full SHA
    46a7327 View commit details

Commits on Aug 22, 2020

  1. First code cleanup.

    pleroy committed Aug 22, 2020
    Copy the full SHA
    27ecd23 View commit details
  2. Second part of the cleanup.

    pleroy committed Aug 22, 2020
    Copy the full SHA
    beaf4e9 View commit details
  3. Work around a Clang bug.

    pleroy committed Aug 22, 2020
    Copy the full SHA
    86056ff View commit details

Commits on Aug 23, 2020

  1. Merge pull request #2690 from pleroy/SetOrigin

    Support for changing the origin of a polynomial
    pleroy authored Aug 23, 2020
    Copy the full SHA
    069b36e View commit details
Showing with 119 additions and 1 deletion.
  1. +3 −0 numerics/polynomial.hpp
  2. +105 −1 numerics/polynomial_body.hpp
  3. +11 −0 numerics/polynomial_test.cpp
3 changes: 3 additions & 0 deletions numerics/polynomial.hpp
Original file line number Diff line number Diff line change
@@ -218,6 +218,9 @@ class PolynomialInMonomialBasis<Value, Point<Argument>, degree_, Evaluator>

Point<Argument> const& origin() const;

// Returns a copy of this polynomial adjusted to the given origin.
PolynomialInMonomialBasis AtOrigin(Point<Argument> const& origin) const;

template<int order = 1>
PolynomialInMonomialBasis<
Derivative<Value, Argument, order>, Point<Argument>, degree_ - order,
106 changes: 105 additions & 1 deletion numerics/polynomial_body.hpp
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
#pragma once
#pragma once

#include "numerics/polynomial.hpp"

@@ -31,6 +31,99 @@ using geometry::cartesian_product::operator*;
using geometry::cartesian_product::operator/;
using geometry::polynomial_ring::operator*;
using quantities::Apply;
using quantities::Exponentiation;
using quantities::Pow;
using quantities::Time;

// A helper for changing the origin of a monomial (x - x₁)ⁿ. It computes the
// coefficients of the same monomial as a polynomial of (x - x₂), i.e.:
// cₙ(x - x₁)ⁿ = cₙ((x - x₂) + (x₁ - x₂))ⁿ =
// Σ cₙ(n k)(x - x₂)ᵏ(x₁ - x₂)ⁿ⁻ᵏ
// where (n k) is the binomial coefficient. The coefficients are for a
// polynomial of the given degree, with zeros for the unneeded high-degree
// terms.
template<typename Value, typename Argument, int degree, int n,
template<typename, typename, int> class Evaluator,
typename = std::make_index_sequence<degree + 1>>
struct MonomialAtOrigin;

template<typename Value, typename Argument, int degree, int n,
template<typename, typename, int> class Evaluator,
std::size_t... k>
struct MonomialAtOrigin<Value, Argument, degree, n,
Evaluator,
std::index_sequence<k...>> {
using Coefficients =
typename
PolynomialInMonomialBasis<Value, Point<Argument>, degree, Evaluator>::
Coefficients;

// The parameter coefficient is the coefficient of the monomial. The
// parameter shift is x₁ - x₂, computed only once by the caller.
static Coefficients MakeCoefficients(
std::tuple_element_t<n, Coefficients> const& coefficient,
Argument const& shift);
};

template<typename Value, typename Argument, int degree, int n,
template<typename, typename, int> class Evaluator,
std::size_t... k>
auto MonomialAtOrigin<Value, Argument, degree, n,
Evaluator,
std::index_sequence<k...>>::MakeCoefficients(
std::tuple_element_t<n, Coefficients> const& coefficient,
Argument const& shift) -> Coefficients {
return {(k <= n ? coefficient * Binomial(n, k) *
Pow<static_cast<int>(n - k)>(shift)
: std::tuple_element_t<k, Coefficients>{})...};
}

// A helper for changing the origin of an entire polynomial, by repeatedly
// using MonomialAtOrigin. We need two helpers because changing the origin is
// a quadratic operation in terms of the degree.
template<typename Value, typename Argument, int degree_,
template<typename, typename, int> class Evaluator,
typename = std::make_index_sequence<degree_ + 1>>
struct PolynomialAtOrigin;

template<typename Value, typename Argument, int degree,
template<typename, typename, int> class Evaluator,
std::size_t... indices>
struct PolynomialAtOrigin<Value, Argument, degree, Evaluator,
std::index_sequence<indices...>> {
using Polynomial =
PolynomialInMonomialBasis<Value, Point<Argument>, degree, Evaluator>;

static Polynomial MakePolynomial(
typename Polynomial::Coefficients const& coefficients,
Point<Argument> const& from_origin,
Point<Argument> const& to_origin);
};

template<typename Value, typename Argument, int degree,
template<typename, typename, int> class Evaluator,
std::size_t ...indices>
auto PolynomialAtOrigin<Value, Argument, degree,
Evaluator,
std::index_sequence<indices...>>::
MakePolynomial(typename Polynomial::Coefficients const& coefficients,
Point<Argument> const& from_origin,
Point<Argument> const& to_origin) -> Polynomial {
Argument const shift = to_origin - from_origin;
std::array<typename Polynomial::Coefficients, degree + 1> const
all_coefficients{
MonomialAtOrigin<Value, Argument, degree, indices, Evaluator>::
MakeCoefficients(std::get<indices>(coefficients), shift)...};

// It would be nicer to compute the sum using a fold expression, but Clang
// refuses to find the operator + in that context. Fold expressions, the
// final frontier...
typename Polynomial::Coefficients sum_coefficients;
for (auto const& coefficients : all_coefficients) {
sum_coefficients = sum_coefficients + coefficients;
}
return Polynomial(sum_coefficients, to_origin);
}

// Index-by-index assignment of RTuple to LTuple, which must have at least as
// many elements (and the types must match).
@@ -416,6 +509,17 @@ origin() const {
return origin_;
}

template<typename Value, typename Argument, int degree_,
template<typename, typename, int> class Evaluator>
PolynomialInMonomialBasis<Value, Point<Argument>, degree_, Evaluator>
PolynomialInMonomialBasis<Value, Point<Argument>, degree_, Evaluator>::AtOrigin(
Point<Argument> const& origin) const {
return PolynomialAtOrigin<Value, Argument, degree_, Evaluator>::
MakePolynomial(coefficients_,
/*from_origin=*/origin_,
/*to_origin=*/origin);
}

template<typename Value, typename Argument, int degree_,
template<typename, typename, int> class Evaluator>
template<int order>
11 changes: 11 additions & 0 deletions numerics/polynomial_test.cpp
Original file line number Diff line number Diff line change
@@ -236,6 +236,17 @@ TEST_F(PolynomialTest, Ring) {
}
}

TEST_F(PolynomialTest, AtOrigin) {
Instant const t0 = Instant() + 3 * Second;
P2A const p(coefficients_, t0);
P2A const q = p.AtOrigin(Instant() - 2 * Second);
for (Instant t = Instant() - 10 * Second;
t < Instant() + 10 * Second;
t += 0.3 * Second) {
EXPECT_THAT(q.Evaluate(t), AlmostEquals(p.Evaluate(t), 0, 942));
}
}

TEST_F(PolynomialTest, Derivative) {
using P2 = PolynomialInMonomialBasis<Temperature, Time, 2, HornerEvaluator>;
using P3 = PolynomialInMonomialBasis<Current, Time, 3, HornerEvaluator>;