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: 8135436f32d6
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: 4348aeeb8350
Choose a head ref
  • 8 commits
  • 5 files changed
  • 1 contributor

Commits on Aug 29, 2020

  1. Funky product of polynomials.

    # Conflicts:
    #	numerics/frequency_analysis_test.cpp
    pleroy committed Aug 29, 2020
    Copy the full SHA
    faff45f View commit details
  2. Copy the full SHA
    4d85188 View commit details
  3. Copy the full SHA
    c2db7c0 View commit details
  4. Remove the operators.

    pleroy committed Aug 29, 2020
    Copy the full SHA
    d3ba933 View commit details
  5. Copy the full SHA
    015c6b9 View commit details
  6. Formatting.

    pleroy committed Aug 29, 2020
    Copy the full SHA
    55cb9e3 View commit details
  7. Copy the full SHA
    34080a7 View commit details

Commits on Aug 30, 2020

  1. Merge pull request #2699 from pleroy/PointwiseInnerProduct

    Pointwise inner product of polynomials
    pleroy authored Aug 30, 2020
    Copy the full SHA
    4348aee View commit details
Showing with 228 additions and 23 deletions.
  1. +10 −0 geometry/cartesian_product.hpp
  2. +121 −23 geometry/cartesian_product_body.hpp
  3. +32 −0 numerics/polynomial.hpp
  4. +28 −0 numerics/polynomial_body.hpp
  5. +37 −0 numerics/polynomial_test.cpp
10 changes: 10 additions & 0 deletions geometry/cartesian_product.hpp
Original file line number Diff line number Diff line change
@@ -52,6 +52,16 @@ template<typename LTuple, typename RTuple,
constexpr auto operator*(LTuple const& left, RTuple const& right);

} // namespace polynomial_ring

namespace pointwise_inner_product {

template<typename LTuple, typename RTuple,
typename = std::enable_if_t<quantities::is_tuple_v<LTuple>>,
typename = std::enable_if_t<quantities::is_tuple_v<RTuple>>>
constexpr auto PointwiseInnerProduct(LTuple const& left, RTuple const& right);

} // namespace pointwise_inner_product

} // namespace geometry
} // namespace principia

144 changes: 121 additions & 23 deletions geometry/cartesian_product_body.hpp
Original file line number Diff line number Diff line change
@@ -6,6 +6,7 @@
#include <type_traits>

#include "base/macros.hpp"
#include "geometry/hilbert.hpp"
#include "quantities/named_quantities.hpp"
#include "quantities/tuples.hpp"

@@ -175,6 +176,12 @@ constexpr auto CartesianProductVectorSpace<
return {std::get<indices>(left) / right...};
}

} // namespace internal_cartesian_product

namespace internal_polynomial_ring {

using internal_cartesian_product::CartesianProductAdditiveGroup;

// A helper for prepending an element to a tuple. Somewhat similar to
// std::tuple_cat but conveniently exports the type of the result.
template<typename Element, typename Tuple,
@@ -210,7 +217,12 @@ Tail(Tuple const& tuple) -> Type {
return std::make_tuple(std::get<tail_indices>(tuple)...);
}

// Technically, this is only a ring if |CartesianProductMultiplicativeSpace| is
// the bona fide |CartesianProductVectorSpace|. If it is not, this implements a
// multiplication-like operation which has the expected properties with respect
// to the cartesian product additive group.
template<typename LTuple, typename RTuple,
template<typename, typename> class CartesianProductMultiplicativeSpace,
int lsize_ = std::tuple_size_v<LTuple>,
int rsize_ = std::tuple_size_v<RTuple>>
class PolynomialRing {
@@ -223,15 +235,16 @@ class PolynomialRing {
// insert a zero for the lowest degree (because of the valuation 1). This is
// the type of that zero.
using LHead = std::tuple_element_t<0, LTuple>;
using Zero = quantities::Product<LHead, RHead>;
using Zero = typename Hilbert<LHead, RHead>::InnerProductType;

// The hard part: generating the type of the result.
using LTupleRHeadProduct =
decltype(CartesianProductVectorSpace<RHead, LTuple>::Multiply(
decltype(CartesianProductMultiplicativeSpace<RHead, LTuple>::Multiply(
std::declval<LTuple>(),
std::declval<RHead>()));
using LTupleRTailProduct =
decltype(PolynomialRing<LTuple, RTail>::Multiply(
decltype(PolynomialRing<LTuple, RTail,
CartesianProductMultiplicativeSpace>::Multiply(
std::declval<LTuple>(),
std::declval<RTail>()));
using ZeroLTupleRTailProduct =
@@ -248,45 +261,115 @@ class PolynomialRing {
Result Multiply(LTuple const& left, RTuple const& right);
};

template<typename LTuple, typename RTuple, int lsize_>
class PolynomialRing<LTuple, RTuple, lsize_, 1> {
template<typename LTuple, typename RTuple,
template<typename, typename> class CartesianProductMultiplicativeSpace,
int lsize_>
class PolynomialRing<LTuple, RTuple,
CartesianProductMultiplicativeSpace,
lsize_, 1> {
using RHead = std::tuple_element_t<0, RTuple>;
using Result = decltype(CartesianProductVectorSpace<RHead, LTuple>::Multiply(
std::declval<LTuple>(),
std::declval<RHead>()));
using Result = decltype(
CartesianProductMultiplicativeSpace<RHead, LTuple>::Multiply(
std::declval<LTuple>(),
std::declval<RHead>()));

public:
FORCE_INLINE(static constexpr)
Result Multiply(LTuple const& left, RTuple const& right);
};

template<typename LTuple, typename RTuple, int lsize_, int rsize_>
constexpr auto PolynomialRing<LTuple, RTuple, lsize_, rsize_>::Multiply(
template<typename LTuple, typename RTuple,
template<typename, typename> class CartesianProductMultiplicativeSpace,
int lsize_, int rsize_>
constexpr auto PolynomialRing<LTuple, RTuple,
CartesianProductMultiplicativeSpace,
lsize_, rsize_>::Multiply(
LTuple const& left,
RTuple const& right) -> Result {
using cartesian_product::operator+;
using cartesian_product::operator*;
using polynomial_ring::operator*;

auto const right_head = std::get<0>(right);
auto const right_tail = TailGenerator<RTuple>::Tail(right);

return left * right_head +
ConsGenerator<Zero, LTupleRTailProduct>::Cons(
Zero{}, left * right_tail);
using RightHeadType = std::remove_const_t<decltype(right_head)>;
using RightTailType = std::remove_const_t<decltype(right_tail)>;

auto const left_times_right_head =
CartesianProductMultiplicativeSpace<RightHeadType, LTuple>::
Multiply(left, right_head);
auto const left_times_right_tail = PolynomialRing<
LTuple, RightTailType,
CartesianProductMultiplicativeSpace>::Multiply(left, right_tail);

return left_times_right_head +
ConsGenerator<Zero, LTupleRTailProduct>::Cons(Zero{},
left_times_right_tail);
}

template<typename LTuple, typename RTuple, int lsize_>
constexpr auto PolynomialRing<LTuple, RTuple, lsize_, 1>::Multiply(
template<typename LTuple, typename RTuple,
template<typename, typename> class CartesianProductMultiplicativeSpace,
int lsize_>
constexpr auto PolynomialRing<LTuple, RTuple,
CartesianProductMultiplicativeSpace,
lsize_, 1>::Multiply(
LTuple const& left,
RTuple const& right) -> Result {
using cartesian_product::operator*;

auto const right_head = std::get<0>(right);
return left * right_head;
return CartesianProductMultiplicativeSpace<
decltype(right_head), LTuple>::Multiply(left, right_head);
}

} // namespace internal_cartesian_product
} // namespace internal_polynomial_ring

namespace internal_pointwise_inner_product {

using geometry::Hilbert;
using quantities::Apply;

using internal_cartesian_product::CartesianProductAdditiveGroup;

template<typename Scalar, typename Tuple,
typename = std::make_index_sequence<std::tuple_size_v<Tuple>>>
class CartesianProductPointwiseMultiplicativeSpace;

template<typename Scalar, typename Tuple, std::size_t... indices>
class CartesianProductPointwiseMultiplicativeSpace<
Scalar, Tuple, std::index_sequence<indices...>> {
template<typename T>
using ScalarLeftProduct = typename Hilbert<Scalar, T>::InnerProductType;
template<typename T>
using ScalarRightProduct = typename Hilbert<T, Scalar>::InnerProductType;

public:
FORCE_INLINE(static constexpr) Apply<ScalarLeftProduct, Tuple> Multiply(
Scalar const& left,
Tuple const& right);
FORCE_INLINE(static constexpr) Apply<ScalarRightProduct, Tuple> Multiply(
Tuple const& left,
Scalar const& right);
};

template<typename Scalar, typename Tuple, std::size_t... indices>
constexpr auto CartesianProductPointwiseMultiplicativeSpace<
Scalar, Tuple,
std::index_sequence<indices...>>::Multiply(Scalar const& left,
Tuple const& right)
-> Apply<ScalarLeftProduct, Tuple> {
return {Hilbert<Scalar, std::tuple_element_t<indices, Tuple>>::InnerProduct(
left, std::get<indices>(right))...};
}

template<typename Scalar, typename Tuple, std::size_t... indices>
constexpr auto CartesianProductPointwiseMultiplicativeSpace<
Scalar, Tuple,
std::index_sequence<indices...>>::Multiply(Tuple const& left,
Scalar const& right)
-> Apply<ScalarRightProduct, Tuple> {
return {Hilbert<std::tuple_element_t<indices, Tuple>, Scalar>::InnerProduct(
std::get<indices>(left), right)...};
}

} // namespace internal_pointwise_inner_product

namespace cartesian_product {

@@ -343,11 +426,26 @@ namespace polynomial_ring {
template<typename LTuple, typename RTuple, typename, typename>
FORCE_INLINE(constexpr)
auto operator*(LTuple const& left, RTuple const& right) {
return internal_cartesian_product::
PolynomialRing<LTuple, RTuple>::Multiply(left, right);
return internal_polynomial_ring::PolynomialRing<
LTuple, RTuple,
internal_cartesian_product::CartesianProductVectorSpace>::Multiply(
left, right);
}

} // namespace polynomial_ring

namespace pointwise_inner_product {

template<typename LTuple, typename RTuple,
typename, typename>
constexpr auto PointwiseInnerProduct(LTuple const& left, RTuple const& right) {
return internal_polynomial_ring::PolynomialRing<
LTuple, RTuple,
internal_pointwise_inner_product::
CartesianProductPointwiseMultiplicativeSpace>::Multiply(left, right);
}

} // namespace pointwise_inner_product

} // namespace geometry
} // namespace principia
32 changes: 32 additions & 0 deletions numerics/polynomial.hpp
Original file line number Diff line number Diff line change
@@ -10,6 +10,7 @@

#include "base/not_null.hpp"
#include "base/traits.hpp"
#include "geometry/hilbert.hpp"
#include "geometry/point.hpp"
#include "quantities/named_quantities.hpp"
#include "quantities/tuples.hpp"
@@ -43,6 +44,7 @@ namespace internal_polynomial {
using base::is_instance_of_v;
using base::not_constructible;
using base::not_null;
using geometry::Hilbert;
using geometry::Point;
using quantities::Derivative;
using quantities::Derivatives;
@@ -177,6 +179,14 @@ class PolynomialInMonomialBasis : public Polynomial<Value, Argument> {
friend operator*(
PolynomialInMonomialBasis<L, A, 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> class E>
constexpr PolynomialInMonomialBasis<
typename Hilbert<L, R>::InnerProductType, A, l + r, E>
friend PointwiseInnerProduct(
PolynomialInMonomialBasis<L, A, l, E> const& left,
PolynomialInMonomialBasis<R, A, r, E> const& right);
template<typename V, typename A, int d,
template<typename, typename, int> class E>
friend std::ostream& operator<<(
@@ -290,6 +300,14 @@ class PolynomialInMonomialBasis<Value, Point<Argument>, degree_, Evaluator>
friend operator*(
PolynomialInMonomialBasis<L, A, 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> class E>
constexpr PolynomialInMonomialBasis<
typename Hilbert<L, R>::InnerProductType, A, l + r, E>
friend PointwiseInnerProduct(
PolynomialInMonomialBasis<L, A, l, E> const& left,
PolynomialInMonomialBasis<R, A, r, E> const& right);
template<typename V, typename A, int d,
template<typename, typename, int> class E>
friend std::ostream& operator<<(
@@ -375,6 +393,20 @@ operator*(
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,
typename Argument, int ldegree_, int rdegree_,
template<typename, typename, int> class Evaluator>
constexpr PolynomialInMonomialBasis<
typename Hilbert<LValue, RValue>::InnerProductType, Argument,
ldegree_ + rdegree_, Evaluator>
PointwiseInnerProduct(
PolynomialInMonomialBasis<LValue, Argument, ldegree_, Evaluator> const&
left,
PolynomialInMonomialBasis<RValue, Argument, rdegree_, Evaluator> const&
right);

// Output.

template<typename Value, typename Argument, int degree_,
28 changes: 28 additions & 0 deletions numerics/polynomial_body.hpp
Original file line number Diff line number Diff line change
@@ -30,6 +30,7 @@ using geometry::cartesian_product::operator+;
using geometry::cartesian_product::operator-;
using geometry::cartesian_product::operator*;
using geometry::cartesian_product::operator/;
using geometry::pointwise_inner_product::PointwiseInnerProduct;
using geometry::polynomial_ring::operator*;
using quantities::Apply;
using quantities::Exponentiation;
@@ -757,6 +758,33 @@ operator*(
}
}

template<typename LValue, typename RValue,
typename Argument, int ldegree_, int rdegree_,
template<typename, typename, int> class Evaluator>
FORCE_INLINE(constexpr)
PolynomialInMonomialBasis<
typename Hilbert<LValue, RValue>::InnerProductType, Argument,
ldegree_ + rdegree_, Evaluator>
PointwiseInnerProduct(
PolynomialInMonomialBasis<LValue, Argument, ldegree_, Evaluator> const&
left,
PolynomialInMonomialBasis<RValue, Argument, rdegree_, Evaluator> const&
right) {
if constexpr (is_instance_of_v<Point, Argument>) {
CONSTEXPR_CHECK(left.origin_ == right.origin_);
return PolynomialInMonomialBasis<
typename Hilbert<LValue, RValue>::InnerProductType, Argument,
ldegree_ + rdegree_, Evaluator>(
PointwiseInnerProduct(left.coefficients_, right.coefficients_),
left.origin_);
} else {
return PolynomialInMonomialBasis<
typename Hilbert<LValue, RValue>::InnerProductType, Argument,
ldegree_ + rdegree_, Evaluator>(
PointwiseInnerProduct(left.coefficients_, right.coefficients_));
}
}

template<typename Value, typename Argument, int degree_,
template<typename, typename, int> class Evaluator>
std::ostream& operator<<(
37 changes: 37 additions & 0 deletions numerics/polynomial_test.cpp
Original file line number Diff line number Diff line change
@@ -266,6 +266,43 @@ TEST_F(PolynomialTest, Ring) {
}
}

TEST_F(PolynomialTest, PointwiseInnerProduct) {
P2V::Coefficients const coefficients({
Displacement<World>({0 * Metre,
2 * Metre,
3 * Metre}),
Velocity<World>({-1 * Metre / Second,
1 * Metre / Second,
0 * Metre / Second}),
Vector<Acceleration, World>({1 * Metre / Second / Second,
1 * Metre / Second / Second,
-2 * Metre / Second / Second})});
P2V const p2va(coefficients_);
P2V const p2vb(coefficients);

auto const p = PointwiseInnerProduct(p2va, p2vb);
{
auto const actual = p.Evaluate(0 * Second);
EXPECT_THAT(actual, AlmostEquals(3 * Metre * Metre, 0));
}
{
auto const actual = p.Evaluate(1 * Second);
EXPECT_THAT(actual, AlmostEquals(5 * Metre * Metre, 0));
}
{
auto const actual = p.Evaluate(-1 * Second);
EXPECT_THAT(actual, AlmostEquals(1 * Metre * Metre, 0));
}
{
auto const actual = p.Evaluate(2 * Second);
EXPECT_THAT(actual, AlmostEquals(19 * Metre * Metre, 0));
}
{
auto const actual = p.Evaluate(-2 * Second);
EXPECT_THAT(actual, AlmostEquals(11 * Metre * Metre, 0));
}
}

TEST_F(PolynomialTest, AtOrigin) {
Instant const t0 = Instant() + 3 * Second;
P2A const p(coefficients_, t0);