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

Adapted complex exp() function so that it can handle inf and nan arguments as well #104

Merged
merged 9 commits into from Jun 17, 2022

Conversation

JorisDeRidder
Copy link
Contributor

Why this PR?

The current version of the complex exp() function is not able to handle arguments that contain +/- inf or NaN in their real or imaginary part. This impacts other complex functions that use the exp() function.

For example, for the most widely used implementation of the complex Faddeeva function w(), the current exp() implementation leads to w(1e160 - 1e159*i) = NaN + NaN *i , while the correct value is -5.586035480670854e-162 + 5.5860354806708545e-161 * i. The underlying reason is that the current exp() implementation erroneously returns exp(-inf + inf *i) = NaN + Nan *i instead of the correct 0 + 0*i.

Cf also issue #103.

Contents of this PR

  • I propose a modified complex exp() function that does deal with inf and nan. The added logic was strongly inspired by the one implemented in the <complex> C++ header that comes with clang++.
  • I added extra unit tests. The relevant values were taken from this page.
  • For the unit tests I also implemented close_naninf() and close_naninf_to_tol() as the existing functions close() and close_to_tol() are not able to deal with inf and nan.

src/lib.rs Outdated Show resolved Hide resolved
@cuviper
Copy link
Member

cuviper commented May 9, 2022

  • I added extra unit tests. The relevant values were taken from this page.

Could you please add that reference link as a comment in the test?

@JorisDeRidder
Copy link
Contributor Author

  • I added extra unit tests. The relevant values were taken from this page.

Could you please add that reference link as a comment in the test?

OK. Done.

@JorisDeRidder
Copy link
Contributor Author

A gentle reminder...

src/lib.rs Outdated
// formula: e^(a + bi) = e^a (cos(b) + i*sin(b)) = from_polar(e^a, b)

// Treat the corner cases +∞, -∞, and NaN
let mut 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.

I think this may be easier to follow if we unpack both parts, and then use them directly throughout.

Suggested change
let mut im = self.im;
let Complex { re, mut im } = self;

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Thanks for the suggestion, I just implemented, tested, and git-pushed it.

src/lib.rs Outdated
if self.re.is_infinite() {
if self.re < T::zero() {
if !self.im.is_finite() {
im = T::one();
Copy link
Member

Choose a reason for hiding this comment

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

Where does this 1.0 come from -- is that an arbitrary choice?

I think the result will be 0+0i, so why not just return that?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

It's indeed an arbitrary choice (following the clang C++ implementation). I agree that it would be simpler to just return 0+0i, so I just implemented, tested, and git-pushed that.

src/lib.rs Outdated
// Compare the imaginary parts
if a.im.is_finite() {
if b.im.is_finite() {
close = (a.im == b.im) || (a.im - b.im).abs() < tol;
Copy link
Member

Choose a reason for hiding this comment

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

We shouldn't lose the re comparison.

Suggested change
close = (a.im == b.im) || (a.im - b.im).abs() < tol;
close &= (a.im == b.im) || (a.im - b.im).abs() < tol;

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Aargh, sorry about that... None of the unit tests caught this. I just fixed the bug.

Copy link
Member

Choose a reason for hiding this comment

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

Yeah, "Who tests the tester?" It could just return true, and the unit tests would pass. :)

@cuviper
Copy link
Member

cuviper commented Jun 17, 2022

bors r+

@bors
Copy link
Contributor

bors bot commented Jun 17, 2022

@bors bors bot merged commit 96db82d into rust-num:master Jun 17, 2022
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.

None yet

2 participants