Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add support for angles in [-π/4, π/4] #4110

Closed
wants to merge 3 commits into from
Closed

Add support for angles in [-π/4, π/4] #4110

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
83ff055
Comment.
pleroy Oct 6, 2024
8dad7da
Add support for negative arguments.
pleroy Oct 7, 2024
1f149db
Merge branch 'master' into SinCos
pleroy Oct 7, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 4 additions & 5 deletions functions/sin_cos_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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) {
Expand Down Expand Up @@ -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
Expand Down
12 changes: 8 additions & 4 deletions numerics/double_precision.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,25 +64,29 @@ DoublePrecision<Product<T, U>> Scale(T const& scale,
template<FMAPolicy fma_policy = FMAPolicy::Auto, typename T, typename U>
DoublePrecision<Product<T, U>> 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<FMAPolicy fma_policy = FMAPolicy::Auto, typename T, typename U>
DoublePrecision<Product<T, U>> TwoProductAdd(T const& a,
U const& b,
Product<T, U> 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<FMAPolicy fma_policy = FMAPolicy::Auto, typename T, typename U>
DoublePrecision<Product<T, U>> TwoProductSubtract(T const& a,
U const& b,
Product<T, U> 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<FMAPolicy fma_policy = FMAPolicy::Auto, typename T, typename U>
DoublePrecision<Product<T, U>> TwoProductNegatedAdd(T const& a,
U const& b,
Product<T, U> 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<FMAPolicy fma_policy = FMAPolicy::Auto, typename T, typename U>
DoublePrecision<Product<T, U>>
TwoProductNegatedSubtract(T const& a,
Expand Down
21 changes: 12 additions & 9 deletions numerics/sin_cos.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,38 +52,41 @@ Value CosPolynomial(Argument const x) {
template<FMAPolicy fma_policy>
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<fma_policy>(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<double> const sin_x₀_plus_h_cos_x₀ =
TwoProductAdd<fma_policy>(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<fma_policy>(h²) +
cos_x₀ * h³ * SinPolynomial<fma_policy>(h²)) +
sin_x₀_plus_h_cos_x₀.error);
return sign * sin_x₀_plus_h_cos_x₀.value +
((sign * sin_x₀ * h² * CosPolynomial<fma_policy>(h²) +
sign * cos_x₀ * h³ * SinPolynomial<fma_policy>(h²)) +
sign * sin_x₀_plus_h_cos_x₀.error);
}
}

template<FMAPolicy fma_policy>
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<double> const cos_x₀_minus_h_sin_x₀ =
TwoProductNegatedAdd<fma_policy>(sin_x₀, h, cos_x₀);
Expand Down