src/numerical/slatec/zbiry.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.6149266274460007)
  22       (c2 0.4482883573538264)
  23       (coef 0.5773502691896257)
  24       (pi$ 3.141592653589793)
  25       (coner 1.0)
  26       (conei 0.0))
  27   (declare (type (double-float) tth c1 c2 coef pi$ coner conei))
  28   (defun zbiry (zr zi id kode bir bii ierr)
  29     (declare (type (f2cl-lib:integer4) ierr kode id)
  30              (type (double-float) bii bir zi zr))
  31     (prog ((cyr (make-array 2 :element-type 'double-float))
  32            (cyi (make-array 2 :element-type 'double-float)) (k 0) (k1 0) (k2 0)
  33            (nz 0) (aa 0.0) (ad 0.0) (ak 0.0) (alim 0.0) (atrm 0.0) (az 0.0)
  34            (az3 0.0) (bb 0.0) (bk 0.0) (cc 0.0) (ck 0.0) (csqi 0.0) (csqr 0.0)
  35            (dig 0.0) (dk 0.0) (d1 0.0) (d2 0.0) (eaa 0.0) (elim 0.0) (fid 0.0)
  36            (fmr 0.0) (fnu 0.0) (fnul 0.0) (rl 0.0) (r1m5 0.0) (sfac 0.0)
  37            (sti 0.0) (str 0.0) (s1i 0.0) (s1r 0.0) (s2i 0.0) (s2r 0.0)
  38            (tol 0.0) (trm1i 0.0) (trm1r 0.0) (trm2i 0.0) (trm2r 0.0) (ztai 0.0)
  39            (ztar 0.0) (z3i 0.0) (z3r 0.0))
  40       (declare (type (simple-array double-float (2)) cyr cyi)
  41                (type (double-float) z3r z3i ztar ztai trm2r trm2i trm1r trm1i
  42                                     tol s2r s2i s1r s1i str sti sfac r1m5 rl
  43                                     fnul fnu fmr fid elim eaa d2 d1 dk dig csqr
  44                                     csqi ck cc bk bb az3 az atrm alim ak ad aa)
  45                (type (f2cl-lib:integer4) nz k2 k1 k))
  46       (setf ierr 0)
  47       (setf nz 0)
  48       (if (or (< id 0) (> id 1)) (setf ierr 1))
  49       (if (or (< kode 1) (> kode 2)) (setf ierr 1))
  50       (if (/= ierr 0) (go end_label))
  51       (setf az (coerce (realpart (zabs zr zi)) 'double-float))
  52       (setf tol (max (f2cl-lib:d1mach 4) 1.0e-18))
  53       (setf fid (coerce (the f2cl-lib:integer4 id) 'double-float))
  54       (if (> az 1.0f0) (go label70))
  55       (setf s1r coner)
  56       (setf s1i conei)
  57       (setf s2r coner)
  58       (setf s2i conei)
  59       (if (< az tol) (go label130))
  60       (setf aa (* az az))
  61       (if (< aa (/ tol az)) (go label40))
  62       (setf trm1r coner)
  63       (setf trm1i conei)
  64       (setf trm2r coner)
  65       (setf trm2i conei)
  66       (setf atrm 1.0)
  67       (setf str (- (* zr zr) (* zi zi)))
  68       (setf sti (+ (* zr zi) (* zi zr)))
  69       (setf z3r (- (* str zr) (* sti zi)))
  70       (setf z3i (+ (* str zi) (* sti zr)))
  71       (setf az3 (* az aa))
  72       (setf ak (+ 2.0 fid))
  73       (setf bk (- 3.0 fid fid))
  74       (setf ck (- 4.0 fid))
  75       (setf dk (+ 3.0 fid fid))
  76       (setf d1 (* ak dk))
  77       (setf d2 (* bk ck))
  78       (setf ad (min d1 d2))
  79       (setf ak (+ 24.0 (* 9.0 fid)))
  80       (setf bk (- 30.0 (* 9.0 fid)))
  81       (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
  82                     ((> k 25) nil)
  83         (tagbody
  84           (setf str (/ (- (* trm1r z3r) (* trm1i z3i)) d1))
  85           (setf trm1i (/ (+ (* trm1r z3i) (* trm1i z3r)) d1))
  86           (setf trm1r str)
  87           (setf s1r (+ s1r trm1r))
  88           (setf s1i (+ s1i trm1i))
  89           (setf str (/ (- (* trm2r z3r) (* trm2i z3i)) d2))
  90           (setf trm2i (/ (+ (* trm2r z3i) (* trm2i z3r)) d2))
  91           (setf trm2r str)
  92           (setf s2r (+ s2r trm2r))
  93           (setf s2i (+ s2i trm2i))
  94           (setf atrm (/ (* atrm az3) ad))
  95           (setf d1 (+ d1 ak))
  96           (setf d2 (+ d2 bk))
  97           (setf ad (min d1 d2))
  98           (if (< atrm (* tol ad)) (go label40))
  99           (setf ak (+ ak 18.0))
 100           (setf bk (+ bk 18.0))
 101          label30))
 102      label40
 103       (if (= id 1) (go label50))
 104       (setf bir (+ (* c1 s1r) (* c2 (- (* zr s2r) (* zi s2i)))))
 105       (setf bii (+ (* c1 s1i) (* c2 (+ (* zr s2i) (* zi s2r)))))
 106       (if (= kode 1) (go end_label))
 107       (multiple-value-bind (var-0 var-1 var-2 var-3)
 108           (zsqrt$ zr zi str sti)
 109         (declare (ignore var-0 var-1))
 110         (setf str var-2)
 111         (setf sti var-3))
 112       (setf ztar (* tth (- (* zr str) (* zi sti))))
 113       (setf ztai (* tth (+ (* zr sti) (* zi str))))
 114       (setf aa ztar)
 115       (setf aa (- (abs aa)))
 116       (setf eaa (exp aa))
 117       (setf bir (* bir eaa))
 118       (setf bii (* bii eaa))
 119       (go end_label)
 120      label50
 121       (setf bir (* s2r c2))
 122       (setf bii (* s2i c2))
 123       (if (<= az tol) (go label60))
 124       (setf cc (/ c1 (+ 1.0 fid)))
 125       (setf str (- (* s1r zr) (* s1i zi)))
 126       (setf sti (+ (* s1r zi) (* s1i zr)))
 127       (setf bir (+ bir (* cc (- (* str zr) (* sti zi)))))
 128       (setf bii (+ bii (* cc (+ (* str zi) (* sti zr)))))
 129      label60
 130       (if (= kode 1) (go end_label))
 131       (multiple-value-bind (var-0 var-1 var-2 var-3)
 132           (zsqrt$ zr zi str sti)
 133         (declare (ignore var-0 var-1))
 134         (setf str var-2)
 135         (setf sti var-3))
 136       (setf ztar (* tth (- (* zr str) (* zi sti))))
 137       (setf ztai (* tth (+ (* zr sti) (* zi str))))
 138       (setf aa ztar)
 139       (setf aa (- (abs aa)))
 140       (setf eaa (exp aa))
 141       (setf bir (* bir eaa))
 142       (setf bii (* bii eaa))
 143       (go end_label)
 144      label70
 145       (setf fnu (/ (+ 1.0 fid) 3.0))
 146       (setf k1 (f2cl-lib:i1mach 15))
 147       (setf k2 (f2cl-lib:i1mach 16))
 148       (setf r1m5 (f2cl-lib:d1mach 5))
 149       (setf k
 150               (min (the f2cl-lib:integer4 (abs k1))
 151                    (the f2cl-lib:integer4 (abs k2))))
 152       (setf elim (* 2.303 (- (* k r1m5) 3.0)))
 153       (setf k1 (f2cl-lib:int-sub (f2cl-lib:i1mach 14) 1))
 154       (setf aa (* r1m5 k1))
 155       (setf dig (min aa 18.0))
 156       (setf aa (* aa 2.303))
 157       (setf alim (+ elim (max (- aa) -41.45)))
 158       (setf rl (+ (* 1.2 dig) 3.0))
 159       (setf fnul (+ 10.0 (* 6.0 (- dig 3.0))))
 160       (setf aa (/ 0.5 tol))
 161       (setf bb (* (f2cl-lib:i1mach 9) 0.5))
 162       (setf aa (min aa bb))
 163       (setf aa (expt aa tth))
 164       (if (> az aa) (go label260))
 165       (setf aa (f2cl-lib:fsqrt aa))
 166       (if (> az aa) (setf ierr 3))
 167       (multiple-value-bind (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 sfac 1.0)
 175       (setf ak ztai)
 176       (if (>= zr 0.0) (go label80))
 177       (setf bk ztar)
 178       (setf ck (- (abs bk)))
 179       (setf ztar ck)
 180       (setf ztai ak)
 181      label80
 182       (if (or (/= zi 0.0) (> zr 0.0)) (go label90))
 183       (setf ztar 0.0)
 184       (setf ztai ak)
 185      label90
 186       (setf aa ztar)
 187       (if (= kode 2) (go label100))
 188       (setf bb (abs aa))
 189       (if (< bb alim) (go label100))
 190       (setf bb (+ bb (* 0.25 (f2cl-lib:flog az))))
 191       (setf sfac tol)
 192       (if (> bb elim) (go label190))
 193      label100
 194       (setf fmr 0.0)
 195       (if (and (>= aa 0.0) (> zr 0.0)) (go label110))
 196       (setf fmr pi$)
 197       (if (< zi 0.0) (setf fmr (- pi$)))
 198       (setf ztar (- ztar))
 199       (setf ztai (- ztai))
 200      label110
 201       (multiple-value-bind
 202             (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10
 203              var-11 var-12)
 204           (zbinu ztar ztai fnu kode 1 cyr cyi nz rl fnul tol elim alim)
 205         (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-8 var-9
 206                          var-10 var-11 var-12))
 207         (setf nz var-7))
 208       (if (< nz 0) (go label200))
 209       (setf aa (* fmr fnu))
 210       (setf z3r sfac)
 211       (setf str (cos aa))
 212       (setf sti (sin aa))
 213       (setf s1r
 214               (*
 215                (- (* str (f2cl-lib:fref cyr (1) ((1 2))))
 216                   (* sti (f2cl-lib:fref cyi (1) ((1 2)))))
 217                z3r))
 218       (setf s1i
 219               (*
 220                (+ (* str (f2cl-lib:fref cyi (1) ((1 2))))
 221                   (* sti (f2cl-lib:fref cyr (1) ((1 2)))))
 222                z3r))
 223       (setf fnu (/ (- 2.0 fid) 3.0))
 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 var-10
 226              var-11 var-12)
 227           (zbinu ztar ztai fnu kode 2 cyr cyi nz rl fnul 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 var-11 var-12))
 230         (setf nz var-7))
 231       (setf (f2cl-lib:fref cyr (1) ((1 2)))
 232               (* (f2cl-lib:fref cyr (1) ((1 2))) z3r))
 233       (setf (f2cl-lib:fref cyi (1) ((1 2)))
 234               (* (f2cl-lib:fref cyi (1) ((1 2))) z3r))
 235       (setf (f2cl-lib:fref cyr (2) ((1 2)))
 236               (* (f2cl-lib:fref cyr (2) ((1 2))) z3r))
 237       (setf (f2cl-lib:fref cyi (2) ((1 2)))
 238               (* (f2cl-lib:fref cyi (2) ((1 2))) z3r))
 239       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 240           (zdiv (f2cl-lib:fref cyr (1) ((1 2))) (f2cl-lib:fref cyi (1) ((1 2)))
 241            ztar ztai str sti)
 242         (declare (ignore var-0 var-1 var-2 var-3))
 243         (setf str var-4)
 244         (setf sti var-5))
 245       (setf s2r (+ (* (+ fnu fnu) str) (f2cl-lib:fref cyr (2) ((1 2)))))
 246       (setf s2i (+ (* (+ fnu fnu) sti) (f2cl-lib:fref cyi (2) ((1 2)))))
 247       (setf aa (* fmr (- fnu 1.0)))
 248       (setf str (cos aa))
 249       (setf sti (sin aa))
 250       (setf s1r (* coef (- (+ s1r (* s2r str)) (* s2i sti))))
 251       (setf s1i (* coef (+ s1i (* s2r sti) (* s2i str))))
 252       (if (= id 1) (go label120))
 253       (setf str (- (* csqr s1r) (* csqi s1i)))
 254       (setf s1i (+ (* csqr s1i) (* csqi s1r)))
 255       (setf s1r str)
 256       (setf bir (/ s1r sfac))
 257       (setf bii (/ s1i sfac))
 258       (go end_label)
 259      label120
 260       (setf str (- (* zr s1r) (* zi s1i)))
 261       (setf s1i (+ (* zr s1i) (* zi s1r)))
 262       (setf s1r str)
 263       (setf bir (/ s1r sfac))
 264       (setf bii (/ s1i sfac))
 265       (go end_label)
 266      label130
 267       (setf aa (+ (* c1 (- 1.0 fid)) (* fid c2)))
 268       (setf bir aa)
 269       (setf bii 0.0)
 270       (go end_label)
 271      label190
 272       (setf ierr 2)
 273       (setf nz 0)
 274       (go end_label)
 275      label200
 276       (if (= nz -1) (go label190))
 277       (setf nz 0)
 278       (setf ierr 5)
 279       (go end_label)
 280      label260
 281       (setf ierr 4)
 282       (setf nz 0)
 283       (go end_label)
 284      end_label
 285       (return (values nil nil nil nil bir bii ierr)))))
 286
 287 (in-package #:cl-user)
 288 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 289 (eval-when (:load-toplevel :compile-toplevel :execute)
 290   (setf (gethash 'fortran-to-lisp::zbiry fortran-to-lisp::*f2cl-function-info*)
 291           (fortran-to-lisp::make-f2cl-finfo
 292            :arg-types '((double-float) (double-float)
 293                         (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 294                         (double-float) (double-float)
 295                         (fortran-to-lisp::integer4))
 296            :return-values '(nil nil nil nil fortran-to-lisp::bir
 297                             fortran-to-lisp::bii fortran-to-lisp::ierr)
 298            :calls '(fortran-to-lisp::zdiv fortran-to-lisp::zbinu
 299                     fortran-to-lisp::i1mach fortran-to-lisp::zsqrt$
 300                     fortran-to-lisp::d1mach fortran-to-lisp::zabs))))
 301