src/numerical/slatec/dbie.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       (nbip1 0)
  25       (nbip2 0)
  26       (x3sml 0.0)
  27       (x32sml 0.0)
  28       (xbig 0.0)
  29       (bifcs
  30        (make-array 13
  31                    :element-type 'double-float
  32                    :initial-contents '(-0.01673021647198665 0.10252335834249446
  33                                        0.0017083092507381517
  34                                        1.1862545467744682e-5
  35                                        4.4932907017792135e-8
  36                                        1.0698207143387889e-10
  37                                        1.7480643399771825e-13
  38                                        2.081023107176171e-16
  39                                        1.8849814695665417e-19
  40                                        1.3425779173097804e-22
  41                                        7.715959342965888e-26
  42                                        3.653387961747857e-29
  43                                        1.4497565927953065e-32)))
  44       (bigcs
  45        (make-array 13
  46                    :element-type 'double-float
  47                    :initial-contents '(0.022466223248574523 0.03736477545301955
  48                                        4.4476218957212283e-4
  49                                        2.4708075636329383e-6
  50                                        7.919135339514964e-9
  51                                        1.649807985182778e-11
  52                                        2.4119906664835456e-14
  53                                        2.6103736236091437e-17
  54                                        2.1753082977160324e-20
  55                                        1.4386946400390432e-23
  56                                        7.734912561208347e-27
  57                                        3.4469292033849e-30
  58                                        1.2938919273216e-33)))
  59       (bif2cs
  60        (make-array 15
  61                    :element-type 'double-float
  62                    :initial-contents '(0.0998457269381604 0.47862497786300556
  63                                        0.02515521196043301 5.820693885232646e-4
  64                                        7.499765964437787e-6
  65                                        6.134602870349384e-8
  66                                        3.462753885148063e-10
  67                                        1.4288910080270254e-12
  68                                        4.496270429833464e-15
  69                                        1.1142323065833012e-17
  70                                        2.2304791066175003e-20
  71                                        3.6815778736393144e-23
  72                                        5.096086844933826e-26
  73                                        6.000338692628856e-29
  74                                        6.082749744657067e-32)))
  75       (big2cs
  76        (make-array 15
  77                    :element-type 'double-float
  78                    :initial-contents '(0.03330566214551434 0.16130921512319707
  79                                        0.006319007309613428
  80                                        1.1879045681625174e-4
  81                                        1.3045345886200265e-6
  82                                        9.374125995535217e-9
  83                                        4.745801886747251e-11
  84                                        1.783107265094814e-13
  85                                        5.167591927849581e-16
  86                                        1.1900450838682712e-18
  87                                        2.229828806664035e-21
  88                                        3.465519230276894e-24
  89                                        4.539263363205045e-27
  90                                        5.078849965135223e-30
  91                                        4.910206746965333e-33)))
  92       (bip1cs
  93        (make-array 47
  94                    :element-type 'double-float
  95                    :initial-contents '(-0.08322047477943448
  96                                        0.011461189273711743
  97                                        4.289644071891151e-4
  98                                        -1.4906639379950513e-4
  99                                        -1.307659726787629e-5
 100                                        6.327598396103035e-6
 101                                        -4.2226696982681925e-7
 102                                        -1.914718629865469e-7
 103                                        6.453106284558317e-8
 104                                        -7.844854677139772e-9
 105                                        -9.607721662378508e-10
 106                                        7.000471331644396e-10
 107                                        -1.7731789132814931e-10
 108                                        2.2720894783465238e-11
 109                                        1.654045631397205e-12
 110                                        -1.8517125559292317e-12
 111                                        5.957631247711729e-13
 112                                        -1.2194348147346564e-13
 113                                        1.3347869253513048e-14
 114                                        1.7278311524339746e-15
 115                                        -1.459073201301672e-15
 116                                        4.901031992711582e-16
 117                                        -1.1556545519261548e-16
 118                                        1.909880736707241e-17
 119                                        -1.176896685449218e-18
 120                                        -6.327192514953006e-19
 121                                        3.386183888071536e-19
 122                                        -1.0725825321758626e-19
 123                                        2.599570960561717e-20
 124                                        -4.8477583571081194e-21
 125                                        5.529891398212162e-22
 126                                        4.942166082606947e-23
 127                                        -5.51621219241457e-23
 128                                        2.143756041763255e-23
 129                                        -6.19103133876556e-24
 130                                        1.4629362707391247e-24
 131                                        -2.7918484471059006e-25
 132                                        3.645570316857025e-26
 133                                        5.851182190618871e-28
 134                                        -2.4946950487566512e-27
 135                                        1.0979323980338381e-27
 136                                        -3.4743388345961113e-28
 137                                        9.13734026353497e-29
 138                                        -2.0510352728210628e-29
 139                                        3.797698569854646e-30
 140                                        -4.847945849775557e-31
 141                                        -1.0558306941230714e-32)))
 142       (bip2cs
 143        (make-array 88
 144                    :element-type 'double-float
 145                    :initial-contents '(-0.11359673758598868 0.00413814739478816
 146                                        1.353470622119333e-4
 147                                        1.0427316653015353e-5
 148                                        1.3474954767849909e-6
 149                                        1.6965374054383983e-7
 150                                        -1.0096500865641625e-8
 151                                        -1.6729119493778474e-8
 152                                        -4.5815364485068385e-9
 153                                        3.736681366565548e-10
 154                                        5.766930320145245e-10
 155                                        6.218126508785033e-11
 156                                        -6.329412028274307e-11
 157                                        -1.4915047908598768e-11
 158                                        7.889621394248677e-12
 159                                        2.4960513721857797e-12
 160                                        -1.213007528729166e-12
 161                                        -3.740493910872728e-13
 162                                        2.2377278140321477e-13
 163                                        4.7490296312192465e-14
 164                                        -4.5261607991821224e-14
 165                                        -3.0172271841986073e-15
 166                                        9.105860355875405e-15
 167                                        -9.814923803380705e-16
 168                                        -1.6429400647889466e-15
 169                                        5.533483421427422e-16
 170                                        2.1750479864482656e-16
 171                                        -1.7379236200220656e-16
 172                                        -1.0470023471443715e-18
 173                                        3.9219145986056385e-17
 174                                        -1.1621293686345197e-17
 175                                        -5.402747449175425e-18
 176                                        4.544158212388461e-18
 177                                        -2.8775599625221075e-19
 178                                        -1.0017340927225342e-18
 179                                        4.482393121506837e-19
 180                                        7.613596865490894e-20
 181                                        -1.4448324094881347e-19
 182                                        4.046085944920536e-20
 183                                        2.0321085700338447e-20
 184                                        -1.9602795471446798e-20
 185                                        3.4273038443944824e-21
 186                                        3.7023705853905135e-21
 187                                        -2.687959517204159e-21
 188                                        2.812167846353171e-22
 189                                        6.09339636361778e-22
 190                                        -3.8666621897150843e-22
 191                                        2.5989331253566943e-23
 192                                        9.71943936229385e-23
 193                                        -5.93928178343751e-23
 194                                        3.8864949977113015e-24
 195                                        1.5334307393617272e-23
 196                                        -9.751355520976262e-24
 197                                        9.634064444048946e-25
 198                                        2.384199940020888e-24
 199                                        -1.6896986315019705e-24
 200                                        2.7352715888928363e-25
 201                                        3.566001618540958e-25
 202                                        -3.0234026608258827e-25
 203                                        7.500204160597394e-26
 204                                        4.840328757585139e-26
 205                                        -5.436413765444789e-26
 206                                        1.9281214470820962e-26
 207                                        5.0116355020532654e-27
 208                                        -9.504074458269326e-27
 209                                        4.637264615710198e-27
 210                                        2.1177170704466955e-29
 211                                        -1.5404850268168594e-27
 212                                        1.0387944293201214e-27
 213                                        -1.9890078156915416e-28
 214                                        -2.1022173878658494e-28
 215                                        2.1353099724525795e-28
 216                                        -7.904081074796134e-29
 217                                        -1.6575359960435586e-29
 218                                        3.886834285012411e-29
 219                                        -2.2309237330896867e-29
 220                                        2.777724442017626e-30
 221                                        5.707854347265773e-30
 222                                        -5.1743084445303856e-30
 223                                        1.841328075109584e-30
 224                                        4.4422562390957095e-31
 225                                        -9.85041426396298e-31
 226                                        5.88572013535851e-31
 227                                        -9.763607544042979e-32
 228                                        -1.3581011996074696e-31
 229                                        1.3999743518492413e-31
 230                                        -5.975490454524848e-32
 231                                        -4.039165387542831e-33)))
 232       (atr 8.750690570848434)
 233       (btr -2.0938363213560542)
 234       (first$ nil))
 235   (declare (type (f2cl-lib:integer4) nbif nbig nbif2 nbig2 nbip1 nbip2)
 236            (type (double-float) x3sml x32sml xbig atr btr)
 237            (type (simple-array double-float (13)) bifcs bigcs)
 238            (type (simple-array double-float (15)) bif2cs big2cs)
 239            (type (simple-array double-float (47)) bip1cs)
 240            (type (simple-array double-float (88)) bip2cs)
 241            (type f2cl-lib:logical first$))
 242   (setq first$ f2cl-lib:%true%)
 243   (defun dbie (x)
 244     (declare (type (double-float) x))
 245     (prog ((sqrtx 0.0) (theta 0.0) (xm 0.0) (z 0.0) (dbie 0.0) (eta 0.0f0))
 246       (declare (type (single-float) eta)
 247                (type (double-float) dbie z xm theta sqrtx))
 248       (cond
 249         (first$
 250          (setf eta (* 0.1f0 (f2cl-lib:freal (f2cl-lib:d1mach 3))))
 251          (setf nbif (initds bifcs 13 eta))
 252          (setf nbig (initds bigcs 13 eta))
 253          (setf nbif2 (initds bif2cs 15 eta))
 254          (setf nbig2 (initds big2cs 15 eta))
 255          (setf nbip1 (initds bip1cs 47 eta))
 256          (setf nbip2 (initds bip2cs 88 eta))
 257          (setf x3sml (coerce (expt eta 0.3333f0) 'double-float))
 258          (setf x32sml (* 1.3104 (expt x3sml 2)))
 259          (setf xbig (expt (f2cl-lib:d1mach 2) 0.6666))))
 260       (setf first$ f2cl-lib:%false%)
 261       (if (>= x -1.0) (go label20))
 262       (multiple-value-bind (var-0 var-1 var-2)
 263           (d9aimp x xm theta)
 264         (declare (ignore var-0))
 265         (setf xm var-1)
 266         (setf theta var-2))
 267       (setf dbie (* xm (sin theta)))
 268       (go end_label)
 269      label20
 270       (if (> x 1.0) (go label30))
 271       (setf z 0.0)
 272       (if (> (abs x) x3sml) (setf z (expt x 3)))
 273       (setf dbie
 274               (+ 0.625
 275                  (dcsevl z bifcs nbif)
 276                  (* x (+ 0.4375 (dcsevl z bigcs nbig)))))
 277       (if (> x x32sml)
 278           (setf dbie (* dbie (exp (/ (* -2.0 x (f2cl-lib:fsqrt x)) 3.0)))))
 279       (go end_label)
 280      label30
 281       (if (> x 2.0) (go label40))
 282       (setf z (/ (- (* 2.0 (expt x 3)) 9.0) 7.0))
 283       (setf dbie
 284               (* (exp (/ (* -2.0 x (f2cl-lib:fsqrt x)) 3.0))
 285                  (+ 1.125
 286                     (dcsevl z bif2cs nbif2)
 287                     (* x (+ 0.625 (dcsevl z big2cs nbig2))))))
 288       (go end_label)
 289      label40
 290       (if (> x 4.0) (go label50))
 291       (setf sqrtx (f2cl-lib:fsqrt x))
 292       (setf z (+ (/ atr (* x sqrtx)) btr))
 293       (setf dbie (/ (+ 0.625 (dcsevl z bip1cs nbip1)) (f2cl-lib:fsqrt sqrtx)))
 294       (go end_label)
 295      label50
 296       (setf sqrtx (f2cl-lib:fsqrt x))
 297       (setf z -1.0)
 298       (if (< x xbig) (setf z (- (/ 16.0 (* x sqrtx)) 1.0)))
 299       (setf dbie (/ (+ 0.625 (dcsevl z bip2cs nbip2)) (f2cl-lib:fsqrt sqrtx)))
 300       (go end_label)
 301      end_label
 302       (return (values dbie nil)))))
 303
 304 (in-package #:cl-user)
 305 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 306 (eval-when (:load-toplevel :compile-toplevel :execute)
 307   (setf (gethash 'fortran-to-lisp::dbie fortran-to-lisp::*f2cl-function-info*)
 308           (fortran-to-lisp::make-f2cl-finfo :arg-types '((double-float))
 309                                             :return-values '(nil)
 310                                             :calls '(fortran-to-lisp::dcsevl
 311                                                      fortran-to-lisp::d9aimp
 312                                                      fortran-to-lisp::initds
 313                                                      fortran-to-lisp::d1mach))))
 314