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.

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

A cube root. #1802

Merged
merged 13 commits into from
May 2, 2018
Merged

A cube root. #1802

merged 13 commits into from
May 2, 2018

Conversation

eggrobin
Copy link
Member

This is part of the effort on #1760.


constexpr std::uint64_t C = 0x2A9F7893782DA1CE;
static const __m128d sign_bit =
_mm_castsi128_pd(_mm_cvtsi64_si128(0x8000'0000'0000'0000));
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In other places we condition the use of intrinsics on some preprocessor symbol and use a slow path otherwise. Should we do the same in this function?

namespace principia {
namespace numerics {

constexpr std::uint64_t C = 0x2A9F7893782DA1CE;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No idea what this constant does and how it was computed. Could we have references to some literature in this code?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

double const denominator =
(7 * x³ + 42 * abs_y) * x⁶ + (30 * x³ + 2 * abs_y) * y²;
return _mm_cvtsd_f64(
_mm_or_pd(_mm_set_sd(x - numerator / denominator), sign));
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this faster than multiplying by sign?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, but we could gain a cycle by moving this or away from the critical path.


class CubeRootTest : public ::testing::Test {
protected:
static std::uint64_t bits(double const x) {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this be Bits?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done.


// Computes ∛y with a maximal error in [0.50005, 0.50022] ULPs; the result is
// incorrectly rounded for approximately 5 inputs per million.
double cbrt(double y);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not Cbrt?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done.

Eq(-std::numeric_limits<double>::infinity()));
EXPECT_THAT(cbrt(std::numeric_limits<double>::quiet_NaN()),
Truly(&std::isnan<double>));
EXPECT_THAT(cbrt(-std::numeric_limits<double>::quiet_NaN()),
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What about signalling NaNs?

}

TEST_F(CubeRootTest, ParticularlyBadRounding) {
EXPECT_THAT(cbrt(0x1.14E35E87EA5DFp0), Eq(0x1.06C80FCCA8E18p0));
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How do we know that it's bad? Could we log some indication of error?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done.

@pleroy
Copy link
Member

pleroy commented Apr 30, 2018

retest this please

@pleroy pleroy added the LGTM label May 2, 2018
@eggrobin eggrobin merged commit 483bc93 into mockingbirdnest:master May 2, 2018
@eggrobin eggrobin mentioned this pull request May 2, 2018
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants