rust-num / num-complex

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

powc of zero returns nan #114

Closed domna closed 1 year ago

domna commented 1 year ago

When I execute the following line

Complex64::from(0.).powc(Complex64::from(1.))

I get nan as the output.

Shouldn't this be zero instead of nan?

cuviper commented 1 year ago

Ouch, that's a big oversight. This probably comes from the implementation's r.ln() when r (the magnitude) is zero. We should probably add a special case for this... I'm not sure about how to handle the minutia of float +/-0 though. Maybe there's precedent we can draw from libraries in other languages.

domna commented 1 year ago

Yes I was also suspecting the ln as well. It seems that C does just a check for zero and uses a similar formula: https://github.com/eblot/newlib/blob/2a63fa0fd26ffb6603f69d9e369e944fe449c246/newlib/libm/complex/cpowf.c#L47

I'm currently using a workaround for this:

trait ComplexPower {
    fn pc(self, exp: Complex64) -> Complex64;
}

impl ComplexPower for Complex64 {
    fn pc(self, exp: Complex64) -> Complex64 {
        if self.re == 0. && self.im == 0. {
            return Complex64::from(0.);
        }
        self.powc(exp)
    }
}

which does exactly this and it seems to be fine for me (but should rather use the abs of the complex number, though). If I'm thinking correctly, I think for floats which are very small but not exactly zero it is fine as the ln will be still be properly defined but negative and very large (which we want here). Negative residuals will not occur as it uses the amplitude in the implementation. I could prepare a PR for this.