src/numerical/slatec/dbesk.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 ((nulim
  21        (make-array 2
  22                    :element-type 'f2cl-lib:integer4
  23                    :initial-contents '(35 70))))
  24   (declare (type (simple-array f2cl-lib:integer4 (2)) nulim))
  25   (defun dbesk (x fnu kode n y nz)
  26     (declare (type (simple-array double-float (*)) y)
  27              (type (f2cl-lib:integer4) nz n kode)
  28              (type (double-float) fnu x))
  29     (prog ((w (make-array 2 :element-type 'double-float)) (cn 0.0) (dnu 0.0)
  30            (elim 0.0) (etx 0.0) (flgik 0.0) (fn 0.0) (fnn 0.0) (gln 0.0)
  31            (gnu 0.0) (rtz 0.0) (s 0.0) (s1 0.0) (s2 0.0) (t$ 0.0) (tm 0.0)
  32            (trx 0.0) (xlim 0.0) (zn 0.0) (i 0) (j 0) (k 0) (mz 0) (nb 0) (nd 0)
  33            (nn 0) (nud 0))
  34       (declare (type (f2cl-lib:integer4) nud nn nd nb mz k j i)
  35                (type (simple-array double-float (2)) w)
  36                (type (double-float) zn xlim trx tm t$ s2 s1 s rtz gnu gln fnn
  37                                     fn flgik etx elim dnu cn))
  38       (setf nn (f2cl-lib:int-sub (f2cl-lib:i1mach 15)))
  39       (setf elim (* 2.303 (- (* nn (f2cl-lib:d1mach 5)) 3.0)))
  40       (setf xlim (* (f2cl-lib:d1mach 1) 1000.0))
  41       (if (or (< kode 1) (> kode 2)) (go label280))
  42       (if (< fnu 0.0) (go label290))
  43       (if (<= x 0.0) (go label300))
  44       (if (< x xlim) (go label320))
  45       (if (< n 1) (go label310))
  46       (setf etx
  47               (coerce (the f2cl-lib:integer4 (f2cl-lib:int-sub kode 1))
  48                       'double-float))
  49       (setf nd n)
  50       (setf nz 0)
  51       (setf nud (f2cl-lib:int fnu))
  52       (setf dnu (- fnu nud))
  53       (setf gnu fnu)
  54       (setf nn (min (the f2cl-lib:integer4 2) (the f2cl-lib:integer4 nd)))
  55       (setf fn (- (+ fnu n) 1))
  56       (setf fnn fn)
  57       (if (< fn 2.0) (go label150))
  58       (setf zn (/ x fn))
  59       (if (= zn 0.0) (go label320))
  60       (setf rtz (f2cl-lib:fsqrt (+ 1.0 (* zn zn))))
  61       (setf gln (f2cl-lib:flog (/ (+ 1.0 rtz) zn)))
  62       (setf t$ (+ (* rtz (- 1.0 etx)) (/ etx (+ zn rtz))))
  63       (setf cn (* (- fn) (- t$ gln)))
  64       (if (> cn elim) (go label320))
  65       (if (< nud (f2cl-lib:fref nulim (nn) ((1 2)))) (go label30))
  66       (if (= nn 1) (go label20))
  67      label10
  68       (setf fn gnu)
  69       (setf zn (/ x fn))
  70       (setf rtz (f2cl-lib:fsqrt (+ 1.0 (* zn zn))))
  71       (setf gln (f2cl-lib:flog (/ (+ 1.0 rtz) zn)))
  72       (setf t$ (+ (* rtz (- 1.0 etx)) (/ etx (+ zn rtz))))
  73       (setf cn (* (- fn) (- t$ gln)))
  74      label20
  75       (if (< cn (- elim)) (go label230))
  76       (setf flgik -1.0)
  77       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7)
  78           (dasyik x gnu kode flgik rtz cn nn y)
  79         (declare (ignore var-0 var-1 var-2 var-3 var-6 var-7))
  80         (setf rtz var-4)
  81         (setf cn var-5))
  82       (if (= nn 1) (go label240))
  83       (setf trx (/ 2.0 x))
  84       (setf tm (/ (+ gnu gnu 2.0) x))
  85       (go label130)
  86      label30
  87       (if (= kode 2) (go label40))
  88       (if (> x elim) (go label230))
  89      label40
  90       (if (/= dnu 0.0) (go label80))
  91       (if (= kode 2) (go label50))
  92       (setf s1 (dbesk0 x))
  93       (go label60)
  94      label50
  95       (setf s1 (dbsk0e x))
  96      label60
  97       (if (and (= nud 0) (= nd 1)) (go label120))
  98       (if (= kode 2) (go label70))
  99       (setf s2 (dbesk1 x))
 100       (go label90)
 101      label70
 102       (setf s2 (dbsk1e x))
 103       (go label90)
 104      label80
 105       (setf nb 2)
 106       (if (and (= nud 0) (= nd 1)) (setf nb 1))
 107       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 108           (dbsknu x dnu kode nb w nz)
 109         (declare (ignore var-0 var-1 var-2 var-3 var-4))
 110         (setf nz var-5))
 111       (setf s1 (f2cl-lib:fref w (1) ((1 2))))
 112       (if (= nb 1) (go label120))
 113       (setf s2 (f2cl-lib:fref w (2) ((1 2))))
 114      label90
 115       (setf trx (/ 2.0 x))
 116       (setf tm (/ (+ dnu dnu 2.0) x))
 117       (if (= nd 1) (setf nud (f2cl-lib:int-sub nud 1)))
 118       (if (> nud 0) (go label100))
 119       (if (> nd 1) (go label120))
 120       (setf s1 s2)
 121       (go label120)
 122      label100
 123       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 124                     ((> i nud) nil)
 125         (tagbody
 126           (setf s s2)
 127           (setf s2 (+ (* tm s2) s1))
 128           (setf s1 s)
 129           (setf tm (+ tm trx))
 130          label110))
 131       (if (= nd 1) (setf s1 s2))
 132      label120
 133       (setf (f2cl-lib:fref y (1) ((1 *))) s1)
 134       (if (= nd 1) (go label240))
 135       (setf (f2cl-lib:fref y (2) ((1 *))) s2)
 136      label130
 137       (if (= nd 2) (go label240))
 138       (f2cl-lib:fdo (i 3 (f2cl-lib:int-add i 1))
 139                     ((> i nd) nil)
 140         (tagbody
 141           (setf (f2cl-lib:fref y (i) ((1 *)))
 142                   (+ (* tm (f2cl-lib:fref y ((f2cl-lib:int-sub i 1)) ((1 *))))
 143                      (f2cl-lib:fref y ((f2cl-lib:int-sub i 2)) ((1 *)))))
 144           (setf tm (+ tm trx))
 145          label140))
 146       (go label240)
 147      label150
 148       (if (= kode 2) (go label160))
 149       (if (> x elim) (go label230))
 150      label160
 151       (if (<= fn 1.0) (go label170))
 152       (if (> (* (- fn) (- (f2cl-lib:flog x) 0.693)) elim) (go label320))
 153      label170
 154       (if (= dnu 0.0) (go label180))
 155       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 156           (dbsknu x fnu kode nd y mz)
 157         (declare (ignore var-0 var-1 var-2 var-3 var-4))
 158         (setf mz var-5))
 159       (go label240)
 160      label180
 161       (setf j nud)
 162       (if (= j 1) (go label210))
 163       (setf j (f2cl-lib:int-add j 1))
 164       (if (= kode 2) (go label190))
 165       (setf (f2cl-lib:fref y (j) ((1 *))) (dbesk0 x))
 166       (go label200)
 167      label190
 168       (setf (f2cl-lib:fref y (j) ((1 *))) (dbsk0e x))
 169      label200
 170       (if (= nd 1) (go label240))
 171       (setf j (f2cl-lib:int-add j 1))
 172      label210
 173       (if (= kode 2) (go label220))
 174       (setf (f2cl-lib:fref y (j) ((1 *))) (dbesk1 x))
 175       (go label240)
 176      label220
 177       (setf (f2cl-lib:fref y (j) ((1 *))) (dbsk1e x))
 178       (go label240)
 179      label230
 180       (setf nud (f2cl-lib:int-add nud 1))
 181       (setf nd (f2cl-lib:int-sub nd 1))
 182       (if (= nd 0) (go label240))
 183       (setf nn (min (the f2cl-lib:integer4 2) (the f2cl-lib:integer4 nd)))
 184       (setf gnu (+ gnu 1.0))
 185       (if (< fnn 2.0) (go label230))
 186       (if (< nud (f2cl-lib:fref nulim (nn) ((1 2)))) (go label230))
 187       (go label10)
 188      label240
 189       (setf nz (f2cl-lib:int-sub n nd))
 190       (if (= nz 0) (go end_label))
 191       (if (= nd 0) (go label260))
 192       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 193                     ((> i nd) nil)
 194         (tagbody
 195           (setf j (f2cl-lib:int-add (f2cl-lib:int-sub n i) 1))
 196           (setf k (f2cl-lib:int-add (f2cl-lib:int-sub nd i) 1))
 197           (setf (f2cl-lib:fref y (j) ((1 *))) (f2cl-lib:fref y (k) ((1 *))))
 198          label250))
 199      label260
 200       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 201                     ((> i nz) nil)
 202         (tagbody (setf (f2cl-lib:fref y (i) ((1 *))) 0.0) label270))
 203       (go end_label)
 204      label280
 205       (xermsg "SLATEC" "DBESK" "SCALING OPTION, KODE, NOT 1 OR 2" 2 1)
 206       (go end_label)
 207      label290
 208       (xermsg "SLATEC" "DBESK" "ORDER, FNU, LESS THAN ZERO" 2 1)
 209       (go end_label)
 210      label300
 211       (xermsg "SLATEC" "DBESK" "X LESS THAN OR EQUAL TO ZERO" 2 1)
 212       (go end_label)
 213      label310
 214       (xermsg "SLATEC" "DBESK" "N LESS THAN ONE" 2 1)
 215       (go end_label)
 216      label320
 217       (xermsg "SLATEC" "DBESK" "OVERFLOW, FNU OR N TOO LARGE OR X TOO SMALL" 6
 218        1)
 219       (go end_label)
 220      end_label
 221       (return (values nil nil nil nil nil nz)))))
 222
 223 (in-package #:cl-user)
 224 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 225 (eval-when (:load-toplevel :compile-toplevel :execute)
 226   (setf (gethash 'fortran-to-lisp::dbesk fortran-to-lisp::*f2cl-function-info*)
 227           (fortran-to-lisp::make-f2cl-finfo
 228            :arg-types '((double-float) (double-float)
 229                         (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 230                         (simple-array double-float (*))
 231                         (fortran-to-lisp::integer4))
 232            :return-values '(nil nil nil nil nil fortran-to-lisp::nz)
 233            :calls '(fortran-to-lisp::xermsg fortran-to-lisp::dbsknu
 234                     fortran-to-lisp::dbsk1e fortran-to-lisp::dbesk1
 235                     fortran-to-lisp::dbsk0e fortran-to-lisp::dbesk0
 236                     fortran-to-lisp::dasyik fortran-to-lisp::d1mach
 237                     fortran-to-lisp::i1mach))))
 238