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

Commits on Dec 18, 2020

  1. Copy the full SHA
    77a37eb View commit details
  2. Copy the full SHA
    7b7dd22 View commit details
  3. Fix all compilation errors.

    pleroy committed Dec 18, 2020
    Copy the full SHA
    6c53d63 View commit details
  4. A more complete test.

    pleroy committed Dec 18, 2020
    Copy the full SHA
    95b043a View commit details

Commits on Dec 19, 2020

  1. Test for the affine case.

    pleroy committed Dec 19, 2020
    Copy the full SHA
    69e9d80 View commit details
  2. Merge pull request #2821 from pleroy/PolynomialCompose

    Composition of polynomials
    pleroy authored Dec 19, 2020

    Verified

    This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
    Copy the full SHA
    2ed9ed8 View commit details
Showing with 151 additions and 0 deletions.
  1. +3 −0 geometry/cartesian_product.hpp
  2. +12 −0 geometry/cartesian_product_body.hpp
  3. +25 −0 numerics/polynomial.hpp
  4. +59 −0 numerics/polynomial_body.hpp
  5. +52 −0 numerics/polynomial_test.cpp
3 changes: 3 additions & 0 deletions geometry/cartesian_product.hpp
Original file line number Diff line number Diff line change
@@ -51,6 +51,9 @@ template<typename LTuple, typename RTuple,
typename = std::enable_if_t<quantities::is_tuple_v<RTuple>>>
constexpr auto operator*(LTuple const& left, RTuple const& right);

template<int exponent, typename Tuple>
constexpr auto Pow(Tuple const& tuple);

} // namespace polynomial_ring

namespace pointwise_inner_product {
12 changes: 12 additions & 0 deletions geometry/cartesian_product_body.hpp
Original file line number Diff line number Diff line change
@@ -450,6 +450,18 @@ auto operator*(LTuple const& left, RTuple const& right) {
left, right);
}

template<int exponent, typename Tuple>
constexpr auto Pow(Tuple const& tuple) {
static_assert(exponent > 0, "Cannot raise a tuple to the zero-th power");
if constexpr (exponent == 1) {
return tuple;
} else if constexpr (exponent % 2 == 0) {
return Pow<exponent / 2>(tuple * tuple);
} else {
return Pow<exponent / 2>(tuple * tuple) * tuple;
}
}

} // namespace polynomial_ring

namespace pointwise_inner_product {
25 changes: 25 additions & 0 deletions numerics/polynomial.hpp
Original file line number Diff line number Diff line change
@@ -189,6 +189,12 @@ class PolynomialInMonomialBasis : public Polynomial<Value_, Argument_> {
template<typename L, typename R, typename A,
int l, int r,
template<typename, typename, int> typename E>
constexpr PolynomialInMonomialBasis<L, A, l * r, E>
friend Compose(PolynomialInMonomialBasis<L, R, l, E> const& left,
PolynomialInMonomialBasis<R, A, r, E> const& right);
template<typename L, typename R, typename A,
int l, int r,
template<typename, typename, int> typename E>
constexpr PolynomialInMonomialBasis<
typename Hilbert<L, R>::InnerProductType, A, l + r, E>
friend PointwiseInnerProduct(
@@ -313,6 +319,12 @@ class PolynomialInMonomialBasis<Value_, Point<Argument_>, degree_, Evaluator>
template<typename L, typename R, typename A,
int l, int r,
template<typename, typename, int> typename E>
constexpr PolynomialInMonomialBasis<L, A, l * r, E>
friend Compose(PolynomialInMonomialBasis<L, R, l, E> const& left,
PolynomialInMonomialBasis<R, A, r, E> const& right);
template<typename L, typename R, typename A,
int l, int r,
template<typename, typename, int> typename E>
constexpr PolynomialInMonomialBasis<
typename Hilbert<L, R>::InnerProductType, A, l + r, E>
friend PointwiseInnerProduct(
@@ -405,6 +417,19 @@ operator*(
PolynomialInMonomialBasis<RValue, Argument, rdegree_, Evaluator> const&
right);

// Application monoid.

template<typename LValue, typename RValue,
typename Argument, int ldegree_, int rdegree_,
template<typename, typename, int> typename Evaluator>
constexpr PolynomialInMonomialBasis<LValue, Argument,
ldegree_ * rdegree_, Evaluator>
Compose(
PolynomialInMonomialBasis<LValue, RValue, ldegree_, Evaluator> const&
left,
PolynomialInMonomialBasis<RValue, Argument, rdegree_, Evaluator> const&
right);

// Returns a scalar polynomial obtained by pointwise inner product of two
// vector-valued polynomials.
template<typename LValue, typename RValue,
59 changes: 59 additions & 0 deletions numerics/polynomial_body.hpp
Original file line number Diff line number Diff line change
@@ -156,6 +156,34 @@ void TupleAssigner<LTuple, RTuple, std::index_sequence<indices...>>::Assign(
((std::get<indices>(left_tuple) = std::get<indices>(right_tuple)), ...);
}


// - 1 in the second type is ultimately to avoid evaluating Pow<0> as generating
// a one is hard.
template<typename LTuple, typename RTuple,
typename = std::make_index_sequence<std::tuple_size_v<LTuple> - 1>>
struct TupleComposition;

template<typename LTuple, typename RTuple, std::size_t... left_indices>
struct TupleComposition<LTuple, RTuple, std::index_sequence<left_indices...>> {
static constexpr auto Compose(LTuple const& left_tuple,
RTuple const& right_tuple);
};

template<typename LTuple, typename RTuple, std::size_t... left_indices>
constexpr auto
TupleComposition<LTuple, RTuple, std::index_sequence<left_indices...>>::Compose(
LTuple const& left_tuple,
RTuple const& right_tuple) {
// The + 1 in the expressions below match the - 1 in the primary declaration
// of TupleComposition.
return std::tuple(std::get<0>(left_tuple)) +
((std::get<left_indices + 1>(left_tuple) *
geometry::polynomial_ring::Pow<left_indices + 1>(right_tuple)) +
...);
}



template<typename Tuple, int order,
typename = std::make_index_sequence<std::tuple_size_v<Tuple> - order>>
struct TupleDerivation;
@@ -768,6 +796,37 @@ operator*(
}
}

template<typename LValue, typename RValue,
typename Argument, int ldegree_, int rdegree_,
template<typename, typename, int> typename Evaluator>
constexpr PolynomialInMonomialBasis<LValue, Argument,
ldegree_ * rdegree_, Evaluator>
Compose(
PolynomialInMonomialBasis<LValue, RValue, ldegree_, Evaluator> const&
left,
PolynomialInMonomialBasis<RValue, Argument, rdegree_, Evaluator> const&
right) {
using LCoefficients =
typename PolynomialInMonomialBasis<LValue, RValue, ldegree_,
Evaluator>::Coefficients;
using RCoefficients =
typename PolynomialInMonomialBasis<RValue, Argument, rdegree_,
Evaluator>::Coefficients;
if constexpr (is_instance_of_v<Point, Argument>) {
return PolynomialInMonomialBasis<LValue, Argument,
ldegree_ * rdegree_,
Evaluator>(
TupleComposition<LCoefficients, RCoefficients>::Compose(
left.coefficients_, right.coefficients_),
right.origin_);
} else {
return PolynomialInMonomialBasis<LValue, Argument,
ldegree_ * rdegree_, Evaluator>(
TupleComposition<LCoefficients, RCoefficients>::Compose(
left.coefficients_, right.coefficients_));
}
}

template<typename LValue, typename RValue,
typename Argument, int ldegree_, int rdegree_,
template<typename, typename, int> typename Evaluator>
52 changes: 52 additions & 0 deletions numerics/polynomial_test.cpp
Original file line number Diff line number Diff line change
@@ -266,6 +266,58 @@ TEST_F(PolynomialTest, Ring) {
}
}

TEST_F(PolynomialTest, Monoid) {
using P2A =
PolynomialInMonomialBasis<Temperature, Instant, 2, HornerEvaluator>;
using P2V =
PolynomialInMonomialBasis<Temperature, Time, 2, HornerEvaluator>;
using P3 =
PolynomialInMonomialBasis<Current, Temperature, 3, HornerEvaluator>;
Instant const t0;
P2A const p2a({1 * Kelvin,
3 * Kelvin / Second,
-8 * Kelvin / Second / Second}, t0);
P2V const p2v({1 * Kelvin,
3 * Kelvin / Second,
-8 * Kelvin / Second / Second});
P3 const p3({2 * Ampere,
-4 * Ampere / Kelvin,
3 * Ampere / Kelvin / Kelvin,
1 * Ampere / Kelvin / Kelvin / Kelvin});
auto const pa = Compose(p3, p2a);
auto const pv = Compose(p3, p2v);
{
auto const actual_a = pa(t0 + 0 * Second);
auto const actual_v = pv(0 * Second);
EXPECT_THAT(actual_a, AlmostEquals(2 * Ampere, 0));
EXPECT_THAT(actual_v, AlmostEquals(2 * Ampere, 0));
}
{
auto const actual_a = pa(t0 + 1 * Second);
auto const actual_v = pv(1 * Second);
EXPECT_THAT(actual_a, AlmostEquals(2 * Ampere, 0));
EXPECT_THAT(actual_v, AlmostEquals(2 * Ampere, 0));
}
{
auto const actual_a = pa(t0 - 1 * Second);
auto const actual_v = pv(-1 * Second);
EXPECT_THAT(actual_a, AlmostEquals(-658 * Ampere, 0));
EXPECT_THAT(actual_v, AlmostEquals(-658 * Ampere, 0));
}
{
auto const actual_a = pa(t0 + 2 * Second);
auto const actual_v = pv(2 * Second);
EXPECT_THAT(actual_a, AlmostEquals(-13648 * Ampere, 0));
EXPECT_THAT(actual_v, AlmostEquals(-13648 * Ampere, 0));
}
{
auto const actual_a = pa(t0 - 2 * Second);
auto const actual_v = pv(-2 * Second);
EXPECT_THAT(actual_a, AlmostEquals(-46396 * Ampere, 0));
EXPECT_THAT(actual_v, AlmostEquals(-46396 * Ampere, 0));
}
}

TEST_F(PolynomialTest, PointwiseInnerProduct) {
P2V::Coefficients const coefficients({
Displacement<World>({0 * Metre,