src/numerical/slatec/dbesy.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 '(70 100))))
  24   (declare (type (simple-array f2cl-lib:integer4 (2)) nulim))
  25   (defun dbesy (x fnu n y)
  26     (declare (type (simple-array double-float (*)) y)
  27              (type (f2cl-lib:integer4) n)
  28              (type (double-float) fnu x))
  29     (prog ((w (make-array 2 :element-type 'double-float))
  30            (wk (make-array 7 :element-type 'double-float)) (azn 0.0) (cn 0.0)
  31            (dnu 0.0) (elim 0.0) (flgjy 0.0) (fn 0.0) (ran 0.0) (s 0.0) (s1 0.0)
  32            (s2 0.0) (tm 0.0) (trx 0.0) (w2n 0.0) (xlim 0.0) (xxn 0.0) (i 0)
  33            (iflw 0) (j 0) (nb 0) (nd 0) (nn 0) (nud 0))
  34       (declare (type (f2cl-lib:integer4) nud nn nd nb j iflw i)
  35                (type (simple-array double-float (7)) wk)
  36                (type (simple-array double-float (2)) w)
  37                (type (double-float) xxn xlim w2n trx tm s2 s1 s ran fn flgjy
  38                                     elim dnu cn azn))
  39       (setf nn (f2cl-lib:int-sub (f2cl-lib:i1mach 15)))
  40       (setf elim (* 2.303 (- (* nn (f2cl-lib:d1mach 5)) 3.0)))
  41       (setf xlim (* (f2cl-lib:d1mach 1) 1000.0))
  42       (if (< fnu 0.0) (go label140))
  43       (if (<= x 0.0) (go label150))
  44       (if (< x xlim) (go label170))
  45       (if (< n 1) (go label160))
  46       (setf nd n)
  47       (setf nud (f2cl-lib:int fnu))
  48       (setf dnu (- fnu nud))
  49       (setf nn (min (the f2cl-lib:integer4 2) (the f2cl-lib:integer4 nd)))
  50       (setf fn (- (+ fnu n) 1))
  51       (if (< fn 2.0) (go label100))
  52       (setf xxn (/ x fn))
  53       (setf w2n (- 1.0 (* xxn xxn)))
  54       (if (<= w2n 0.0) (go label10))
  55       (setf ran (f2cl-lib:fsqrt w2n))
  56       (setf azn (- (f2cl-lib:flog (/ (+ 1.0 ran) xxn)) ran))
  57       (setf cn (* fn azn))
  58       (if (> cn elim) (go label170))
  59      label10
  60       (if (< nud (f2cl-lib:fref nulim (nn) ((1 2)))) (go label20))
  61       (setf flgjy -1.0)
  62       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7)
  63           (dasyjy #'dyairy x fnu flgjy nn y wk iflw)
  64         (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6))
  65         (setf iflw var-7))
  66       (if (/= iflw 0) (go label170))
  67       (if (= nn 1) (go end_label))
  68       (setf trx (/ 2.0 x))
  69       (setf tm (/ (+ fnu fnu 2.0) x))
  70       (go label80)
  71      label20
  72       (if (/= dnu 0.0) (go label30))
  73       (setf s1 (dbesy0 x))
  74       (if (and (= nud 0) (= nd 1)) (go label70))
  75       (setf s2 (dbesy1 x))
  76       (go label40)
  77      label30
  78       (setf nb 2)
  79       (if (and (= nud 0) (= nd 1)) (setf nb 1))
  80       (dbsynu x dnu nb w)
  81       (setf s1 (f2cl-lib:fref w (1) ((1 2))))
  82       (if (= nb 1) (go label70))
  83       (setf s2 (f2cl-lib:fref w (2) ((1 2))))
  84      label40
  85       (setf trx (/ 2.0 x))
  86       (setf tm (/ (+ dnu dnu 2.0) x))
  87       (if (= nd 1) (setf nud (f2cl-lib:int-sub nud 1)))
  88       (if (> nud 0) (go label50))
  89       (if (> nd 1) (go label70))
  90       (setf s1 s2)
  91       (go label70)
  92      label50
  93       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  94                     ((> i nud) nil)
  95         (tagbody
  96           (setf s s2)
  97           (setf s2 (- (* tm s2) s1))
  98           (setf s1 s)
  99           (setf tm (+ tm trx))
 100          label60))
 101       (if (= nd 1) (setf s1 s2))
 102      label70
 103       (setf (f2cl-lib:fref y (1) ((1 *))) s1)
 104       (if (= nd 1) (go end_label))
 105       (setf (f2cl-lib:fref y (2) ((1 *))) s2)
 106      label80
 107       (if (= nd 2) (go end_label))
 108       (f2cl-lib:fdo (i 3 (f2cl-lib:int-add i 1))
 109                     ((> i nd) nil)
 110         (tagbody
 111           (setf (f2cl-lib:fref y (i) ((1 *)))
 112                   (- (* tm (f2cl-lib:fref y ((f2cl-lib:int-sub i 1)) ((1 *))))
 113                      (f2cl-lib:fref y ((f2cl-lib:int-sub i 2)) ((1 *)))))
 114           (setf tm (+ tm trx))
 115          label90))
 116       (go end_label)
 117      label100
 118       (if (<= fn 1.0) (go label110))
 119       (if (> (* (- fn) (- (f2cl-lib:flog x) 0.693)) elim) (go label170))
 120      label110
 121       (if (= dnu 0.0) (go label120))
 122       (dbsynu x fnu nd y)
 123       (go end_label)
 124      label120
 125       (setf j nud)
 126       (if (= j 1) (go label130))
 127       (setf j (f2cl-lib:int-add j 1))
 128       (setf (f2cl-lib:fref y (j) ((1 *))) (dbesy0 x))
 129       (if (= nd 1) (go end_label))
 130       (setf j (f2cl-lib:int-add j 1))
 131      label130
 132       (setf (f2cl-lib:fref y (j) ((1 *))) (dbesy1 x))
 133       (if (= nd 1) (go end_label))
 134       (setf trx (/ 2.0 x))
 135       (setf tm trx)
 136       (go label80)
 137      label140
 138       (xermsg "SLATEC" "DBESY" "ORDER, FNU, LESS THAN ZERO" 2 1)
 139       (go end_label)
 140      label150
 141       (xermsg "SLATEC" "DBESY" "X LESS THAN OR EQUAL TO ZERO" 2 1)
 142       (go end_label)
 143      label160
 144       (xermsg "SLATEC" "DBESY" "N LESS THAN ONE" 2 1)
 145       (go end_label)
 146      label170
 147       (xermsg "SLATEC" "DBESY" "OVERFLOW, FNU OR N TOO LARGE OR X TOO SMALL" 6
 148        1)
 149       (go end_label)
 150      end_label
 151       (return (values nil nil nil nil)))))
 152
 153 (in-package #:cl-user)
 154 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 155 (eval-when (:load-toplevel :compile-toplevel :execute)
 156   (setf (gethash 'fortran-to-lisp::dbesy fortran-to-lisp::*f2cl-function-info*)
 157           (fortran-to-lisp::make-f2cl-finfo
 158            :arg-types '((double-float) (double-float)
 159                         (fortran-to-lisp::integer4)
 160                         (simple-array double-float (*)))
 161            :return-values '(nil nil nil nil)
 162            :calls '(fortran-to-lisp::xermsg fortran-to-lisp::dbsynu
 163                     fortran-to-lisp::dbesy1 fortran-to-lisp::dbesy0
 164                     fortran-to-lisp::dasyjy fortran-to-lisp::d1mach
 165                     fortran-to-lisp::i1mach))))
 166