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

revision 1.8 by rtoy, Wed May 23 13:16:33 2007 UTC revision 1.9 by rtoy, Wed May 23 16:48:50 2007 UTC
# Line 1680  Z may be any number, but the result is a Line 1680  Z may be any number, but the result is a
1680  (defun dd-complex-atanh (z)  (defun dd-complex-atanh (z)
1681    "Compute atanh z = (log(1+z) - log(1-z))/2"    "Compute atanh z = (log(1+z) - log(1-z))/2"
1682    (declare (number z))    (declare (number z))
1683    (if (and (realp z) (< z -1))    (cond ((realp z)
1684        ;; atanh is continuous in quadrant III in this case.           ;; Look at the definition:
1685        (dd-complex-atanh (complex z -0f0))           ;;
1686        (let* ( ;; Constants           ;; atanh(z) = 1/2*(log(1+z)-log(1-z))
1687               (theta (/ (sqrt most-positive-double-float) 4.0w0))           ;;
1688               (rho (/ 4.0w0 (sqrt most-positive-double-float)))           (cond ((> z 1)
1689               (half-pi dd-pi/2)                  ;; Let x = z, x > 1.  Then
1690               (rp (float (realpart z) 1.0w0))                  ;;
1691               (beta (float-sign rp 1.0w0))                  ;;   atanh(x) = 1/2*(log(1+x)-log(1-x))
1692               (x (* beta rp))                  ;;
1693               (y (* beta (- (float (imagpart z) 1.0w0))))                  ;; Only the term log(1-x) requires care since the
1694               (eta 0.0w0)                  ;; other term is purely real.  The CLHS says atanh for
1695               (nu 0.0w0))                  ;; x > 1 is continuous with quadrant I.  Assume x is
1696          ;; Shouldn't need this declare.                  ;; really x0 + i*eps, where eps > 0.  Then
1697          (declare (double-double-float x y))                  ;;
1698          (locally                  ;;   log(1-x) = log(x0-1) - i*pi/2
1699              (declare (optimize (speed 3)))                  ;;
1700            (cond ((or (> x theta)                  ;; because arg(1-x) = arg(1-x0-i*eps) = -pi
1701                       (> (abs y) theta))                  ;;
1702                   ;; To avoid overflow...                  ;; Thus
1703                   (setf nu (float-sign y half-pi))                  ;;
1704                   ;; eta is real part of 1/(x + iy).  This is x/(x^2+y^2),                  ;;   atanh(x) = 1/2*log((x+1)/(x-1)) + i*pi/2
1705                   ;; which can cause overflow.  Arrange this computation so                  ;;            = 1/2*log(1+2/(x-1)) + i*pi/2
1706                   ;; that it won't overflow.                  (complex (* 0.5w0 (dd-%log1p (/ 2 (- z 1))))
1707                   (setf eta (let* ((x-bigger (> x (abs y)))                           dd-pi/2))
1708                                    (r (if x-bigger (/ y x) (/ x y)))                 (t
1709                                    (d (+ 1.0d0 (* r r))))                  ;; As above, but z = -x, x > 1.  Then
1710                               (if x-bigger                  ;;
1711                                   (/ (/ x) d)                  ;;   atanh(z) = 1/2*(log(1-x)-log(1+x))
1712                                   (/ (/ r y) d)))))                  ;;
1713                  ((= x 1.0w0)                  ;; And log(1-x) is the interesting term.  The CLHS
1714                   ;; Should this be changed so that if y is zero, eta is set                  ;; says in this case atanh is continuous with quadrant
1715                   ;; to +infinity instead of approx 176?  In any case                  ;; III.  Let x = x0-i*eps.  Then
1716                   ;; tanh(176) is 1.0d0 within working precision.                  ;;
1717                   (let ((t1 (+ 4w0 (square y)))                  ;;   log(1-x) = log(x0-1) + i*pi/2
1718                         (t2 (+ (abs y) rho)))                  ;;
1719                     (setf eta (dd-%log (/ (sqrt (sqrt t1))                  ;; because arg(1-x) = arg(1-x0-i*eps) = pi.  Thus
1720                                           (sqrt t2))))                  ;;
1721                     (setf nu (* 0.5d0                  ;;   atanh(z) = 1/2*log((x-1)/(x+1)) - i*pi/2
1722                                 (float-sign y                  ;;            = -1/2*log((x+1)/(x-1)) - i*pi/2
1723                                             (+ half-pi (dd-%atan (* 0.5d0 t2))))))))                  (complex (* -0.5w0 (dd-%log1p (/ 2 (- (abs z) 1))))
1724                  (t                           (- dd-pi/2)))))
1725                   (let ((t1 (+ (abs y) rho)))          (t
1726                     ;; Normal case using log1p(x) = log(1 + x)           (let* ( ;; Constants
1727                     (setf eta (* 0.25d0                  (theta (/ (sqrt most-positive-double-float) 4.0w0))
1728                                  (dd-%log1p (/ (* 4.0d0 x)                  (rho (/ 4.0w0 (sqrt most-positive-double-float)))
1729                                                (+ (square (- 1.0d0 x))                  (half-pi dd-pi/2)
1730                                                   (square t1))))))                  (rp (float (realpart z) 1.0w0))
1731                     (setf nu (* 0.5d0                  (beta (float-sign rp 1.0w0))
1732                                 (dd-%atan2 (* 2.0d0 y)                  (x (* beta rp))
1733                                            (- (* (- 1.0d0 x)                  (y (* beta (- (float (imagpart z) 1.0w0))))
1734                                                  (+ 1.0d0 x))                  (eta 0.0w0)
1735                                               (square t1))))))))                  (nu 0.0w0))
1736            (complex (* beta eta)             ;; Shouldn't need this declare.
1737                     (- (* beta nu)))))))             (declare (double-double-float x y))
1738               (locally
1739                   (declare (optimize (speed 3)))
1740                 (cond ((or (> x theta)
1741                            (> (abs y) theta))
1742                        ;; To avoid overflow...
1743                        (setf nu (float-sign y half-pi))
1744                        ;; eta is real part of 1/(x + iy).  This is x/(x^2+y^2),
1745                        ;; which can cause overflow.  Arrange this computation so
1746                        ;; that it won't overflow.
1747                        (setf eta (let* ((x-bigger (> x (abs y)))
1748                                         (r (if x-bigger (/ y x) (/ x y)))
1749                                         (d (+ 1.0d0 (* r r))))
1750                                    (if x-bigger
1751                                        (/ (/ x) d)
1752                                        (/ (/ r y) d)))))
1753                       ((= x 1.0w0)
1754                        ;; Should this be changed so that if y is zero, eta is set
1755                        ;; to +infinity instead of approx 176?  In any case
1756                        ;; tanh(176) is 1.0d0 within working precision.
1757                        (let ((t1 (+ 4w0 (square y)))
1758                              (t2 (+ (abs y) rho)))
1759                          (setf eta (dd-%log (/ (sqrt (sqrt t1))
1760                                                (sqrt t2))))
1761                          (setf nu (* 0.5d0
1762                                      (float-sign y
1763                                                  (+ half-pi (dd-%atan (* 0.5d0 t2))))))))
1764                       (t
1765                        (let ((t1 (+ (abs y) rho)))
1766                          ;; Normal case using log1p(x) = log(1 + x)
1767                          (setf eta (* 0.25d0
1768                                       (dd-%log1p (/ (* 4.0d0 x)
1769                                                     (+ (square (- 1.0d0 x))
1770                                                        (square t1))))))
1771                          (setf nu (* 0.5d0
1772                                      (dd-%atan2 (* 2.0d0 y)
1773                                                 (- (* (- 1.0d0 x)
1774                                                       (+ 1.0d0 x))
1775                                                    (square t1))))))))
1776                 (complex (* beta eta)
1777                          (- (* beta nu))))))))
1778
1779  (defun dd-complex-tanh (z)  (defun dd-complex-tanh (z)
1780    "Compute tanh z = sinh z / cosh z"    "Compute tanh z = sinh z / cosh z"

Legend:
 Removed from v.1.8 changed lines Added in v.1.9