|Free42 complex ASIN/ASINH/ACOS/ACOSH redux|
Message #24 Posted by Thomas Okken on 18 Sept 2011, 9:10 p.m.,
in response to message #1 by Thomas Okken
To summarize the whole Free42 complex ASIN saga:
Earlier this year, 聲gel Martin reported to me that Free42 returned an incorrect result for arcsin(sin(200+200i)).
I thought I'd fixed that bug by changing the implementation of complex arcsin from -i*arcsinh(i*x) to conj(-i*arcsinh(conj(-i*x))).
I'm kicking myself now for ever making that change. The original algorithm and the revised one are equivalent, except for the edge case of |Re(x)|>1 and Im(x)=0, and the problem that 聲gel reported had nothing to do with that edge case. The 1.4.67 change didn't fix the problem with arcsin(sin(200+200i)); it merely moved it to a different part of the complex plane: instead of returning wrong results for 200+200i and -200+200i, it now returned wrong results for 200-200i and -200-200i.
The problem with arcsin(2.404+0i) that two people reported to me via email, and that also started this thread (or, rather, the original thread that got moved to the archive a couple of days ago!), was caused by my non-fix to the arcsin(sin(200+200i)) bug. Undoing that non-fix made the arcsin(2.404+0i) bug go away, but of course it didn't do anything about the arcsin(sin(200+200i)) bug, other than moving it back to its original spot.
The problem with arcsin(sin(200+200i)) turned out to be that my naive implementation of arcsinh(x)=ln(x+sqrt(x^2+1)) fails when |x| is large while Re(x)<0. In that case, the sqrt(x^2+1) part becomes -x+epsilon, and so x+sqrt(x^2+1) becomes noise.
Changing the arcsinh algorithm to -arcsinh(-x) when Re(x)<0 causes the sqrt(x^2+1) part to approach x instead of -x for large |x|, and this avoids the loss of precision in x+sqrt(x^2+1).
So, I fixed the arcsin(2.404+0i) bug I introduced in 1.4.67 by simply reverting back to the 1.4.66 code, and I fixed the arcsin(sin(200+200i)) bug by changing the arcsinh algorithm so as to avoid the loss of precision in x+sqrt(x^2+1) for Re(x)<0.
Note that the reason that Free42 1.4.(67,68,69) did return the correct result for arcsin(2.404) but not for arcsin(2.404+0i) is that the case for real X in ASIN is handled by completely different code than the case for complex X. The "real X" code was right all along.
Also note that there was a bug in ACOS and ACOSH, caused by the same loss of precision issue, in the arccosh code, in the x+sqrt(x^2-1) subexpression for large |x| and Re(x)<0. Nobody reported that bug, and it wasn't affected by the 1.4.67 changes, but that code should be OK now, too.
Thank you to Franz Huber and Miguel Toro-Alvarez for contacting me about the arcsin(2.404+0i) bug, and to 聲gel Martin for bringing the arcsin(sin(200+200i)) case to my attention, and to everyone who weighed in on the HP Forum!
Edited: 18 Sept 2011, 10:24 p.m.