rust-num / num-complex

Complex numbers for Rust
Apache License 2.0
231 stars 48 forks source link

more accurate sqrt function #129

Open pascalkuthe opened 4 months ago

pascalkuthe commented 4 months ago

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.

pascalkuthe commented 4 months ago

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 commented 3 months ago

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