src/numerical/slatec/daie.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 ((naif 0)
  21       (naig 0)
  22       (naip1 0)
  23       (naip2 0)
  24       (x3sml 0.0)
  25       (x32sml 0.0)
  26       (xbig 0.0)
  27       (aifcs
  28        (make-array 13
  29                    :element-type 'double-float
  30                    :initial-contents '(-0.03797135849667 0.05919188853726364
  31                                        9.862928057727998e-4
  32                                        6.848843819076567e-6
  33                                        2.5942025962194713e-8
  34                                        6.176612774081375e-11
  35                                        1.0092454172466118e-13
  36                                        1.2014792511179938e-16
  37                                        1.0882945588716992e-19
  38                                        7.751377219668488e-23
  39                                        4.4548112037175636e-26
  40                                        2.1092845231692343e-29
  41                                        8.370173591074134e-33)))
  42       (aigcs
  43        (make-array 13
  44                    :element-type 'double-float
  45                    :initial-contents '(0.018152365581161272
  46                                        0.021572563166010757
  47                                        2.567835698748325e-4
  48                                        1.4265214119792405e-6
  49                                        4.572114920018043e-9
  50                                        9.52517084356471e-12
  51                                        1.3925634605771398e-14
  52                                        1.5070999142762378e-17
  53                                        1.2559148312567778e-20
  54                                        8.306307377082133e-24
  55                                        4.465753849371857e-27
  56                                        1.9900855034518868e-30
  57                                        7.4702885256533335e-34)))
  58       (aip1cs
  59        (make-array 57
  60                    :element-type 'double-float
  61                    :initial-contents '(-0.021469518589105386
  62                                        -0.0075353825350433015
  63                                        5.971527949026381e-4
  64                                        -7.28325125420761e-5
  65                                        1.1102971307392997e-5
  66                                        -1.950386152284406e-6
  67                                        3.786973885159515e-7
  68                                        -7.929675297350979e-8
  69                                        1.762247638674256e-8
  70                                        -4.110767539667195e-9
  71                                        9.984770057857892e-10
  72                                        -2.5100932513871223e-10
  73                                        6.500501929860696e-11
  74                                        -1.7278184053936166e-11
  75                                        4.6993788428245126e-12
  76                                        -1.304675656297744e-12
  77                                        3.6896984784626787e-13
  78                                        -1.0610872066468062e-13
  79                                        3.0984143848781874e-14
  80                                        -9.17490807982414e-15
  81                                        2.7520491403472108e-15
  82                                        -8.353750115922047e-16
  83                                        2.563931129357935e-16
  84                                        -7.950633762598855e-17
  85                                        2.48928363460307e-17
  86                                        -7.864326933928736e-18
  87                                        2.5056873114399757e-18
  88                                        -8.047420364163909e-19
  89                                        2.604097118952054e-19
  90                                        -8.486954164056412e-20
  91                                        2.784706882142338e-20
  92                                        -9.195858953498614e-21
  93                                        3.055304318374239e-21
  94                                        -1.0210354554794778e-21
  95                                        3.431118190743758e-22
  96                                        -1.1591293417977495e-22
  97                                        3.935772844200256e-23
  98                                        -1.3428809802967176e-23
  99                                        4.6032878835200026e-24
 100                                        -1.5850439270040642e-24
 101                                        5.481275667729676e-25
 102                                        -1.9033493718550473e-25
 103                                        6.635682302374009e-26
 104                                        -2.3223116500263143e-26
 105                                        8.157640113429179e-27
 106                                        -2.8758242406329004e-27
 107                                        1.0173294509429014e-27
 108                                        -3.6108791087422165e-28
 109                                        1.2857885403639934e-28
 110                                        -4.5929010373785476e-29
 111                                        1.6455970338207138e-29
 112                                        -5.913421299843502e-30
 113                                        2.131057006604993e-30
 114                                        -7.701158157787598e-31
 115                                        2.7905333079689304e-31
 116                                        -1.013807715111284e-31
 117                                        3.692580158719624e-32)))
 118       (aip2cs
 119        (make-array 37
 120                    :element-type 'double-float
 121                    :initial-contents '(-0.0017431449692937551
 122                                        -0.0016789385432554166
 123                                        3.5965340335216605e-5
 124                                        -1.3808186027392284e-6
 125                                        7.411228077315053e-8
 126                                        -5.00238203900133e-9
 127                                        4.0069391741718425e-10
 128                                        -3.6733124279590504e-11
 129                                        3.760344395923738e-12
 130                                        -4.2232133271874755e-13
 131                                        5.135094540336571e-14
 132                                        -6.690958503904776e-15
 133                                        9.266675456412906e-16
 134                                        -1.3551438241607058e-16
 135                                        2.0811549631283098e-17
 136                                        -3.3411649915917686e-18
 137                                        5.5857858458592435e-19
 138                                        -9.692190401523652e-20
 139                                        1.740457001288932e-20
 140                                        -3.226409797311304e-21
 141                                        6.160744711066252e-22
 142                                        -1.2093634798249005e-22
 143                                        2.436327633101381e-23
 144                                        -5.029142214974575e-24
 145                                        1.062241755436357e-24
 146                                        -2.2928428489598924e-25
 147                                        5.051817339295037e-26
 148                                        -1.134981237144124e-26
 149                                        2.5976556598560697e-27
 150                                        -6.051246215429395e-28
 151                                        1.4335977796677281e-28
 152                                        -3.4514775706089996e-29
 153                                        8.438751902136468e-30
 154                                        -2.0939614229818816e-30
 155                                        5.270088734789455e-31
 156                                        -1.3445743301455338e-31
 157                                        3.475709645266011e-32)))
 158       (first$ nil))
 159   (declare (type (f2cl-lib:integer4) naif naig naip1 naip2)
 160            (type (double-float) x3sml x32sml xbig)
 161            (type (simple-array double-float (13)) aifcs aigcs)
 162            (type (simple-array double-float (57)) aip1cs)
 163            (type (simple-array double-float (37)) aip2cs)
 164            (type f2cl-lib:logical first$))
 165   (setq first$ f2cl-lib:%true%)
 166   (defun daie (x)
 167     (declare (type (double-float) x))
 168     (prog ((sqrtx 0.0) (theta 0.0) (xm 0.0) (z 0.0) (daie 0.0) (eta 0.0f0))
 169       (declare (type (single-float) eta)
 170                (type (double-float) daie z xm theta sqrtx))
 171       (cond
 172         (first$
 173          (setf eta (* 0.1f0 (f2cl-lib:freal (f2cl-lib:d1mach 3))))
 174          (setf naif (initds aifcs 13 eta))
 175          (setf naig (initds aigcs 13 eta))
 176          (setf naip1 (initds aip1cs 57 eta))
 177          (setf naip2 (initds aip2cs 37 eta))
 178          (setf x3sml (coerce (expt eta 0.3333f0) 'double-float))
 179          (setf x32sml (* 1.3104 (expt x3sml 2)))
 180          (setf xbig (expt (f2cl-lib:d1mach 2) 0.6666))))
 181       (setf first$ f2cl-lib:%false%)
 182       (if (>= x -1.0) (go label20))
 183       (multiple-value-bind (var-0 var-1 var-2)
 184           (d9aimp x xm theta)
 185         (declare (ignore var-0))
 186         (setf xm var-1)
 187         (setf theta var-2))
 188       (setf daie (* xm (cos theta)))
 189       (go end_label)
 190      label20
 191       (if (> x 1.0) (go label30))
 192       (setf z 0.0)
 193       (if (> (abs x) x3sml) (setf z (expt x 3)))
 194       (setf daie
 195               (+ 0.375
 196                  (- (dcsevl z aifcs naif)
 197                     (* x (+ 0.25 (dcsevl z aigcs naig))))))
 198       (if (> x x32sml)
 199           (setf daie (* daie (exp (/ (* 2.0 x (f2cl-lib:fsqrt x)) 3.0)))))
 200       (go end_label)
 201      label30
 202       (if (> x 4.0) (go label40))
 203       (setf sqrtx (f2cl-lib:fsqrt x))
 204       (setf z (/ (- (/ 16.0 (* x sqrtx)) 9.0) 7.0))
 205       (setf daie
 206               (/ (+ 0.28125 (dcsevl z aip1cs naip1)) (f2cl-lib:fsqrt sqrtx)))
 207       (go end_label)
 208      label40
 209       (setf sqrtx (f2cl-lib:fsqrt x))
 210       (setf z -1.0)
 211       (if (< x xbig) (setf z (- (/ 16.0 (* x sqrtx)) 1.0)))
 212       (setf daie
 213               (/ (+ 0.28125 (dcsevl z aip2cs naip2)) (f2cl-lib:fsqrt sqrtx)))
 214       (go end_label)
 215      end_label
 216       (return (values daie nil)))))
 217
 218 (in-package #:cl-user)
 219 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 220 (eval-when (:load-toplevel :compile-toplevel :execute)
 221   (setf (gethash 'fortran-to-lisp::daie fortran-to-lisp::*f2cl-function-info*)
 222           (fortran-to-lisp::make-f2cl-finfo :arg-types '((double-float))
 223                                             :return-values '(nil)
 224                                             :calls '(fortran-to-lisp::dcsevl
 225                                                      fortran-to-lisp::d9aimp
 226                                                      fortran-to-lisp::initds
 227                                                      fortran-to-lisp::d1mach))))
 228