Difficulties in z.log for im(z) >> re (z)
Closed this issue · 8 comments
(From investigating josdejong/mathjs#2503)
Consider the following program using complex.js (v2.0.15):
import Complex from 'complex.js'
const reals = [0,-1,1]
const ims = [1e10, 1e15, 1e17, 1e18, 1e19, 1e30]
for (const r of reals) {
for (const im of ims) {
const c = new Complex(r, im)
console.log(r, im, ':', c.log())
}
}
It produces output:
0 10000000000 : { re: 23.025850929940457, im: 1.5707963267948966 }
0 1000000000000000 : { re: 34.538776394910684, im: 1.5707963267948966 }
0 100000000000000000 : { re: 39.14394658089878, im: 1.5707963267948966 }
0 1000000000000000000 : { re: 41.44653167389282, im: 1.5707963267948966 }
0 10000000000000000000 : { re: 43.74911676688687, im: 1.5707963267948966 }
0 1e+30 : { re: 69.07755278982137, im: 1.5707963267948966 }
-1 10000000000 : { re: 23.025849239078866, im: 1.5707963268948968 }
-1 1000000000000000 : { re: 34.490947945738256, im: 1.5707963267948977 }
-1 100000000000000000 : { re: 36.3662940453822, im: 1.5707963267948968 }
-1 1000000000000000000 : { re: 36.3662940453822, im: 1.5707963267948968 }
-1 10000000000000000000 : { re: NaN, im: 1.5707963267948966 }
-1 1e+30 : { re: NaN, im: 1.5707963267948966 }
1 10000000000 : { re: 23.025850234876927, im: 1.5707963266948965 }
1 1000000000000000 : { re: 34.59069013472681, im: 1.5707963267948957 }
1 100000000000000000 : { re: 37.33185619326892, im: 1.5707963267948966 }
1 1000000000000000000 : { re: 37.33185619326892, im: 1.5707963267948966 }
1 10000000000000000000 : { re: 37.33185619326892, im: 1.5707963267948966 }
1 1e+30 : { re: 37.33185619326892, im: 1.5707963267948966 }
So the re = 0 values look plausible, but the re = 1 values have the real part of the log bog down in a way that leaves them vastly off, and the degradation in the re = -1 values is even worse, eventually yielding NaN. Although of course with floats it's impossible to capture the miniscule differences from tau/4 in the theta of these various numbers, the complex module should be able to capture their magnitude and thereby return .log() values with reasonable real parts. This issue is a blocker for a contributor working on an implementation of lngamma for mathjs, so thank you for any attention you're able to give it.
Looking at the code, the real part of z.log() is computed by an internal function loghypot() which actually in its comments lists four methods of computation to avoid overflow. Looking over them, methods 3 and 4 rely on miniscule differences in the theta of the complex number to encode magnitude when im(z) is much larger than re(z), whereas methods 1 and 2 are symmetric in the magnitudes of re and im. I see that in a commit titled "The endless fight for precision" six years ago, the implementation switched from method 2 to method 3, even though it appears from the comments that the overall error was smallest with method 1. Would it be possible to switch again from method 3 to method 1 to take advantage of the latter's symmetry? Or if there's a reason why not, how about a regime check such that if b/a is larger than some threshold, it switches to method 1 or 2? I'd be happy to file a PR, but I'm reluctant without some guidance from the author as to the preferred approach, and if thresholding b/a seems the way to go, where to put that threshold? (I am definitely not a numerical analyst...)
Thanks for bringing this up!! I looked into it. The problem is the loghypot function, which I tried to make avoid squares, but obviously the trick doesn't play well when the difference between the real and imaginary part is huge. The idea was to define log(hypot(a, b)):= log(a / cos(atan2(b, a))), based on the trig composition identity cos(atan2(y, x)) = x / hypot(x, y).
When you have an idea to make it work with huge differences, would be glad to discuss.
I must have been funny back then with the commit message. Did you check your examples with method 2? I will do so as well now, I think method 1 has the problem with the squares, at least for 1e30.
Maybe the Python implementation can help a little
https://github.com/python/cpython/blob/1b2611eb0283055835e5df632a7a735db8c894b8/Modules/cmathmodule.c#L606
Interesting, for our case, they simply do log(hypot(ax/2., ay/2.)) + M_LN2
If method 1 has a problem for the 1e30 example, a reasonable plan seems to me to just check if abs(b)/abs(a) is bigger than (I am guessing at a reasonable threshold) 1e10, and use method 2 if so. Let me know if there's any work you'd like me to do toward a solution along those lines.
The idea of halving and squaring is pretty good, I checked your test cases against wolframalpha and they now have at least at the last digit a slight rounding error. Other than the python implementation, I avoided the square root and divided the log by two.
I also fixed this in quaternion.js. New releases will be published in a few minutes. Thanks @gwhitney !