Bug in tanh? A solution, maybe for the GCL people?

>>>>> "Richard" == Richard Fateman <fateman at cs.berkeley.edu> writes:

    >> From the comments in float.lisp, maxima source, which I presume I wrote 32
    Richard> years ago.
    Richard> ...

    Richard> ;;
    Richard> ;;; The implementation for the special functions for a complex
    Richard> ;;; argument are mostly taken from W. Kahan, "Branch Cuts for Complex
    Richard> ;;; Elementary Functions or Much Ado About Nothing's Sign Bit", in
    Richard> ;;; Iserles and Powell (eds.) "The State of the Art in Numerical
    Richard> ;;; Analysis", pp 165-211, Clarendon Press, 1987

I think I wrote that (1987 is only 20 years ago)....

    Richard> ....

    Richard> We can compute  d(x):= exp(x) - 1, but do it carefully to preserve precision
    Richard> when
    Richard> |x| is small.  X is a FP number, and a FP number is returned.  

    Richard> Then tanh is   d(2*x)/(2+d(2*x))

Yes, this would work, if you have an accurate d(x).

Here is an alternative, based on Kahan's complex tanh, except modified
to work with just real numbers.  Simple tests near zero and for large
x indicates this this works.  Works on GCL too.

One nice advantage is that it only requires an accurate sinh and
sqrt.  GCL has an accurate sinh as does all other Lisp's I know of. 

I think with some minor modifications, this would be a decent
replacement for GCL's current tanh (which uses the simple formula
sinh(x)/cosh(x)).  Of course, GCL would need a different routine to
handle the complex-valued case.

Ray

(defun maxima-tanh (x)
  (declare (double-float x))
  (cond ((> (abs x)
	    #.(/ (+ (log 2.0d0)
		    (log most-positive-double-float)) 4d0))
	 (float-sign x))
	(t
	 (let* ((s (sinh x))
		(sqr (* s s))
		(rho (sqrt (+ 1d0 sqr)))
		(den (+ 1d0 sqr)))
	   (/ (* rho s)
	      den)))))