[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 147fa2c7aa5e99988a7c3fb35800865d22efbc51
Raymond Toy
rtoy at common-lisp.net
Sat Mar 12 04:30:15 UTC 2011
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 147fa2c7aa5e99988a7c3fb35800865d22efbc51 (commit)
via e2d8d63c0c06c474f32b76ebbb7e4cde44aac736 (commit)
from 5c21f133b0ebb511c664ab9fd967732cca6b76ea (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 147fa2c7aa5e99988a7c3fb35800865d22efbc51
Author: Raymond Toy <toy.raymond at gmail.com>
Date: Fri Mar 11 23:11:32 2011 -0500
Add tests for contagion support.
diff --git a/rt-tests.lisp b/rt-tests.lisp
index 85f9b95..7635611 100644
--- a/rt-tests.lisp
+++ b/rt-tests.lisp
@@ -864,3 +864,68 @@
(true #q1.797210352103388311159883738420485817340818994823477337395512429419599q0))
(check-accuracy 212 rd true))
nil)
+
+;; Test some of the contagion stuff.
+
+(rt:deftest oct.carlson-rf.contagion.1
+ ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1)
+ ;; = 1/4*beta(1/2,1/2)
+ ;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4)
+ (let ((rf (carlson-rf 0 2 1))
+ (true 1.31102877714605990523241979494d0))
+ (check-accuracy 23 rf true))
+ nil)
+
+(rt:deftest oct.carlson-rf.contagion.1d
+ ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1)
+ ;; = 1/4*beta(1/2,1/2)
+ ;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4)
+ (let ((rf (carlson-rf 0d0 2 1))
+ (true 1.31102877714605990523241979494d0))
+ (check-accuracy 53 rf true))
+ nil)
+
+(rt:deftest oct.carlson-rf.contagion.2d
+ ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1)
+ ;; = 1/4*beta(1/2,1/2)
+ ;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4)
+ (let ((rf (carlson-rf 0 2d0 1))
+ (true 1.31102877714605990523241979494d0))
+ (check-accuracy 53 rf true))
+ nil)
+
+(rt:deftest oct.carlson-rf.contagion.3d
+ ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1)
+ ;; = 1/4*beta(1/2,1/2)
+ ;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4)
+ (let ((rf (carlson-rf 0 2 1d0))
+ (true 1.31102877714605990523241979494d0))
+ (check-accuracy 53 rf true))
+ nil)
+
+(rt:deftest oct.carlson-rf.contagion.1q
+ ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1)
+ ;; = 1/4*beta(1/2,1/2)
+ ;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4)
+ (let ((rf (carlson-rf #q0q0 2 1))
+ (true #q1.311028777146059905232419794945559706841377475715811581408410851900395q0))
+ (check-accuracy 212 rf true))
+ nil)
+
+(rt:deftest oct.carlson-rf.contagion.2q
+ ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1)
+ ;; = 1/4*beta(1/2,1/2)
+ ;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4)
+ (let ((rf (carlson-rf 0 #q2q0 1))
+ (true #q1.311028777146059905232419794945559706841377475715811581408410851900395q0))
+ (check-accuracy 212 rf true))
+ nil)
+
+(rt:deftest oct.carlson-rf.contagion.3q
+ ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1)
+ ;; = 1/4*beta(1/2,1/2)
+ ;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4)
+ (let ((rf (carlson-rf 0 2 #q1q0))
+ (true #q1.311028777146059905232419794945559706841377475715811581408410851900395q0))
+ (check-accuracy 212 rf true))
+ nil)
\ No newline at end of file
commit e2d8d63c0c06c474f32b76ebbb7e4cde44aac736
Author: Raymond Toy <toy.raymond at gmail.com>
Date: Fri Mar 11 22:57:40 2011 -0500
Clean up float-contagion stuff; use it in Carlson routines.
o FLOAT-CONTAGION now only returns the real type, not a complex type.
o Add APPLY-CONTAGION to make the specified conversion. This handle
complex numbers and makes the components have the specified
precision.
o Change uses of contagion stuff to use APPLY-CONTAGION.
o Use the contagion stuff in CARLSON-RD and CARLSON-RF.
diff --git a/qd-elliptic.lisp b/qd-elliptic.lisp
index 6c8070a..342290e 100644
--- a/qd-elliptic.lisp
+++ b/qd-elliptic.lisp
@@ -72,12 +72,16 @@
(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)))
-
+ max-type))
+
+(defun apply-contagion (number precision)
+ (etypecase number
+ ((or cl:real qd-real)
+ (coerce number precision))
+ ((or cl:complex qd-complex)
+ (complex (coerce (realpart number) precision)
+ (coerce (imagpart number) precision)))))
+
;;; Jacobian elliptic functions
(defun ascending-transform (u m)
@@ -282,9 +286,10 @@
"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)
+ (let* ((precision (float-contagion x y z))
+ (xn (apply-contagion x precision))
+ (yn (apply-contagion y precision))
+ (zn (apply-contagion z precision))
(a (/ (+ xn yn zn) 3))
(epslon (/ (max (abs (- a xn))
(abs (- a yn))
@@ -333,9 +338,10 @@
"Compute Carlson's Rd function:
Rd(x,y,z) = integrate(3/2*(t+x)^(-1/2)*(t+y)^(-1/2)*(t+z)^(-3/2), t, 0, inf)"
- (let* ((xn x)
- (yn y)
- (zn z)
+ (let* ((precision (float-contagion x y z))
+ (xn (apply-contagion x precision))
+ (yn (apply-contagion y precision))
+ (zn (apply-contagion z precision))
(a (/ (+ xn yn (* 3 zn)) 5))
(epslon (/ (max (abs (- a xn))
(abs (- a yn))
@@ -348,20 +354,20 @@
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 sigma (+ sigma (/ power4
- (* znroot (+ zn lam)))))
- (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))
+ (setf xnroot (sqrt xn))
+ (setf ynroot (sqrt yn))
+ (setf znroot (sqrt zn))
+ (setf lam (+ (* xnroot ynroot)
+ (* xnroot znroot)
+ (* ynroot znroot)))
+ (setf sigma (+ sigma (/ power4
+ (* znroot (+ zn lam)))))
+ (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))
@@ -384,9 +390,9 @@
(* 3/20 ee2 ee4)
(* 45/272 ee2 ee2 ee3)
(* -9/68 (+ (* ee2 ee5) (* ee3 ee4))))))
- (+ (* 3 sigma)
- (/ (* power4 s)
- (expt an 3/2))))))
+ (+ (* 3 sigma)
+ (/ (* power4 s)
+ (expt an 3/2))))))
;; Complete elliptic integral of the first kind. This can be computed
;; from Carlson's Rf function:
@@ -403,7 +409,7 @@
(/ (float +pi+ m) 2))
(t
(let ((precision (float-contagion m)))
- (carlson-rf (coerce 0 precision) (- 1 m) (coerce 1 precision))))))
+ (carlson-rf 0 (- 1 m) 1)))))
;; Elliptic integral of the first kind. This is computed using
;; Carlson's Rf function:
@@ -416,8 +422,8 @@
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)))
+ (x (apply-contagion x precision))
+ (m (apply-contagion m precision)))
(cond ((and (realp m) (realp x))
(cond ((> m 1)
;; A&S 17.4.15
@@ -425,7 +431,7 @@
;; 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))
+ (/ (elliptic-f (asin (* (sqrt m) (sin x))) (/ m))
(sqrt m)))
((< m 0)
;; A&S 17.4.17
@@ -445,7 +451,7 @@
;;
;; F(phi,1) = log(sec(phi)+tan(phi))
;; = log(tan(pi/4+pi/2))
- (log (cl:tan (+ (/ x 2) (/ (float-pi x) 4)))))
+ (log (tan (+ (/ x 2) (/ (float-pi x) 4)))))
((minusp x)
(- (elliptic-f (- x) m)))
((> x (float-pi x))
@@ -462,7 +468,7 @@
(carlson-rf (* cos-x cos-x)
(* (- 1 (* k sin-x))
(+ 1 (* k sin-x)))
- 1.0))))
+ 1))))
((< x (float-pi x))
(+ (* 2 (elliptic-k m))
(elliptic-f (- x (float pi x)) m)))))
@@ -485,8 +491,8 @@
E(phi, m) = integrate(sqrt(1-m*sin(x)^2), x, 0, phi)"
(let* ((precision (float-contagion phi m))
- (phi (coerce phi precision))
- (m (coerce m precision)))
+ (phi (apply-contagion phi precision))
+ (m (apply-contagion m precision)))
(cond ((= m 0)
;; A&S 17.4.23
phi)
@@ -500,10 +506,10 @@ E(phi, m) = integrate(sqrt(1-m*sin(x)^2), x, 0, phi)"
(y (* (- 1 (* k sin-phi))
(+ 1 (* k sin-phi)))))
(- (* sin-phi
- (carlson-rf (* cos-phi cos-phi) y (coerce 1 precision)))
+ (carlson-rf (* cos-phi cos-phi) y 1))
(* (/ m 3)
(expt sin-phi 3)
- (carlson-rd (* cos-phi cos-phi) y (coerce 1 precision)))))))))
+ (carlson-rd (* cos-phi cos-phi) y 1))))))))
;; Complete elliptic integral of second kind.
;;
@@ -513,17 +519,16 @@ E(phi, m) = integrate(sqrt(1-m*sin(x)^2), x, 0, phi)"
"Complete elliptic integral of the second kind:
E(m) = integrate(sqrt(1-m*sin(x)^2), x, 0, %pi/2)"
- (let ((precision (float-contagion m)))
- (cond ((= m 0)
- ;; A&S 17.4.23
- (/ (float-pi m) 2))
- ((= m 1)
- ;; A&S 17.4.25
- (coerce 1 precision))
- (t
- (let* ((k (sqrt m))
- (y (* (- 1 k)
- (+ 1 k))))
- (- (carlson-rf 0.0 y 1.0)
- (* (/ m 3)
- (carlson-rd 0.0 y 1.0))))))))
+ (cond ((= m 0)
+ ;; A&S 17.4.23
+ (/ (float-pi m) 2))
+ ((= m 1)
+ ;; A&S 17.4.25
+ (float 1 m))
+ (t
+ (let* ((k (sqrt m))
+ (y (* (- 1 k)
+ (+ 1 k))))
+ (- (carlson-rf 0 y 1)
+ (* (/ m 3)
+ (carlson-rd 0 y 1)))))))
-----------------------------------------------------------------------
Summary of changes:
qd-elliptic.lisp | 111 ++++++++++++++++++++++++++++--------------------------
rt-tests.lisp | 65 +++++++++++++++++++++++++++++++
2 files changed, 123 insertions(+), 53 deletions(-)
hooks/post-receive
--
OCT: A portable Lisp implementation for quad-double precision floats
More information about the oct-cvs
mailing list