>>>>> "Robert" == Robert Dodier <robert.dodier at gmail.com> writes:
Robert> With GCL, Clisp, and SBCL (all Linux, and didn't test CMUCL),
Robert> rtest_trig reports that test #61 fails. In each case the erroneous result is
CMUCL fails in the same way.
Robert> [lambda([],cabs(bfloat(rectform(atanh(-2)))-atanh(bfloat(-2)))),
Robert> lambda([],cabs(bfloat(rectform(atanh(2)))-atanh(bfloat(2))))];
I am the culprit. I implemented more accurate atanh for bigfloats,
using Kahan's algorithm for atanh. I looked into this, and yes, this
is basically a branch cut problem. Kahan's algorithm matches Common
Lisp's (or vice versa) definition of atanh:
The branch cut for the inverse hyperbolic tangent function is
in two pieces: one along the negative real axis to the left of
-1 (inclusive), continuous with quadrant III, and one along the
positive real axis to the right of 1 (inclusive), continuous
with quadrant I. The points -1 and 1 are excluded from the
domain. The range is that strip of the complex plane containing
numbers whose imaginary part is between -<PI>/2 and <PI>/2. A
number with imaginary part equal to -<PI>/2 is in the range if
and only if its real part is strictly negative; a number with
imaginary part equal to <PI>/2 is in the range if and only if
its imaginary part is strictly positive. Thus the range of the
inverse hyperbolic tangent function is identical to that of the
inverse hyperbolic sine function with the points -<PI>i/2 and
<PI>i/2 excluded.
So atanh(bfloat(2)) returns 0.549 + %pi/2 * %i. But
rectform(atanh(2)) returns log(3)/2-%i*%pi/2, which is not continuous
with Quadrant I. atanh(-2) has the same problem, for the other branch
cut.
I do not have a fix for this. Alternatively, bigfloat atanh could be
changed to match, but then atanh(2.0) will produce different results
from atanh(2) and atanh(2b0).
Robert> With SBCL, rtest_trig reports tests #60 and #55 also fail.
Robert> #60 =>
Robert> [lambda([],cabs(float(rectform(acoth(1/2)))-acoth(float(1/2)))),
Robert> lambda([],cabs(float(rectform(atanh(2)))-atanh(float(2)))),
Robert> lambda([],cabs(float(rectform(asec(1/2)))-asec(float(1/2)))),
Robert> lambda([],cabs(float(rectform(acsc(1/2)))-acsc(float(1/2)))),
Robert> lambda([],cabs(float(rectform(asin(2)))-asin(float(2)))),
Robert> lambda([],cabs(float(rectform(acos(2)))-acos(float(2))))]
Robert> then map (lambda ([f], f()), %);
Robert> => [3.14..., 3.14..., 2.63..., 2.63..., 2.63..., 2.63...]
Hmm. What does sbcl say asin(2.0) is? CMUCL and Clisp both say
1.5707-1.316958*%i. Perhaps sbcl has a positive imaginary part?
Robert> Dunno what the 2.63's are about.
Robert> #55 =>
Robert> lambda([],
Robert> ['cabs(1/0.9-cos(asec(0.9))),
Robert> 'cabs(sqrt(0.9^2-1)/0.9-sin(asec(0.9))),
Robert> 'cabs(sqrt(0.9^2-1)-tan(asec(0.9))),
Robert> 'cabs(0.9/sqrt(0.9^2-1)-csc(asec(0.9))),
Robert> 'cabs(0.9-sec(asec(0.9))),
Robert> 'cabs(1/sqrt(0.9^2-1)-cot(asec(0.9)))])
Maxima says asec(0.9) is 0.46714530810326216*%i. Perhaps sbcl says
the sin of that is different from 0.48432210483785276*%i?
Ray