infusion / Complex.js

Complex.js is a com numbers library written in JavaScript
https://raw.org/article/complex-numbers-in-javascript/
MIT License
232 stars 33 forks source link

real part of expm1 is inaccurate near nonzero multiples of 2*pi*i #45

Open donhatch opened 7 months ago

donhatch commented 7 months ago

The helper function cosm1(x) (where x is real) takes care to avoid catastrophic cancellation when x is near 0 (by switching to taylor series when abs(x)<=pi/4), which is good.

But it fails to do that when x is near any other (i.e. nonzero) multiple of 2*pi. As a consequence, the real part of Complex expm1(z) is accurate when z is near 0, but not when z is near a nonzero multiple of 2*pi*i, where catastrophic cancellation still occurs.

To demonstrate the problem, comparing naively-computed cos(x)-1 with the mathematically-equivalent-but-catastrophic-cancellation-free formula -2*sin(.5*x)**2 as suggested on my old page https://www.plunk.org/~hatch/rightway.html :

$ gnuplot
gnuplot> set samples 1001
gnuplot> plot [2*pi-4e-8:2*pi+4e-8] cos(x)-1, -2*sin(.5*x)**2 

image

One possible fix might be to start the implementation of cosm1 by doing argument reduction modulo 2*pi, to the range [-pi,pi]. But doing argument reduction correctly is tricky (you can't just mod out by 2*Math.PI, since that isn't an exact representation of 2*pi so it would introduce error) and it's not necessary.

I suggest simply using the robust formula -2*sin(.5*x)**2 in all cases-- no need for argument reduction, and you can get rid of the taylor series stuff too.