# Diff of /src/code/irrat-dd.lisp

revision 1.11 by rtoy, Fri May 25 17:11:30 2007 UTC revision 1.12 by rtoy, Fri May 25 20:35:15 2007 UTC
# Line 1699  Z may be any number, but the result is a Line 1699  Z may be any number, but the result is a
1699  (defun dd-complex-atanh (z)  (defun dd-complex-atanh (z)
1700    "Compute atanh z = (log(1+z) - log(1-z))/2"    "Compute atanh z = (log(1+z) - log(1-z))/2"
1701    (declare (number z))    (declare (number z))
1702    (cond ((realp z)    (cond ((and (realp z) (< z -1))
1703           ;; Look at the definition:           ;; ATANH is continuous with quadrant III in this case.
1704           ;;           (dd-complex-atanh (complex z -0d0)))
;; atanh(z) = 1/2*(log(1+z)-log(1-z))
;;
(cond ((> z 1)
;; Let x = z, x > 1.  Then
;;
;;   atanh(x) = 1/2*(log(1+x)-log(1-x))
;;
;; Only the term log(1-x) requires care since the
;; other term is purely real.  The CLHS says atanh for
;; x > 1 is continuous with quadrant I.  Assume x is
;; really x0 + i*eps, where eps > 0.  Then
;;
;;   log(1-x) = log(x0-1) - i*pi/2
;;
;; because arg(1-x) = arg(1-x0-i*eps) = -pi
;;
;; Thus
;;
;;   atanh(x) = 1/2*log((x+1)/(x-1)) + i*pi/2
;;            = 1/2*log(1+2/(x-1)) + i*pi/2
(complex (* 0.5w0 (dd-%log1p (/ 2 (- z 1))))
dd-pi/2))
(t
;; As above, but z = -x, x > 1.  Then
;;
;;   atanh(z) = 1/2*(log(1-x)-log(1+x))
;;
;; And log(1-x) is the interesting term.  The CLHS
;; says in this case atanh is continuous with quadrant
;; III.  Let x = x0-i*eps.  Then
;;
;;   log(1-x) = log(x0-1) + i*pi/2
;;
;; because arg(1-x) = arg(1-x0-i*eps) = pi.  Thus
;;
;;   atanh(z) = 1/2*log((x-1)/(x+1)) - i*pi/2
;;            = -1/2*log((x+1)/(x-1)) - i*pi/2
(complex (* -0.5w0 (dd-%log1p (/ 2 (- (abs z) 1))))
(- dd-pi/2)))))
1705          (t          (t
1706           (flet ((careful-mul (a b)           (flet ((careful-mul (a b)
1707                    ;; Carefully multiply a and b, taking care to handle                    ;; Carefully multiply a and b, taking care to handle

Legend:
 Removed from v.1.11 changed lines Added in v.1.12