src/numerical/slatec/zbesy.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 (defun zbesy (zr zi fnu kode n cyr cyi nz cwrkr cwrki ierr)
  21   (declare (type (simple-array double-float (*)) cwrki cwrkr cyi cyr)
  22            (type (f2cl-lib:integer4) ierr nz n kode)
  23            (type (double-float) fnu zi zr))
  24   (prog ((i 0) (k 0) (k1 0) (k2 0) (nz1 0) (nz2 0) (c1i 0.0) (c1r 0.0)
  25          (c2i 0.0) (c2r 0.0) (elim 0.0) (exi 0.0) (exr 0.0) (ey 0.0) (hcii 0.0)
  26          (sti 0.0) (str 0.0) (tay 0.0) (ascle 0.0) (rtol 0.0) (atol 0.0)
  27          (aa 0.0) (bb 0.0) (tol 0.0) (r1m5 0.0))
  28     (declare (type (double-float) r1m5 tol bb aa atol rtol ascle tay str sti
  29                                   hcii ey exr exi elim c2r c2i c1r c1i)
  30              (type (f2cl-lib:integer4) nz2 nz1 k2 k1 k i))
  31     (setf ierr 0)
  32     (setf nz 0)
  33     (if (and (= zr 0.0) (= zi 0.0)) (setf ierr 1))
  34     (if (< fnu 0.0) (setf ierr 1))
  35     (if (or (< kode 1) (> kode 2)) (setf ierr 1))
  36     (if (< n 1) (setf ierr 1))
  37     (if (/= ierr 0) (go end_label))
  38     (setf hcii 0.5)
  39     (multiple-value-bind
  40           (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9)
  41         (zbesh zr zi fnu kode 1 n cyr cyi nz1 ierr)
  42       (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7))
  43       (setf nz1 var-8)
  44       (setf ierr var-9))
  45     (if (and (/= ierr 0) (/= ierr 3)) (go label170))
  46     (multiple-value-bind
  47           (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9)
  48         (zbesh zr zi fnu kode 2 n cwrkr cwrki nz2 ierr)
  49       (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7))
  50       (setf nz2 var-8)
  51       (setf ierr var-9))
  52     (if (and (/= ierr 0) (/= ierr 3)) (go label170))
  53     (setf nz (min (the f2cl-lib:integer4 nz1) (the f2cl-lib:integer4 nz2)))
  54     (if (= kode 2) (go label60))
  55     (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  56                   ((> i n) nil)
  57       (tagbody
  58         (setf str
  59                 (- (f2cl-lib:fref cwrkr (i) ((1 n)))
  60                    (f2cl-lib:fref cyr (i) ((1 n)))))
  61         (setf sti
  62                 (- (f2cl-lib:fref cwrki (i) ((1 n)))
  63                    (f2cl-lib:fref cyi (i) ((1 n)))))
  64         (setf (f2cl-lib:fref cyr (i) ((1 n))) (* (- sti) hcii))
  65         (setf (f2cl-lib:fref cyi (i) ((1 n))) (* str hcii))
  66        label50))
  67     (go end_label)
  68    label60
  69     (setf tol (max (f2cl-lib:d1mach 4) 1.0e-18))
  70     (setf k1 (f2cl-lib:i1mach 15))
  71     (setf k2 (f2cl-lib:i1mach 16))
  72     (setf k
  73             (min (the f2cl-lib:integer4 (abs k1))
  74                  (the f2cl-lib:integer4 (abs k2))))
  75     (setf r1m5 (f2cl-lib:d1mach 5))
  76     (setf elim (* 2.303 (- (* k r1m5) 3.0)))
  77     (setf exr (cos zr))
  78     (setf exi (sin zr))
  79     (setf ey 0.0)
  80     (setf tay (abs (+ zi zi)))
  81     (if (< tay elim) (setf ey (exp (- tay))))
  82     (if (< zi 0.0) (go label90))
  83     (setf c1r (* exr ey))
  84     (setf c1i (* exi ey))
  85     (setf c2r exr)
  86     (setf c2i (- exi))
  87    label70
  88     (setf nz 0)
  89     (setf rtol (/ 1.0 tol))
  90     (setf ascle (* (f2cl-lib:d1mach 1) rtol 1000.0))
  91     (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  92                   ((> i n) nil)
  93       (tagbody
  94         (setf aa (f2cl-lib:fref cwrkr (i) ((1 n))))
  95         (setf bb (f2cl-lib:fref cwrki (i) ((1 n))))
  96         (setf atol 1.0)
  97         (if (> (max (abs aa) (abs bb)) ascle) (go label75))
  98         (setf aa (* aa rtol))
  99         (setf bb (* bb rtol))
 100         (setf atol tol)
 101        label75
 102         (setf str (* (- (* aa c2r) (* bb c2i)) atol))
 103         (setf sti (* (+ (* aa c2i) (* bb c2r)) atol))
 104         (setf aa (f2cl-lib:fref cyr (i) ((1 n))))
 105         (setf bb (f2cl-lib:fref cyi (i) ((1 n))))
 106         (setf atol 1.0)
 107         (if (> (max (abs aa) (abs bb)) ascle) (go label85))
 108         (setf aa (* aa rtol))
 109         (setf bb (* bb rtol))
 110         (setf atol tol)
 111        label85
 112         (setf str (- str (* (+ (* aa c1r) (* -1 bb c1i)) atol)))
 113         (setf sti (- sti (* (+ (* aa c1i) (* bb c1r)) atol)))
 114         (setf (f2cl-lib:fref cyr (i) ((1 n))) (* (- sti) hcii))
 115         (setf (f2cl-lib:fref cyi (i) ((1 n))) (* str hcii))
 116         (if (and (= str 0.0) (= sti 0.0) (= ey 0.0))
 117             (setf nz (f2cl-lib:int-add nz 1)))
 118        label80))
 119     (go end_label)
 120    label90
 121     (setf c1r exr)
 122     (setf c1i exi)
 123     (setf c2r (* exr ey))
 124     (setf c2i (* (- exi) ey))
 125     (go label70)
 126    label170
 127     (setf nz 0)
 128     (go end_label)
 129    end_label
 130     (return (values nil nil nil nil nil nil nil nz nil nil ierr))))
 131
 132 (in-package #:cl-user)
 133 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 134 (eval-when (:load-toplevel :compile-toplevel :execute)
 135   (setf (gethash 'fortran-to-lisp::zbesy fortran-to-lisp::*f2cl-function-info*)
 136           (fortran-to-lisp::make-f2cl-finfo
 137            :arg-types '((double-float) (double-float) (double-float)
 138                         (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 139                         (simple-array double-float (*))
 140                         (simple-array double-float (*))
 141                         (fortran-to-lisp::integer4)
 142                         (simple-array double-float (*))
 143                         (simple-array double-float (*))
 144                         (fortran-to-lisp::integer4))
 145            :return-values '(nil nil nil nil nil nil nil fortran-to-lisp::nz nil
 146                             nil fortran-to-lisp::ierr)
 147            :calls '(fortran-to-lisp::i1mach fortran-to-lisp::d1mach
 148                     fortran-to-lisp::zbesh))))
 149