-
Notifications
You must be signed in to change notification settings - Fork 69
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 correctly-rounded cube root. #3062
Conversation
std::array<double, n> s{}; | ||
int k = 0; | ||
s[0] = f[0]; | ||
for (int j = 1; j < n; ++j) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Use a new-style loop.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
No, following [Nie04] and thus using indices. I guess we could loop for j in some range object, but if we want to go there we need to change a lot of loops.
} else { | ||
s[l + 1] = dʹ; | ||
} | ||
l = l - 1; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
--l
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
No, following [Nie04].
s[k] = c; | ||
if (d != 0) { | ||
int l = k - 1; | ||
k = k + 1; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
++k
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
No, following [Nie04].
auto const [cʹ, dʹ] = TwoSum(s[l], s[l + 1]); | ||
s[l] = cʹ; | ||
if (dʹ == 0) { | ||
k = k - 1; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
--k
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
No, following [Nie04].
static_assert((C - 296) * 0x1p-52 == (2 * 1023 - Γᴋ) / 3); | ||
std::uint64_t const Y = _mm_cvtsi128_si64(_mm_castpd_si128(Y_0)); | ||
std::uint64_t const Q = C + Y / 3; | ||
double const q = _mm_cvtsd_f64(_mm_castsi128_pd(_mm_cvtsi64_si128(Q))); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The code up to this point is identical in both methods, except for the constants. Could we factor it out by passing the constants as template parameters to an auxiliary function? It should be straightforward to inline.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
There isn’t that much code there ; a lot of it is comments on the constants, which differ.
Factoring things out would be very messy since we have a lot of local variables that we want to use throughout.
double const y² = y * y; | ||
double const q² = q * q; | ||
double const q³ = q² * q; | ||
double const d = |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I have no idea where these formulas come from.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Appendix D, it says so in the comment above.
@@ -38,6 +43,77 @@ class CubeRootTest : public ::testing::Test { | |||
: quiet_dead_beef_(FromBits(0x7FF8'0000'DEAD'BEEF)), | |||
signaling_dead_beef_(FromBits(0x7FF0'0000'DEAD'BEEF)) {} | |||
|
|||
struct RoundedReal { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
"Real" is not a word in C++. I'd prefer RoundedDouble
.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Real is a real word in the real world though. We round a real, that is, an element of ℝ (which in this case is the cube root of an element of binary64) to binary64 in any of three ways.
Rounding a double
makes little sense, it is in binary64 already.
Also corrected a typo in the documentation,
|
retest this please |
Towards #1760; follow-up on #2973.