Changeset e2d8d63c0c06c474f32b76ebbb7e4cde44aac736

Timestamp:

03/11/11 19:57:40 (2 years ago)

Author:

Raymond Toy <toy.raymond@…>

Children:

147fa2c7aa5e99988a7c3fb35800865d22efbc51

Parents:

5c21f133b0ebb511c664ab9fd967732cca6b76ea

git-committer:

Raymond Toy <toy.raymond@…> (03/11/11 19:57:40)

Message:

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.

Files:

• qd-elliptic.lisp (modified) (13 diffs)

qd-elliptic.lisp

@@ -73 +73 @@
            (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
 
@@ -283 +287 @@
 
   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))
@@ -334 +339 @@
 
   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))
@@ -349 +355 @@
     (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))
@@ -385 +391 @@
                  (* 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
@@ -404 +410 @@
         (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
@@ -417 +423 @@
   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)
@@ -426 +432 @@
                   ;;
                   ;; 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)
@@ -446 +452 @@
                   ;; 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)))
@@ -463 +469 @@
                                    (* (- 1 (* k sin-x))
                                       (+ 1 (* k sin-x)))
-                                   1.0))))
+                                   1))))
                  ((< x (float-pi x))
                   (+ (* 2 (elliptic-k m))
@@ -486 +492 @@
 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
@@ -501 +507 @@
                         (+ 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.
@@ -514 +520 @@
 
 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)))))))