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

Fixes nan value for powc of zero #116

Merged
merged 8 commits into from
Aug 13, 2023
Merged
Changes from 2 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
28 changes: 12 additions & 16 deletions src/lib.rs
Original file line number Diff line number Diff line change
Expand Up @@ -361,22 +361,8 @@ impl<T: Float> Complex<T> {
/// Raises `self` to a complex power.
#[inline]
pub fn powc(self, exp: Self) -> Self {
// formula: x^y = (a + i b)^(c + i d)
// = (ρ e^(i θ))^c (ρ e^(i θ))^(i d)
// where ρ=|x| and θ=arg(x)
// = ρ^c e^(−d θ) e^(i c θ) ρ^(i d)
// = p^c e^(−d θ) (cos(c θ)
// + i sin(c θ)) (cos(d ln(ρ)) + i sin(d ln(ρ)))
// = p^c e^(−d θ) (
// cos(c θ) cos(d ln(ρ)) − sin(c θ) sin(d ln(ρ))
// + i(cos(c θ) sin(d ln(ρ)) + sin(c θ) cos(d ln(ρ))))
// = p^c e^(−d θ) (cos(c θ + d ln(ρ)) + i sin(c θ + d ln(ρ)))
// = from_polar(p^c e^(−d θ), c θ + d ln(ρ))
let (r, theta) = self.to_polar();
Self::from_polar(
r.powf(exp.re) * (-exp.im * theta).exp(),
exp.re * theta + exp.im * r.ln(),
)
// formula: x^y = exp(y * ln(x))
(exp * self.ln()).exp()
}

/// Raises a floating point number to the complex power `self`.
Expand Down Expand Up @@ -1715,6 +1701,9 @@ pub(crate) mod test {

#[cfg(any(feature = "std", feature = "libm"))]
pub(crate) mod float {

use core::f64::INFINITY;

use super::*;
use num_traits::{Float, Pow};

Expand Down Expand Up @@ -1908,6 +1897,13 @@ pub(crate) mod test {
Complex::new(1.65826, -0.33502),
1e-5
));
let z = Complex::new(0.0, 0.0);
assert!(close(z.powc(b), z));
assert!(z.powc(Complex64::new(0., 0.)).is_nan());
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 result should be one -- in fact the standard f64::powf(x, 0.0) == 1.0 for any base, even NaN! I'm going to fix both powf and powc here...

Choose a reason for hiding this comment

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

That's not mathematically correct though, 0^0 is indeterminate.

Copy link
Member

Choose a reason for hiding this comment

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

See the prior discussions in rust-num/num-traits#78 and rust-num/num-traits#79. It's a pragmatic choice, especially to be consistent with std.

Copy link
Member

Choose a reason for hiding this comment

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

That wikipedia article also says:

Another example is the expression 0^0. Whether this expression is left undefined, or is defined to equal 1, depends on the field of application and may vary between authors. For more, see the article Zero to the power of zero.

assert!(z.powc(Complex64::new(0., INFINITY)).is_nan());
assert!(z.powc(Complex64::new(10., INFINITY)).is_nan());
assert!(z.powc(Complex64::new(INFINITY, INFINITY)).is_nan());
assert!(close(z.powc(Complex64::new(INFINITY, 0.)), z));
domna marked this conversation as resolved.
Show resolved Hide resolved
domna marked this conversation as resolved.
Show resolved Hide resolved
}

#[test]
Expand Down