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 5bd5df9360268fae90fc41fcdf86b728f8a54e86 (commit) from 390a7483f5658fe802d5d239070cdaa573adf4a5 (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 5bd5df9360268fae90fc41fcdf86b728f8a54e86 Author: Raymond Toy <toy.raymond@gmail.com> Date: Sat Mar 12 18:21:55 2011 -0500 Fix possible bug in elliptic-pi; add comments. qd-elliptic.lisp: o Add some comments o Fix a possible bug if n is a complex number or a negative number. rt-tests.lisp: o Remove one broken test. o Fix the other tests for elliptic-pi and adjust required precision down a bit so the tests can pass. diff --git a/qd-elliptic.lisp b/qd-elliptic.lisp index 1b80f13..eafca11 100644 --- a/qd-elliptic.lisp +++ b/qd-elliptic.lisp @@ -694,12 +694,18 @@ E(m) = integrate(sqrt(1-m*sin(x)^2), x, 0, %pi/2)" ;; ;; PI(n; phi|m) = integrate(1/sqrt(1-m*sin(x)^2)/(1-n*sin(x)^2), x, 0, phi) ;; +;; +;; Carlson writes +;; +;; P(phi,k,n) = integrate((1+n*sin(t)^2)^(-1)*(1-k^2*sin(t)^2)^(-1/2), t, 0, phi) +;; = sin(phi)*Rf(cos(phi)^2, 1-k^2*sin(phi)^2, 1) +;; - n/3*sin(phi)^3*Rj(cos(phi)^2, 1-k^2*sin(phi)^2, 1, 1+n*sin(phi)^2) +;; +;; Note that this definition as a different sign for the n parameter from A&S! (defun elliptic-pi (n phi m) "Compute elliptic integral of the third kind: PI(n; phi|m) = integrate(1/sqrt(1-m*sin(x)^2)/(1-n*sin(x)^2), x, 0, phi)" - ;; Note: Carlson's DRJ has n defined as the negative of the n given - ;; in A&S. (let* ((precision (float-contagion n phi m)) (n (apply-contagion n precision)) (phi (apply-contagion phi precision)) @@ -707,10 +713,8 @@ E(m) = integrate(sqrt(1-m*sin(x)^2), x, 0, %pi/2)" (nn (- n)) (sin-phi (sin phi)) (cos-phi (cos phi)) - (k (sqrt m)) - (k2sin (* (- 1 (* k sin-phi)) - (+ 1 (* k sin-phi))))) - (- (* sin-phi (carlson-rf (expt cos-phi 2) k2sin 1)) + (m-sin2 (- 1 (* m sin-phi sin-phi))) + (- (* sin-phi (carlson-rf (expt cos-phi 2) m-sin2 1)) (* (/ nn 3) (expt sin-phi 3) - (carlson-rj (expt cos-phi 2) k2sin 1 + (carlson-rj (expt cos-phi 2) m-sin2 1 (+ 1 (* nn (expt sin-phi 2)))))))) diff --git a/rt-tests.lisp b/rt-tests.lisp index 80d5c57..ab81f6a 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -978,32 +978,26 @@ append (list (list (list k n) result))) nil) -#+nil -(rt:deftest oct.elliptic-pi.19.6.2.d - (loop for k from 0 below 100 - for n = (random 1d0) - for epi = (elliptic-pi (- n) (/ (float-pi n) 2) n) - for true = (+ (/ (float-pi n) 4 (sqrt (+ 1 (sqrt n)))) - (/ (elliptic-k n) 2)) - for result = (check-accuracy 53 epi true) - when result - append (list (list (list k n) result))) - nil) - - -#|| ;; elliptic-pi(n, phi, 0) = -;; atanh(sqrt(1-n)*tan(phi))/sqrt(1-n) n < 1 +;; atan(sqrt(1-n)*tan(phi))/sqrt(1-n) n < 1 ;; atanh(sqrt(n-1)*tan(phi))/sqrt(n-1) n > 1 ;; tan(phi) n = 1 +;; +;; These are easy to derive if you look at the integral: +;; +;; ellipti-pi(n, phi, 0) = integrate(1/(1-n*sin(t)^2), t, 0, phi) +;; +;; and this can be easily integrated to give the above expressions for +;; the different values of n. (rt:deftest oct.elliptic-pi.n0.d + ;; Tests for random values for phi in [0, pi/2] and n in [0, 1] (loop for k from 0 below 100 for phi = (random (/ pi 2)) for n = (random 1d0) for epi = (elliptic-pi n phi 0) - for true = (/ (atanh (* (tan phi) (sqrt (- 1 n)))) + for true = (/ (atan (* (tan phi) (sqrt (- 1 n)))) (sqrt (- 1 n))) - for result = (check-accuracy 53 epi true) + for result = (check-accuracy 50 epi true) unless (eq nil result) append (list (list (list k n phi) result))) nil) @@ -1011,9 +1005,9 @@ (rt:deftest oct.elliptic-pi.n1.d (loop for k from 0 below 100 for phi = (random (/ pi 2)) - for epi = (elliptic-pi 0 phi 0) + for epi = (elliptic-pi 1 phi 0) for true = (tan phi) - for result = (check-accuracy 53 epi true) + for result = (check-accuracy 43 epi true) unless (eq nil result) append (list (list (list k phi) result))) nil) @@ -1025,7 +1019,7 @@ for epi = (elliptic-pi n phi 0) for true = (/ (atanh (* (tan phi) (sqrt (- n 1)))) (sqrt (- n 1))) - for result = (check-accuracy 52 epi true) + for result = (check-accuracy 49 epi true) ;; Not sure if this formula holds when atanh gives a complex ;; result. Wolfram doesn't say when (and (not (complexp true)) result) @@ -1033,13 +1027,14 @@ nil) (rt:deftest oct.elliptic-pi.n0.q + ;; Tests for random values for phi in [0, pi/2] and n in [0, 1] (loop for k from 0 below 100 for phi = (random (/ +pi+ 2)) for n = (random #q1) for epi = (elliptic-pi n phi 0) - for true = (/ (atanh (* (tan phi) (sqrt (- 1 n)))) + for true = (/ (atan (* (tan phi) (sqrt (- 1 n)))) (sqrt (- 1 n))) - for result = (check-accuracy 212 epi true) + for result = (check-accuracy 208 epi true) unless (eq nil result) append (list (list (list k n phi) result))) nil) @@ -1047,9 +1042,9 @@ (rt:deftest oct.elliptic-pi.n1.q (loop for k from 0 below 100 for phi = (random (/ +pi+ 2)) - for epi = (elliptic-pi 0 phi 0) + for epi = (elliptic-pi 1 phi 0) for true = (tan phi) - for result = (check-accuracy 212 epi true) + for result = (check-accuracy 205 epi true) unless (eq nil result) append (list (list (list k phi) result))) nil) @@ -1061,10 +1056,9 @@ for epi = (elliptic-pi n phi 0) for true = (/ (atanh (* (tan phi) (sqrt (- n 1)))) (sqrt (- n 1))) - for result = (check-accuracy 209 epi true) + for result = (check-accuracy 208 epi true) ;; Not sure if this formula holds when atanh gives a complex ;; result. Wolfram doesn't say when (and (not (complexp true)) result) append (list (list (list k n phi) result))) nil) -||# ----------------------------------------------------------------------- Summary of changes: qd-elliptic.lisp | 18 +++++++++++------- rt-tests.lisp | 46 ++++++++++++++++++++-------------------------- 2 files changed, 31 insertions(+), 33 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats