home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / numerical / slatec / zairy.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  10.0 KB  |  297 lines
  1. ;;; Compiled by f2cl version 2.0 beta 2002-05-06
  2. ;;; 
  3. ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
  4. ;;;           (:coerce-assigns :as-needed) (:array-type ':simple-array)
  5. ;;;           (:array-slicing nil) (:declare-common nil)
  6. ;;;           (:float-format double-float))
  7.  
  8. (in-package "SLATEC")
  9.  
  10.  
  11. (let ((tth 0.6666666666666667)
  12.       (c1 0.3550280538878172)
  13.       (c2 0.2588194037928068)
  14.       (coef 0.1837762984739307)
  15.       (zeror 0.0)
  16.       (zeroi 0.0)
  17.       (coner 1.0)
  18.       (conei 0.0))
  19.   (declare (type double-float conei coner zeroi zeror coef c2 c1 tth))
  20.   (defun zairy (zr zi id kode air aii nz ierr)
  21.     (declare (type double-float zr zi air aii)
  22.              (type f2cl-lib:integer4 id kode nz ierr))
  23.     (prog ((cyr (make-array 1 :element-type 'double-float))
  24.            (cyi (make-array 1 :element-type 'double-float)) (iflag 0) (k 0)
  25.            (k1 0) (k2 0) (mr 0) (nn 0) (aa 0.0) (ad 0.0) (ak 0.0) (alim 0.0)
  26.            (atrm 0.0) (az 0.0) (az3 0.0) (bk 0.0) (cc 0.0) (ck 0.0) (csqi 0.0)
  27.            (csqr 0.0) (dig 0.0) (dk 0.0) (d1 0.0) (d2 0.0) (elim 0.0) (fid 0.0)
  28.            (fnu 0.0) (ptr 0.0) (rl 0.0) (r1m5 0.0) (sfac 0.0) (sti 0.0)
  29.            (str 0.0) (s1i 0.0) (s1r 0.0) (s2i 0.0) (s2r 0.0) (tol 0.0)
  30.            (trm1i 0.0) (trm1r 0.0) (trm2i 0.0) (trm2r 0.0) (ztai 0.0)
  31.            (ztar 0.0) (z3i 0.0) (z3r 0.0) (alaz 0.0) (bb 0.0))
  32.       (declare (type (simple-array double-float (1)) cyr cyi)
  33.                (type double-float bb alaz z3r z3i ztar ztai trm2r trm2i trm1r
  34.                 trm1i tol s2r s2i s1r s1i str sti sfac r1m5 rl ptr fnu fid elim
  35.                 d2 d1 dk dig csqr csqi ck cc bk az3 az atrm alim ak ad aa)
  36.                (type f2cl-lib:integer4 nn mr k2 k1 k iflag))
  37.       (setf ierr 0)
  38.       (setf nz 0)
  39.       (if (or (< id 0) (> id 1)) (setf ierr 1))
  40.       (if (or (< kode 1) (> kode 2)) (setf ierr 1))
  41.       (if (/= ierr 0) (go end_label))
  42.       (setf az (zabs zr zi))
  43.       (setf tol (max (f2cl-lib:d1mach 4) 1.0e-18))
  44.       (setf fid (coerce (the f2cl-lib:integer4 id) 'double-float))
  45.       (if (> az 1.0) (go label70))
  46.       (setf s1r coner)
  47.       (setf s1i conei)
  48.       (setf s2r coner)
  49.       (setf s2i conei)
  50.       (if (< az tol) (go label170))
  51.       (setf aa (* az az))
  52.       (if (< aa (/ tol az)) (go label40))
  53.       (setf trm1r coner)
  54.       (setf trm1i conei)
  55.       (setf trm2r coner)
  56.       (setf trm2i conei)
  57.       (setf atrm 1.0)
  58.       (setf str (- (* zr zr) (* zi zi)))
  59.       (setf sti (+ (* zr zi) (* zi zr)))
  60.       (setf z3r (- (* str zr) (* sti zi)))
  61.       (setf z3i (+ (* str zi) (* sti zr)))
  62.       (setf az3 (* az aa))
  63.       (setf ak (+ 2.0 fid))
  64.       (setf bk (- 3.0 fid fid))
  65.       (setf ck (- 4.0 fid))
  66.       (setf dk (+ 3.0 fid fid))
  67.       (setf d1 (* ak dk))
  68.       (setf d2 (* bk ck))
  69.       (setf ad (min d1 d2))
  70.       (setf ak (+ 24.0 (* 9.0 fid)))
  71.       (setf bk (- 30.0 (* 9.0 fid)))
  72.       (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
  73.                     ((> k 25) nil)
  74.         (tagbody
  75.           (setf str (/ (- (* trm1r z3r) (* trm1i z3i)) d1))
  76.           (setf trm1i (/ (+ (* trm1r z3i) (* trm1i z3r)) d1))
  77.           (setf trm1r str)
  78.           (setf s1r (+ s1r trm1r))
  79.           (setf s1i (+ s1i trm1i))
  80.           (setf str (/ (- (* trm2r z3r) (* trm2i z3i)) d2))
  81.           (setf trm2i (/ (+ (* trm2r z3i) (* trm2i z3r)) d2))
  82.           (setf trm2r str)
  83.           (setf s2r (+ s2r trm2r))
  84.           (setf s2i (+ s2i trm2i))
  85.           (setf atrm (/ (* atrm az3) ad))
  86.           (setf d1 (+ d1 ak))
  87.           (setf d2 (+ d2 bk))
  88.           (setf ad (min d1 d2))
  89.           (if (< atrm (* tol ad)) (go label40))
  90.           (setf ak (+ ak 18.0))
  91.           (setf bk (+ bk 18.0))
  92.          label30))
  93.      label40
  94.       (if (= id 1) (go label50))
  95.       (setf air (- (* s1r c1) (* c2 (- (* zr s2r) (* zi s2i)))))
  96.       (setf aii (- (* s1i c1) (* c2 (+ (* zr s2i) (* zi s2r)))))
  97.       (if (= kode 1) (go end_label))
  98.       (multiple-value-bind
  99.           (var-0 var-1 var-2 var-3)
  100.           (zsqrt zr zi str sti)
  101.         (declare (ignore var-0 var-1))
  102.         (setf str var-2)
  103.         (setf sti var-3))
  104.       (setf ztar (* tth (- (* zr str) (* zi sti))))
  105.       (setf ztai (* tth (+ (* zr sti) (* zi str))))
  106.       (multiple-value-bind
  107.           (var-0 var-1 var-2 var-3)
  108.           (zexp ztar ztai str sti)
  109.         (declare (ignore var-0 var-1))
  110.         (setf str var-2)
  111.         (setf sti var-3))
  112.       (setf ptr (- (* air str) (* aii sti)))
  113.       (setf aii (+ (* air sti) (* aii str)))
  114.       (setf air ptr)
  115.       (go end_label)
  116.      label50
  117.       (setf air (* (- s2r) c2))
  118.       (setf aii (* (- s2i) c2))
  119.       (if (<= az tol) (go label60))
  120.       (setf str (- (* zr s1r) (* zi s1i)))
  121.       (setf sti (+ (* zr s1i) (* zi s1r)))
  122.       (setf cc (/ c1 (+ 1.0 fid)))
  123.       (setf air (+ air (* cc (- (* str zr) (* sti zi)))))
  124.       (setf aii (+ aii (* cc (+ (* str zi) (* sti zr)))))
  125.      label60
  126.       (if (= kode 1) (go end_label))
  127.       (multiple-value-bind
  128.           (var-0 var-1 var-2 var-3)
  129.           (zsqrt zr zi str sti)
  130.         (declare (ignore var-0 var-1))
  131.         (setf str var-2)
  132.         (setf sti var-3))
  133.       (setf ztar (* tth (- (* zr str) (* zi sti))))
  134.       (setf ztai (* tth (+ (* zr sti) (* zi str))))
  135.       (multiple-value-bind
  136.           (var-0 var-1 var-2 var-3)
  137.           (zexp ztar ztai str sti)
  138.         (declare (ignore var-0 var-1))
  139.         (setf str var-2)
  140.         (setf sti var-3))
  141.       (setf ptr (- (* str air) (* sti aii)))
  142.       (setf aii (+ (* str aii) (* sti air)))
  143.       (setf air ptr)
  144.       (go end_label)
  145.      label70
  146.       (setf fnu (/ (+ 1.0 fid) 3.0))
  147.       (setf k1 (f2cl-lib:i1mach 15))
  148.       (setf k2 (f2cl-lib:i1mach 16))
  149.       (setf r1m5 (f2cl-lib:d1mach 5))
  150.       (setf k (f2cl-lib:int (min (abs k1) (abs k2))))
  151.       (setf elim (* 2.303 (- (* k r1m5) 3.0)))
  152.       (setf k1 (f2cl-lib:int-sub (f2cl-lib:i1mach 14) 1))
  153.       (setf aa (* r1m5 k1))
  154.       (setf dig (min aa 18.0))
  155.       (setf aa (* aa 2.303))
  156.       (setf alim (+ elim (max (- aa) -41.45)))
  157.       (setf rl (+ (* 1.2 dig) 3.0))
  158.       (setf alaz (f2cl-lib:flog az))
  159.       (setf aa (/ 0.5 tol))
  160.       (setf bb (* (f2cl-lib:i1mach 9) 0.5))
  161.       (setf aa (min aa bb))
  162.       (setf aa (expt aa tth))
  163.       (if (> az aa) (go label260))
  164.       (setf aa (f2cl-lib:fsqrt aa))
  165.       (if (> az aa) (setf ierr 3))
  166.       (multiple-value-bind
  167.           (var-0 var-1 var-2 var-3)
  168.           (zsqrt zr zi csqr csqi)
  169.         (declare (ignore var-0 var-1))
  170.         (setf csqr var-2)
  171.         (setf csqi var-3))
  172.       (setf ztar (* tth (- (* zr csqr) (* zi csqi))))
  173.       (setf ztai (* tth (+ (* zr csqi) (* zi csqr))))
  174.       (setf iflag 0)
  175.       (setf sfac 1.0)
  176.       (setf ak ztai)
  177.       (if (>= zr 0.0) (go label80))
  178.       (setf bk ztar)
  179.       (setf ck (coerce (- (abs bk)) 'double-float))
  180.       (setf ztar ck)
  181.       (setf ztai ak)
  182.      label80
  183.       (if (/= zi 0.0) (go label90))
  184.       (if (> zr 0.0) (go label90))
  185.       (setf ztar 0.0)
  186.       (setf ztai ak)
  187.      label90
  188.       (setf aa ztar)
  189.       (if (and (>= aa 0.0) (> zr 0.0)) (go label110))
  190.       (if (= kode 2) (go label100))
  191.       (if (> aa (- alim)) (go label100))
  192.       (setf aa (- (* 0.25 alaz) aa))
  193.       (setf iflag 1)
  194.       (setf sfac tol)
  195.       (if (> aa elim) (go label270))
  196.      label100
  197.       (setf mr 1)
  198.       (if (< zi 0.0) (setf mr -1))
  199.       (multiple-value-bind
  200.           (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10
  201.            var-11 var-12)
  202.           (zacai ztar ztai fnu kode mr 1 cyr cyi nn rl tol elim alim)
  203.         (declare
  204.          (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-9 var-10
  205.           var-11 var-12))
  206.         (setf nn var-8))
  207.       (if (< nn 0) (go label280))
  208.       (setf nz (f2cl-lib:int-add nz nn))
  209.       (go label130)
  210.      label110
  211.       (if (= kode 2) (go label120))
  212.       (if (< aa alim) (go label120))
  213.       (setf aa (- (* -0.25 alaz) aa))
  214.       (setf iflag 2)
  215.       (setf sfac (/ 1.0 tol))
  216.       (if (< aa (- elim)) (go label210))
  217.      label120
  218.       (multiple-value-bind
  219.           (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10)
  220.           (zbknu ztar ztai fnu kode 1 cyr cyi nz tol elim alim)
  221.         (declare
  222.          (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-8 var-9 var-10))
  223.         (setf nz var-7))
  224.      label130
  225.       (setf s1r (* (f2cl-lib:fref cyr (1) ((1 1))) coef))
  226.       (setf s1i (* (f2cl-lib:fref cyi (1) ((1 1))) coef))
  227.       (if (/= iflag 0) (go label150))
  228.       (if (= id 1) (go label140))
  229.       (setf air (- (* csqr s1r) (* csqi s1i)))
  230.       (setf aii (+ (* csqr s1i) (* csqi s1r)))
  231.       (go end_label)
  232.      label140
  233.       (setf air (- (- (* zr s1r) (* zi s1i))))
  234.       (setf aii (- (+ (* zr s1i) (* zi s1r))))
  235.       (go end_label)
  236.      label150
  237.       (setf s1r (* s1r sfac))
  238.       (setf s1i (* s1i sfac))
  239.       (if (= id 1) (go label160))
  240.       (setf str (- (* s1r csqr) (* s1i csqi)))
  241.       (setf s1i (+ (* s1r csqi) (* s1i csqr)))
  242.       (setf s1r str)
  243.       (setf air (/ s1r sfac))
  244.       (setf aii (/ s1i sfac))
  245.       (go end_label)
  246.      label160
  247.       (setf str (- (- (* s1r zr) (* s1i zi))))
  248.       (setf s1i (- (+ (* s1r zi) (* s1i zr))))
  249.       (setf s1r str)
  250.       (setf air (/ s1r sfac))
  251.       (setf aii (/ s1i sfac))
  252.       (go end_label)
  253.      label170
  254.       (setf aa (* 1000.0 (f2cl-lib:d1mach 1)))
  255.       (setf s1r zeror)
  256.       (setf s1i zeroi)
  257.       (if (= id 1) (go label190))
  258.       (if (<= az aa) (go label180))
  259.       (setf s1r (* c2 zr))
  260.       (setf s1i (* c2 zi))
  261.      label180
  262.       (setf air (- c1 s1r))
  263.       (setf aii (- s1i))
  264.       (go end_label)
  265.      label190
  266.       (setf air (- c2))
  267.       (setf aii 0.0)
  268.       (setf aa (f2cl-lib:fsqrt aa))
  269.       (if (<= az aa) (go label200))
  270.       (setf s1r (* 0.5 (- (* zr zr) (* zi zi))))
  271.       (setf s1i (* zr zi))
  272.      label200
  273.       (setf air (+ air (* c1 s1r)))
  274.       (setf aii (+ aii (* c1 s1i)))
  275.       (go end_label)
  276.      label210
  277.       (setf nz 1)
  278.       (setf air zeror)
  279.       (setf aii zeroi)
  280.       (go end_label)
  281.      label270
  282.       (setf nz 0)
  283.       (setf ierr 2)
  284.       (go end_label)
  285.      label280
  286.       (if (= nn -1) (go label270))
  287.       (setf nz 0)
  288.       (setf ierr 5)
  289.       (go end_label)
  290.      label260
  291.       (setf ierr 4)
  292.       (setf nz 0)
  293.       (go end_label)
  294.      end_label
  295.       (return (values nil nil nil nil air aii nz ierr)))))
  296.  
  297.