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

Commits on Dec 27, 2020

  1. Verified

    This commit was signed with the committer’s verified signature. The key has expired.
    LnL7 Daiderd Jordan
    Copy the full SHA
    31b4b42 View commit details
  2. Verified

    This commit was signed with the committer’s verified signature. The key has expired.
    LnL7 Daiderd Jordan
    Copy the full SHA
    f6176f4 View commit details
  3. Merge pull request #2822 from pleroy/DoublePrecision

    New double precision operators
    pleroy authored Dec 27, 2020

    Verified

    This commit was signed with the committer’s verified signature. The key has expired.
    LnL7 Daiderd Jordan
    Copy the full SHA
    c85c131 View commit details
Showing with 65 additions and 0 deletions.
  1. +9 −0 numerics/double_precision.hpp
  2. +25 −0 numerics/double_precision_body.hpp
  3. +31 −0 numerics/double_precision_test.cpp
9 changes: 9 additions & 0 deletions numerics/double_precision.hpp
Original file line number Diff line number Diff line change
@@ -14,6 +14,7 @@ namespace internal_double_precision {
using base::not_null;
using quantities::Difference;
using quantities::Product;
using quantities::Quotient;
using quantities::Sum;

// A simple container for accumulating a value using double precision. The
@@ -103,6 +104,14 @@ template<typename T, typename U>
DoublePrecision<Difference<T, U>> operator-(DoublePrecision<T> const& left,
DoublePrecision<U> const& right);

template<typename T, typename U>
DoublePrecision<Product<T, U>> operator*(DoublePrecision<T> const& left,
DoublePrecision<U> const& right);

template<typename T, typename U>
DoublePrecision<Quotient<T, U>> operator/(DoublePrecision<T> const& left,
DoublePrecision<U> const& right);

template<typename T>
std::string DebugString(DoublePrecision<T> const& double_precision);

25 changes: 25 additions & 0 deletions numerics/double_precision_body.hpp
Original file line number Diff line number Diff line change
@@ -255,6 +255,31 @@ DoublePrecision<Difference<T, U>> operator-(DoublePrecision<T> const& left,
return QuickTwoSum(sum.value, (sum.error + left.error) - right.error);
}

template<typename T, typename U>
DoublePrecision<Product<T, U>> operator*(DoublePrecision<T> const& left,
DoublePrecision<U> const& right) {
// Linnainmaa (1981), Software for Doubled-Precision Floating-Point
// Computations, algorithm longmul.
auto product = TwoProduct(left.value, right.value);
product.error +=
(left.value + left.error) * right.error + left.error * right.value;
return QuickTwoSum(product.value, product.error);
}

template<typename T, typename U>
DoublePrecision<Quotient<T, U>> operator/(DoublePrecision<T> const& left,
DoublePrecision<U> const& right) {
// Linnainmaa (1981), Software for Doubled-Precision Floating-Point
// Computations, algorithm longdiv.
auto const z = left.value / right.value;
auto const product = TwoProduct(right.value, z);
auto const zz =
((((left.value - product.value) - product.error) + left.error) -
z * right.error) /
(right.value + right.error);
return QuickTwoSum(z, zz);
}

template<typename T>
std::string DebugString(DoublePrecision<T> const& double_precision) {
// We use |DebugString| to get all digits when |T| is |double|. In that case
31 changes: 31 additions & 0 deletions numerics/double_precision_test.cpp
Original file line number Diff line number Diff line change
@@ -30,6 +30,7 @@ using quantities::Mass;
using quantities::Momentum;
using quantities::Speed;
using quantities::Sqrt;
using quantities::Square;
using quantities::Tan;
using quantities::si::Kilogram;
using quantities::si::Metre;
@@ -288,6 +289,36 @@ TEST_F(DoublePrecisionTest, Product) {
0));
}

TEST_F(DoublePrecisionTest, LongProduct) {
DoublePrecision<Length> a(3 * Metre);
a.Increment(474 * ε * Metre);
DoublePrecision<Length> b(7 * Metre);
b.Increment(-517 * ε * Metre);
DoublePrecision<Square<Length>> const c = a * b;
// The numbers below were obtained using Mathematica.
EXPECT_THAT(
c.value,
AlmostEquals(5910974510923886.0 * std::pow(0.5, 48) * Metre * Metre,
0));
EXPECT_THAT(
c.error,
AlmostEquals(7881299347837104.0 * std::pow(0.5, 102) * Metre * Metre,
0));
}

TEST_F(DoublePrecisionTest, LongQuotient) {
DoublePrecision<Length> a(3 * Metre);
a.Increment(474 * ε * Metre);
DoublePrecision<Length> b(7 * Metre);
b.Increment(-517 * ε * Metre);
DoublePrecision<double> const c = a / b;
// The numbers below were obtained using Mathematica.
EXPECT_THAT(c.value,
AlmostEquals(7720456504064105.0 * std::pow(0.5, 54), 0));
EXPECT_THAT(c.error,
AlmostEquals(-7352815717686216.0 * std::pow(0.5, 110), 0));
}

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