src/numerical/slatec/zuni1.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) (zeroi 0.0) (coner 1.0))
  21   (declare (type (double-float) zeror zeroi coner))
  22   (defun zuni1 (zr zi fnu kode n yr yi nz nlast fnul tol elim alim)
  23     (declare (type (simple-array double-float (*)) yi yr)
  24              (type (f2cl-lib:integer4) nlast nz n kode)
  25              (type (double-float) alim elim tol fnul fnu zi zr))
  26     (prog ((bry (make-array 3 :element-type 'double-float))
  27            (cwrkr (make-array 16 :element-type 'double-float))
  28            (cwrki (make-array 16 :element-type 'double-float))
  29            (cssr (make-array 3 :element-type 'double-float))
  30            (csrr (make-array 3 :element-type 'double-float))
  31            (cyr (make-array 2 :element-type 'double-float))
  32            (cyi (make-array 2 :element-type 'double-float)) (i 0) (iflag 0)
  33            (init 0) (k 0) (m 0) (nd 0) (nn 0) (nuf 0) (nw 0) (aphi 0.0)
  34            (ascle 0.0) (crsc 0.0) (cscl 0.0) (c1r 0.0) (c2i 0.0) (c2m 0.0)
  35            (c2r 0.0) (fn 0.0) (phii 0.0) (phir 0.0) (rast 0.0) (rs1 0.0)
  36            (rzi 0.0) (rzr 0.0) (sti 0.0) (str 0.0) (sumi 0.0) (sumr 0.0)
  37            (s1i 0.0) (s1r 0.0) (s2i 0.0) (s2r 0.0) (zeta1i 0.0) (zeta1r 0.0)
  38            (zeta2i 0.0) (zeta2r 0.0))
  39       (declare (type (simple-array double-float (2)) cyi cyr)
  40                (type (simple-array double-float (16)) cwrkr cwrki)
  41                (type (simple-array double-float (3)) cssr csrr bry)
  42                (type (double-float) zeta2r zeta2i zeta1r zeta1i s2r s2i s1r s1i
  43                                     sumr sumi str sti rzr rzi rs1 rast phir
  44                                     phii fn c2r c2m c2i c1r cscl crsc ascle
  45                                     aphi)
  46                (type (f2cl-lib:integer4) nw nuf nn nd m k init iflag i))
  47       (setf nz 0)
  48       (setf nd n)
  49       (setf nlast 0)
  50       (setf cscl (/ 1.0 tol))
  51       (setf crsc tol)
  52       (setf (f2cl-lib:fref cssr (1) ((1 3))) cscl)
  53       (setf (f2cl-lib:fref cssr (2) ((1 3))) coner)
  54       (setf (f2cl-lib:fref cssr (3) ((1 3))) crsc)
  55       (setf (f2cl-lib:fref csrr (1) ((1 3))) crsc)
  56       (setf (f2cl-lib:fref csrr (2) ((1 3))) coner)
  57       (setf (f2cl-lib:fref csrr (3) ((1 3))) cscl)
  58       (setf (f2cl-lib:fref bry (1) ((1 3)))
  59               (/ (* 1000.0 (f2cl-lib:d1mach 1)) tol))
  60       (setf fn (max fnu 1.0))
  61       (setf init 0)
  62       (multiple-value-bind
  63             (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10
  64              var-11 var-12 var-13 var-14 var-15 var-16)
  65           (zunik zr zi fn 1 1 tol init phir phii zeta1r zeta1i zeta2r zeta2i
  66            sumr sumi cwrkr cwrki)
  67         (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-15 var-16))
  68         (setf init var-6)
  69         (setf phir var-7)
  70         (setf phii var-8)
  71         (setf zeta1r var-9)
  72         (setf zeta1i var-10)
  73         (setf zeta2r var-11)
  74         (setf zeta2i var-12)
  75         (setf sumr var-13)
  76         (setf sumi var-14))
  77       (if (= kode 1) (go label10))
  78       (setf str (+ zr zeta2r))
  79       (setf sti (+ zi zeta2i))
  80       (setf rast (coerce (realpart (/ fn (zabs str sti))) 'double-float))
  81       (setf str (* str rast rast))
  82       (setf sti (* (- sti) rast rast))
  83       (setf s1r (- str zeta1r))
  84       (setf s1i (- sti zeta1i))
  85       (go label20)
  86      label10
  87       (setf s1r (- zeta2r zeta1r))
  88       (setf s1i (- zeta2i zeta1i))
  89      label20
  90       (setf rs1 s1r)
  91       (if (> (abs rs1) elim) (go label130))
  92      label30
  93       (setf nn (min (the f2cl-lib:integer4 2) (the f2cl-lib:integer4 nd)))
  94       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  95                     ((> i nn) nil)
  96         (tagbody
  97           (setf fn (+ fnu (f2cl-lib:int-sub nd i)))
  98           (setf init 0)
  99           (multiple-value-bind
 100                 (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9
 101                  var-10 var-11 var-12 var-13 var-14 var-15 var-16)
 102               (zunik zr zi fn 1 0 tol init phir phii zeta1r zeta1i zeta2r
 103                zeta2i sumr sumi cwrkr cwrki)
 104             (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-15
 105                              var-16))
 106             (setf init var-6)
 107             (setf phir var-7)
 108             (setf phii var-8)
 109             (setf zeta1r var-9)
 110             (setf zeta1i var-10)
 111             (setf zeta2r var-11)
 112             (setf zeta2i var-12)
 113             (setf sumr var-13)
 114             (setf sumi var-14))
 115           (if (= kode 1) (go label40))
 116           (setf str (+ zr zeta2r))
 117           (setf sti (+ zi zeta2i))
 118           (setf rast (coerce (realpart (/ fn (zabs str sti))) 'double-float))
 119           (setf str (* str rast rast))
 120           (setf sti (* (- sti) rast rast))
 121           (setf s1r (- str zeta1r))
 122           (setf s1i (+ (- sti zeta1i) zi))
 123           (go label50)
 124          label40
 125           (setf s1r (- zeta2r zeta1r))
 126           (setf s1i (- zeta2i zeta1i))
 127          label50
 128           (setf rs1 s1r)
 129           (if (> (abs rs1) elim) (go label110))
 130           (if (= i 1) (setf iflag 2))
 131           (if (< (abs rs1) alim) (go label60))
 132           (setf aphi (coerce (realpart (zabs phir phii)) 'double-float))
 133           (setf rs1 (+ rs1 (f2cl-lib:flog aphi)))
 134           (if (> (abs rs1) elim) (go label110))
 135           (if (= i 1) (setf iflag 1))
 136           (if (< rs1 0.0) (go label60))
 137           (if (= i 1) (setf iflag 3))
 138          label60
 139           (setf s2r (- (* phir sumr) (* phii sumi)))
 140           (setf s2i (+ (* phir sumi) (* phii sumr)))
 141           (setf str (* (exp s1r) (f2cl-lib:fref cssr (iflag) ((1 3)))))
 142           (setf s1r (* str (cos s1i)))
 143           (setf s1i (* str (sin s1i)))
 144           (setf str (- (* s2r s1r) (* s2i s1i)))
 145           (setf s2i (+ (* s2r s1i) (* s2i s1r)))
 146           (setf s2r str)
 147           (if (/= iflag 1) (go label70))
 148           (multiple-value-bind (var-0 var-1 var-2 var-3 var-4)
 149               (zuchk s2r s2i nw (f2cl-lib:fref bry (1) ((1 3))) tol)
 150             (declare (ignore var-0 var-1 var-3 var-4))
 151             (setf nw var-2))
 152           (if (/= nw 0) (go label110))
 153          label70
 154           (setf (f2cl-lib:fref cyr (i) ((1 2))) s2r)
 155           (setf (f2cl-lib:fref cyi (i) ((1 2))) s2i)
 156           (setf m (f2cl-lib:int-add (f2cl-lib:int-sub nd i) 1))
 157           (setf (f2cl-lib:fref yr (m) ((1 n)))
 158                   (* s2r (f2cl-lib:fref csrr (iflag) ((1 3)))))
 159           (setf (f2cl-lib:fref yi (m) ((1 n)))
 160                   (* s2i (f2cl-lib:fref csrr (iflag) ((1 3)))))
 161          label80))
 162       (if (<= nd 2) (go label100))
 163       (setf rast (coerce (realpart (/ 1.0 (zabs zr zi))) 'double-float))
 164       (setf str (* zr rast))
 165       (setf sti (* (- zi) rast))
 166       (setf rzr (* (+ str str) rast))
 167       (setf rzi (* (+ sti sti) rast))
 168       (setf (f2cl-lib:fref bry (2) ((1 3)))
 169               (/ 1.0 (f2cl-lib:fref bry (1) ((1 3)))))
 170       (setf (f2cl-lib:fref bry (3) ((1 3))) (f2cl-lib:d1mach 2))
 171       (setf s1r (f2cl-lib:fref cyr (1) ((1 2))))
 172       (setf s1i (f2cl-lib:fref cyi (1) ((1 2))))
 173       (setf s2r (f2cl-lib:fref cyr (2) ((1 2))))
 174       (setf s2i (f2cl-lib:fref cyi (2) ((1 2))))
 175       (setf c1r (f2cl-lib:fref csrr (iflag) ((1 3))))
 176       (setf ascle (f2cl-lib:fref bry (iflag) ((1 3))))
 177       (setf k (f2cl-lib:int-sub nd 2))
 178       (setf fn (coerce (the f2cl-lib:integer4 k) 'double-float))
 179       (f2cl-lib:fdo (i 3 (f2cl-lib:int-add i 1))
 180                     ((> i nd) nil)
 181         (tagbody
 182           (setf c2r s2r)
 183           (setf c2i s2i)
 184           (setf s2r (+ s1r (* (+ fnu fn) (- (* rzr c2r) (* rzi c2i)))))
 185           (setf s2i (+ s1i (* (+ fnu fn) (+ (* rzr c2i) (* rzi c2r)))))
 186           (setf s1r c2r)
 187           (setf s1i c2i)
 188           (setf c2r (* s2r c1r))
 189           (setf c2i (* s2i c1r))
 190           (setf (f2cl-lib:fref yr (k) ((1 n))) c2r)
 191           (setf (f2cl-lib:fref yi (k) ((1 n))) c2i)
 192           (setf k (f2cl-lib:int-sub k 1))
 193           (setf fn (- fn 1.0))
 194           (if (>= iflag 3) (go label90))
 195           (setf str (abs c2r))
 196           (setf sti (abs c2i))
 197           (setf c2m (max str sti))
 198           (if (<= c2m ascle) (go label90))
 199           (setf iflag (f2cl-lib:int-add iflag 1))
 200           (setf ascle (f2cl-lib:fref bry (iflag) ((1 3))))
 201           (setf s1r (* s1r c1r))
 202           (setf s1i (* s1i c1r))
 203           (setf s2r c2r)
 204           (setf s2i c2i)
 205           (setf s1r (* s1r (f2cl-lib:fref cssr (iflag) ((1 3)))))
 206           (setf s1i (* s1i (f2cl-lib:fref cssr (iflag) ((1 3)))))
 207           (setf s2r (* s2r (f2cl-lib:fref cssr (iflag) ((1 3)))))
 208           (setf s2i (* s2i (f2cl-lib:fref cssr (iflag) ((1 3)))))
 209           (setf c1r (f2cl-lib:fref csrr (iflag) ((1 3))))
 210          label90))
 211      label100
 212       (go end_label)
 213      label110
 214       (if (> rs1 0.0) (go label120))
 215       (setf (f2cl-lib:fref yr (nd) ((1 n))) zeror)
 216       (setf (f2cl-lib:fref yi (nd) ((1 n))) zeroi)
 217       (setf nz (f2cl-lib:int-add nz 1))
 218       (setf nd (f2cl-lib:int-sub nd 1))
 219       (if (= nd 0) (go label100))
 220       (multiple-value-bind
 221             (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10
 222              var-11)
 223           (zuoik zr zi fnu kode 1 nd yr yi nuf tol elim alim)
 224         (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-9
 225                          var-10 var-11))
 226         (setf nuf var-8))
 227       (if (< nuf 0) (go label120))
 228       (setf nd (f2cl-lib:int-sub nd nuf))
 229       (setf nz (f2cl-lib:int-add nz nuf))
 230       (if (= nd 0) (go label100))
 231       (setf fn (+ fnu (f2cl-lib:int-sub nd 1)))
 232       (if (>= fn fnul) (go label30))
 233       (setf nlast nd)
 234       (go end_label)
 235      label120
 236       (setf nz -1)
 237       (go end_label)
 238      label130
 239       (if (> rs1 0.0) (go label120))
 240       (setf nz n)
 241       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 242                     ((> i n) nil)
 243         (tagbody
 244           (setf (f2cl-lib:fref yr (i) ((1 n))) zeror)
 245           (setf (f2cl-lib:fref yi (i) ((1 n))) zeroi)
 246          label140))
 247       (go end_label)
 248      end_label
 249       (return (values nil nil nil nil nil nil nil nz nlast nil nil nil nil)))))
 250
 251 (in-package #:cl-user)
 252 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 253 (eval-when (:load-toplevel :compile-toplevel :execute)
 254   (setf (gethash 'fortran-to-lisp::zuni1 fortran-to-lisp::*f2cl-function-info*)
 255           (fortran-to-lisp::make-f2cl-finfo
 256            :arg-types '((double-float) (double-float) (double-float)
 257                         (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 258                         (simple-array double-float (*))
 259                         (simple-array double-float (*))
 260                         (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 261                         (double-float) (double-float) (double-float)
 262                         (double-float))
 263            :return-values '(nil nil nil nil nil nil nil fortran-to-lisp::nz
 264                             fortran-to-lisp::nlast nil nil nil nil)
 265            :calls '(fortran-to-lisp::zuoik fortran-to-lisp::zuchk
 266                     fortran-to-lisp::zabs fortran-to-lisp::zunik
 267                     fortran-to-lisp::d1mach))))
 268