(tanh #c(200d0 -200d0)) is incorrect
Describe the bug
(tanh #c(200d0 -200d0))
returns #c(1.0d0 -1d0)
. That's incorrect.
To Reproduce
Steps to reproduce the behavior:
- Evaluate
(tanh #c(200d0 -200d0))
- See the result is
#c(1d0 -1d0)
Expected behavior
The correct answer is closer to #c(1d0 -0d0)
. This can be seen by
using the simple definition of tanh(x) = sinh(x)/cosh(x)
. Using
this formula returns #c(1d0 0d0)
.
Desktop (please complete the following information):
- OS: All
- Version: [c275fc68] (2001-04-16) and later
Additional context
This is caused by this bit of code in complex-tanh
in
irrat.lisp:
(defun complex-tanh (z)
"Compute tanh z = sinh z / cosh z"
(declare (number z))
#+double-double
(when (typep z '(or double-double-float (complex double-double-float)))
(return-from complex-tanh (dd-complex-tanh z)))
(let ((x (float (realpart z) 1.0d0))
(y (float (imagpart z) 1.0d0)))
(locally
;; space 0 to get maybe-inline functions inlined
(declare (optimize (speed 3) (space 0)))
(cond ((> (abs x)
#-(or linux hpux) #.(/ (%asinh most-positive-double-float) 4d0)
;; This is more accurate under linux.
#+(or linux hpux) #.(/ (+ (%log 2.0d0)
(%log most-positive-double-float)) 4d0))
(coerce-to-complex-type (float-sign x)
(float-sign y) z))
The bug is in the sexp (float-sign y)
. It should be (float-sign y 0d0)
because we want to copy the sign of y
to 0. See Kahan's
Branch Cuts for Complex Elementary
Functions. It says
the answer should be copysign(1, ξ) + i*copysign(0, η)
where ξ + ηi
is the arg of tanh
.
We messed up the copysign
of the imaginary part. This is an ancient
bug!
While we're at it, there's a comment about using an alternative method
to compute (/ (asinh most-positive-double-float) 4)
. Now that we
use fdlibm to compute asinh
, these two expressions are the same on
Linux. Presumably it would be true on hpux as well.
Note: This was found in Maxima bug 4301.