Changeset e2d8d63c0c06c474f32b76ebbb7e4cde44aac736
 Timestamp:
 03/11/11 19:57:40 (3 years ago)
 Author:
 Raymond Toy <toy.raymond@…>
 Children:
 147fa2c7aa5e99988a7c3fb35800865d22efbc51
 Parents:
 5c21f133b0ebb511c664ab9fd967732cca6b76ea
 gitcommitter:
 Raymond Toy <toy.raymond@…> (03/11/11 19:57:40)
 Message:

Clean up floatcontagion stuff; use it in Carlson routines.
o FLOATCONTAGION now only returns the real type, not a complex type.
o Add APPLYCONTAGION 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 APPLYCONTAGION.
o Use the contagion stuff in CARLSONRD and CARLSONRF.
 Files:

Legend:
 Unmodified
 Added
 Removed

r5f6c62

re2d8d6


73  73  (doublefloat 'doublefloat) 
74  74  (qdreal 'qdreal)))) 
75   (if complexp 
76   (if (eq maxtype 'qdreal) 
77   'qdcomplex 
78   `(cl:complex ,maxtype)) 
79   maxtype))) 
80   
 75  maxtype)) 
 76  
 77  (defun applycontagion (number precision) 
 78  (etypecase number 
 79  ((or cl:real qdreal) 
 80  (coerce number precision)) 
 81  ((or cl:complex qdcomplex) 
 82  (complex (coerce (realpart number) precision) 
 83  (coerce (imagpart number) precision))))) 
 84  
81  85  ;;; Jacobian elliptic functions 
82  86  
… 
… 

283  287  
284  288  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 (floatcontagion x y z)) 
 290  (xn (applycontagion x precision)) 
 291  (yn (applycontagion y precision)) 
 292  (zn (applycontagion z precision)) 
288  293  (a (/ (+ xn yn zn) 3)) 
289  294  (epslon (/ (max (abs ( a xn)) 
… 
… 

334  339  
335  340  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 (floatcontagion x y z)) 
 342  (xn (applycontagion x precision)) 
 343  (yn (applycontagion y precision)) 
 344  (zn (applycontagion z precision)) 
339  345  (a (/ (+ xn yn (* 3 zn)) 5)) 
340  346  (epslon (/ (max (abs ( a xn)) 
… 
… 

349  355  (loop while (> (* power4 epslon) (abs an)) 
350  356  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)) 
365  371  ;; c1=3/14,c2=1/6,c3=9/88,c4=9/22,c5=3/22,c6=9/52,c7=3/26 
366  372  (let* ((xndev (/ (* ( a x) power4) an)) 
… 
… 

385  391  (* 45/272 ee2 ee2 ee3) 
386  392  (* 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)))))) 
390  396  
391  397  ;; Complete elliptic integral of the first kind. This can be computed 
… 
… 

404  410  (t 
405  411  (let ((precision (floatcontagion m))) 
406   (carlsonrf (coerce 0 precision) ( 1 m) (coerce 1 precision)))))) 
 412  (carlsonrf 0 ( 1 m) 1))))) 
407  413  
408  414  ;; Elliptic integral of the first kind. This is computed using 
… 
… 

417  423  Note for the complete elliptic integral, you can use elliptick" 
418  424  (let* ((precision (floatcontagion x m)) 
419   (x (coerce x precision)) 
420   (m (coerce m precision))) 
 425  (x (applycontagion x precision)) 
 426  (m (applycontagion m precision))) 
421  427  (cond ((and (realp m) (realp x)) 
422  428  (cond ((> m 1) 
… 
… 

426  432  ;; 
427  433  ;; with sin(theta) = sqrt(m)*sin(phi) 
428   (/ (ellipticf (cl:asin (* (sqrt m) (sin x))) (/ m)) 
 434  (/ (ellipticf (asin (* (sqrt m) (sin x))) (/ m)) 
429  435  (sqrt m))) 
430  436  ((< m 0) 
… 
… 

446  452  ;; F(phi,1) = log(sec(phi)+tan(phi)) 
447  453  ;; = log(tan(pi/4+pi/2)) 
448   (log (cl:tan (+ (/ x 2) (/ (floatpi x) 4))))) 
 454  (log (tan (+ (/ x 2) (/ (floatpi x) 4))))) 
449  455  ((minusp x) 
450  456  ( (ellipticf ( x) m))) 
… 
… 

463  469  (* ( 1 (* k sinx)) 
464  470  (+ 1 (* k sinx))) 
465   1.0)))) 
 471  1)))) 
466  472  ((< x (floatpi x)) 
467  473  (+ (* 2 (elliptick m)) 
… 
… 

486  492  E(phi, m) = integrate(sqrt(1m*sin(x)^2), x, 0, phi)" 
487  493  (let* ((precision (floatcontagion phi m)) 
488   (phi (coerce phi precision)) 
489   (m (coerce m precision))) 
 494  (phi (applycontagion phi precision)) 
 495  (m (applycontagion m precision))) 
490  496  (cond ((= m 0) 
491  497  ;; A&S 17.4.23 
… 
… 

501  507  (+ 1 (* k sinphi))))) 
502  508  ( (* sinphi 
503   (carlsonrf (* cosphi cosphi) y (coerce 1 precision))) 
 509  (carlsonrf (* cosphi cosphi) y 1)) 
504  510  (* (/ m 3) 
505  511  (expt sinphi 3) 
506   (carlsonrd (* cosphi cosphi) y (coerce 1 precision))))))))) 
 512  (carlsonrd (* cosphi cosphi) y 1)))))))) 
507  513  
508  514  ;; Complete elliptic integral of second kind. 
… 
… 

514  520  
515  521  E(m) = integrate(sqrt(1m*sin(x)^2), x, 0, %pi/2)" 
516   (let ((precision (floatcontagion m))) 
517   (cond ((= m 0) 
518   ;; A&S 17.4.23 
519   (/ (floatpi 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   ( (carlsonrf 0.0 y 1.0) 
528   (* (/ m 3) 
529   (carlsonrd 0.0 y 1.0)))))))) 
 522  (cond ((= m 0) 
 523  ;; A&S 17.4.23 
 524  (/ (floatpi 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  ( (carlsonrf 0 y 1) 
 533  (* (/ m 3) 
 534  (carlsonrd 0 y 1))))))) 