From 83ff055d2ff3e8f78655139a082329a0ddcb14eb Mon Sep 17 00:00:00 2001 From: pleroy Date: Sun, 6 Oct 2024 23:24:49 +0200 Subject: [PATCH 1/2] Comment. --- numerics/double_precision.hpp | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/numerics/double_precision.hpp b/numerics/double_precision.hpp index 8a03a181f8..dadb8a4e7b 100644 --- a/numerics/double_precision.hpp +++ b/numerics/double_precision.hpp @@ -64,25 +64,29 @@ DoublePrecision> Scale(T const& scale, template DoublePrecision> TwoProduct(T const& a, U const& b); -// Returns the exact value of `a * b + c`. +// Returns the exact value of `a * b + c` if `|a * b|` is small compared to +// `|c|`. See [SZ05], section 2.1. template DoublePrecision> TwoProductAdd(T const& a, U const& b, Product const& c); -// Returns the exact value of `a * b - c`. +// Returns the exact value of `a * b - c` if `|a * b|` is small compared to +// `|c|`. See [SZ05], section 2.1. template DoublePrecision> TwoProductSubtract(T const& a, U const& b, Product const& c); -// Returns the exact value of `-a * b + c`. +// Returns the exact value of `-a * b + c` if `|a * b|` is small compared to +// `|c|`. See [SZ05], section 2.1. template DoublePrecision> TwoProductNegatedAdd(T const& a, U const& b, Product const& c); -// Returns the exact value of `-a * b - c`. +// Returns the exact value of `-a * b - c` if `|a * b|` is small compared to +// `|c|`. See [SZ05], section 2.1. template DoublePrecision> TwoProductNegatedSubtract(T const& a, From 8dad7da7bae372af2f0a485b43dd465736da0080 Mon Sep 17 00:00:00 2001 From: pleroy Date: Mon, 7 Oct 2024 08:50:29 +0200 Subject: [PATCH 2/2] Add support for negative arguments. --- functions/sin_cos_test.cpp | 9 ++++----- numerics/sin_cos.cpp | 21 ++++++++++++--------- 2 files changed, 16 insertions(+), 14 deletions(-) diff --git a/functions/sin_cos_test.cpp b/functions/sin_cos_test.cpp index 5c2b93f57c..659913577b 100644 --- a/functions/sin_cos_test.cpp +++ b/functions/sin_cos_test.cpp @@ -23,8 +23,7 @@ class SinCosTest : public ::testing::Test {}; TEST_F(SinCosTest, Random) { std::mt19937_64 random(42); - // TODO(phl): Negative angles. - std::uniform_real_distribution<> uniformly_at(0, π / 4); + std::uniform_real_distribution<> uniformly_at(-π / 4, π / 4); cpp_bin_float_50 max_sin_ulps_error = 0; cpp_bin_float_50 max_cos_ulps_error = 0; @@ -34,7 +33,7 @@ TEST_F(SinCosTest, Random) { #if _DEBUG static constexpr std::int64_t iterations = 100; #else - static constexpr std::int64_t iterations = 400'000; + static constexpr std::int64_t iterations = 300'000; #endif for (std::int64_t i = 0; i < iterations; ++i) { @@ -69,8 +68,8 @@ TEST_F(SinCosTest, Random) { } // This implementation is not quite correctly rounded, but not far from it. - EXPECT_LE(max_sin_ulps_error, 0.500002); - EXPECT_LE(max_cos_ulps_error, 0.500002); + EXPECT_LE(max_sin_ulps_error, 0.500003); + EXPECT_LE(max_cos_ulps_error, 0.500001); LOG(ERROR) << "Sin error: " << max_sin_ulps_error << std::setprecision(25) << " ulps for argument: " << worst_sin_argument diff --git a/numerics/sin_cos.cpp b/numerics/sin_cos.cpp index 232673d634..f4a2c8caee 100644 --- a/numerics/sin_cos.cpp +++ b/numerics/sin_cos.cpp @@ -52,38 +52,41 @@ Value CosPolynomial(Argument const x) { template FORCE_INLINE(inline) Value SinImplementation(Argument const x) { - if (x < sin_near_zero_cutoff) { + double const abs_x = std::abs(x); + double const sign = std::copysign(1, x); + if (abs_x < sin_near_zero_cutoff) { double const x² = x * x; double const x³ = x² * x; return x + x³ * SinPolynomialNearZero(x²); } else { - std::int64_t const i = std::lround(x * table_spacing_reciprocal); + std::int64_t const i = std::lround(abs_x * table_spacing_reciprocal); auto const& accurate_values = SinCosAccurateTable[i]; double const& x₀ = accurate_values.x; double const& sin_x₀ = accurate_values.sin_x; double const& cos_x₀ = accurate_values.cos_x; - double const h = x - x₀; + double const h = abs_x - x₀; DoublePrecision const sin_x₀_plus_h_cos_x₀ = TwoProductAdd(cos_x₀, h, sin_x₀); double const h² = h * h; double const h³ = h² * h; - return sin_x₀_plus_h_cos_x₀.value + - ((sin_x₀ * h² * CosPolynomial(h²) + - cos_x₀ * h³ * SinPolynomial(h²)) + - sin_x₀_plus_h_cos_x₀.error); + return sign * sin_x₀_plus_h_cos_x₀.value + + ((sign * sin_x₀ * h² * CosPolynomial(h²) + + sign * cos_x₀ * h³ * SinPolynomial(h²)) + + sign * sin_x₀_plus_h_cos_x₀.error); } } template FORCE_INLINE(inline) Value CosImplementation(Argument const x) { - std::int64_t const i = std::lround(x * table_spacing_reciprocal); + double const abs_x = std::abs(x); + std::int64_t const i = std::lround(abs_x * table_spacing_reciprocal); auto const& accurate_values = SinCosAccurateTable[i]; double const& x₀ = accurate_values.x; double const& sin_x₀ = accurate_values.sin_x; double const& cos_x₀ = accurate_values.cos_x; - double const h = x - x₀; + double const h = abs_x - x₀; DoublePrecision const cos_x₀_minus_h_sin_x₀ = TwoProductNegatedAdd(sin_x₀, h, cos_x₀);