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

more accurate sqrt function #129

Open
wants to merge 1 commit into
base: master
Choose a base branch
from

Conversation

pascalkuthe
Copy link

I noticed that the square root implementation in num-complex uses conversion trough polar coordinates to compute complex squre roots. Usually the algorithm from https://dl.acm.org/doi/abs/10.1145/363717.363780 is used to compute the complex square root instead.

The algorithm uses hypot/norm and square root to compute csqrt. This approach should be faster since less transcendental calls are needed. Hypot and sqrt also tend to be faster compared to exp/ln/atan/cos/sin. Both hypot and sqart also have much higher precision. Most implementations guarantee that these two functions return the correctly rounded infinite accuracy result.

For prior art you can look at the glibc and musl implementation (https://git.musl-libc.org/cgit/musl/tree/src/complex/csqrt.c). The glibc implementation is a lot more complicated/hard to read because they ensure that underflow floating point exceptions are triggered correctly. I don't think that is something num-complex needs to do. There is also some accuracy loss for subnormal numbers that I didn't handle yet. I left a comment about it, but it is very minor.

Copy link
Member

@cuviper cuviper left a comment

Choose a reason for hiding this comment

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

Thanks! I didn't know about this algorithm.

// Complex64::new(2.4421097261308304e-162, 1.0115549693666347e-162)
// );

if self.re.is_zero() && self.im.is_zero() {
Copy link
Member

Choose a reason for hiding this comment

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

Can you add a source for all these special cases? e.g.
https://en.cppreference.com/w/c/numeric/complex/csqrt
(and make sure all those are covered)

Copy link
Author

Choose a reason for hiding this comment

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

I added more test to test_nan to make sure all of these are covered by theses and added a comment

src/lib.rs Outdated
}
if self.re.is_nan() {
// nan + nan i
return Self::new(self.re, (self.im - self.im) / (self.im - self.im));
Copy link
Member

Choose a reason for hiding this comment

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

We have a direct NaN -- not sure if this should also copysign though.

Suggested change
return Self::new(self.re, (self.im - self.im) / (self.im - self.im));
return Self::new(self.re, T::nan());

Copy link
Author

Choose a reason for hiding this comment

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

I wasn't 100% sure about the sign of nan here and 100% matched the code just to be sure but it seems T::nan() is indeed sufficient and the sign doesn't change.

// √(inf +/- x i) = inf +/- 0 i
// √(-inf +/- NaN i) = NaN +/- inf i
// √(-inf +/- x i) = 0 +/- inf i

Copy link
Member

Choose a reason for hiding this comment

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

Maybe add a variable to make this clearer:

Suggested change
#[allow(clippy::eq_op)]
let zero_or_nan = self.im - self.im;

Copy link
Author

Choose a reason for hiding this comment

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

that is indeed more readable, I also added a comments. good point

if scale {
self = self / four;
}
if self.re.is_sign_negative() {
Copy link
Member

Choose a reason for hiding this comment

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

We could also use a citation and link in a comment for the algorithm you mentioned.

Copy link
Author

Choose a reason for hiding this comment

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

I added a citation to the algorithm and the musl libc implementation as well as provide some additional background in a comement

@pascalkuthe
Copy link
Author

Thanks for the fast review! Apologies for not getting back to this yet.

I had to unexpectetly travel for work so I will not be able to get back to this before the weekend.

@pascalkuthe
Copy link
Author

This should be ready for another round of review, apologies again for the delay

Copy link

@ralphtandetzky ralphtandetzky left a comment

Choose a reason for hiding this comment

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

LGTM.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants