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