/[cmucl]/src/code/irrat-dd.lisp
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

  ViewVC Help
Powered by ViewVC 1.1.5