cmucl-cvs April 2015

cmucl-cvs@common-lisp.net

3 participants
318 discussions

[git] CMU Common Lisp branch master updated. snapshot-2013-12-a-44-ga00d8a9
by rtoy＠common-lisp.net 08 Apr '15

08 Apr '15

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(a)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 -- CMU Common Lisp

1 0

[git] CMU Common Lisp branch master updated. snapshot-2013-05-3-g7889e98
by rtoy＠common-lisp.net 08 Apr '15

08 Apr '15

1 0

[git] CMU Common Lisp branch master updated. snapshot-2013-05-9-gb3b0725
by rtoy＠common-lisp.net 08 Apr '15

08 Apr '15

1 0

[git] CMU Common Lisp branch master updated. snapshot-2014-06-81-g0a2163e
by rtoy＠common-lisp.net 08 Apr '15

08 Apr '15

1 0

[git] CMU Common Lisp annotated tag snapshot-2014-01 created. snapshot-2014-01
by rtoy＠common-lisp.net 08 Apr '15

08 Apr '15

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 annotated tag, snapshot-2014-01 has been created at d45720c2c896beada8be3dec56a2e1f95e1dd78d (tag) tagging 7bb20e4cfe7f12567dcbffa0a584f9d01fad5e64 (commit) replaces snapshot-2013-12-a tagged by Raymond Toy on Sat Jan 4 10:12:01 2014 -0800 - Log ----------------------------------------------------------------- Snapshot 2014-01 -----BEGIN PGP SIGNATURE----- Version: GnuPG v1.4.13 (Darwin) iEYEABECAAYFAlLITvgACgkQJ5IjUmgZO7IAtgCgyNzKZvbTBi+9onW1g89S/1cv SP4AoIPOqjKJXUOaGcNYJU/Wl+QnTkmI =bO3U -----END PGP SIGNATURE----- Raymond Toy (49): Optimizations to many vops to allow a descriptor for the second arg. Simplify dd-%%sin, dd-%%cos, and dd-%%tan. Add RT. Implement trig functions in Lisp Implement sincos using the new Lisp trig routines. This can now be Add test for sincos(-0d0). Small cleanups. Correct bug in tan. Need -ffloat-store when compiling e_rem/k_rem. Merge branch 'master' into rtoy-lisp-trig Convert to using lisp-unit. Remove RT. Add some comments and remove unneeded package qualifiers. Need to build e_rem_pio2.c on sparc now. Move ieee754_rem_pio2 to os-common.c. Remove the version in from Add more tests for exceptional values for cos and tan. Add a few Make all top-level disassembly functions accept base, case, and radix Regenerated. Merge branch 'master' into rtoy-simp-dd-trig Update for ppc to get pot files generated correctly. Merge branch 'master' into rtoy-simp-dd-trig For dd-%%sin, return x if x is small enough. (Makes sin(-0w0) be Add tests for double-double trig functions. Add dd-%sincos and use it as needed instead of calling sin and cos Add tests for dd-%sincos. Increase accuracy of tan near multiples of pi/4. Add tests for branch cuts. Document the value of atan better. The CLHS is confusing. Move src/tests/trig.lisp to tests/trig.lisp. Add some selected tests from trac bug reports. Fix trac.50 to use princ-to-string. Add more tests from trac. Add test file needed by a trac test. Use the formulas from the CLHS in the docstrings. More documentation. Fix typo. Fix test trac.69: Move the setting of the logical pathname Regenerated. Simple test script to run all of the tests in the tests directory. Float printer tests from reading the logs for print.lisp. Add tests for the branch cuts for asinh and acosh, and some fixes. Add tests for the branch cut for atanh. Not clear that this is Fix ticket:90 Produce correct results on branch cuts for atanh. Make interactive use of tests a little easier. Update according to commit logs. Regenerate. Don't run trac.15 on x87. We expect this to fail there and won't fix Fix typo in cis deftransform for double-doubles. ----------------------------------------------------------------------- hooks/post-receive -- CMU Common Lisp

1 0

[git] CMU Common Lisp branch master updated. snapshot-2013-05-8-g9be6530
by rtoy＠common-lisp.net 08 Apr '15

08 Apr '15

1 0

[git] CMU Common Lisp branch master updated. snapshot-2014-06-31-g3b8bd83
by rtoy＠common-lisp.net 08 Apr '15

08 Apr '15

1 0

[git] CMU Common Lisp branch master updated. snapshot-2014-01-5-g7b313ca
by rtoy＠common-lisp.net 08 Apr '15

08 Apr '15

1 0

[git] CMU Common Lisp branch master updated. snapshot-2014-03-3-g179e05d
by rtoy＠common-lisp.net 08 Apr '15

08 Apr '15

1 0

[git] CMU Common Lisp branch master updated. begin-x87-removal-17-g7fe70d3
by rtoy＠common-lisp.net 08 Apr '15

08 Apr '15

1 0

← Newer
1
...
14
15
16
17
18
19
20
...
32
Older →

Jump to page:

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

cmucl-cvs April 2015