Changeset e2d8d63c0c06c474f32b76ebbb7e4cde44aac736

Show
Ignore:
Timestamp:
03/11/11 19:57:40 (3 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:
1 modified

Legend:

Unmodified
Added
Removed
  • qd-elliptic.lisp

    r5f6c62 re2d8d6  
    7373           (double-float 'double-float) 
    7474           (qd-real 'qd-real)))) 
    75     (if complexp 
    76         (if (eq max-type 'qd-real) 
    77             'qd-complex 
    78             `(cl:complex ,max-type)) 
    79         max-type))) 
    80    
     75    max-type)) 
     76 
     77(defun apply-contagion (number precision) 
     78  (etypecase number 
     79    ((or cl:real qd-real) 
     80     (coerce number precision)) 
     81    ((or cl:complex qd-complex) 
     82     (complex (coerce (realpart number) precision) 
     83              (coerce (imagpart number) precision))))) 
     84 
    8185;;; Jacobian elliptic functions 
    8286 
     
    283287 
    284288  Rf(x, y, z) = 1/2*integrate((t+x)^(-1/2)*(t+y)^(-1/2)*(t+z)^(-1/2), t, 0, inf)" 
    285   (let* ((xn x) 
    286          (yn y) 
    287          (zn z) 
     289  (let* ((precision (float-contagion x y z)) 
     290         (xn (apply-contagion x precision)) 
     291         (yn (apply-contagion y precision)) 
     292         (zn (apply-contagion z precision)) 
    288293         (a (/ (+ xn yn zn) 3)) 
    289294         (epslon (/ (max (abs (- a xn)) 
     
    334339 
    335340  Rd(x,y,z) = integrate(3/2*(t+x)^(-1/2)*(t+y)^(-1/2)*(t+z)^(-3/2), t, 0, inf)" 
    336   (let* ((xn x) 
    337          (yn y) 
    338          (zn z) 
     341  (let* ((precision (float-contagion x y z)) 
     342         (xn (apply-contagion x precision)) 
     343         (yn (apply-contagion y precision)) 
     344         (zn (apply-contagion z precision)) 
    339345         (a (/ (+ xn yn (* 3 zn)) 5)) 
    340346         (epslon (/ (max (abs (- a xn)) 
     
    349355    (loop while (> (* power4 epslon) (abs an)) 
    350356       do 
    351         (setf xnroot (sqrt xn)) 
    352         (setf ynroot (sqrt yn)) 
    353         (setf znroot (sqrt zn)) 
    354         (setf lam (+ (* xnroot ynroot) 
    355                       (* xnroot znroot) 
    356                       (* ynroot znroot))) 
    357         (setf sigma (+ sigma (/ power4 
    358                                 (* znroot (+ zn lam))))) 
    359         (setf power4 (* power4 1/4)) 
    360         (setf xn (* (+ xn lam) 1/4)) 
    361         (setf yn (* (+ yn lam) 1/4)) 
    362         (setf zn (* (+ zn lam) 1/4)) 
    363         (setf an (* (+ an lam) 1/4)) 
    364         (incf n)) 
     357      (setf xnroot (sqrt xn)) 
     358      (setf ynroot (sqrt yn)) 
     359      (setf znroot (sqrt zn)) 
     360      (setf lam (+ (* xnroot ynroot) 
     361                    (* xnroot znroot) 
     362                    (* ynroot znroot))) 
     363      (setf sigma (+ sigma (/ power4 
     364                              (* znroot (+ zn lam))))) 
     365      (setf power4 (* power4 1/4)) 
     366      (setf xn (* (+ xn lam) 1/4)) 
     367      (setf yn (* (+ yn lam) 1/4)) 
     368      (setf zn (* (+ zn lam) 1/4)) 
     369      (setf an (* (+ an lam) 1/4)) 
     370      (incf n)) 
    365371    ;; c1=-3/14,c2=1/6,c3=9/88,c4=9/22,c5=-3/22,c6=-9/52,c7=3/26 
    366372    (let* ((xndev (/ (* (- a x) power4) an)) 
     
    385391                 (* 45/272 ee2 ee2 ee3) 
    386392                 (* -9/68 (+ (* ee2 ee5) (* ee3 ee4)))))) 
    387     (+ (* 3 sigma) 
    388       (/ (* power4 s) 
    389           (expt an 3/2)))))) 
     393      (+ (* 3 sigma) 
     394        (/ (* power4 s) 
     395            (expt an 3/2)))))) 
    390396 
    391397;; Complete elliptic integral of the first kind.  This can be computed 
     
    404410        (t 
    405411         (let ((precision (float-contagion m))) 
    406            (carlson-rf (coerce 0 precision) (- 1 m) (coerce 1 precision)))))) 
     412           (carlson-rf 0 (- 1 m) 1))))) 
    407413 
    408414;; Elliptic integral of the first kind.  This is computed using 
     
    417423  Note for the complete elliptic integral, you can use elliptic-k" 
    418424  (let* ((precision (float-contagion x m)) 
    419          (x (coerce x precision)) 
    420          (m (coerce m precision))) 
     425         (x (apply-contagion x precision)) 
     426         (m (apply-contagion m precision))) 
    421427    (cond ((and (realp m) (realp x)) 
    422428           (cond ((> m 1) 
     
    426432                  ;; 
    427433                  ;; with sin(theta) = sqrt(m)*sin(phi) 
    428                   (/ (elliptic-f (cl:asin (* (sqrt m) (sin x))) (/ m)) 
     434                  (/ (elliptic-f (asin (* (sqrt m) (sin x))) (/ m)) 
    429435                     (sqrt m))) 
    430436                 ((< m 0) 
     
    446452                  ;; F(phi,1) = log(sec(phi)+tan(phi)) 
    447453                  ;;          = log(tan(pi/4+pi/2)) 
    448                   (log (cl:tan (+ (/ x 2) (/ (float-pi x) 4))))) 
     454                  (log (tan (+ (/ x 2) (/ (float-pi x) 4))))) 
    449455                 ((minusp x) 
    450456                  (- (elliptic-f (- x) m))) 
     
    463469                                   (* (- 1 (* k sin-x)) 
    464470                                      (+ 1 (* k sin-x))) 
    465                                    1.0)))) 
     471                                   1)))) 
    466472                 ((< x (float-pi x)) 
    467473                  (+ (* 2 (elliptic-k m)) 
     
    486492E(phi, m) = integrate(sqrt(1-m*sin(x)^2), x, 0, phi)" 
    487493  (let* ((precision (float-contagion phi m)) 
    488          (phi (coerce phi precision)) 
    489          (m (coerce m precision))) 
     494         (phi (apply-contagion phi precision)) 
     495         (m (apply-contagion m precision))) 
    490496    (cond ((= m 0) 
    491497           ;; A&S 17.4.23 
     
    501507                        (+ 1 (* k sin-phi))))) 
    502508             (- (* sin-phi 
    503                    (carlson-rf (* cos-phi cos-phi) y (coerce 1 precision))) 
     509                   (carlson-rf (* cos-phi cos-phi) y 1)) 
    504510                (* (/ m 3) 
    505511                   (expt sin-phi 3) 
    506                    (carlson-rd (* cos-phi cos-phi) y (coerce 1 precision))))))))) 
     512                   (carlson-rd (* cos-phi cos-phi) y 1)))))))) 
    507513 
    508514;; Complete elliptic integral of second kind. 
     
    514520 
    515521E(m) = integrate(sqrt(1-m*sin(x)^2), x, 0, %pi/2)" 
    516   (let ((precision (float-contagion m))) 
    517     (cond ((= m 0) 
    518            ;; A&S 17.4.23 
    519            (/ (float-pi m) 2)) 
    520           ((= m 1) 
    521            ;; A&S 17.4.25 
    522            (coerce 1 precision)) 
    523           (t 
    524            (let* ((k (sqrt m)) 
    525                   (y (* (- 1 k) 
    526                         (+ 1 k)))) 
    527              (- (carlson-rf 0.0 y 1.0) 
    528                 (* (/ m 3) 
    529                    (carlson-rd 0.0 y 1.0)))))))) 
     522  (cond ((= m 0) 
     523         ;; A&S 17.4.23 
     524         (/ (float-pi m) 2)) 
     525        ((= m 1) 
     526         ;; A&S 17.4.25 
     527         (float 1 m)) 
     528        (t 
     529         (let* ((k (sqrt m)) 
     530                (y (* (- 1 k) 
     531                      (+ 1 k)))) 
     532           (- (carlson-rf 0 y 1) 
     533              (* (/ m 3) 
     534                 (carlson-rd 0 y 1)))))))