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

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

11 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 c239a8514ff5f1fc96d33a4c581cbaec834f6641 (commit) via 46359140b7aa3c0c4af6b3aabc7bcbb4cf248301 (commit) via 85f3e2b22f917ff35d080180ab1e2dda4476e881 (commit) from 68432cb1855605c346216788f1aa64517a96a808 (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 c239a8514ff5f1fc96d33a4c581cbaec834f6641 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Thu Mar 10 22:23:50 2011 -0500 Implement elliptic integrals of the first kind. o Add FLOAT-CONTAGION to determine the max precision of the given arguments so we can do appopriate contagion in the routines. o Add some docstrings and other documentation of the algorithms. o Add implmentation of ELLIPTIC-K and ELLIPTIC-F for the complete and incomplete elliptic integrals of the first kind, respectively. diff --git a/qd-elliptic.lisp b/qd-elliptic.lisp index 2fadecd..63c3687 100644 --- a/qd-elliptic.lisp +++ b/qd-elliptic.lisp @@ -29,6 +29,57 @@ (declaim (inline descending-transform ascending-transform)) +;; Determine which of x and y has the higher precision and return the +;; value of the higher precision number. If both x and y are +;; rationals, just return 1f0, for a single-float value. +(defun float-contagion-2 (x y) + (etypecase x + (cl:rational + (etypecase y + (cl:rational + 1f0) + (cl:float + y) + (qd-real + y))) + (single-float + (etypecase y + ((or cl:rational single-float) + x) + ((or double-float qd-real) + y))) + (double-float + (etypecase y + ((or cl:rational single-float double-float) + x) + (qd-real + y))) + (qd-real + x))) + +;; Return a floating point (or complex) type of the highest precision +;; among all of the given arguments. +(defun float-contagion (&rest args) + ;; It would be easy if we could just add the args together and let + ;; normal contagion do the work, but we could easily introduce + ;; overflows or other errors that way. So look at each argument and + ;; determine the precision and choose the highest precision. + (let ((complexp (some #'complexp args)) + (max-type + (etypecase (reduce #'float-contagion-2 (mapcar #'realpart (if (cdr args) + args + (list (car args) 0)))) + (single-float 'single-float) + (double-float 'double-float) + (qd-real 'qd-real)))) + (if complexp + (if (eq max-type 'qd-real) + 'qd-complex + `(cl:complex ,max-type)) + max-type))) + +;;; Jacobian elliptic functions + (defun ascending-transform (u m) ;; A&S 16.14.1 ;; @@ -153,12 +204,14 @@ new-dn))))) (defun jacobi-sn (u m) + "Compute Jacobian sn for argument u and parameter m" (let ((s (elliptic-sn-descending u m))) (if (and (realp u) (realp m)) (realpart s) s))) (defun jacobi-dn (u m) + "Compute Jacobi dn for argument u and parameter m" ;; Use the Gauss transformation from ;; http://functions.wolfram.com/09.29.16.0013.01: ;; @@ -192,6 +245,7 @@ (+ 1 p)))) (defun jacobi-cn (u m) + "Compute Jacobi cn for argument u and parameter m" ;; Use the ascending Landen transformation, A&S 16.14.3. ;; ;; cn(u,m) = (1+sqrt(mu1))/mu * (dn(v,mu)^2-sqrt(mu1))/dn(v,mu) @@ -202,17 +256,32 @@ (/ (- (* d d) root-mu1) d))))) +;;; Elliptic Integrals +;;; +;; Translation of Jim FitzSimons' bigfloat implementation of elliptic +;; integrals from http://www.getnet.com/~cherry/elliptbf3.mac. +;; +;; The algorithms are based on B.C. Carlson's "Numerical Computation +;; of Real or Complex Elliptic Integrals". These are updated to the +;; algorithms in Journal of Computational and Applied Mathematics 118 +;; (2000) 71-85 "Reduction Theorems for Elliptic Integrands with the +;; Square Root of two quadritic factors" +;; + (defun errtol (&rest args) - ;; Compute error tolerance as sqrt(2^(-fpprec)). Not sure this is - ;; quite right, but it makes the routines more accurate as fpprec - ;; increases. + ;; Compute error tolerance as sqrt(<float-precision>). Not sure + ;; this is quite right, but it makes the routines more accurate as + ;; precision increases increases. (sqrt (reduce #'min (mapcar #'(lambda (x) (if (rationalp x) - double-float-epsilon + single-float-epsilon (epsilon x))) args)))) (defun carlson-rf (x y z) + "Compute Carlson's Rf function: + + Rf(x, y, z) = 1/2*integrate((t+x)^(-1/2)*(t+y)^(-1/2)*(t+z)^(-1/2), t, 0, inf)" (let* ((xn x) (yn y) (zn z) @@ -252,8 +321,90 @@ (* -3/44 ee2 ee3)))) (/ s (sqrt an))))) +;; Complete elliptic integral of the first kind. This can be computed +;; from Carlson's Rf function: +;; +;; K(m) = Rf(0, 1 - m, 1) (defun elliptic-k (m) + "Complete elliptic integral of the first kind K for parameter m + + K(m) = integrate(1/sqrt(1-m*sin(x)^2), x, 0, %pi/2). + + Note: K(m) = F(%pi/2, m), where F is the (incomplete) elliptic + integral of the first kind." (cond ((= m 0) (/ (float +pi+ m) 2)) (t - (carlson-rf 0 (- 1 m) 1)))) + (let ((precision (float-contagion m))) + (carlson-rf (coerce 0 precision) (- 1 m) (coerce 1 precision)))))) + +;; Elliptic integral of the first kind. This is computed using +;; Carlson's Rf function: +;; +;; F(phi, m) = sin(phi) * Rf(cos(phi)^2, 1 - m*sin(phi)^2, 1) +(defun elliptic-f (x m) + "Incomplete Elliptic integral of the first kind: + + F(x, m) = integrate(1/sqrt(1-m*sin(phi)^2), phi, 0, x) + + Note for the complete elliptic integral, you can use elliptic-k" + (let* ((precision (float-contagion x m)) + (x (coerce x precision)) + (m (coerce m precision))) + (cond ((and (realp m) (realp x)) + (cond ((> m 1) + ;; A&S 17.4.15 + ;; + ;; F(phi|m) = 1/sqrt(m)*F(theta|1/m) + ;; + ;; with sin(theta) = sqrt(m)*sin(phi) + (/ (elliptic-f (cl:asin (* (sqrt m) (sin x))) (/ m)) + (sqrt m))) + ((< m 0) + ;; A&S 17.4.17 + (let* ((m (- m)) + (m+1 (+ 1 m)) + (root (sqrt m+1)) + (m/m+1 (/ m m+1))) + (- (/ (elliptic-f (/ (float-pi m) 2) m/m+1) + root) + (/ (elliptic-f (- (/ (float-pi x) 2) x) m/m+1) + root)))) + ((= m 0) + ;; A&S 17.4.19 + x) + ((= m 1) + ;; A&S 17.4.21 + ;; + ;; F(phi,1) = log(sec(phi)+tan(phi)) + ;; = log(tan(pi/4+pi/2)) + (log (cl:tan (+ (/ x 2) (/ (float-pi x) 4))))) + ((minusp x) + (- (elliptic-f (- x) m))) + ((> x (float-pi x)) + ;; A&S 17.4.3 + (multiple-value-bind (s x-rem) + (truncate x (float-pi x)) + (+ (* 2 s (elliptic-k m)) + (elliptic-f x-rem m)))) + ((<= x (/ (float-pi x) 2)) + (let ((sin-x (sin x)) + (cos-x (cos x)) + (k (sqrt m))) + (* sin-x + (carlson-rf (* cos-x cos-x) + (* (- 1 (* k sin-x)) + (+ 1 (* k sin-x))) + 1.0)))) + ((< x (float-pi x)) + (+ (* 2 (elliptic-k m)) + (elliptic-f (- x (float pi x)) m))))) + (t + (let ((sin-x (sin x)) + (cos-x (cos x)) + (k (sqrt m))) + (* sin-x + (carlson-rf (* cos-x cos-x) + (* (- 1 (* k sin-x)) + (+ 1 (* k sin-x))) + 1))))))) commit 46359140b7aa3c0c4af6b3aabc7bcbb4cf248301 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Thu Mar 10 22:21:26 2011 -0500 Add documentation for EPSILON and add FLOAT-PI. FLOAT-PI returns a value of pi that matches the precision of the argument. diff --git a/qd-methods.lisp b/qd-methods.lisp index 6f41857..966d85f 100644 --- a/qd-methods.lisp +++ b/qd-methods.lisp @@ -1063,7 +1063,13 @@ underlying floating-point format" (frob two-arg-- cl:- sub-qd sub-d-qd sub-qd-d) (frob two-arg-* cl:* mul-qd mul-d-qd mul-qd-d) (frob two-arg-/ cl:/ div-qd nil nil)) - + +(defgeneric epsilon (m) + (:documentation +"Return an epsilon value of the same precision as the argument. It is +the smallest number x such that 1+x /= x. The argument can be +complex")) + (defmethod epsilon ((m cl:float)) (etypecase m (single-float single-float-epsilon) @@ -1082,3 +1088,23 @@ underlying floating-point format" (defmethod epsilon ((m qd-complex)) (epsilon (realpart m))) + +(defgeneric float-pi (x) + (:documentation +"Return a floating-point value of the mathematical constant pi that is +the same precision as the argument. The argument can be complex.")) + +(defmethod float-pi ((x cl:rational)) + (float pi 1f0)) + +(defmethod float-pi ((x cl:float)) + (float pi x)) + +(defmethod float-pi ((x qd-real)) + +pi+) + +(defmethod float-pi ((z cl:complex)) + (float pi (realpart z))) + +(defmethod float-pi ((z qd-complex)) + +pi+) \ No newline at end of file commit 85f3e2b22f917ff35d080180ab1e2dda4476e881 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Thu Mar 10 22:20:27 2011 -0500 Move inline function NEG-QD-T before first use. diff --git a/qd.lisp b/qd.lisp index af5a779..fb3e4e1 100644 --- a/qd.lisp +++ b/qd.lisp @@ -476,11 +476,6 @@ If TARGET is given, TARGET is destructively modified to contain the result." (%store-qd-d target (+ a0 b0) 0d0 0d0 0d0) (%store-qd-d target s0 s1 s2 s3))))))))))) -(defun neg-qd (a &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) - "Return the negative of the %QUAD-DOUBLE number A. -If TARGET is given, TARGET is destructively modified to contain the result." - (neg-qd-t a target)) - (defun neg-qd-t (a target) (declare (type %quad-double a #+oct-array target) #+(and cmu (not oct-array)) (ignore target)) @@ -489,6 +484,13 @@ If TARGET is given, TARGET is destructively modified to contain the result." (declare (double-float a0 a1 a2 a3)) (%store-qd-d target (cl:- a0) (cl:- a1) (cl:- a2) (cl:- a3)))) +(defun neg-qd (a &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) + "Return the negative of the %QUAD-DOUBLE number A. +If TARGET is given, TARGET is destructively modified to contain the result." + (neg-qd-t a target)) + + + (defun sub-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) "Return the difference between the %QUAD-DOUBLE numbers A and B. If TARGET is given, TARGET is destructively modified to contain the result." ----------------------------------------------------------------------- Summary of changes: qd-elliptic.lisp | 161 ++++++++++++++++++++++++++++++++++++++++++++++++++++-- qd-methods.lisp | 28 +++++++++- qd.lisp | 12 ++-- 3 files changed, 190 insertions(+), 11 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. 68432cb1855605c346216788f1aa64517a96a808
by Raymond Toy 08 Mar '11

08 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 68432cb1855605c346216788f1aa64517a96a808 (commit) via 3d37e3b9346092b58e820723ae20f114dc912d44 (commit) from 6c91244701c5867328a58c94d2118db6e09310d1 (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 68432cb1855605c346216788f1aa64517a96a808 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Tue Mar 8 09:02:42 2011 -0500 Oops. Remove the #q reader functions from here. They've been moved to qd-reader.lisp. diff --git a/qd-methods.lisp b/qd-methods.lisp index 4f8065f..6f41857 100644 --- a/qd-methods.lisp +++ b/qd-methods.lisp @@ -1064,36 +1064,6 @@ underlying floating-point format" (frob two-arg-* cl:* mul-qd mul-d-qd mul-qd-d) (frob two-arg-/ cl:/ div-qd nil nil)) - -(defun read-qd-real-or-complex (stream) - (let ((c (peek-char t stream))) - (cond ((char= c #\() - ;; Read a QD complex - (read-char stream) ; Skip the paren - (let ((real (read stream t nil t)) - (imag (read stream t nil t))) - (unless (char= (peek-char t stream) #\)) - (error "Illegal QD-COMPLEX number format")) - ;; Read closing paren - (read-char stream) - (make-instance 'qd-complex - :real (qd-value (float real +qd-real-one+)) - :imag (qd-value (float imag +qd-real-one+))))) - (t - (make-instance 'qd-real :value (read-qd stream)))))) - -(defun qd-class-reader (stream subchar arg) - (declare (ignore subchar)) - (when arg - (warn "Numeric argument ignored in #~DQ" arg)) - (read-qd-real-or-complex stream)) - -;; Yow! We redefine the #q reader that is in qd-io.lisp to read in -;; 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) commit 3d37e3b9346092b58e820723ae20f114dc912d44 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Tue Mar 8 09:01:06 2011 -0500 Move #Q reader functions to new file qd-reader.lisp. qd-reader.lisp: o New file containing #Q reader functions. o Create a new readtable containing our reader functions so we don't destructively modify the default *readtable*. oct.asd: o Add qd-reader.lisp. qd-complex.lisp: qd-elliptic.lisp: qd-format.lisp: o Set the *readtable* to *oct-readtable* before compiling these files. diff --git a/oct.asd b/oct.asd index 153706c..77828aa 100644 --- a/oct.asd +++ b/oct.asd @@ -54,12 +54,14 @@ :depends-on ("qd-fun")) (:file "qd-methods" :depends-on ("qd-class")) - (:file "qd-format" + (:file "qd-reader" :depends-on ("qd-methods")) + (:file "qd-format" + :depends-on ("qd-methods" "qd-reader")) (:file "qd-complex" - :depends-on ("qd-methods")) + :depends-on ("qd-methods" "qd-reader")) (:file "qd-elliptic" - :depends-on ("qd-methods")) + :depends-on ("qd-methods" "qd-reader")) )) (defmethod perform ((op test-op) (c (eql (find-system :oct)))) diff --git a/qd-complex.lisp b/qd-complex.lisp index 2a2aaf7..409e94f 100644 --- a/qd-complex.lisp +++ b/qd-complex.lisp @@ -25,6 +25,9 @@ (in-package #:oct) +(eval-when (:compile-toplevel :load-toplevel :execute) + (setf *readtable* *oct-readtable*)) + (defmethod add1 ((a qd-complex)) (make-instance 'qd-complex :real (add-qd-d (qd-value (realpart a)) 1d0) diff --git a/qd-elliptic.lisp b/qd-elliptic.lisp index 3336444..2fadecd 100644 --- a/qd-elliptic.lisp +++ b/qd-elliptic.lisp @@ -24,6 +24,9 @@ (in-package #:oct) +(eval-when (:compile-toplevel :load-toplevel :execute) + (setf *readtable* *oct-readtable*)) + (declaim (inline descending-transform ascending-transform)) (defun ascending-transform (u m) diff --git a/qd-format.lisp b/qd-format.lisp index a63d814..5ae5289 100644 --- a/qd-format.lisp +++ b/qd-format.lisp @@ -25,6 +25,9 @@ (in-package #:oct) +(eval-when (:compile-toplevel :load-toplevel :execute) + (setf *readtable* *oct-readtable*)) + (defun qd-scale-exponent (original-x) (let* ((x original-x)) (multiple-value-bind (sig exponent) diff --git a/qd-reader.lisp b/qd-reader.lisp new file mode 100644 index 0000000..d2159ea --- /dev/null +++ b/qd-reader.lisp @@ -0,0 +1,58 @@ +;;;; -*- Mode: lisp -*- +;;;; +;;;; Copyright (c) 2007, 2008, 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) + + + +(defun read-qd-real-or-complex (stream) + (let ((c (peek-char t stream))) + (cond ((char= c #\() + ;; Read a QD complex + (read-char stream) ; Skip the paren + (let ((real (read stream t nil t)) + (imag (read stream t nil t))) + (unless (char= (peek-char t stream) #\)) + (error "Illegal QD-COMPLEX number format")) + ;; Read closing paren + (read-char stream) + (make-instance 'qd-complex + :real (qd-value (float real +qd-real-one+)) + :imag (qd-value (float imag +qd-real-one+))))) + (t + (make-instance 'qd-real :value (read-qd stream)))))) + +(defun qd-class-reader (stream subchar arg) + (declare (ignore subchar)) + (when arg + (warn "Numeric argument ignored in #~DQ" arg)) + (read-qd-real-or-complex stream)) + +(defvar *oct-readtable* + (let ((rt (copy-readtable nil))) + (set-dispatch-macro-character #\# #\Q #'qd-class-reader rt) + rt) + "Readtable that extends the standard readtable to include #q for + reading QD-REAL and QD-COMPLEX numbers") ----------------------------------------------------------------------- Summary of changes: oct.asd | 8 +++-- qd-complex.lisp | 3 ++ qd-elliptic.lisp | 3 ++ qd-format.lisp | 3 ++ qd-methods.lisp | 30 ---------------------- oct-test.system => qd-reader.lisp | 50 +++++++++++++++++++++++------------- 6 files changed, 46 insertions(+), 51 deletions(-) copy oct-test.system => qd-reader.lisp (51%) 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. 6c91244701c5867328a58c94d2118db6e09310d1
by Raymond Toy 08 Mar '11

08 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 6c91244701c5867328a58c94d2118db6e09310d1 (commit) from f86bc7e7f3c5b617a141ac5d0f32b5fc6b5e77d5 (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 6c91244701c5867328a58c94d2118db6e09310d1 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Tue Mar 8 08:56:13 2011 -0500 For clisp, disable floating-point underflow when running tests. diff --git a/rt-tests.lisp b/rt-tests.lisp index f4e06c1..82ca0b7 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -25,6 +25,11 @@ (in-package #:oct) +;; For the tests, we need to turn off underflow for clisp. +#+clisp +(ext:without-package-lock () + (setq sys::*inhibit-floating-point-underflow* t)) + ;; Compute how many bits are the same for two numbers EST and TRUE. ;; Return T if they are identical. (defun bit-accuracy (est true) ----------------------------------------------------------------------- Summary of changes: rt-tests.lisp | 5 +++++ 1 files changed, 5 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. f86bc7e7f3c5b617a141ac5d0f32b5fc6b5e77d5
by Raymond Toy 08 Mar '11

08 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 f86bc7e7f3c5b617a141ac5d0f32b5fc6b5e77d5 (commit) via 1a38c9e3dd557620c053919bbccde06741118a36 (commit) from 8481ee30bafaec66ccb99a330e1aee392c4d894b (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 f86bc7e7f3c5b617a141ac5d0f32b5fc6b5e77d5 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Mon Mar 7 19:36:35 2011 -0500 Add elliptic K function and tests for it and Jacobi functions. qd-elliptic.lisp: o Add support for the complete elliptic integral K using Carlson's Rf function. rt-tests.lisp: o Fix indentation for oct.atan.5 o Add tests for elliptic K o Add tests for Jacobi sn, cn, and dn functions. diff --git a/qd-elliptic.lisp b/qd-elliptic.lisp index 9240eb0..3336444 100644 --- a/qd-elliptic.lisp +++ b/qd-elliptic.lisp @@ -198,3 +198,59 @@ (* (/ (+ 1 root-mu1) mu) (/ (- (* d d) root-mu1) d))))) + +(defun errtol (&rest args) + ;; Compute error tolerance as sqrt(2^(-fpprec)). Not sure this is + ;; quite right, but it makes the routines more accurate as fpprec + ;; increases. + (sqrt (reduce #'min (mapcar #'(lambda (x) + (if (rationalp x) + double-float-epsilon + (epsilon x))) + args)))) + +(defun carlson-rf (x y z) + (let* ((xn x) + (yn y) + (zn z) + (a (/ (+ xn yn zn) 3)) + (epslon (/ (max (abs (- a xn)) + (abs (- a yn)) + (abs (- a zn))) + (errtol x y z))) + (an a) + (power4 1) + (n 0) + xnroot ynroot znroot lam) + (loop while (> (* power4 epslon) (abs an)) + do + (setf xnroot (sqrt xn)) + (setf ynroot (sqrt yn)) + (setf znroot (sqrt zn)) + (setf lam (+ (* xnroot ynroot) + (* xnroot znroot) + (* ynroot znroot))) + (setf power4 (* power4 1/4)) + (setf xn (* (+ xn lam) 1/4)) + (setf yn (* (+ yn lam) 1/4)) + (setf zn (* (+ zn lam) 1/4)) + (setf an (* (+ an lam) 1/4)) + (incf n)) + ;; c1=-3/14,c2=1/6,c3=9/88,c4=9/22,c5=-3/22,c6=-9/52,c7=3/26 + (let* ((xndev (/ (* (- a x) power4) an)) + (yndev (/ (* (- a y) power4) an)) + (zndev (- (+ xndev yndev))) + (ee2 (- (* xndev yndev) (* 6 zndev zndev))) + (ee3 (* xndev yndev zndev)) + (s (+ 1 + (* -1/10 ee2) + (* 1/14 ee3) + (* 1/24 ee2 ee2) + (* -3/44 ee2 ee3)))) + (/ s (sqrt an))))) + +(defun elliptic-k (m) + (cond ((= m 0) + (/ (float +pi+ m) 2)) + (t + (carlson-rf 0 (- 1 m) 1)))) diff --git a/rt-tests.lisp b/rt-tests.lisp index b9e5867..f4e06c1 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -118,8 +118,8 @@ (rt:deftest oct.atan.5 (let* ((arg #q-1q100) - (y (/ (atan arg) +pi+)) - (true #q-.5)) + (y (/ (atan arg) +pi+)) + (true #q-.5)) (check-accuracy 212 y true)) nil) @@ -761,3 +761,70 @@ (let ((true (cl:atanh #c(2d0 1d-20)))) (check-signs #'atanh #q2 true)) t) + +;; elliptic_k(-1) = gamma(1/4)^2/2^(5/2)/sqrt(%pi) +(rt:deftest oct.elliptic-k.1d + (let* ((val (elliptic-k -1d0)) + (true #q1.311028777146059905232419794945559706841377475715811581408410851900395293535207125115147766480714547q0)) + (check-accuracy 53 val true)) + nil) + +(rt:deftest oct.elliptic-k.1q + (let* ((val (elliptic-k #q-1q0)) + (true #q1.311028777146059905232419794945559706841377475715811581408410851900395293535207125115147766480714547q0)) + (check-accuracy 210 val true)) + nil) + +;; elliptic_k(1/2) = %pi^(3/2)/2/gamma(3/4)^2 +(rt:deftest oct.elliptic-k.2d + (let* ((val (elliptic-k 0.5d0)) + (true #q1.854074677301371918433850347195260046217598823521766905585928045056021776838119978357271861650371897q0)) + (check-accuracy 53 val true)) + nil) + +(rt:deftest oct.elliptic-k.2q + (let* ((val (elliptic-k #q.5)) + (true #q1.854074677301371918433850347195260046217598823521766905585928045056021776838119978357271861650371897q0)) + (check-accuracy 210 val true)) + nil) + +;; jacobi_sn(K,1/2) = 1, where K = elliptic_k(1/2) +(rt:deftest oct.jacobi-sn.1d + (let* ((ek (elliptic-k .5d0)) + (val (jacobi-sn ek .5d0))) + (check-accuracy 54 val 1d0)) + nil) + +(rt:deftest oct.jacobi-sn.1q + (let* ((ek (elliptic-k #q.5)) + (val (jacobi-sn ek #q.5))) + (check-accuracy 212 val #q1)) + nil) + +;; jacobi_cn(K,1/2) = 0 +(rt:deftest oct.jacobi-cn.1d + (let* ((ek (elliptic-k .5d0)) + (val (jacobi-cn ek .5d0))) + (check-accuracy 50 val 0d0)) + nil) + +(rt:deftest oct.jacobi-sn.1q + (let* ((ek (elliptic-k #q.5)) + (val (jacobi-cn ek #q.5))) + (check-accuracy 210 val #q0)) + nil) + +;; jacobi-dn(K, 1/2) = sqrt(1/2) +(rt:deftest oct.jacobi-dn.1d + (let* ((ek (elliptic-k .5d0)) + (true (sqrt .5d0)) + (val (jacobi-dn ek .5d0))) + (check-accuracy 52 val true)) + nil) + +(rt:deftest oct.jacobi-dn.1q + (let* ((ek (elliptic-k #q.5)) + (true (sqrt #q.5)) + (val (jacobi-dn ek #q.5))) + (check-accuracy 212 val true)) + nil) commit 1a38c9e3dd557620c053919bbccde06741118a36 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Mon Mar 7 19:33:34 2011 -0500 Fix bug in qd-scale-float. MAKE-QD-D returns a %quad-double, not qd-real. Use +QD-REAL-ONE+ instead. diff --git a/qd-methods.lisp b/qd-methods.lisp index 037bcc3..4f8065f 100644 --- a/qd-methods.lisp +++ b/qd-methods.lisp @@ -1108,7 +1108,7 @@ underlying floating-point format" ;; 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 (make-qd-d 1d0) -212)) + (scale-float +qd-real-one+ -212)) (defmethod epsilon ((m qd-complex)) (epsilon (realpart m))) ----------------------------------------------------------------------- Summary of changes: qd-elliptic.lisp | 56 ++++++++++++++++++++++++++++++++++++++++++ qd-methods.lisp | 2 +- rt-tests.lisp | 71 ++++++++++++++++++++++++++++++++++++++++++++++++++++- 3 files changed, 126 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. 8481ee30bafaec66ccb99a330e1aee392c4d894b
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 8481ee30bafaec66ccb99a330e1aee392c4d894b (commit) from 25a15317f76cba7a1c38341edf5a418020e9b477 (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 8481ee30bafaec66ccb99a330e1aee392c4d894b Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Mon Mar 7 17:55:50 2011 -0500 Oops. Use jacobi-dn, not dn. diff --git a/qd-elliptic.lisp b/qd-elliptic.lisp index c29d9ce..9240eb0 100644 --- a/qd-elliptic.lisp +++ b/qd-elliptic.lisp @@ -194,7 +194,7 @@ ;; cn(u,m) = (1+sqrt(mu1))/mu * (dn(v,mu)^2-sqrt(mu1))/dn(v,mu) (multiple-value-bind (v mu root-mu1) (ascending-transform u m) - (let ((d (dn v mu))) + (let ((d (jacobi-dn v mu))) (* (/ (+ 1 root-mu1) mu) (/ (- (* d d) root-mu1) - d))))) \ No newline at end of file + d))))) ----------------------------------------------------------------------- Summary of changes: qd-elliptic.lisp | 4 ++-- 1 files changed, 2 insertions(+), 2 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. 25a15317f76cba7a1c38341edf5a418020e9b477
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 25a15317f76cba7a1c38341edf5a418020e9b477 (commit) from 3d2184f115538262f60b359249d76fb1c4c7048d (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 25a15317f76cba7a1c38341edf5a418020e9b477 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Mon Mar 7 15:48:05 2011 -0500 Clean up printing and fix printing bug in non-cmucl. Clisp was printing #q(1 2) as #q(#q#q1 #q#q2). This is caused by PRINT-OBJECT for QD-COMPLEX printing a #q in addition to a #q printed by QD-FORMAT. qd-class.lisp: o Use the same PRINT-OBJECT method for qd-real for all Lisps. o Remove extraneous #q from PRINT-OBJECT method for qd-complex. qd-io.lisp: o QD-FORMAT for cmucl needs to output #q. diff --git a/qd-class.lisp b/qd-class.lisp index 5cc3cb4..2134ef7 100644 --- a/qd-class.lisp +++ b/qd-class.lisp @@ -47,21 +47,9 @@ :type %quad-double)) (:documentation "Complex number consisting of QUAD-DOUBLE components")) -#-cmu (defmethod print-object ((qd qd-real) stream) (format stream "~/octi::qd-format/" (qd-value qd))) -#+cmu -(defun print-qd (q stream) - (declare (type %quad-double q)) - (if (or (float-infinity-p (qd-0 q)) - (float-nan-p (qd-0 q))) - (format stream "~/octi::qd-format/" q) - (format stream "#q~/octi::qd-format/" q))) -#+cmu -(defmethod print-object ((qd qd-real) stream) - (print-qd (qd-value qd) stream)) - (defmethod make-qd ((x cl:float)) (make-instance 'qd-real :value (make-qd-d (float x 1d0)))) @@ -69,7 +57,7 @@ (make-instance 'qd-real :value (qd-value x))) (defmethod print-object ((qd qd-complex) stream) - (format stream "#q(~<#q~/octi::qd-format/ #q~/octi::qd-format/~:@>)" + (format stream "#q(~<~/octi::qd-format/ ~/octi::qd-format/~:@>)" (list (qd-real qd) (qd-imag qd)))) diff --git a/qd-io.lisp b/qd-io.lisp index 1df662e..40d7c90 100644 --- a/qd-io.lisp +++ b/qd-io.lisp @@ -303,6 +303,7 @@ ((ext:float-nan-p (qd-0 arg)) (qd-output-nan arg stream)) (t + (write-string "#q" stream) (qd-output-aux arg stream)))) ;; This version has problems with roundoff. ----------------------------------------------------------------------- Summary of changes: qd-class.lisp | 14 +------------- qd-io.lisp | 1 + 2 files changed, 2 insertions(+), 13 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. 3d2184f115538262f60b359249d76fb1c4c7048d
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 3d2184f115538262f60b359249d76fb1c4c7048d (commit) via 1e343ca0be49b9e56938a8063b4d2e679e338965 (commit) via 935f7e0f335c1ba01a17c6c6c39e0421f5f62b50 (commit) from 1421224af506145ab18f9f4b412f300c97c11751 (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 3d2184f115538262f60b359249d76fb1c4c7048d Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Mon Mar 7 14:59:00 2011 -0500 Support testing with Lisps without signed floating-point zeros. CHECK-SIGN just returns T if Lisp thinks -0d0 and 0d0 are the same number. This is only true if Lisp doesn't support signed zeroes. diff --git a/rt-tests.lisp b/rt-tests.lisp index b87100e..b9e5867 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -583,8 +583,11 @@ (let* ((z (funcall fun arg)) (x (realpart z)) (y (imagpart z))) - (if (and (= (float-sign x) (float-sign (realpart expected))) - (= (float-sign y) (float-sign (imagpart expected)))) + ;; If the Lisp doesn't support signed zeroes, then this test + ;; should always pass. + (if (or (eql -0d0 0d0) + (and (= (float-sign x) (float-sign (realpart expected))) + (= (float-sign y) (float-sign (imagpart expected))))) t (list z expected fun arg)))) commit 1e343ca0be49b9e56938a8063b4d2e679e338965 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Mon Mar 7 14:50:44 2011 -0500 Fix issue with float-sign calling qfloat-sign. If the optional arg to FLOAT-SIGN is not given, don't call qfloat-sign with a second arg of NIL. This breaks things. Call qfloat-sign with the same number of arguments as float-sign. diff --git a/qd-methods.lisp b/qd-methods.lisp index 205c9fe..037bcc3 100644 --- a/qd-methods.lisp +++ b/qd-methods.lisp @@ -717,6 +717,7 @@ underlying floating-point format" (defmethod qfloat-sign ((a real) &optional (f (float 1 a))) (cl:float-sign a f)) + (defmethod qfloat-sign ((a qd-real) &optional f) (if f (make-instance 'qd-real @@ -725,8 +726,10 @@ underlying floating-point format" (make-instance 'qd-real :value (make-qd-d (cl:float-sign (qd-0 (qd-value a))))))) (declaim (inline float-sign)) -(defun float-sign (n &optional float2) - (qfloat-sign n float2)) +(defun float-sign (n &optional (float2 nil float2p)) + (if float2p + (qfloat-sign n float2) + (qfloat-sign n))) (defun max (number &rest more-numbers) "Returns the greatest of its arguments." commit 935f7e0f335c1ba01a17c6c6c39e0421f5f62b50 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Mon Mar 7 14:47:00 2011 -0500 Ignore *.fas files. diff --git a/.gitignore b/.gitignore index 39ef52b..62c3cc0 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,7 @@ *.sse2f *.x86f *.sparcf +*.fas *.fasl *.err *~ ----------------------------------------------------------------------- Summary of changes: .gitignore | 1 + qd-methods.lisp | 7 +++++-- rt-tests.lisp | 7 +++++-- 3 files changed, 11 insertions(+), 4 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. 1421224af506145ab18f9f4b412f300c97c11751
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 1421224af506145ab18f9f4b412f300c97c11751 (commit) from c014c5f4e4b941aacf16a897a2cf689ed49b37ef (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 1421224af506145ab18f9f4b412f300c97c11751 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Sun Mar 6 22:38:09 2011 -0500 Document algorithms better. diff --git a/qd-elliptic.lisp b/qd-elliptic.lisp index 9c56fcd..c29d9ce 100644 --- a/qd-elliptic.lisp +++ b/qd-elliptic.lisp @@ -34,6 +34,12 @@ ;; (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). + ;; + ;; mu = 4*sqrt(m)/(1+sqrt(m))^2 + ;; sqrt(mu1) = (1-sqrt(m))/(1+sqrt(m)) + ;; v = u/(1+sqrt(mu1)) + ;; + ;; Return v, mu, sqrt(mu1) (let* ((root-m (sqrt m)) (mu (/ (* 4 root-m) (expt (1+ root-m) 2))) @@ -42,10 +48,19 @@ (values v mu root-mu1))) (defun descending-transform (u m) + ;; A&S 16.12.1 + ;; ;; 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. + ;; + ;; sqrt(mu) = (1-sqrt(m1))/(1+sqrt(m1)) + ;; v = u/(1+sqrt(mu)) + ;; + ;; where m1 = 1-m + ;; + ;; Return v, mu, sqrt(mu) (let* ((root-m1 (sqrt (- 1 m))) (root-mu (/ (- 1 root-m1) (+ 1 root-m1))) (mu (* root-mu root-mu)) @@ -108,6 +123,9 @@ ;; A&S 16.6.1 (sin u)) (t + ;; A&S 16.12.2 + ;; + ;; sn(u|m) = (1 + sqrt(mu))*sn(v|u)/(1 + sqrt(mu)*sn(v|mu)^2) (multiple-value-bind (v mu root-mu) (descending-transform u m) (let* ((new-sn (elliptic-sn-descending v mu))) @@ -172,6 +190,8 @@ (defun jacobi-cn (u m) ;; Use the ascending Landen transformation, A&S 16.14.3. + ;; + ;; cn(u,m) = (1+sqrt(mu1))/mu * (dn(v,mu)^2-sqrt(mu1))/dn(v,mu) (multiple-value-bind (v mu root-mu1) (ascending-transform u m) (let ((d (dn v mu))) ----------------------------------------------------------------------- Summary of changes: qd-elliptic.lisp | 20 ++++++++++++++++++++ 1 files changed, 20 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. c014c5f4e4b941aacf16a897a2cf689ed49b37ef
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 c014c5f4e4b941aacf16a897a2cf689ed49b37ef (commit) from 535d2509c2974c113e6c1bfa27ece5d5b8f99bd1 (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 c014c5f4e4b941aacf16a897a2cf689ed49b37ef Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Sun Mar 6 21:47:26 2011 -0500 Oops. Can't use #q yet,so change #q1 to (make-qd-d 1d0). diff --git a/qd-methods.lisp b/qd-methods.lisp index 3af7bd0..205c9fe 100644 --- a/qd-methods.lisp +++ b/qd-methods.lisp @@ -1105,7 +1105,7 @@ underlying floating-point format" ;; 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)) + (scale-float (make-qd-d 1d0) -212)) (defmethod epsilon ((m qd-complex)) (epsilon (realpart m))) ----------------------------------------------------------------------- Summary of changes: qd-methods.lisp | 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. 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