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@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@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