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

Commits on Dec 28, 2020

  1. Reduction mod 2 pi.

    pleroy committed Dec 28, 2020
    Copy the full SHA
    8d39616 View commit details
  2. After egg's review.

    pleroy committed Dec 28, 2020
    Copy the full SHA
    c200c60 View commit details

Commits on Dec 29, 2020

  1. Merge pull request #2825 from pleroy/Mod2Pi

    Reduction mod 2π
    pleroy authored Dec 29, 2020
    Copy the full SHA
    8fde9a6 View commit details
Showing with 28 additions and 0 deletions.
  1. +4 −0 numerics/double_precision.hpp
  2. +10 −0 numerics/double_precision_body.hpp
  3. +14 −0 numerics/double_precision_test.cpp
4 changes: 4 additions & 0 deletions numerics/double_precision.hpp
Original file line number Diff line number Diff line change
@@ -12,6 +12,7 @@ namespace numerics {
namespace internal_double_precision {

using base::not_null;
using quantities::Angle;
using quantities::Difference;
using quantities::Product;
using quantities::Quotient;
@@ -80,6 +81,8 @@ DoublePrecision<Difference<T, U>> TwoDifference(T const& a, U const& b);
template<typename T, typename U, typename = Difference<Difference<T, U>, T>>
DoublePrecision<Difference<T, U>> TwoDifference(T const& a, U const& b);

DoublePrecision<Angle> Mod2π(DoublePrecision<Angle> const& θ);

template<typename T>
bool operator==(DoublePrecision<T> const& left,
DoublePrecision<T> const& right);
@@ -122,6 +125,7 @@ std::ostream& operator<<(std::ostream& os,
} // namespace internal_double_precision

using internal_double_precision::DoublePrecision;
using internal_double_precision::Mod2π;
using internal_double_precision::TwoDifference;
using internal_double_precision::TwoProduct;
using internal_double_precision::TwoSum;
10 changes: 10 additions & 0 deletions numerics/double_precision_body.hpp
Original file line number Diff line number Diff line change
@@ -21,6 +21,7 @@ using geometry::DoubleOrQuantityOrMultivectorSerializer;
using quantities::Abs;
using quantities::FusedMultiplyAdd;
using quantities::Quantity;
using quantities::si::Radian;
namespace si = quantities::si;

// Assumes that |T| and |U| have a memory representation that is a sequence of
@@ -204,6 +205,15 @@ DoublePrecision<Difference<T, U>> TwoDifference(T const& a, U const& b) {
return TwoSum(a, -b);
}

inline DoublePrecision<Angle> Mod2π(DoublePrecision<Angle> const& θ) {
static DoublePrecision<Angle> const two_π = []() {
return QuickTwoSum(0x1.921FB54442D18p1 * Radian,
0x1.1A62633145C07p-53 * Radian);
}();
auto const θ_over_2π = θ / two_π;
return θ - two_π * DoublePrecision<double>(static_cast<int>(θ_over_2π.value));
}

template<typename T>
bool operator==(DoublePrecision<T> const& left,
DoublePrecision<T> const& right) {
14 changes: 14 additions & 0 deletions numerics/double_precision_test.cpp
Original file line number Diff line number Diff line change
@@ -319,6 +319,20 @@ TEST_F(DoublePrecisionTest, LongQuotient) {
AlmostEquals(-7352815717686216.0 * std::pow(0.5, 110), 0));
}

TEST_F(DoublePrecisionTest, Mod2π) {
// Slightly above 2000 π.
DoublePrecision<Angle> a((2e3 * 103'993 + 2) / 33'102 * Radian);
auto const c = Mod2π(a);
// HexLiteral[
// CorrectlyRound[Mod[CorrectlyRound[(2000 103993 + 2)/33102], 2 π]]]
EXPECT_THAT(c.value + c.error,
AlmostEquals(0x1.F12375A5877D6p-15 * Radian, 0));
// HexLiteral[CorrectlyRound[(2000 103993 + 2)/33102] -
// CorrectlyRound[2000 CorrectlyRound[π]]]
EXPECT_THAT(a.value - 2000 * π * Radian,
AlmostEquals(0x1.F123760000000p-15 * Radian, 0));
}

} // namespace internal_double_precision
} // namespace numerics
} // namespace principia