March 2011 - oct-cvs - mailman3.common-lisp.net

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 535d2509c2974c113e6c1bfa27ece5d5b8f99bd1
by Raymond Toy 07 Mar '11

07 Mar '11

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 "OCT: A portable Lisp implementation for quad-double precision floats". The branch, master has been updated via 535d2509c2974c113e6c1bfa27ece5d5b8f99bd1 (commit) from 258cd829289826458d6075caa279b253b8848749 (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 535d2509c2974c113e6c1bfa27ece5d5b8f99bd1 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Sun Mar 6 21:30:15 2011 -0500 Add support for Jacobi elliptic functions qd-elliptic.lisp: o Implementation for Jacobi sn, cn, and dn functions. qd-methods.lisp: o Add EPSILON method to return the floating-point epsilon value for the given float. oct.asd: o Build qd-elliptic.lisp. diff --git a/oct.asd b/oct.asd index 947b047..153706c 100644 --- a/oct.asd +++ b/oct.asd @@ -58,6 +58,8 @@ :depends-on ("qd-methods")) (:file "qd-complex" :depends-on ("qd-methods")) + (:file "qd-elliptic" + :depends-on ("qd-methods")) )) (defmethod perform ((op test-op) (c (eql (find-system :oct)))) diff --git a/qd-elliptic.lisp b/qd-elliptic.lisp new file mode 100644 index 0000000..9c56fcd --- /dev/null +++ b/qd-elliptic.lisp @@ -0,0 +1,180 @@ +;;;; -*- Mode: lisp -*- +;;;; +;;;; Copyright (c) 2011 Raymond Toy +;;;; Permission is hereby granted, free of charge, to any person +;;;; obtaining a copy of this software and associated documentation +;;;; files (the "Software"), to deal in the Software without +;;;; restriction, including without limitation the rights to use, +;;;; copy, modify, merge, publish, distribute, sublicense, and/or sell +;;;; copies of the Software, and to permit persons to whom the +;;;; Software is furnished to do so, subject to the following +;;;; conditions: +;;;; +;;;; The above copyright notice and this permission notice shall be +;;;; included in all copies or substantial portions of the Software. +;;;; +;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +;;;; OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +;;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +;;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +;;;; OTHER DEALINGS IN THE SOFTWARE. + +(in-package #:oct) + +(declaim (inline descending-transform ascending-transform)) + +(defun ascending-transform (u m) + ;; A&S 16.14.1 + ;; + ;; Take care in computing this transform. For the case where + ;; m is complex, we should compute sqrt(mu1) first as + ;; (1-sqrt(m))/(1+sqrt(m)), and then square this to get mu1. + ;; If not, we may choose the wrong branch when computing + ;; sqrt(mu1). + (let* ((root-m (sqrt m)) + (mu (/ (* 4 root-m) + (expt (1+ root-m) 2))) + (root-mu1 (/ (- 1 root-m) (+ 1 root-m))) + (v (/ u (1+ root-mu1)))) + (values v mu root-mu1))) + +(defun descending-transform (u m) + ;; Note: Don't calculate mu first, as given in 16.12.1. We + ;; should calculate sqrt(mu) = (1-sqrt(m1)/(1+sqrt(m1)), and + ;; then compute mu = sqrt(mu)^2. If we calculate mu first, + ;; sqrt(mu) loses information when m or m1 is complex. + (let* ((root-m1 (sqrt (- 1 m))) + (root-mu (/ (- 1 root-m1) (+ 1 root-m1))) + (mu (* root-mu root-mu)) + (v (/ u (1+ root-mu)))) + (values v mu root-mu))) + + +;; Could use the descending transform, but some of my tests show +;; that it has problems with roundoff errors. + +;; WARNING: This doesn't work very well for u > 1000 or so. For +;; example (elliptic-dn-ascending 1000b0 .5b0) -> 3.228b324, but dn <= 1. +#+nil +(defun elliptic-dn-ascending (u m) + (cond ((zerop m) + ;; A&S 16.6.3 + 1.0) + ((< (abs (- 1 m)) (* 4 (epsilon u))) + ;; A&S 16.6.3 + (/ (cosh u))) + (t + (multiple-value-bind (v mu root-mu1) + (ascending-transform u m) + ;; A&S 16.14.4 + (let* ((new-dn (elliptic-dn-ascending v mu))) + (* (/ (- 1 root-mu1) mu) + (/ (+ root-mu1 (* new-dn new-dn)) + new-dn))))))) + +;; Don't use the descending version because it requires cn, dn, and +;; sn. +;; +;; WARNING: This doesn't work very well for large u. +;; (elliptic-cn-ascending 1000b0 .5b0) -> 4.565b324. But |cn| <= 1. +#+nil +(defun elliptic-cn-ascending (u m) + (cond ((zerop m) + ;; A&S 16.6.2 + (cos u)) + ((< (abs (- 1 m)) (* 4 (epsilon u))) + ;; A&S 16.6.2 + (/ (cl:cosh u))) + (t + (multiple-value-bind (v mu root-mu1) + (ascending-transform u m) + ;; A&S 16.14.3 + (let* ((new-dn (elliptic-dn-ascending v mu))) + (* (/ (+ 1 root-mu1) mu) + (/ (- (* new-dn new-dn) root-mu1) + new-dn))))))) + +;; +;; This appears to work quite well for both real and complex values +;; of u. +(defun elliptic-sn-descending (u m) + (cond ((= m 1) + ;; A&S 16.6.1 + (tanh u)) + ((< (abs m) (epsilon u)) + ;; A&S 16.6.1 + (sin u)) + (t + (multiple-value-bind (v mu root-mu) + (descending-transform u m) + (let* ((new-sn (elliptic-sn-descending v mu))) + (/ (* (1+ root-mu) new-sn) + (1+ (* root-mu new-sn new-sn)))))))) + +;; We don't use the ascending transform here because it requires +;; evaluating sn, cn, and dn. The ascending transform only needs +;; sn. +#+nil +(defun elliptic-sn-ascending (u m) + (if (< (abs (- 1 m)) (* 4 flonum-epsilon)) + ;; A&S 16.6.1 + (tanh u) + (multiple-value-bind (v mu root-mu1) + (ascending-transform u m) + ;; A&S 16.14.2 + (let* ((new-cn (elliptic-cn-ascending v mu)) + (new-dn (elliptic-dn-ascending v mu)) + (new-sn (elliptic-sn-ascending v mu))) + (/ (* (+ 1 root-mu1) new-sn new-cn) + new-dn))))) + +(defun jacobi-sn (u m) + (let ((s (elliptic-sn-descending u m))) + (if (and (realp u) (realp m)) + (realpart s) + s))) + +(defun jacobi-dn (u m) + ;; Use the Gauss transformation from + ;; http://functions.wolfram.com/09.29.16.0013.01: + ;; + ;; + ;; dn((1+sqrt(m))*z, 4*sqrt(m)/(1+sqrt(m))^2) + ;; = (1-sqrt(m)*sn(z, m)^2)/(1+sqrt(m)*sn(z,m)^2) + ;; + ;; So + ;; + ;; dn(y, mu) = (1-sqrt(m)*sn(z, m)^2)/(1+sqrt(m)*sn(z,m)^2) + ;; + ;; where z = y/(1+sqrt(m)) and mu=4*sqrt(m)/(1+sqrt(m))^2. + ;; + ;; Solve for m, and we get + ;; + ;; sqrt(m) = -(mu+2*sqrt(1-mu)-2)/mu or (-mu+2*sqrt(1-mu)+2)/mu. + ;; + ;; I don't think it matters which sqrt we use, so I (rtoy) + ;; arbitrarily choose the first one above. + ;; + ;; Note that (1-sqrt(1-mu))/(1+sqrt(1-mu)) is the same as + ;; -(mu+2*sqrt(1-mu)-2)/mu. Also, the former is more + ;; accurate for small mu. + (let* ((root (let ((root-1-m (sqrt (- 1 m)))) + (/ (- 1 root-1-m) + (+ 1 root-1-m)))) + (z (/ u (+ 1 root))) + (s (elliptic-sn-descending z (* root root))) + (p (* root s s ))) + (/ (- 1 p) + (+ 1 p)))) + +(defun jacobi-cn (u m) + ;; Use the ascending Landen transformation, A&S 16.14.3. + (multiple-value-bind (v mu root-mu1) + (ascending-transform u m) + (let ((d (dn v mu))) + (* (/ (+ 1 root-mu1) mu) + (/ (- (* d d) root-mu1) + d))))) \ No newline at end of file diff --git a/qd-methods.lisp b/qd-methods.lisp index b6a83cf..3af7bd0 100644 --- a/qd-methods.lisp +++ b/qd-methods.lisp @@ -1089,3 +1089,23 @@ underlying floating-point format" ;; and make a real qd-real float, instead of the hackish ;; %qd-real. (set-dispatch-macro-character #\# #\Q #'qd-class-reader) + + +(defmethod epsilon ((m cl:float)) + (etypecase m + (single-float single-float-epsilon) + (double-float double-float-epsilon))) + +(defmethod epsilon ((m cl:complex)) + (epsilon (realpart m))) + +(defmethod epsilon ((m qd-real)) + ;; What is the epsilon value for a quad-double? This is complicated + ;; by the fact that things like (+ #q1 #q1q-100) is representable as + ;; a quad-double. For most purposes we want epsilon to be close to + ;; the 212 bits of precision (4*53 bits) that we normally have with + ;; a quad-double. + (scale-float #q1 -212)) + +(defmethod epsilon ((m qd-complex)) + (epsilon (realpart m))) ----------------------------------------------------------------------- Summary of changes: oct.asd | 2 + qd-elliptic.lisp | 180 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ qd-methods.lisp | 20 ++++++ 3 files changed, 202 insertions(+), 0 deletions(-) create mode 100644 qd-elliptic.lisp hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 258cd829289826458d6075caa279b253b8848749
by Raymond Toy 07 Mar '11

07 Mar '11

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 "OCT: A portable Lisp implementation for quad-double precision floats". The branch, master has been updated via 258cd829289826458d6075caa279b253b8848749 (commit) from 20558955b7953055a77ef1cbe690dcb1b1f54291 (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 258cd829289826458d6075caa279b253b8848749 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Sun Mar 6 21:27:16 2011 -0500 Remove extra blank lines. Fix #1. diff --git a/qd.lisp b/qd.lisp index 2b72769..af5a779 100644 --- a/qd.lisp +++ b/qd.lisp @@ -1346,11 +1346,6 @@ is a fixnum." (cl:* (cl:* q2 s1) s2) (cl:* (cl:* q3 s1) s2))))) - - - - - (defun decode-float-qd (q) "Like DECODE-FLOAT, but for %QUAD-DOUBLE numbers. Returns three values: 1) a %QUAD-DOUBLE number representing the significand. This is always ----------------------------------------------------------------------- Summary of changes: qd.lisp | 5 ----- 1 files changed, 0 insertions(+), 5 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 20558955b7953055a77ef1cbe690dcb1b1f54291
by Raymond Toy 04 Mar '11

04 Mar '11

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 "OCT: A portable Lisp implementation for quad-double precision floats". The branch, master has been updated via 20558955b7953055a77ef1cbe690dcb1b1f54291 (commit) from 0597a1124076fd119c55c7ecb57b765e1fb6e2ad (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 20558955b7953055a77ef1cbe690dcb1b1f54291 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Fri Mar 4 16:30:41 2011 -0500 Note this defsystem is deprecated. Fix #1. diff --git a/oct.system b/oct.system index 91d5e58..d815a15 100644 --- a/oct.system +++ b/oct.system @@ -23,6 +23,10 @@ ;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR ;;;; OTHER DEALINGS IN THE SOFTWARE. +;; This is no longer used, but we leave it here for reference and in +;; case someone wants to build oct with mk:defsys. Oct now uses asdf +;; to build. + (in-package #:cl-user) (mk:defsystem oct :source-pathname (make-pathname :directory (pathname-directory *load-pathname*)) ----------------------------------------------------------------------- Summary of changes: oct.system | 4 ++++ 1 files changed, 4 insertions(+), 0 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 0597a1124076fd119c55c7ecb57b765e1fb6e2ad
by Raymond Toy 04 Mar '11

04 Mar '11

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 "OCT: A portable Lisp implementation for quad-double precision floats". The branch, master has been updated via 0597a1124076fd119c55c7ecb57b765e1fb6e2ad (commit) from d28bce81063ce796763a7b7de995ee0157603ed9 (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 0597a1124076fd119c55c7ecb57b765e1fb6e2ad Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Fri Mar 4 16:24:11 2011 -0500 Remove trailing empty lines. Fix #1. diff --git a/branch-test.lisp b/branch-test.lisp index ee40765..1ae66fb 100644 --- a/branch-test.lisp +++ b/branch-test.lisp @@ -138,6 +138,3 @@ #+cmu (check-signs #'atanh 2w0 tr ti) (check-signs #'atanh #q2 tr ti))) - - - ----------------------------------------------------------------------- Summary of changes: branch-test.lisp | 3 --- 1 files changed, 0 insertions(+), 3 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. d28bce81063ce796763a7b7de995ee0157603ed9
by Raymond Toy 04 Mar '11

04 Mar '11

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 "OCT: A portable Lisp implementation for quad-double precision floats". The branch, master has been updated via d28bce81063ce796763a7b7de995ee0157603ed9 (commit) from 7310f8e9c4f8d905d60121a1bc89665b9f620689 (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 d28bce81063ce796763a7b7de995ee0157603ed9 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Fri Mar 4 16:08:27 2011 -0500 Add newline to end of file. Fix #1. diff --git a/LICENSE b/LICENSE index d4bf97c..90cf16a 100644 --- a/LICENSE +++ b/LICENSE @@ -16,4 +16,4 @@ FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN -THE SOFTWARE. \ No newline at end of file +THE SOFTWARE. ----------------------------------------------------------------------- Summary of changes: LICENSE | 2 +- 1 files changed, 1 insertions(+), 1 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 7310f8e9c4f8d905d60121a1bc89665b9f620689
by Raymond Toy 04 Mar '11

04 Mar '11

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 "OCT: A portable Lisp implementation for quad-double precision floats". The branch, master has been updated via 7310f8e9c4f8d905d60121a1bc89665b9f620689 (commit) via cc9eb34a6fa1e8971a09dddc91e93848ba360dbf (commit) from 6ab0812b935e4fa1c63ba9f5d58f7e66eb6f375b (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 7310f8e9c4f8d905d60121a1bc89665b9f620689 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Fri Mar 4 16:03:17 2011 -0500 Add newline to end of file. Fix #1. diff --git a/oct.asd b/oct.asd index 3cef6d7..947b047 100644 --- a/oct.asd +++ b/oct.asd @@ -75,4 +75,4 @@ (defmethod perform ((op test-op) (c (eql (find-system :oct-tests)))) (or (funcall (intern "DO-TESTS" (find-package "RT"))) - (error "TEST-OP failed for OCT-TESTS"))) \ No newline at end of file + (error "TEST-OP failed for OCT-TESTS"))) commit cc9eb34a6fa1e8971a09dddc91e93848ba360dbf Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Sun Feb 13 18:41:29 2011 -0500 Add test-op for oct to run the tests. diff --git a/oct.asd b/oct.asd index ba4d6b6..3cef6d7 100644 --- a/oct.asd +++ b/oct.asd @@ -27,7 +27,7 @@ ;;; so it might be out of date. Use at your own risk. (defpackage #:oct-system - (:use #:cl)) + (:use #:cl #:asdf)) (in-package #:oct-system) @@ -59,3 +59,20 @@ (:file "qd-complex" :depends-on ("qd-methods")) )) + +(defmethod perform ((op test-op) (c (eql (find-system :oct)))) + (oos 'test-op 'oct-tests)) + +(asdf:defsystem oct-tests + :depends-on (oct) + :version "2011-02-09" ; Just use the date + :in-order-to ((compile-op (load-op :rt)) + (test-op (load-op :rt :oct))) + :components + ((:file "qd-extra") + (:file "qd-test") + (:file "rt-tests"))) + +(defmethod perform ((op test-op) (c (eql (find-system :oct-tests)))) + (or (funcall (intern "DO-TESTS" (find-package "RT"))) + (error "TEST-OP failed for OCT-TESTS"))) \ No newline at end of file ----------------------------------------------------------------------- Summary of changes: oct.asd | 19 ++++++++++++++++++- 1 files changed, 18 insertions(+), 1 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0