The value of eps in trigi.lisp is too small -- this gives poor numerical
accuracy for some trigonometric functions. An example
/* Homebrew acosh function. */
(C1) xacosh(x) := log(x + sqrt(x * x - 1))$
(C2) xacosh(4100.0d0);
(D2) 9.011889418380244
(C3) acosh(4100.0d0);
(D3) 9.011889433252344
(C4)
Each digit of xacosh(4100.0d0) is correct. [I choose 4100.0 because it
just exceeds 1 / sqrt(eps).]
We can fix this by changing eps to a value appropriate for double floats.
To do this, change (in trigi.lisp) to something like
(machine-mantissa-precision is now declared a constant in simp.lisp)
#.(SETQ EPS #+PDP10 (FSC 1.0 -26.)
#+cl ;(ASH 1.0 #+3600 -24. #-3600 -31.)
;(scale-float 1.0 -24)
(scale-float 1.0d0 (- machine-mantissa-precision))
#-(or PDP10 Cl) 1.4E-8)
alternatively, we could turn the numerical work over to Lisp. For example
(defun lisp-complex-to-maxima (z)
(cond ((complexp z)
`((mplus add) ,(realpart z) ((mtimes simp) $%i ,(imagpart z))))
(t
z)))
(defun acosh (x)
(lisp-complex-to-maxima (lisp:acosh (float x))))
The trig functions that aren't Lisp functions (csc, acsc, ...) would need
to
be handled differently.
I favor the policy of turning numerical work over to Lisp when ever
possible. Other than possible bugs in the
numerical functions in Lisp, what are arguments
against doing this?
Barton
P.S. gcl 2.4.0 (and 2.4.3) has a bug in its inverse hyperbolic cosine
function.
I'll report it.
[barton@localhost maxima]$ gcl
GCL (GNU Common Lisp) Version(2.4.0) Fri Jun 22 06:13:39 CDT 2001
Licensed under GNU Library General Public License
Contains Enhancements by W. Schelter
;; See A&S 4.6.12 page 87. The function id should evaluate to zero.
>(defun id (z)
(+ (acosh (* -1 z)) (acosh z) (* pi #c(0 -1))))
ID
;; For real values between -1 and 1, it seems that id evaluates to zero.
>(id -0.2)
#C(0.0 0.0)
>(id 0.2)
#C(0.0 0.0)
>(id -0.9)
#C(0.0 0.0)
>(id 0.9)
#C(0.0 0.0)
;; Moving into the complex plane, we've got trouble.
>(id #c(1 1))
#C(2.1225501238100715 -4.4740715185748234)
>(id -42)
#C(8.8613500906990961 0.0)