src/numerical/slatec/zbesj.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 zbesj (zr zi fnu kode n cyr cyi nz ierr)
  23     (declare (type (simple-array double-float (*)) cyi cyr)
  24              (type (f2cl-lib:integer4) ierr nz n kode)
  25              (type (double-float) fnu zi zr))
  26     (prog ((i 0) (inu 0) (inuh 0) (ir 0) (k 0) (k1 0) (k2 0) (nl 0) (aa 0.0)
  27            (alim 0.0) (arg 0.0) (cii 0.0) (csgni 0.0) (csgnr 0.0) (dig 0.0)
  28            (elim 0.0) (fnul 0.0) (rl 0.0) (r1m5 0.0) (str 0.0) (tol 0.0)
  29            (zni 0.0) (znr 0.0) (bb 0.0) (fn 0.0) (az 0.0) (ascle 0.0)
  30            (rtol 0.0) (atol 0.0) (sti 0.0))
  31       (declare (type (double-float) sti atol rtol ascle az fn bb znr zni tol
  32                                     str r1m5 rl fnul elim dig csgnr csgni cii
  33                                     arg alim aa)
  34                (type (f2cl-lib:integer4) nl k2 k1 k ir inuh inu i))
  35       (setf ierr 0)
  36       (setf nz 0)
  37       (if (< fnu 0.0) (setf ierr 1))
  38       (if (or (< kode 1) (> kode 2)) (setf ierr 1))
  39       (if (< n 1) (setf ierr 1))
  40       (if (/= ierr 0) (go end_label))
  41       (setf tol (max (f2cl-lib:d1mach 4) 1.0e-18))
  42       (setf k1 (f2cl-lib:i1mach 15))
  43       (setf k2 (f2cl-lib:i1mach 16))
  44       (setf r1m5 (f2cl-lib:d1mach 5))
  45       (setf k
  46               (min (the f2cl-lib:integer4 (abs k1))
  47                    (the f2cl-lib:integer4 (abs k2))))
  48       (setf elim (* 2.303 (- (* k r1m5) 3.0)))
  49       (setf k1 (f2cl-lib:int-sub (f2cl-lib:i1mach 14) 1))
  50       (setf aa (* r1m5 k1))
  51       (setf dig (min aa 18.0))
  52       (setf aa (* aa 2.303))
  53       (setf alim (+ elim (max (- aa) -41.45)))
  54       (setf rl (+ (* 1.2 dig) 3.0))
  55       (setf fnul (+ 10.0 (* 6.0 (- dig 3.0))))
  56       (setf az (coerce (realpart (zabs zr zi)) 'double-float))
  57       (setf fn (+ fnu (f2cl-lib:int-sub n 1)))
  58       (setf aa (/ 0.5 tol))
  59       (setf bb (* (f2cl-lib:i1mach 9) 0.5))
  60       (setf aa (min aa bb))
  61       (if (> az aa) (go label260))
  62       (if (> fn aa) (go label260))
  63       (setf aa (f2cl-lib:fsqrt aa))
  64       (if (> az aa) (setf ierr 3))
  65       (if (> fn aa) (setf ierr 3))
  66       (setf cii 1.0)
  67       (setf inu (f2cl-lib:int fnu))
  68       (setf inuh (the f2cl-lib:integer4 (truncate inu 2)))
  69       (setf ir (f2cl-lib:int-sub inu (f2cl-lib:int-mul 2 inuh)))
  70       (setf arg (* (- fnu (f2cl-lib:int-sub inu ir)) hpi))
  71       (setf csgnr (cos arg))
  72       (setf csgni (sin arg))
  73       (if (= (mod inuh 2) 0) (go label40))
  74       (setf csgnr (- csgnr))
  75       (setf csgni (- csgni))
  76      label40
  77       (setf znr zi)
  78       (setf zni (- zr))
  79       (if (>= zi 0.0) (go label50))
  80       (setf znr (- znr))
  81       (setf zni (- zni))
  82       (setf csgni (- csgni))
  83       (setf cii (- cii))
  84      label50
  85       (multiple-value-bind
  86             (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10
  87              var-11 var-12)
  88           (zbinu znr zni fnu kode n cyr cyi nz rl fnul tol elim alim)
  89         (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-8 var-9
  90                          var-10 var-11 var-12))
  91         (setf nz var-7))
  92       (if (< nz 0) (go label130))
  93       (setf nl (f2cl-lib:int-sub n nz))
  94       (if (= nl 0) (go end_label))
  95       (setf rtol (/ 1.0 tol))
  96       (setf ascle (* (f2cl-lib:d1mach 1) rtol 1000.0))
  97       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  98                     ((> i nl) nil)
  99         (tagbody
 100           (setf aa (f2cl-lib:fref cyr (i) ((1 n))))
 101           (setf bb (f2cl-lib:fref cyi (i) ((1 n))))
 102           (setf atol 1.0)
 103           (if (> (max (abs aa) (abs bb)) ascle) (go label55))
 104           (setf aa (* aa rtol))
 105           (setf bb (* bb rtol))
 106           (setf atol tol)
 107          label55
 108           (setf str (- (* aa csgnr) (* bb csgni)))
 109           (setf sti (+ (* aa csgni) (* bb csgnr)))
 110           (setf (f2cl-lib:fref cyr (i) ((1 n))) (* str atol))
 111           (setf (f2cl-lib:fref cyi (i) ((1 n))) (* sti atol))
 112           (setf str (* (- csgni) cii))
 113           (setf csgni (* csgnr cii))
 114           (setf csgnr str)
 115          label60))
 116       (go end_label)
 117      label130
 118       (if (= nz -2) (go label140))
 119       (setf nz 0)
 120       (setf ierr 2)
 121       (go end_label)
 122      label140
 123       (setf nz 0)
 124       (setf ierr 5)
 125       (go end_label)
 126      label260
 127       (setf nz 0)
 128       (setf ierr 4)
 129       (go end_label)
 130      end_label
 131       (return (values nil nil nil nil nil nil nil nz ierr)))))
 132
 133 (in-package #:cl-user)
 134 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 135 (eval-when (:load-toplevel :compile-toplevel :execute)
 136   (setf (gethash 'fortran-to-lisp::zbesj fortran-to-lisp::*f2cl-function-info*)
 137           (fortran-to-lisp::make-f2cl-finfo
 138            :arg-types '((double-float) (double-float) (double-float)
 139                         (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 140                         (simple-array double-float (*))
 141                         (simple-array double-float (*))
 142                         (fortran-to-lisp::integer4)
 143                         (fortran-to-lisp::integer4))
 144            :return-values '(nil nil nil nil nil nil nil fortran-to-lisp::nz
 145                             fortran-to-lisp::ierr)
 146            :calls '(fortran-to-lisp::zbinu fortran-to-lisp::zabs
 147                     fortran-to-lisp::i1mach fortran-to-lisp::d1mach))))
 148