Imprecision of complex natural logarithm for arguments close to 1

[Expander]

I found that in version 0.4.5 for complex numbers close to 1 the implementation of the complex logarithm ln looses many digits of precision:

use num::complex::Complex;

fn main() {
    let z = Complex::new(0.999995168714454, -0.004396919500211628);
    println!("z.ln() = {}", z.ln());
    // result:    0.0000048351532916926595 - 0.0043969124079359465i
    // expected result: 4.8351532916397e-6 - 0.0043969124079359465i
}

Link to rust playground: https://play.rust-lang.org/?version=stable&mode=debug&edition=2021&gist=d300b26da9fc5bed4304f52b99989386

To my knowledge there are algorithms to avoid such problems, for example by scaling the magnitude of the complex number accordingly.

[cuviper]

What is the basis for your expected answer? e.g. Wolfram Alpha says something even more different, although I'm sure that isn't accounting for the error in the floating point representation we're working with.

4.835153291547937069316406051483811016080566141758251053204... × 10^-6 -
0.004396912407935946425834708152509808769009396062847666986031... i

Still the current implementation here is pretty naive:

num-complex/src/lib.rs

Lines 239 to 243 in 0eb3e90

	pub fn ln(self) -> Self {
	// formula: ln(z) = ln|z| + i*arg(z)
	let (r, theta) = self.to_polar();
	Self::new(r.ln(), theta)
	}

Improvements are welcome!

[Expander]

My reference value is from Mathematica 13.3.0. I've written the input value z as a rational and evaluated Log[z] with up to 17 decimal digits:

z = 499997584357227045897/500000000000000000000 - (1099229875052907*I)/250000000000000000;
N[z, 17]   (* yields: 0.99999516871445409 - 0.00439691950021163 I *)
Log[z, 17] (* yields: 4.8351532916397 10^(-6) - 0.0043969124079359460 I *)

I am pretty confident that Mathematica's result is correct, because one obtains the same result with a Taylor expansion, evaluated in terms of rational numbers (increasing the order of the Taylor polynomial does not change the result):

Normal[Series[Log[1+y], {y,0,7}]] /. y -> (z-1) // N[#,17]&
(* yields: 4.8351532916397 10^(-6) - 0.0043969124079359460 I *)

I've cross-checked with julia, which even gives a different answer:

z = 0.999995168714454 - 0.004396919500211628im;
log(z) # yields: 4.835153291589251e-6 - 0.0043969124079359465im

[cuviper]

That ratio is not the input we have here, nor the 17-digit decimal approximation, but rather a pair of f64 -- which are the “binary64” type defined in IEEE 754-2008. Here's a tool for exploring those: real, imaginary. So it's more like this:

Log[9007155738389423 / 2^53 - (5069303045819176 / 2^60) I]

4.835153291589251684832824070075951847538991697458660122994... × 10^-6 -
0.004396912407935946439796111089250739704261084637130222808251... i

Julia's answer looks really good!

[Expander]

Thank a lot for this excellent comment! You are of course right! :)

Just for my understanding: My reference value for Log[z] was not correct, because the rational number that I used to approximate z did not quite match the decimal f64 number that I used as input in my original comment, right?

Anyway, I think it is possible to mimic the julia implementation (or any other appropriate one) to obtain a more precise value for the complex logarithm.

[cuviper]

Just for my understanding: My reference value for Log[z] was not correct, because the rational number that I used to approximate z did not quite match the decimal f64 number that I used as input in my original comment, right?

Even further, the decimal value that you wrote in your source does not match the binary f64 value that the compiler will actually use. The compiler picks the closest representation it can, but it's often not exact.