Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
donhatch opened this issue Feb 2, 2024 · 0 comments
Open

Comments

@donhatch
Copy link

donhatch commented Feb 2, 2024

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant