src/numerical/slatec/dbi.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 ((nbif 0)
  21       (nbig 0)
  22       (nbif2 0)
  23       (nbig2 0)
  24       (x3sml 0.0)
  25       (xmax 0.0)
  26       (bifcs
  27        (make-array 13
  28                    :element-type 'double-float
  29                    :initial-contents '(-0.01673021647198665 0.10252335834249446
  30                                        0.0017083092507381517
  31                                        1.1862545467744682e-5
  32                                        4.4932907017792135e-8
  33                                        1.0698207143387889e-10
  34                                        1.7480643399771825e-13
  35                                        2.081023107176171e-16
  36                                        1.8849814695665417e-19
  37                                        1.3425779173097804e-22
  38                                        7.715959342965888e-26
  39                                        3.653387961747857e-29
  40                                        1.4497565927953065e-32)))
  41       (bigcs
  42        (make-array 13
  43                    :element-type 'double-float
  44                    :initial-contents '(0.022466223248574523 0.03736477545301955
  45                                        4.4476218957212283e-4
  46                                        2.4708075636329383e-6
  47                                        7.919135339514964e-9
  48                                        1.649807985182778e-11
  49                                        2.4119906664835456e-14
  50                                        2.6103736236091437e-17
  51                                        2.1753082977160324e-20
  52                                        1.4386946400390432e-23
  53                                        7.734912561208347e-27
  54                                        3.4469292033849e-30
  55                                        1.2938919273216e-33)))
  56       (bif2cs
  57        (make-array 15
  58                    :element-type 'double-float
  59                    :initial-contents '(0.0998457269381604 0.47862497786300556
  60                                        0.02515521196043301 5.820693885232646e-4
  61                                        7.499765964437787e-6
  62                                        6.134602870349384e-8
  63                                        3.462753885148063e-10
  64                                        1.4288910080270254e-12
  65                                        4.496270429833464e-15
  66                                        1.1142323065833012e-17
  67                                        2.2304791066175003e-20
  68                                        3.6815778736393144e-23
  69                                        5.096086844933826e-26
  70                                        6.000338692628856e-29
  71                                        6.082749744657067e-32)))
  72       (big2cs
  73        (make-array 15
  74                    :element-type 'double-float
  75                    :initial-contents '(0.03330566214551434 0.16130921512319707
  76                                        0.006319007309613428
  77                                        1.1879045681625174e-4
  78                                        1.3045345886200265e-6
  79                                        9.374125995535217e-9
  80                                        4.745801886747251e-11
  81                                        1.783107265094814e-13
  82                                        5.167591927849581e-16
  83                                        1.1900450838682712e-18
  84                                        2.229828806664035e-21
  85                                        3.465519230276894e-24
  86                                        4.539263363205045e-27
  87                                        5.078849965135223e-30
  88                                        4.910206746965333e-33)))
  89       (first$ nil))
  90   (declare (type (f2cl-lib:integer4) nbif nbig nbif2 nbig2)
  91            (type (double-float) x3sml xmax)
  92            (type (simple-array double-float (13)) bifcs bigcs)
  93            (type (simple-array double-float (15)) bif2cs big2cs)
  94            (type f2cl-lib:logical first$))
  95   (setq first$ f2cl-lib:%true%)
  96   (defun dbi (x)
  97     (declare (type (double-float) x))
  98     (prog ((theta 0.0) (xm 0.0) (z 0.0) (dbi 0.0) (eta 0.0f0))
  99       (declare (type (single-float) eta) (type (double-float) dbi z xm theta))
 100       (cond
 101         (first$
 102          (setf eta (* 0.1f0 (f2cl-lib:freal (f2cl-lib:d1mach 3))))
 103          (setf nbif (initds bifcs 13 eta))
 104          (setf nbig (initds bigcs 13 eta))
 105          (setf nbif2 (initds bif2cs 15 eta))
 106          (setf nbig2 (initds big2cs 15 eta))
 107          (setf x3sml (coerce (expt eta 0.3333f0) 'double-float))
 108          (setf xmax
 109                  (expt (* 1.5f0 (f2cl-lib:flog (f2cl-lib:d1mach 2))) 0.6666))))
 110       (setf first$ f2cl-lib:%false%)
 111       (if (>= x -1.0) (go label20))
 112       (multiple-value-bind (var-0 var-1 var-2)
 113           (d9aimp x xm theta)
 114         (declare (ignore var-0))
 115         (setf xm var-1)
 116         (setf theta var-2))
 117       (setf dbi (* xm (sin theta)))
 118       (go end_label)
 119      label20
 120       (if (> x 1.0) (go label30))
 121       (setf z 0.0)
 122       (if (> (abs x) x3sml) (setf z (expt x 3)))
 123       (setf dbi
 124               (+ 0.625
 125                  (dcsevl z bifcs nbif)
 126                  (* x (+ 0.4375 (dcsevl z bigcs nbig)))))
 127       (go end_label)
 128      label30
 129       (if (> x 2.0) (go label40))
 130       (setf z (/ (- (* 2.0 (expt x 3)) 9.0) 7.0))
 131       (setf dbi
 132               (+ 1.125
 133                  (dcsevl z bif2cs nbif2)
 134                  (* x (+ 0.625 (dcsevl z big2cs nbig2)))))
 135       (go end_label)
 136      label40
 137       (if (> x xmax) (xermsg "SLATEC" "DBI" "X SO BIG THAT BI OVERFLOWS" 1 2))
 138       (setf dbi (* (dbie x) (exp (/ (* 2.0 x (f2cl-lib:fsqrt x)) 3.0))))
 139       (go end_label)
 140      end_label
 141       (return (values dbi nil)))))
 142
 143 (in-package #:cl-user)
 144 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 145 (eval-when (:load-toplevel :compile-toplevel :execute)
 146   (setf (gethash 'fortran-to-lisp::dbi fortran-to-lisp::*f2cl-function-info*)
 147           (fortran-to-lisp::make-f2cl-finfo :arg-types '((double-float))
 148                                             :return-values '(nil)
 149                                             :calls '(fortran-to-lisp::dbie
 150                                                      fortran-to-lisp::xermsg
 151                                                      fortran-to-lisp::dcsevl
 152                                                      fortran-to-lisp::d9aimp
 153                                                      fortran-to-lisp::initds
 154                                                      fortran-to-lisp::d1mach))))
 155