more accurate sqrt function #129
Conversation
| // Complex64::new(2.4421097261308304e-162, 1.0115549693666347e-162) | ||
| // ); | ||
|
|
||
| if self.re.is_zero() && self.im.is_zero() { |
There was a problem hiding this comment.
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)
There was a problem hiding this comment.
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)); |
There was a problem hiding this comment.
We have a direct NaN -- not sure if this should also copysign though.
| return Self::new(self.re, (self.im - self.im) / (self.im - self.im)); | |
| return Self::new(self.re, T::nan()); |
There was a problem hiding this comment.
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 | ||
|
|
There was a problem hiding this comment.
Maybe add a variable to make this clearer:
| #[allow(clippy::eq_op)] | |
| let zero_or_nan = self.im - self.im; |
There was a problem hiding this comment.
that is indeed more readable, I also added a comments. good point
| if scale { | ||
| self = self / four; | ||
| } | ||
| if self.re.is_sign_negative() { |
There was a problem hiding this comment.
We could also use a citation and link in a comment for the algorithm you mentioned.
There was a problem hiding this comment.
I added a citation to the algorithm and the musl libc implementation as well as provide some additional background in a comement
|
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. |
|
This should be ready for another round of review, apologies again for the delay |
|
@cuviper still quite interested in this. Anything I can do to push this over the finish line? |
|
How about tests that demonstrate the improved accuracy? |
|
@cuviper testing for accuracy is tricky. I used numpy (which generates the same results as glibc/musl) to generate some reference numbers for some cases that round differently with the current implementation (but the same with the new one). This is a bit of a chicken and egg problem with showing that the rounding is truely more accurate but since this is a well known algorithm (and very well known implementations) I think it's a fair approach. I choose some arbitrary examples to demonstrate accuracy. The current implementation rounds incorrectly for a ton of numbers so it's not practical to test exhaustively |
|
Does the accuracy show itself by comparing the result squared to the original value? You may be able to use simple inputs without having to hard-code expected results. Hopefully that round trip is better now than it was before, at least for some inputs. (and hopefully not any worse in general, but I agree we can't really be exhaustive about it) |
|
I did find a few where the roundtrip was better by an ULP or two. Some of the accuracy gets lost (particularly since it's multiple operations for complex) during the multiplication so its not perfect test (makes it a bit harder to find cases) but a good addition, thanks! During some trial and error I also never had any case where accuracy was worse. It was quite a while since I worked on this but I think I also semi exhaustively (somewhat sampled) tested against glibc at some point IIRC. With this algorithm I am pretty certain that existing implementations get this right so perfectly matching those is the test that gave me the most confidence in the implementation |
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.