src/numerical/slatec/derfc.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 ((nterf 0)
  21       (nterfc 0)
  22       (nterc2 0)
  23       (xsml 0.0)
  24       (xmax 0.0)
  25       (sqeps 0.0)
  26       (erfcs
  27        (make-array 21
  28                    :element-type 'double-float
  29                    :initial-contents '(-0.049046121234691806
  30                                        -0.14226120510371365
  31                                        0.010035582187599796
  32                                        -5.768764699767485e-4
  33                                        2.741993125219606e-5
  34                                        -1.1043175507344507e-6
  35                                        3.8488755420345036e-8
  36                                        -1.1808582533875466e-9
  37                                        3.2334215826050907e-11
  38                                        -7.991015947004549e-13
  39                                        1.7990725113961456e-14
  40                                        -3.718635487818693e-16
  41                                        7.103599003714253e-18
  42                                        -1.2612455119155226e-19
  43                                        2.0916406941769294e-21
  44                                        -3.2539731029314073e-23
  45                                        4.766867209797675e-25
  46                                        -6.598012078285134e-27
  47                                        8.655011469963763e-29
  48                                        -1.0788925177498064e-30
  49                                        1.2811883993017003e-32)))
  50       (erc2cs
  51        (make-array 49
  52                    :element-type 'double-float
  53                    :initial-contents '(-0.0696013466023095 -0.04110133936262089
  54                                        0.003914495866689627
  55                                        -4.906395650548979e-4
  56                                        7.157479001377036e-5
  57                                        -1.1530716341312328e-5
  58                                        1.9946705902019974e-6
  59                                        -3.642666471599223e-7
  60                                        6.944372610005012e-8
  61                                        -1.371220902104366e-8
  62                                        2.7883896610071373e-9
  63                                        -5.814164724331161e-10
  64                                        1.2389204917527532e-10
  65                                        -2.6906391453067435e-11
  66                                        5.942614350847911e-12
  67                                        -1.3323867357581197e-12
  68                                        3.0280468061771323e-13
  69                                        -6.966648814941033e-14
  70                                        1.620854541053923e-14
  71                                        -3.809934465250492e-15
  72                                        9.040487815978831e-16
  73                                        -2.1640061950896072e-16
  74                                        5.222102233995855e-17
  75                                        -1.2697296023645554e-17
  76                                        3.1091455042761977e-18
  77                                        -7.663762920320386e-19
  78                                        1.9008192513627452e-19
  79                                        -4.7422072790690393e-20
  80                                        1.1896492000765284e-20
  81                                        -3.0000355903257804e-21
  82                                        7.602993453043246e-22
  83                                        -1.9359094476068728e-22
  84                                        4.951399124773338e-23
  85                                        -1.2718074813363719e-23
  86                                        3.2800496004695132e-24
  87                                        -8.492320176822897e-25
  88                                        2.2069178928075603e-25
  89                                        -5.755617245696529e-26
  90                                        1.5061915336392342e-26
  91                                        -3.954502959018797e-27
  92                                        1.0415297041515009e-27
  93                                        -2.751487795278765e-28
  94                                        7.290058205497557e-29
  95                                        -1.936939645915948e-29
  96                                        5.1603571120514875e-30
  97                                        -1.3784193221930942e-30
  98                                        3.691326793107069e-31
  99                                        -9.909389590624365e-32
 100                                        2.666491705195388e-32)))
 101       (erfccs
 102        (make-array 59
 103                    :element-type 'double-float
 104                    :initial-contents '(0.07151793102029248
 105                                        -0.026532434337606717
 106                                        0.0017111539779208558
 107                                        -1.6375166345851787e-4
 108                                        1.9871293500552038e-5
 109                                        -2.843712412766555e-6
 110                                        4.6061613089631305e-7
 111                                        -8.227753025879209e-8
 112                                        1.5921418727709012e-8
 113                                        -3.295071362252843e-9
 114                                        7.223439760400556e-10
 115                                        -1.6648558133987297e-10
 116                                        4.010392588237665e-11
 117                                        -1.004816214425731e-11
 118                                        2.608275913300334e-12
 119                                        -6.991110560404025e-13
 120                                        1.9294923332617072e-13
 121                                        -5.470131188754331e-14
 122                                        1.5896633097626975e-14
 123                                        -4.726893980197555e-15
 124                                        1.4358733767849847e-15
 125                                        -4.449510561817358e-16
 126                                        1.4048108847682335e-16
 127                                        -4.5138183877642106e-17
 128                                        1.474521541045133e-17
 129                                        -4.8926214069457765e-18
 130                                        1.6476121414106467e-18
 131                                        -5.626817176329408e-19
 132                                        1.9474433822320786e-19
 133                                        -6.82630564294842e-20
 134                                        2.4219888872986492e-20
 135                                        -8.693414133503071e-21
 136                                        3.1551803462280855e-21
 137                                        -1.1573723240496087e-21
 138                                        4.288947161605654e-22
 139                                        -1.6050307420576167e-22
 140                                        6.063298757453803e-23
 141                                        -2.3114042516979585e-23
 142                                        8.888778540661885e-24
 143                                        -3.447260576651376e-24
 144                                        1.347865460206965e-24
 145                                        -5.311794071125021e-25
 146                                        2.109341058619783e-25
 147                                        -8.438365587923789e-26
 148                                        3.399982524945209e-26
 149                                        -1.3794523880732422e-26
 150                                        5.6344903118332525e-27
 151                                        -2.3164904344770655e-27
 152                                        9.58446284460181e-28
 153                                        -3.9907228803301096e-28
 154                                        1.6721292259444773e-28
 155                                        -7.045991522766014e-29
 156                                        2.9797684028642063e-29
 157                                        -1.2625224664606192e-29
 158                                        5.395438704542488e-30
 159                                        -2.380992882531459e-30
 160                                        1.0990528301027615e-30
 161                                        -4.867713741644966e-31
 162                                        1.5258772641103575e-31)))
 163       (sqrtpi 1.772453850905516)
 164       (first$ nil))
 165   (declare (type (f2cl-lib:integer4) nterf nterfc nterc2)
 166            (type (double-float) xsml xmax sqeps sqrtpi)
 167            (type (simple-array double-float (21)) erfcs)
 168            (type (simple-array double-float (49)) erc2cs)
 169            (type (simple-array double-float (59)) erfccs)
 170            (type f2cl-lib:logical first$))
 171   (setq first$ f2cl-lib:%true%)
 172   (defun derfc (x)
 173     (declare (type (double-float) x))
 174     (prog ((txmax 0.0) (y 0.0) (derfc 0.0) (eta 0.0f0))
 175       (declare (type (single-float) eta) (type (double-float) derfc y txmax))
 176       (cond
 177         (first$
 178          (setf eta (* 0.1f0 (f2cl-lib:freal (f2cl-lib:d1mach 3))))
 179          (setf nterf (initds erfcs 21 eta))
 180          (setf nterfc (initds erfccs 59 eta))
 181          (setf nterc2 (initds erc2cs 49 eta))
 182          (setf xsml
 183                  (-
 184                   (f2cl-lib:fsqrt
 185                    (- (f2cl-lib:flog (* sqrtpi (f2cl-lib:d1mach 3)))))))
 186          (setf txmax
 187                  (f2cl-lib:fsqrt
 188                   (- (f2cl-lib:flog (* sqrtpi (f2cl-lib:d1mach 1))))))
 189          (setf xmax (- (+ txmax (/ (* -0.5 (f2cl-lib:flog txmax)) txmax)) 0.01))
 190          (setf sqeps (f2cl-lib:fsqrt (* 2.0 (f2cl-lib:d1mach 3))))))
 191       (setf first$ f2cl-lib:%false%)
 192       (if (> x xsml) (go label20))
 193       (setf derfc 2.0)
 194       (go end_label)
 195      label20
 196       (if (> x xmax) (go label40))
 197       (setf y (abs x))
 198       (if (> y 1.0) (go label30))
 199       (if (< y sqeps) (setf derfc (+ 1.0 (/ (* -2.0 x) sqrtpi))))
 200       (if (>= y sqeps)
 201           (setf derfc
 202                   (- 1.0
 203                      (* x (+ 1.0 (dcsevl (- (* 2.0 x x) 1.0) erfcs nterf))))))
 204       (go end_label)
 205      label30
 206       (setf y (* y y))
 207       (if (<= y 4.0)
 208           (setf derfc
 209                   (* (/ (exp (- y)) (abs x))
 210                      (+ 0.5
 211                         (dcsevl (/ (- (/ 8.0 y) 5.0) 3.0) erc2cs nterc2)))))
 212       (if (> y 4.0)
 213           (setf derfc
 214                   (* (/ (exp (- y)) (abs x))
 215                      (+ 0.5 (dcsevl (- (/ 8.0 y) 1.0) erfccs nterfc)))))
 216       (if (< x 0.0) (setf derfc (- 2.0 derfc)))
 217       (go end_label)
 218      label40
 219       (xermsg "SLATEC" "DERFC" "X SO BIG ERFC UNDERFLOWS" 1 1)
 220       (setf derfc 0.0)
 221       (go end_label)
 222      end_label
 223       (return (values derfc nil)))))
 224
 225 (in-package #:cl-user)
 226 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 227 (eval-when (:load-toplevel :compile-toplevel :execute)
 228   (setf (gethash 'fortran-to-lisp::derfc fortran-to-lisp::*f2cl-function-info*)
 229           (fortran-to-lisp::make-f2cl-finfo :arg-types '((double-float))
 230                                             :return-values '(nil)
 231                                             :calls '(fortran-to-lisp::xermsg
 232                                                      fortran-to-lisp::dcsevl
 233                                                      fortran-to-lisp::initds
 234                                                      fortran-to-lisp::d1mach))))
 235