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 ca2bf8c29d22aaf44caf31f1d7b6ba77ab418be5 (commit) from 894af6c4aaa3b83f3c13d2e59735c33f79abdc20 (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 ca2bf8c29d22aaf44caf31f1d7b6ba77ab418be5 Author: Raymond Toy toy.raymond@gmail.com Date: Tue Dec 24 12:57:45 2013 -0800

Fix ticket:90

src/code/irrat.lisp: src/code/irrat-dd.lisp:

* Remove the special case that made atanh continuous with quadrant III for x < -1 on the branch cut.

tests/trig.lisp:

* Update tests for atanh * Rename rel-or-abs-error to close-to.

diff --git a/src/code/irrat-dd.lisp b/src/code/irrat-dd.lisp index 024ec93..dbe9159 100644 --- a/src/code/irrat-dd.lisp +++ b/src/code/irrat-dd.lisp @@ -1663,63 +1663,59 @@ 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)) - (cond ((and (realp z) (< z -1)) - ;; ATANH is continuous with quadrant III in this case. - (dd-complex-atanh (complex z -0d0))) - (t - (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)))))))) + (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 ff23a5e..e4e744d 100644 --- a/src/code/irrat.lisp +++ b/src/code/irrat.lisp @@ -1741,62 +1741,59 @@ Z may be any number, but the result is always a complex." (when (typep z '(or double-double-float (complex double-double-float))) (return-from complex-atanh (dd-complex-atanh z)))

- (if (and (realp z) (< z -1)) - ;; atanh is continuous in quadrant III in this case. - (complex-atanh (complex z -0f0)) - (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))))) + (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 05437e5..df76231 100644 --- a/tests/trig.lisp +++ b/tests/trig.lisp @@ -215,18 +215,20 @@ (assert-eql nil (sincos-test (scale-float 1d0 1023) 1000)))

-;; Compute the relative error between actual and expected if expected -;; is not zero. Otherwise, return absolute error between actual and -;; expected. If the error is less than the threshold, return T. -;; Otherwise return the actual (relative or absolute) error. -(defun rel-or-abs-error (actual expected &optional (threshold double-float-epsilon)) +(defun close-to (actual expected &optional (threshold double-float-epsilon)) + "Determine if Actual is close to Expected. If Expected is not zero, + then close-to returns t if |Actual - Expected|/|Expected| <= + Threshold. If Expected is 0, then close-to returns T if |Actual - + Expected| <= threshold. In either of these conditions does not + hold, then a list of the actual error (relative or absolute), the + actual value and the expected value is returned." (let ((err (if (zerop expected) (abs (- actual expected)) (/ (abs (- actual expected)) (abs expected))))) (if (<= err threshold) t - err))) + (list err actual expected))))

;;; Tests for double-double-floats @@ -239,11 +241,11 @@ (define-test dd-sin.no-reduction "Test sin for small args without reduction" (:tag :sin :double-double) - (assert-eq t (rel-or-abs-error + (assert-eq t (close-to (sin .5w0) 4.794255386042030002732879352155713880818033679406006751886166131w-1 1w-32)) - (assert-eq t (rel-or-abs-error + (assert-eq t (close-to (sin -0.5w0) -4.794255386042030002732879352155713880818033679406006751886166131w-1 1w-32))) @@ -251,7 +253,7 @@ (define-test dd-sin.pi/2 "Test for arg near pi/2" (:tag :sin :double-double) - (assert-eq t (rel-or-abs-error + (assert-eq t (close-to (sin (/ kernel:dd-pi 2)) 1w0 1w-50))) @@ -265,27 +267,27 @@ "Test for sin with arg reduction" (:tag :sin :double-double) ;; Test for argument reduction with n mod 4 = 0 - (assert-eq t (rel-or-abs-error + (assert-eq t (close-to (sin (* 7/4 kernel:dd-pi)) -7.07106781186547524400844362104849691328261037289050238659653433w-1 0w0)) ;; Test for argument reduction with n mod 4 = 1 - (assert-eq t (rel-or-abs-error + (assert-eq t (close-to (sin (* 9/4 kernel:dd-pi)) 7.07106781186547524400844362104858161816423215627023442400880643w-1 0w0)) ;; Test for argument reduction with n mod 4 = 2 - (assert-eq t (rel-or-abs-error + (assert-eq t (close-to (sin (* 11/4 kernel:dd-pi)) 7.071067811865475244008443621048998682901731241858306822215522497w-1 8.716w-33)) ;; Test for argument reduction with n mod 4 = 3 - (assert-eq t (rel-or-abs-error + (assert-eq t (close-to (sin (* 13/4 kernel:dd-pi)) -7.071067811865475244008443621048777109664479707052746581685893187w-1 8.716w-33)) ;; Test for argument reduction, big value - (assert-eq t (rel-or-abs-error + (assert-eq t (close-to (sin (scale-float 1w0 120)) 3.778201093607520226555484700569229919605866976512306642257987199w-1 8.156w-33))) @@ -299,11 +301,11 @@ (define-test dd-cos.no-reduction "Test cos for small args without reduction" (:tag :cos :double-double) - (assert-eq t (rel-or-abs-error + (assert-eq t (close-to (cos .5w0) 8.775825618903727161162815826038296519916451971097440529976108683w-1 0w0)) - (assert-eq t (rel-or-abs-error + (assert-eq t (close-to (cos -0.5w0) 8.775825618903727161162815826038296519916451971097440529976108683w-1 0w0))) @@ -311,7 +313,7 @@ (define-test dd-cos.pi/2 "Test for arg near pi/2" (:tag :cos :double-double) - (assert-eq t (rel-or-abs-error + (assert-eq t (close-to (cos (/ kernel:dd-pi 2)) -1.497384904859169777320797133937725094986669701841027904483071358w-33 0w0))) @@ -320,27 +322,27 @@ "Test for cos with arg reduction" (:tag :cos :double-double) ;; Test for argument reduction with n mod 4 = 0 - (assert-eq t (rel-or-abs-error + (assert-eq t (close-to (cos (* 7/4 kernel:dd-pi)) 7.07106781186547524400844362104849691328261037289050238659653433w-1 0w0)) ;; Test for argument reduction with n mod 4 = 1 - (assert-eq t (rel-or-abs-error + (assert-eq t (close-to (cos (* 9/4 kernel:dd-pi)) 7.07106781186547524400844362104858161816423215627023442400880643w-1 3.487w-32)) ;; Test for argument reduction with n mod 4 = 2 - (assert-eq t (rel-or-abs-error + (assert-eq t (close-to (cos (* 11/4 kernel:dd-pi)) -7.071067811865475244008443621048998682901731241858306822215522497w-1 1.482w-31)) ;; Test for argument reduction with n mod 4 = 3 - (assert-eq t (rel-or-abs-error + (assert-eq t (close-to (cos (* 13/4 kernel:dd-pi)) -7.071067811865475244008443621048777109664479707052746581685893187w-1 7.845w-32)) ;; Test for argument reduction, big value - (assert-eq t (rel-or-abs-error + (assert-eq t (close-to (cos (scale-float 1w0 120)) -9.258790228548378673038617641074149467308332099286564602360493726w-1 0w0))) @@ -354,11 +356,11 @@ (define-test dd-tan.no-reduction "Test tan for small args without reduction" (:tag :tan :double-double) - (assert-eq t (rel-or-abs-error + (assert-eq t (close-to (tan .5w0) 5.463024898437905132551794657802853832975517201797912461640913859w-1 0w0)) - (assert-eq t (rel-or-abs-error + (assert-eq t (close-to (tan -0.5w0) -5.463024898437905132551794657802853832975517201797912461640913859w-1 0w0))) @@ -366,7 +368,7 @@ (define-test dd-tan.pi/2 "Test for arg near pi/2" (:tag :tan :double-double) - (assert-eq t (rel-or-abs-error + (assert-eq t (close-to (tan (/ kernel:dd-pi 2)) -6.67830961000672557834948096545679895621313886078988606234681001w32 0w0))) @@ -375,17 +377,17 @@ "Test for tan with arg reduction" (:tag :tan :double-double) ;; Test for argument reduction with n even - (assert-eq t (rel-or-abs-error + (assert-eq t (close-to (tan (* 7/4 kernel:dd-pi)) -1.000000000000000000000000000000001844257310064121018312678894979w0 3.422w-49)) ;; Test for argument reduction with n odd - (assert-eq t (rel-or-abs-error + (assert-eq t (close-to (tan (* 9/4 kernel:dd-pi)) 1.000000000000000000000000000000025802415787810837455445433037983w0 0w0)) ;; Test for argument reduction, big value - (assert-eq t (rel-or-abs-error + (assert-eq t (close-to (tan (scale-float 1w0 120)) -4.080663888418042385451434945255951177650840227682488471558860153w-1 1.888w-33))) @@ -753,7 +755,8 @@ ;; Thus, atanh(-2) is continuous with Quadrant II, NOT continuous with ;; Quadrant III! ;; -;; What do we do? +;; The formula, however, is clear. We will use the formula and ignore +;; the text in the CLHS. (defun atanh-def (z) (r*z 1/2 (- (log (1+z z)) @@ -761,19 +764,19 @@

(define-test branch-cut.atanh (:tag :atanh :branch-cuts) - ;; Test for x < -1, which is continuous with Quadrant III. Use the - ;; the value at #c(-2d0 -1d-20) as the reference. + ;; Test for x < -1, which is continuous with Quadrant II. Use the + ;; the value at #c(-2d0 +1d-20) as the reference. (multiple-value-bind (tr ti) - (get-signs (atanh-def #c(-2d0 -1d-20))) + (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))) + (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

-----------------------------------------------------------------------

Summary of changes: src/code/irrat-dd.lisp | 110 +++++++++++++++++++++++------------------------- src/code/irrat.lisp | 109 +++++++++++++++++++++++------------------------ tests/trig.lisp | 77 +++++++++++++++++---------------- 3 files changed, 146 insertions(+), 150 deletions(-)

hooks/post-receive