src/numerical/slatec/dbsk1e.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 ((ntk1 0)
  21       (ntak1 0)
  22       (ntak12 0)
  23       (xmin 0.0)
  24       (xsml 0.0)
  25       (bk1cs
  26        (make-array 16
  27                    :element-type 'double-float
  28                    :initial-contents '(0.02530022733894777 -0.3531559607765449
  29                                        -0.12261118082265715
  30                                        -0.006975723859639864
  31                                        -1.730288957513052e-4
  32                                        -2.4334061415659684e-6
  33                                        -2.213387630734726e-8
  34                                        -1.4114883926335278e-10
  35                                        -6.666901694199329e-13
  36                                        -2.427449850519366e-15
  37                                        -7.023863479386288e-18
  38                                        -1.6543275155100994e-20
  39                                        -3.233834745994449e-23
  40                                        -5.331275052926527e-26
  41                                        -7.513040716215723e-29
  42                                        -9.155085717654187e-32)))
  43       (ak1cs
  44        (make-array 38
  45                    :element-type 'double-float
  46                    :initial-contents '(0.2744313406973883 0.07571989953199368
  47                                        -0.0014410515564754062
  48                                        6.650116955125748e-5
  49                                        -4.369984709520141e-6
  50                                        3.5402774997630525e-7
  51                                        -3.311163779293292e-8
  52                                        3.4459775819010535e-9
  53                                        -3.898932347475427e-10
  54                                        4.720819750465836e-11
  55                                        -6.047835662875356e-12
  56                                        8.128494874865875e-13
  57                                        -1.138694574714789e-13
  58                                        1.654035840846228e-14
  59                                        -2.4809025677068848e-15
  60                                        3.8292378907024097e-16
  61                                        -6.064734104001242e-17
  62                                        9.832425623264862e-18
  63                                        -1.628416873828438e-18
  64                                        2.750153649675262e-19
  65                                        -4.728966646395325e-20
  66                                        8.268150002810994e-21
  67                                        -1.4681405136624957e-21
  68                                        2.6447639269208245e-22
  69                                        -4.829015756485639e-23
  70                                        8.929302074361012e-24
  71                                        -1.6708397168972516e-24
  72                                        3.1616456034040695e-25
  73                                        -6.046205531227498e-26
  74                                        1.1678798942042733e-26
  75                                        -2.2773741582653997e-27
  76                                        4.481109730077368e-28
  77                                        -8.893288476902019e-29
  78                                        1.7794680018850274e-29
  79                                        -3.58845559673291e-30
  80                                        7.290629049269426e-31
  81                                        -1.4918449845546228e-31
  82                                        3.0736573872934276e-32)))
  83       (ak12cs
  84        (make-array 33
  85                    :element-type 'double-float
  86                    :initial-contents '(0.06379308343739001 0.02832887813049721
  87                                        -2.4753706739052506e-4
  88                                        5.771972451607249e-6
  89                                        -2.0689392195365484e-7
  90                                        9.739983441381804e-9
  91                                        -5.585336140380625e-10
  92                                        3.7329966340461855e-11
  93                                        -2.8250519610232256e-12
  94                                        2.372019002484144e-13
  95                                        -2.176677387991754e-14
  96                                        2.1579141616160325e-15
  97                                        -2.290196930718269e-16
  98                                        2.582885729823275e-17
  99                                        -3.076752641268463e-18
 100                                        3.8514877212804914e-19
 101                                        -5.044794897641529e-20
 102                                        6.888673850418544e-21
 103                                        -9.775041541950119e-22
 104                                        1.4374162185238365e-22
 105                                        -2.1850594973443474e-23
 106                                        3.426245621809221e-24
 107                                        -5.531064394246408e-25
 108                                        9.176601505685995e-26
 109                                        -1.562287203618025e-26
 110                                        2.725419375484333e-27
 111                                        -4.865674910074828e-28
 112                                        8.879388552723502e-29
 113                                        -1.6545859180392576e-29
 114                                        3.1451113213578485e-30
 115                                        -6.092998312193127e-31
 116                                        1.2020219393698158e-31
 117                                        -2.412930801459409e-32)))
 118       (first$ nil))
 119   (declare (type (f2cl-lib:integer4) ntk1 ntak1 ntak12)
 120            (type (double-float) xmin xsml)
 121            (type (simple-array double-float (16)) bk1cs)
 122            (type (simple-array double-float (38)) ak1cs)
 123            (type (simple-array double-float (33)) ak12cs)
 124            (type f2cl-lib:logical first$))
 125   (setq first$ f2cl-lib:%true%)
 126   (defun dbsk1e (x)
 127     (declare (type (double-float) x))
 128     (prog ((y 0.0) (dbsk1e 0.0) (eta 0.0f0))
 129       (declare (type (single-float) eta) (type (double-float) dbsk1e y))
 130       (cond
 131         (first$
 132          (setf eta (* 0.1f0 (f2cl-lib:freal (f2cl-lib:d1mach 3))))
 133          (setf ntk1 (initds bk1cs 16 eta))
 134          (setf ntak1 (initds ak1cs 38 eta))
 135          (setf ntak12 (initds ak12cs 33 eta))
 136          (setf xmin
 137                  (exp
 138                   (+
 139                    (max (f2cl-lib:flog (f2cl-lib:d1mach 1))
 140                         (- (f2cl-lib:flog (f2cl-lib:d1mach 2))))
 141                    0.01)))
 142          (setf xsml (f2cl-lib:fsqrt (* 4.0 (f2cl-lib:d1mach 3))))))
 143       (setf first$ f2cl-lib:%false%)
 144       (if (<= x 0.0) (xermsg "SLATEC" "DBSK1E" "X IS ZERO OR NEGATIVE" 2 2))
 145       (if (> x 2.0) (go label20))
 146       (if (< x xmin) (xermsg "SLATEC" "DBSK1E" "X SO SMALL K1 OVERFLOWS" 3 2))
 147       (setf y 0.0)
 148       (if (> x xsml) (setf y (* x x)))
 149       (setf dbsk1e
 150               (* (exp x)
 151                  (+ (* (f2cl-lib:flog (* 0.5 x)) (dbesi1 x))
 152                     (/ (+ 0.75 (dcsevl (- (* 0.5 y) 1.0) bk1cs ntk1)) x))))
 153       (go end_label)
 154      label20
 155       (if (<= x 8.0)
 156           (setf dbsk1e
 157                   (/ (+ 1.25 (dcsevl (/ (- (/ 16.0 x) 5.0) 3.0) ak1cs ntak1))
 158                      (f2cl-lib:fsqrt x))))
 159       (if (> x 8.0)
 160           (setf dbsk1e
 161                   (/ (+ 1.25 (dcsevl (- (/ 16.0 x) 1.0) ak12cs ntak12))
 162                      (f2cl-lib:fsqrt x))))
 163       (go end_label)
 164      end_label
 165       (return (values dbsk1e nil)))))
 166
 167 (in-package #:cl-user)
 168 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 169 (eval-when (:load-toplevel :compile-toplevel :execute)
 170   (setf (gethash 'fortran-to-lisp::dbsk1e
 171                  fortran-to-lisp::*f2cl-function-info*)
 172           (fortran-to-lisp::make-f2cl-finfo :arg-types '((double-float))
 173                                             :return-values '(nil)
 174                                             :calls '(fortran-to-lisp::dcsevl
 175                                                      fortran-to-lisp::dbesi1
 176                                                      fortran-to-lisp::xermsg
 177                                                      fortran-to-lisp::initds
 178                                                      fortran-to-lisp::d1mach))))
 179