Complex.js
Complex.js copied to clipboard
real part of expm1 is inaccurate near nonzero multiples of 2*pi*i
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
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.