Changeset e2d8d6


Ignore:
Timestamp:
03/12/11 03:57:40 (4 years ago)
Author:
Raymond Toy <toy.raymond@…>
Branches:
master
Children:
147fa2
Parents:
5c21f1
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.

File:
1 edited

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
     
    284288
    285289  Rf(x, y, z) = 1/2*integrate((t+x)^(-1/2)*(t+y)^(-1/2)*(t+z)^(-1/2), t, 0, inf)"
    286   (let* ((xn x)
    287          (yn y)
    288          (zn z)
     290  (let* ((precision (float-contagion x y z))
     291         (xn (apply-contagion x precision))
     292         (yn (apply-contagion y precision))
     293         (zn (apply-contagion z precision))
    289294         (a (/ (+ xn yn zn) 3))
    290295         (epslon (/ (max (abs (- a xn))
     
    335340
    336341  Rd(x,y,z) = integrate(3/2*(t+x)^(-1/2)*(t+y)^(-1/2)*(t+z)^(-3/2), t, 0, inf)"
    337   (let* ((xn x)
    338          (yn y)
    339          (zn z)
     342  (let* ((precision (float-contagion x y z))
     343         (xn (apply-contagion x precision))
     344         (yn (apply-contagion y precision))
     345         (zn (apply-contagion z precision))
    340346         (a (/ (+ xn yn (* 3 zn)) 5))
    341347         (epslon (/ (max (abs (- a xn))
     
    350356    (loop while (> (* power4 epslon) (abs an))
    351357       do
    352         (setf xnroot (sqrt xn))
    353         (setf ynroot (sqrt yn))
    354         (setf znroot (sqrt zn))
    355         (setf lam (+ (* xnroot ynroot)
    356                       (* xnroot znroot)
    357                       (* ynroot znroot)))
    358         (setf sigma (+ sigma (/ power4
    359                                 (* znroot (+ zn lam)))))
    360         (setf power4 (* power4 1/4))
    361         (setf xn (* (+ xn lam) 1/4))
    362         (setf yn (* (+ yn lam) 1/4))
    363         (setf zn (* (+ zn lam) 1/4))
    364         (setf an (* (+ an lam) 1/4))
    365         (incf n))
     358      (setf xnroot (sqrt xn))
     359      (setf ynroot (sqrt yn))
     360      (setf znroot (sqrt zn))
     361      (setf lam (+ (* xnroot ynroot)
     362                    (* xnroot znroot)
     363                    (* ynroot znroot)))
     364      (setf sigma (+ sigma (/ power4
     365                              (* znroot (+ zn lam)))))
     366      (setf power4 (* power4 1/4))
     367      (setf xn (* (+ xn lam) 1/4))
     368      (setf yn (* (+ yn lam) 1/4))
     369      (setf zn (* (+ zn lam) 1/4))
     370      (setf an (* (+ an lam) 1/4))
     371      (incf n))
    366372    ;; c1=-3/14,c2=1/6,c3=9/88,c4=9/22,c5=-3/22,c6=-9/52,c7=3/26
    367373    (let* ((xndev (/ (* (- a x) power4) an))
     
    386392                 (* 45/272 ee2 ee2 ee3)
    387393                 (* -9/68 (+ (* ee2 ee5) (* ee3 ee4))))))
    388     (+ (* 3 sigma)
    389       (/ (* power4 s)
    390           (expt an 3/2))))))
     394      (+ (* 3 sigma)
     395        (/ (* power4 s)
     396            (expt an 3/2))))))
    391397
    392398;; Complete elliptic integral of the first kind.  This can be computed
     
    405411        (t
    406412         (let ((precision (float-contagion m)))
    407            (carlson-rf (coerce 0 precision) (- 1 m) (coerce 1 precision))))))
     413           (carlson-rf 0 (- 1 m) 1)))))
    408414
    409415;; Elliptic integral of the first kind.  This is computed using
     
    418424  Note for the complete elliptic integral, you can use elliptic-k"
    419425  (let* ((precision (float-contagion x m))
    420          (x (coerce x precision))
    421          (m (coerce m precision)))
     426         (x (apply-contagion x precision))
     427         (m (apply-contagion m precision)))
    422428    (cond ((and (realp m) (realp x))
    423429           (cond ((> m 1)
     
    427433                  ;;
    428434                  ;; with sin(theta) = sqrt(m)*sin(phi)
    429                   (/ (elliptic-f (cl:asin (* (sqrt m) (sin x))) (/ m))
     435                  (/ (elliptic-f (asin (* (sqrt m) (sin x))) (/ m))
    430436                     (sqrt m)))
    431437                 ((< m 0)
     
    447453                  ;; F(phi,1) = log(sec(phi)+tan(phi))
    448454                  ;;          = log(tan(pi/4+pi/2))
    449                   (log (cl:tan (+ (/ x 2) (/ (float-pi x) 4)))))
     455                  (log (tan (+ (/ x 2) (/ (float-pi x) 4)))))
    450456                 ((minusp x)
    451457                  (- (elliptic-f (- x) m)))
     
    464470                                   (* (- 1 (* k sin-x))
    465471                                      (+ 1 (* k sin-x)))
    466                                    1.0))))
     472                                   1))))
    467473                 ((< x (float-pi x))
    468474                  (+ (* 2 (elliptic-k m))
     
    487493E(phi, m) = integrate(sqrt(1-m*sin(x)^2), x, 0, phi)"
    488494  (let* ((precision (float-contagion phi m))
    489          (phi (coerce phi precision))
    490          (m (coerce m precision)))
     495         (phi (apply-contagion phi precision))
     496         (m (apply-contagion m precision)))
    491497    (cond ((= m 0)
    492498           ;; A&S 17.4.23
     
    502508                        (+ 1 (* k sin-phi)))))
    503509             (- (* sin-phi
    504                    (carlson-rf (* cos-phi cos-phi) y (coerce 1 precision)))
     510                   (carlson-rf (* cos-phi cos-phi) y 1))
    505511                (* (/ m 3)
    506512                   (expt sin-phi 3)
    507                    (carlson-rd (* cos-phi cos-phi) y (coerce 1 precision)))))))))
     513                   (carlson-rd (* cos-phi cos-phi) y 1))))))))
    508514
    509515;; Complete elliptic integral of second kind.
     
    515521
    516522E(m) = integrate(sqrt(1-m*sin(x)^2), x, 0, %pi/2)"
    517   (let ((precision (float-contagion m)))
    518     (cond ((= m 0)
    519            ;; A&S 17.4.23
    520            (/ (float-pi m) 2))
    521           ((= m 1)
    522            ;; A&S 17.4.25
    523            (coerce 1 precision))
    524           (t
    525            (let* ((k (sqrt m))
    526                   (y (* (- 1 k)
    527                         (+ 1 k))))
    528              (- (carlson-rf 0.0 y 1.0)
    529                 (* (/ m 3)
    530                    (carlson-rd 0.0 y 1.0))))))))
     523  (cond ((= m 0)
     524         ;; A&S 17.4.23
     525         (/ (float-pi m) 2))
     526        ((= m 1)
     527         ;; A&S 17.4.25
     528         (float 1 m))
     529        (t
     530         (let* ((k (sqrt m))
     531                (y (* (- 1 k)
     532                      (+ 1 k))))
     533           (- (carlson-rf 0 y 1)
     534              (* (/ m 3)
     535                 (carlson-rd 0 y 1)))))))
Note: See TracChangeset for help on using the changeset viewer.