src/numerical/slatec/zunik.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 ((zeror 0.0)
  21       (zeroi 0.0)
  22       (coner 1.0)
  23       (conei 0.0)
  24       (con
  25        (make-array 2
  26                    :element-type 'double-float
  27                    :initial-contents '(0.3989422804014327 1.2533141373155003)))
  28       (c
  29        (make-array 120
  30                    :element-type 'double-float
  31                    :initial-contents '(1.0 -0.20833333333333334 0.125
  32                                        0.3342013888888889 -0.4010416666666667
  33                                        0.0703125 -1.0258125964506173
  34                                        1.8464626736111112 -0.8912109375
  35                                        0.0732421875 4.669584423426247
  36                                        -11.207002616222994 8.78912353515625
  37                                        -2.3640869140625 0.112152099609375
  38                                        -28.212072558200244 84.63621767460073
  39                                        -91.81824154324002 42.53499874538846
  40                                        -7.368794359479632 0.22710800170898438
  41                                        212.57013003921713 -765.2524681411817
  42                                        1059.9904525279999 -699.5796273761325
  43                                        218.1905117442116 -26.491430486951554
  44                                        0.5725014209747314 -1919.457662318407
  45                                        8061.722181737309 -13586.550006434138
  46                                        11655.393336864534 -5305.646978613403
  47                                        1200.9029132163525 -108.09091978839466
  48                                        1.7277275025844574 20204.29133096615
  49                                        -96980.59838863752 192547.00123253153
  50                                        -203400.17728041555 122200.46498301746
  51                                        -41192.65496889755 7109.514302489364
  52                                        -493.915304773088 6.074042001273483
  53                                        -242919.18790055133 1311763.6146629772
  54                                        -2998015.9185381066 3763271.297656404
  55                                        -2813563.226586534 1268365.2733216248
  56                                        -331645.1724845636 45218.76898136273
  57                                        -2499.8304818112097 24.380529699556064
  58                                        3284469.853072038 -1.9706819118432228e7
  59                                        5.095260249266464e7 -7.410514821153265e7
  60                                        6.634451227472903e7 -3.756717666076335e7
  61                                        1.3288767166421818e7 -2785618.1280864547
  62                                        308186.4046126624 -13886.08975371704
  63                                        110.01714026924674 -4.932925366450996e7
  64                                        3.2557307418576574e8
  65                                        -9.394623596815784e8 1.55359689957058e9
  66                                        -1.6210805521083372e9
  67                                        1.1068428168230145e9
  68                                        -4.958897842750303e8 1.420629077975331e8
  69                                        -2.447406272573873e7 2243768.1779224495
  70                                        -84005.43360302408 551.3358961220206
  71                                        8.147890961183121e8 -5.866481492051847e9
  72                                        1.8688207509295826e10
  73                                        -3.4632043388158775e10
  74                                        4.1280185579753975e10
  75                                        -3.3026599749800724e10
  76                                        1.79542137311556e10 -6.563293792619285e9
  77                                        1.5592798648792574e9
  78                                        -2.2510566188941526e8
  79                                        1.7395107553978164e7 -549842.3275722887
  80                                        3038.090510922384 -1.4679261247695616e10
  81                                        1.144982377320258e11
  82                                        -3.990961752244665e11
  83                                        8.192186695485773e11
  84                                        -1.0983751560812233e12
  85                                        1.0081581068653821e12
  86                                        -6.453648692453765e11
  87                                        2.879006499061506e11
  88                                        -8.786707217802327e10
  89                                        1.763473060683497e10
  90                                        -2.167164983223795e9
  91                                        1.4315787671888897e8 -3871833.442572613
  92                                        18257.755474293175 2.86464035717679e11
  93                                        -2.406297900028504e12
  94                                        9.109341185239898e12
  95                                        -2.0516899410934438e13
  96                                        3.056512551993532e13
  97                                        -3.166708858478516e13
  98                                        2.334836404458184e13
  99                                        -1.2320491305598287e13
 100                                        4.612725780849132e12
 101                                        -1.1965528801961816e12
 102                                        2.0591450323241e11
 103                                        -2.1822927757529224e10
 104                                        1.2470092935127103e9
 105                                        -2.9188388122220814e7
 106                                        118838.42625678325))))
 107   (declare (type (double-float) zeror zeroi coner conei)
 108            (type (simple-array double-float (2)) con)
 109            (type (simple-array double-float (120)) c))
 110   (defun zunik
 111          (zrr zri fnu ikflg ipmtr tol init phir phii zeta1r zeta1i zeta2r
 112           zeta2i sumr sumi cwrkr cwrki)
 113     (declare (type (simple-array double-float (*)) cwrki cwrkr)
 114              (type (f2cl-lib:integer4) init ipmtr ikflg)
 115              (type (double-float) sumi sumr zeta2i zeta2r zeta1i zeta1r phii
 116                                   phir tol fnu zri zrr))
 117     (prog ((i 0) (idum 0) (j 0) (k 0) (l 0) (ac 0.0) (crfni 0.0) (crfnr 0.0)
 118            (rfn 0.0) (si 0.0) (sr 0.0) (sri 0.0) (srr 0.0) (sti 0.0) (str 0.0)
 119            (test 0.0) (ti 0.0) (tr 0.0) (t2i 0.0) (t2r 0.0) (zni 0.0)
 120            (znr 0.0))
 121       (declare (type (double-float) znr zni t2r t2i tr ti test str sti srr sri
 122                                     sr si rfn crfnr crfni ac)
 123                (type (f2cl-lib:integer4) l k j idum i))
 124       (if (/= init 0) (go label40))
 125       (setf rfn (/ 1.0 fnu))
 126       (setf test (* (f2cl-lib:d1mach 1) 1000.0))
 127       (setf ac (* fnu test))
 128       (if (or (> (abs zrr) ac) (> (abs zri) ac)) (go label15))
 129       (setf zeta1r (+ (* 2.0 (abs (f2cl-lib:flog test))) fnu))
 130       (setf zeta1i 0.0)
 131       (setf zeta2r fnu)
 132       (setf zeta2i 0.0)
 133       (setf phir 1.0)
 134       (setf phii 0.0)
 135       (go end_label)
 136      label15
 137       (setf tr (* zrr rfn))
 138       (setf ti (* zri rfn))
 139       (setf sr (+ coner (- (* tr tr) (* ti ti))))
 140       (setf si (+ conei (+ (* tr ti) (* ti tr))))
 141       (multiple-value-bind (var-0 var-1 var-2 var-3)
 142           (zsqrt$ sr si srr sri)
 143         (declare (ignore var-0 var-1))
 144         (setf srr var-2)
 145         (setf sri var-3))
 146       (setf str (+ coner srr))
 147       (setf sti (+ conei sri))
 148       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 149           (zdiv str sti tr ti znr zni)
 150         (declare (ignore var-0 var-1 var-2 var-3))
 151         (setf znr var-4)
 152         (setf zni var-5))
 153       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4)
 154           (zlog znr zni str sti idum)
 155         (declare (ignore var-0 var-1))
 156         (setf str var-2)
 157         (setf sti var-3)
 158         (setf idum var-4))
 159       (setf zeta1r (* fnu str))
 160       (setf zeta1i (* fnu sti))
 161       (setf zeta2r (* fnu srr))
 162       (setf zeta2i (* fnu sri))
 163       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 164           (zdiv coner conei srr sri tr ti)
 165         (declare (ignore var-0 var-1 var-2 var-3))
 166         (setf tr var-4)
 167         (setf ti var-5))
 168       (setf srr (* tr rfn))
 169       (setf sri (* ti rfn))
 170       (multiple-value-bind (var-0 var-1 var-2 var-3)
 171           (zsqrt$ srr sri (f2cl-lib:fref cwrkr (16) ((1 16)))
 172            (f2cl-lib:fref cwrki (16) ((1 16))))
 173         (declare (ignore var-0 var-1))
 174         (setf (f2cl-lib:fref cwrkr (16) ((1 16))) var-2)
 175         (setf (f2cl-lib:fref cwrki (16) ((1 16))) var-3))
 176       (setf phir
 177               (* (f2cl-lib:fref cwrkr (16) ((1 16)))
 178                  (f2cl-lib:fref con (ikflg) ((1 2)))))
 179       (setf phii
 180               (* (f2cl-lib:fref cwrki (16) ((1 16)))
 181                  (f2cl-lib:fref con (ikflg) ((1 2)))))
 182       (if (/= ipmtr 0) (go end_label))
 183       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 184           (zdiv coner conei sr si t2r t2i)
 185         (declare (ignore var-0 var-1 var-2 var-3))
 186         (setf t2r var-4)
 187         (setf t2i var-5))
 188       (setf (f2cl-lib:fref cwrkr (1) ((1 16))) coner)
 189       (setf (f2cl-lib:fref cwrki (1) ((1 16))) conei)
 190       (setf crfnr coner)
 191       (setf crfni conei)
 192       (setf ac 1.0)
 193       (setf l 1)
 194       (f2cl-lib:fdo (k 2 (f2cl-lib:int-add k 1))
 195                     ((> k 15) nil)
 196         (tagbody
 197           (setf sr zeror)
 198           (setf si zeroi)
 199           (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 200                         ((> j k) nil)
 201             (tagbody
 202               (setf l (f2cl-lib:int-add l 1))
 203               (setf str
 204                       (+ (- (* sr t2r) (* si t2i))
 205                          (f2cl-lib:fref c (l) ((1 120)))))
 206               (setf si (+ (* sr t2i) (* si t2r)))
 207               (setf sr str)
 208              label10))
 209           (setf str (- (* crfnr srr) (* crfni sri)))
 210           (setf crfni (+ (* crfnr sri) (* crfni srr)))
 211           (setf crfnr str)
 212           (setf (f2cl-lib:fref cwrkr (k) ((1 16)))
 213                   (- (* crfnr sr) (* crfni si)))
 214           (setf (f2cl-lib:fref cwrki (k) ((1 16)))
 215                   (+ (* crfnr si) (* crfni sr)))
 216           (setf ac (* ac rfn))
 217           (setf test
 218                   (+ (abs (f2cl-lib:fref cwrkr (k) ((1 16))))
 219                      (abs (f2cl-lib:fref cwrki (k) ((1 16))))))
 220           (if (and (< ac tol) (< test tol)) (go label30))
 221          label20))
 222       (setf k 15)
 223      label30
 224       (setf init k)
 225      label40
 226       (if (= ikflg 2) (go label60))
 227       (setf sr zeror)
 228       (setf si zeroi)
 229       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 230                     ((> i init) nil)
 231         (tagbody
 232           (setf sr (+ sr (f2cl-lib:fref cwrkr (i) ((1 16)))))
 233           (setf si (+ si (f2cl-lib:fref cwrki (i) ((1 16)))))
 234          label50))
 235       (setf sumr sr)
 236       (setf sumi si)
 237       (setf phir
 238               (* (f2cl-lib:fref cwrkr (16) ((1 16)))
 239                  (f2cl-lib:fref con (1) ((1 2)))))
 240       (setf phii
 241               (* (f2cl-lib:fref cwrki (16) ((1 16)))
 242                  (f2cl-lib:fref con (1) ((1 2)))))
 243       (go end_label)
 244      label60
 245       (setf sr zeror)
 246       (setf si zeroi)
 247       (setf tr coner)
 248       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 249                     ((> i init) nil)
 250         (tagbody
 251           (setf sr (+ sr (* tr (f2cl-lib:fref cwrkr (i) ((1 16))))))
 252           (setf si (+ si (* tr (f2cl-lib:fref cwrki (i) ((1 16))))))
 253           (setf tr (- tr))
 254          label70))
 255       (setf sumr sr)
 256       (setf sumi si)
 257       (setf phir
 258               (* (f2cl-lib:fref cwrkr (16) ((1 16)))
 259                  (f2cl-lib:fref con (2) ((1 2)))))
 260       (setf phii
 261               (* (f2cl-lib:fref cwrki (16) ((1 16)))
 262                  (f2cl-lib:fref con (2) ((1 2)))))
 263       (go end_label)
 264      end_label
 265       (return
 266        (values nil
 267                nil
 268                nil
 269                nil
 270                nil
 271                nil
 272                init
 273                phir
 274                phii
 275                zeta1r
 276                zeta1i
 277                zeta2r
 278                zeta2i
 279                sumr
 280                sumi
 281                nil
 282                nil)))))
 283
 284 (in-package #:cl-user)
 285 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 286 (eval-when (:load-toplevel :compile-toplevel :execute)
 287   (setf (gethash 'fortran-to-lisp::zunik fortran-to-lisp::*f2cl-function-info*)
 288           (fortran-to-lisp::make-f2cl-finfo
 289            :arg-types '((double-float) (double-float) (double-float)
 290                         (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 291                         (double-float) (fortran-to-lisp::integer4)
 292                         (double-float) (double-float) (double-float)
 293                         (double-float) (double-float) (double-float)
 294                         (double-float) (double-float)
 295                         (simple-array double-float (*))
 296                         (simple-array double-float (*)))
 297            :return-values '(nil nil nil nil nil nil fortran-to-lisp::init
 298                             fortran-to-lisp::phir fortran-to-lisp::phii
 299                             fortran-to-lisp::zeta1r fortran-to-lisp::zeta1i
 300                             fortran-to-lisp::zeta2r fortran-to-lisp::zeta2i
 301                             fortran-to-lisp::sumr fortran-to-lisp::sumi nil
 302                             nil)
 303            :calls '(fortran-to-lisp::zlog fortran-to-lisp::zdiv
 304                     fortran-to-lisp::zsqrt$ fortran-to-lisp::d1mach))))
 305