src/numerical/slatec/zairy.lisp

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