This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "CMU Common Lisp".
The branch, master has been updated via a00d8a9a746933ad02f39fb22149a11299ea6a9b (commit) from ca2bf8c29d22aaf44caf31f1d7b6ba77ab418be5 (commit)
Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below.
- Log ----------------------------------------------------------------- commit a00d8a9a746933ad02f39fb22149a11299ea6a9b Author: Raymond Toy toy.raymond@gmail.com Date: Fri Jan 3 10:13:12 2014 -0800
Produce correct results on branch cuts for atanh.
For (atanh x), where x is a real number on the branch cuts for atanh, make the code return the correct value based on the definition of atanh. This is inconsistent with the description in the CLHS because the values on the branch cut are now continuous with different quadrants. The formula is clear, however, for the values on the branch cut.
* src/code/irrat.lisp * src/code/irrat-dd.lisp * Use the correct branch cut values for atanh. * tests/trig.lisp * Update tests accordingly, with additional comments.
diff --git a/src/code/irrat-dd.lisp b/src/code/irrat-dd.lisp index dbe9159..29627b5 100644 --- a/src/code/irrat-dd.lisp +++ b/src/code/irrat-dd.lisp @@ -1663,59 +1663,62 @@ Z may be any number, but the result is always a complex." (defun dd-complex-atanh (z) _N"Compute atanh z = (log(1+z) - log(1-z))/2" (declare (number z)) - (let* ( ;; Constants - (theta (/ (sqrt most-positive-double-float) 4.0w0)) - (rho (/ 4.0w0 (sqrt most-positive-double-float))) - (half-pi dd-pi/2) - (rp (float (realpart z) 1.0w0)) - (beta (float-sign rp 1.0w0)) - (x (* beta rp)) - (y (* beta (- (float (imagpart z) 1.0w0)))) - (eta 0.0w0) - (nu 0.0w0)) - ;; Shouldn't need this declare. - (declare (double-double-float x y)) - (locally - (declare (optimize (speed 3) - (inhibit-warnings 3))) - (cond ((or (> x theta) - (> (abs y) theta)) - ;; To avoid overflow... - (setf nu (float-sign y half-pi)) - ;; eta is real part of 1/(x + iy). This is x/(x^2+y^2), - ;; which can cause overflow. Arrange this computation so - ;; that it won't overflow. - (setf eta (let* ((x-bigger (> x (abs y))) - (r (if x-bigger (/ y x) (/ x y))) - (d (+ 1.0d0 (* r r)))) - (if x-bigger - (/ (/ x) d) - (/ (/ r y) d))))) - ((= x 1.0w0) - ;; Should this be changed so that if y is zero, eta is set - ;; to +infinity instead of approx 176? In any case - ;; tanh(176) is 1.0d0 within working precision. - (let ((t1 (+ 4w0 (square y))) - (t2 (+ (abs y) rho))) - (setf eta (dd-%log (/ (sqrt (sqrt t1)) - (sqrt t2)))) - (setf nu (* 0.5d0 - (float-sign y - (+ half-pi (dd-%atan (* 0.5d0 t2)))))))) - (t - (let ((t1 (+ (abs y) rho))) - ;; Normal case using log1p(x) = log(1 + x) - (setf eta (* 0.25d0 - (dd-%log1p (/ (* 4.0d0 x) - (+ (square (- 1.0d0 x)) - (square t1)))))) - (setf nu (* 0.5d0 - (dd-%atan2 (* 2.0d0 y) - (- (* (- 1.0d0 x) - (+ 1.0d0 x)) - (square t1)))))))) - (complex (* beta eta) - (- (* beta nu)))))) + (if (realp z) + ;; See complex-atanh for why we do this. + (dd-complex-atanh (complex z (- (* 0.0 z)))) + (let* ( ;; Constants + (theta (/ (sqrt most-positive-double-float) 4.0w0)) + (rho (/ 4.0w0 (sqrt most-positive-double-float))) + (half-pi dd-pi/2) + (rp (float (realpart z) 1.0w0)) + (beta (float-sign rp 1.0w0)) + (x (* beta rp)) + (y (* beta (- (float (imagpart z) 1.0w0)))) + (eta 0.0w0) + (nu 0.0w0)) + ;; Shouldn't need this declare. + (declare (double-double-float x y)) + (locally + (declare (optimize (speed 3) + (inhibit-warnings 3))) + (cond ((or (> x theta) + (> (abs y) theta)) + ;; To avoid overflow... + (setf nu (float-sign y half-pi)) + ;; eta is real part of 1/(x + iy). This is x/(x^2+y^2), + ;; which can cause overflow. Arrange this computation so + ;; that it won't overflow. + (setf eta (let* ((x-bigger (> x (abs y))) + (r (if x-bigger (/ y x) (/ x y))) + (d (+ 1.0d0 (* r r)))) + (if x-bigger + (/ (/ x) d) + (/ (/ r y) d))))) + ((= x 1.0w0) + ;; Should this be changed so that if y is zero, eta is set + ;; to +infinity instead of approx 176? In any case + ;; tanh(176) is 1.0d0 within working precision. + (let ((t1 (+ 4w0 (square y))) + (t2 (+ (abs y) rho))) + (setf eta (dd-%log (/ (sqrt (sqrt t1)) + (sqrt t2)))) + (setf nu (* 0.5d0 + (float-sign y + (+ half-pi (dd-%atan (* 0.5d0 t2)))))))) + (t + (let ((t1 (+ (abs y) rho))) + ;; Normal case using log1p(x) = log(1 + x) + (setf eta (* 0.25d0 + (dd-%log1p (/ (* 4.0d0 x) + (+ (square (- 1.0d0 x)) + (square t1)))))) + (setf nu (* 0.5d0 + (dd-%atan2 (* 2.0d0 y) + (- (* (- 1.0d0 x) + (+ 1.0d0 x)) + (square t1)))))))) + (complex (* beta eta) + (- (* beta nu)))))))
(defun dd-complex-tanh (z) _N"Compute tanh z = sinh z / cosh z" diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp index e4e744d..bcd40c3 100644 --- a/src/code/irrat.lisp +++ b/src/code/irrat.lisp @@ -1740,60 +1740,78 @@ Z may be any number, but the result is always a complex." #+double-double (when (typep z '(or double-double-float (complex double-double-float))) (return-from complex-atanh (dd-complex-atanh z))) - - (let* ( ;; Constants - (theta (/ (sqrt most-positive-double-float) 4.0d0)) - (rho (/ 4.0d0 (sqrt most-positive-double-float))) - (half-pi (/ pi 2.0d0)) - (rp (float (realpart z) 1.0d0)) - (beta (float-sign rp 1.0d0)) - (x (* beta rp)) - (y (* beta (- (float (imagpart z) 1.0d0)))) - (eta 0.0d0) - (nu 0.0d0)) - ;; Shouldn't need this declare. - (declare (double-float x y)) - (locally - (declare (optimize (speed 3))) - (cond ((or (> x theta) - (> (abs y) theta)) - ;; To avoid overflow... - (setf nu (float-sign y half-pi)) - ;; eta is real part of 1/(x + iy). This is x/(x^2+y^2), - ;; which can cause overflow. Arrange this computation so - ;; that it won't overflow. - (setf eta (let* ((x-bigger (> x (abs y))) - (r (if x-bigger (/ y x) (/ x y))) - (d (+ 1.0d0 (* r r)))) - (if x-bigger - (/ (/ x) d) - (/ (/ r y) d))))) - ((= x 1.0d0) - ;; Should this be changed so that if y is zero, eta is set - ;; to +infinity instead of approx 176? In any case - ;; tanh(176) is 1.0d0 within working precision. - (let ((t1 (+ 4d0 (square y))) - (t2 (+ (abs y) rho))) - (setf eta (log (/ (sqrt (sqrt t1)) - (sqrt t2)))) - (setf nu (* 0.5d0 - (float-sign y - (+ half-pi (atan (* 0.5d0 t2)))))))) - (t - (let ((t1 (+ (abs y) rho))) - ;; Normal case using log1p(x) = log(1 + x) - (setf eta (* 0.25d0 - (%log1p (/ (* 4.0d0 x) - (+ (square (- 1.0d0 x)) - (square t1)))))) - (setf nu (* 0.5d0 - (atan (* 2.0d0 y) - (- (* (- 1.0d0 x) - (+ 1.0d0 x)) - (square t1)))))))) - (coerce-to-complex-type (* beta eta) - (- (* beta nu)) - z)))) + + ;; When z is purely real and x > 1, we have + ;; + ;; atanh(z) = 1/2*(log(1+z) - log(1-z)) + ;; = 1/2*(log(1+x) - (log(x-1)+i*pi)) + ;; = 1/2*(log(1+x) - log(x-1) - i*pi) + ;; = 1/2(log(1+x)-log(x-1)) - i*pi/2 + ;; + ;; Thus, the imaginary part is negative. Thus, for z > 1, atanh is + ;; continuous with quadrant IV. To preserve the identity atanh(-z) + ;; = -atanh(z), we must have atanh be continuous with Quadrant II + ;; for z < -1. (The identity is obviously true from the + ;; definition.) + ;; + ;; NOTE: this differs from what the CLHS says for the continuity. + ;; Instead of the text in the CLHS, we choose to use the definition + ;; to derive the correct values. + (if (realp z) + (complex-atanh (complex (float z) (- (* 0 (float z))))) + (let* ( ;; Constants + (theta (/ (sqrt most-positive-double-float) 4.0d0)) + (rho (/ 4.0d0 (sqrt most-positive-double-float))) + (half-pi (/ pi 2.0d0)) + (rp (float (realpart z) 1.0d0)) + (beta (float-sign rp 1.0d0)) + (x (* beta rp)) + (y (* beta (- (float (imagpart z) 1.0d0)))) + (eta 0.0d0) + (nu 0.0d0)) + ;; Shouldn't need this declare. + (declare (double-float x y)) + (locally + (declare (optimize (speed 3))) + (cond ((or (> x theta) + (> (abs y) theta)) + ;; To avoid overflow... + (setf nu (float-sign y half-pi)) + ;; eta is real part of 1/(x + iy). This is x/(x^2+y^2), + ;; which can cause overflow. Arrange this computation so + ;; that it won't overflow. + (setf eta (let* ((x-bigger (> x (abs y))) + (r (if x-bigger (/ y x) (/ x y))) + (d (+ 1.0d0 (* r r)))) + (if x-bigger + (/ (/ x) d) + (/ (/ r y) d))))) + ((= x 1.0d0) + ;; Should this be changed so that if y is zero, eta is set + ;; to +infinity instead of approx 176? In any case + ;; tanh(176) is 1.0d0 within working precision. + (let ((t1 (+ 4d0 (square y))) + (t2 (+ (abs y) rho))) + (setf eta (log (/ (sqrt (sqrt t1)) + (sqrt t2)))) + (setf nu (* 0.5d0 + (float-sign y + (+ half-pi (atan (* 0.5d0 t2)))))))) + (t + (let ((t1 (+ (abs y) rho))) + ;; Normal case using log1p(x) = log(1 + x) + (setf eta (* 0.25d0 + (%log1p (/ (* 4.0d0 x) + (+ (square (- 1.0d0 x)) + (square t1)))))) + (setf nu (* 0.5d0 + (atan (* 2.0d0 y) + (- (* (- 1.0d0 x) + (+ 1.0d0 x)) + (square t1)))))))) + (coerce-to-complex-type (* beta eta) + (- (* beta nu)) + z)))))
(defun complex-tanh (z) "Compute tanh z = sinh z / cosh z" diff --git a/tests/trig.lisp b/tests/trig.lisp index df76231..695c50d 100644 --- a/tests/trig.lisp +++ b/tests/trig.lisp @@ -770,24 +770,24 @@ (get-signs (atanh-def #c(-2d0 1d-20))) (assert-true (check-signs #'atanh -2d0 tr ti)) (assert-true (check-signs #'atanh -2w0 tr ti)) - (assert-true (check-signs #'atanh #c(-2d0 +0d0) tr ti)) - (assert-true (check-signs #'atanh #c(-2w0 +0w0) tr ti))) + (assert-true (check-signs #'atanh #c(-2d0 0d0) tr ti)) + (assert-true (check-signs #'atanh #c(-2w0 0w0) tr ti))) ;; Test the other side of the branch cut for x < -1. (multiple-value-bind (tr ti) (get-signs (atanh-def #c(-2d0 -1d-20))) (assert-true (check-signs #'atanh #c(-2d0 -0d0) tr ti)) (assert-true (check-signs #'atanh #c(-2w0 -0w0) tr ti)))
- ;; Test for x > 1, which is continuous with Quadrant I, using the - ;; value at #c(+2d0 1d-10) as the reference + ;; Test for x > 1, which is continuous with Quadrant IV, using the + ;; value at #c(+2d0 -1d-10) as the reference (multiple-value-bind (tr ti) - (get-signs (atanh-def #c(2d0 1d-20))) + (get-signs (atanh-def #c(2d0 -1d-10))) (assert-true (check-signs #'atanh 2d0 tr ti)) (assert-true (check-signs #'atanh 2w0 tr ti)) - (assert-true (check-signs #'atanh #c(2d0 0) tr ti)) - (assert-true (check-signs #'atanh #c(2w0 0) tr ti))) + (assert-true (check-signs #'atanh #c(2d0 -0d0) tr ti)) + (assert-true (check-signs #'atanh #c(2w0 -0w0) tr ti))) ;; Test the other side of the branch cut for x > 1. (multiple-value-bind (tr ti) - (get-signs (atanh-def #c(2d0 -1d-20))) - (assert-true (check-signs #'atanh #c(2d0 -0d0) tr ti)) - (assert-true (check-signs #'atanh #c(2w0 -0w0) tr ti)))) + (get-signs (atanh-def #c(2d0 +1d-20))) + (assert-true (check-signs #'atanh #c(2d0 0d0) tr ti)) + (assert-true (check-signs #'atanh #c(2w0 0w0) tr ti))))
-----------------------------------------------------------------------
Summary of changes: src/code/irrat-dd.lisp | 109 +++++++++++++++++++++-------------------- src/code/irrat.lisp | 126 +++++++++++++++++++++++++++--------------------- tests/trig.lisp | 20 ++++---- 3 files changed, 138 insertions(+), 117 deletions(-)
hooks/post-receive