src/numerical/slatec/zbesh.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 ((hpi 1.5707963267948966))
  21   (declare (type (double-float) hpi))
  22   (defun zbesh (zr zi fnu kode m n cyr cyi nz ierr)
  23     (declare (type (simple-array double-float (*)) cyi cyr)
  24              (type (f2cl-lib:integer4) ierr nz n m kode)
  25              (type (double-float) fnu zi zr))
  26     (prog ((i 0) (inu 0) (inuh 0) (ir 0) (k 0) (k1 0) (k2 0) (mm 0) (mr 0)
  27            (nn 0) (nuf 0) (nw 0) (aa 0.0) (alim 0.0) (aln 0.0) (arg 0.0)
  28            (az 0.0) (dig 0.0) (elim 0.0) (fmm 0.0) (fn 0.0) (fnul 0.0)
  29            (rhpi 0.0) (rl 0.0) (r1m5 0.0) (sgn 0.0) (str 0.0) (tol 0.0)
  30            (ufl 0.0) (zni 0.0) (znr 0.0) (zti 0.0) (bb 0.0) (ascle 0.0)
  31            (rtol 0.0) (atol 0.0) (sti 0.0) (csgnr 0.0) (csgni 0.0))
  32       (declare (type (double-float) csgni csgnr sti atol rtol ascle bb zti znr
  33                                     zni ufl tol str sgn r1m5 rl rhpi fnul fn
  34                                     fmm elim dig az arg aln alim aa)
  35                (type (f2cl-lib:integer4) nw nuf nn mr mm k2 k1 k ir inuh inu
  36                                          i))
  37       (setf ierr 0)
  38       (setf nz 0)
  39       (if (and (= zr 0.0) (= zi 0.0)) (setf ierr 1))
  40       (if (< fnu 0.0) (setf ierr 1))
  41       (if (or (< m 1) (> m 2)) (setf ierr 1))
  42       (if (or (< kode 1) (> kode 2)) (setf ierr 1))
  43       (if (< n 1) (setf ierr 1))
  44       (if (/= ierr 0) (go end_label))
  45       (setf nn n)
  46       (setf tol (max (f2cl-lib:d1mach 4) 1.0e-18))
  47       (setf k1 (f2cl-lib:i1mach 15))
  48       (setf k2 (f2cl-lib:i1mach 16))
  49       (setf r1m5 (f2cl-lib:d1mach 5))
  50       (setf k
  51               (min (the f2cl-lib:integer4 (abs k1))
  52                    (the f2cl-lib:integer4 (abs k2))))
  53       (setf elim (* 2.303 (- (* k r1m5) 3.0)))
  54       (setf k1 (f2cl-lib:int-sub (f2cl-lib:i1mach 14) 1))
  55       (setf aa (* r1m5 k1))
  56       (setf dig (min aa 18.0))
  57       (setf aa (* aa 2.303))
  58       (setf alim (+ elim (max (- aa) -41.45)))
  59       (setf fnul (+ 10.0 (* 6.0 (- dig 3.0))))
  60       (setf rl (+ (* 1.2 dig) 3.0))
  61       (setf fn (+ fnu (f2cl-lib:int-sub nn 1)))
  62       (setf mm (f2cl-lib:int-sub 3 m m))
  63       (setf fmm (coerce (the f2cl-lib:integer4 mm) 'double-float))
  64       (setf znr (* fmm zi))
  65       (setf zni (* (- fmm) zr))
  66       (setf az (coerce (realpart (zabs zr zi)) 'double-float))
  67       (setf aa (/ 0.5 tol))
  68       (setf bb (* (f2cl-lib:i1mach 9) 0.5))
  69       (setf aa (min aa bb))
  70       (if (> az aa) (go label260))
  71       (if (> fn aa) (go label260))
  72       (setf aa (f2cl-lib:fsqrt aa))
  73       (if (> az aa) (setf ierr 3))
  74       (if (> fn aa) (setf ierr 3))
  75       (setf ufl (* (f2cl-lib:d1mach 1) 1000.0))
  76       (if (< az ufl) (go label230))
  77       (if (> fnu fnul) (go label90))
  78       (if (<= fn 1.0) (go label70))
  79       (if (> fn 2.0) (go label60))
  80       (if (> az tol) (go label70))
  81       (setf arg (* 0.5 az))
  82       (setf aln (* (- fn) (f2cl-lib:flog arg)))
  83       (if (> aln elim) (go label230))
  84       (go label70)
  85      label60
  86       (multiple-value-bind
  87             (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10
  88              var-11)
  89           (zuoik znr zni fnu kode 2 nn cyr cyi nuf tol elim alim)
  90         (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-9
  91                          var-10 var-11))
  92         (setf nuf var-8))
  93       (if (< nuf 0) (go label230))
  94       (setf nz (f2cl-lib:int-add nz nuf))
  95       (setf nn (f2cl-lib:int-sub nn nuf))
  96       (if (= nn 0) (go label140))
  97      label70
  98       (if (or (< znr 0.0) (and (= znr 0.0) (< zni 0.0) (= m 2))) (go label80))
  99       (multiple-value-bind
 100             (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9
 101              var-10)
 102           (zbknu znr zni fnu kode nn cyr cyi nz tol elim alim)
 103         (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-8 var-9
 104                          var-10))
 105         (setf nz var-7))
 106       (go label110)
 107      label80
 108       (setf mr (f2cl-lib:int-sub mm))
 109       (multiple-value-bind
 110             (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10
 111              var-11 var-12 var-13)
 112           (zacon znr zni fnu kode mr nn cyr cyi nw rl fnul tol elim alim)
 113         (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-9
 114                          var-10 var-11 var-12 var-13))
 115         (setf nw var-8))
 116       (if (< nw 0) (go label240))
 117       (setf nz nw)
 118       (go label110)
 119      label90
 120       (setf mr 0)
 121       (if (and (>= znr 0.0) (or (/= znr 0.0) (>= zni 0.0) (/= m 2)))
 122           (go label100))
 123       (setf mr (f2cl-lib:int-sub mm))
 124       (if (or (/= znr 0.0) (>= zni 0.0)) (go label100))
 125       (setf znr (- znr))
 126       (setf zni (- zni))
 127      label100
 128       (multiple-value-bind
 129             (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10
 130              var-11)
 131           (zbunk znr zni fnu kode mr nn cyr cyi nw tol elim alim)
 132         (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-9
 133                          var-10 var-11))
 134         (setf nw var-8))
 135       (if (< nw 0) (go label240))
 136       (setf nz (f2cl-lib:int-add nz nw))
 137      label110
 138       (setf sgn (coerce (f2cl-lib:dsign hpi (- fmm)) 'double-float))
 139       (setf inu (f2cl-lib:int fnu))
 140       (setf inuh (the f2cl-lib:integer4 (truncate inu 2)))
 141       (setf ir (f2cl-lib:int-sub inu (f2cl-lib:int-mul 2 inuh)))
 142       (setf arg (* (- fnu (f2cl-lib:int-sub inu ir)) sgn))
 143       (setf rhpi (/ 1.0 sgn))
 144       (setf csgni (* rhpi (cos arg)))
 145       (setf csgnr (* (- rhpi) (sin arg)))
 146       (if (= (mod inuh 2) 0) (go label120))
 147       (setf csgnr (- csgnr))
 148       (setf csgni (- csgni))
 149      label120
 150       (setf zti (- fmm))
 151       (setf rtol (/ 1.0 tol))
 152       (setf ascle (* ufl rtol))
 153       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 154                     ((> i nn) nil)
 155         (tagbody
 156           (setf aa (f2cl-lib:fref cyr (i) ((1 n))))
 157           (setf bb (f2cl-lib:fref cyi (i) ((1 n))))
 158           (setf atol 1.0)
 159           (if (> (max (abs aa) (abs bb)) ascle) (go label135))
 160           (setf aa (* aa rtol))
 161           (setf bb (* bb rtol))
 162           (setf atol tol)
 163          label135
 164           (setf str (- (* aa csgnr) (* bb csgni)))
 165           (setf sti (+ (* aa csgni) (* bb csgnr)))
 166           (setf (f2cl-lib:fref cyr (i) ((1 n))) (* str atol))
 167           (setf (f2cl-lib:fref cyi (i) ((1 n))) (* sti atol))
 168           (setf str (* (- csgni) zti))
 169           (setf csgni (* csgnr zti))
 170           (setf csgnr str)
 171          label130))
 172       (go end_label)
 173      label140
 174       (if (< znr 0.0) (go label230))
 175       (go end_label)
 176      label230
 177       (setf nz 0)
 178       (setf ierr 2)
 179       (go end_label)
 180      label240
 181       (if (= nw -1) (go label230))
 182       (setf nz 0)
 183       (setf ierr 5)
 184       (go end_label)
 185      label260
 186       (setf nz 0)
 187       (setf ierr 4)
 188       (go end_label)
 189      end_label
 190       (return (values nil nil nil nil nil nil nil nil nz ierr)))))
 191
 192 (in-package #:cl-user)
 193 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 194 (eval-when (:load-toplevel :compile-toplevel :execute)
 195   (setf (gethash 'fortran-to-lisp::zbesh fortran-to-lisp::*f2cl-function-info*)
 196           (fortran-to-lisp::make-f2cl-finfo
 197            :arg-types '((double-float) (double-float) (double-float)
 198                         (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 199                         (fortran-to-lisp::integer4)
 200                         (simple-array double-float (*))
 201                         (simple-array double-float (*))
 202                         (fortran-to-lisp::integer4)
 203                         (fortran-to-lisp::integer4))
 204            :return-values '(nil nil nil nil nil nil nil nil fortran-to-lisp::nz
 205                             fortran-to-lisp::ierr)
 206            :calls '(fortran-to-lisp::zbunk fortran-to-lisp::zacon
 207                     fortran-to-lisp::zbknu fortran-to-lisp::zuoik
 208                     fortran-to-lisp::zabs fortran-to-lisp::i1mach
 209                     fortran-to-lisp::d1mach))))
 210