src/numerical/slatec/zuni2.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       (cipr
  24        (make-array 4
  25                    :element-type 'double-float
  26                    :initial-contents '(1.0 0.0 -1.0 0.0)))
  27       (cipi
  28        (make-array 4
  29                    :element-type 'double-float
  30                    :initial-contents '(0.0 1.0 0.0 -1.0)))
  31       (hpi 1.5707963267948966)
  32       (aic 1.2655121234846454))
  33   (declare (type (double-float) zeror zeroi coner hpi aic)
  34            (type (simple-array double-float (4)) cipr cipi))
  35   (defun zuni2 (zr zi fnu kode n yr yi nz nlast fnul tol elim alim)
  36     (declare (type (simple-array double-float (*)) yi yr)
  37              (type (f2cl-lib:integer4) nlast nz n kode)
  38              (type (double-float) alim elim tol fnul fnu zi zr))
  39     (prog ((bry (make-array 3 :element-type 'double-float))
  40            (cssr (make-array 3 :element-type 'double-float))
  41            (csrr (make-array 3 :element-type 'double-float))
  42            (cyr (make-array 2 :element-type 'double-float))
  43            (cyi (make-array 2 :element-type 'double-float)) (i 0) (iflag 0)
  44            (in 0) (inu 0) (j 0) (k 0) (nai 0) (nd 0) (ndai 0) (nn 0) (nuf 0)
  45            (nw 0) (idum 0) (aarg 0.0) (aii 0.0) (air 0.0) (ang 0.0) (aphi 0.0)
  46            (argi 0.0) (argr 0.0) (ascle 0.0) (asumi 0.0) (asumr 0.0)
  47            (bsumi 0.0) (bsumr 0.0) (cidi 0.0) (crsc 0.0) (cscl 0.0) (c1r 0.0)
  48            (c2i 0.0) (c2m 0.0) (c2r 0.0) (daii 0.0) (dair 0.0) (fn 0.0)
  49            (phii 0.0) (phir 0.0) (rast 0.0) (raz 0.0) (rs1 0.0) (rzi 0.0)
  50            (rzr 0.0) (sti 0.0) (str 0.0) (s1i 0.0) (s1r 0.0) (s2i 0.0)
  51            (s2r 0.0) (zbi 0.0) (zbr 0.0) (zeta1i 0.0) (zeta1r 0.0) (zeta2i 0.0)
  52            (zeta2r 0.0) (zni 0.0) (znr 0.0) (car$ 0.0) (sar 0.0))
  53       (declare (type (simple-array double-float (2)) cyi cyr)
  54                (type (simple-array double-float (3)) cssr csrr bry)
  55                (type (double-float) sar car$ znr zni zeta2r zeta2i zeta1r
  56                                     zeta1i zbr zbi s2r s2i s1r s1i str sti rzr
  57                                     rzi rs1 raz rast phir phii fn dair daii c2r
  58                                     c2m c2i c1r cscl crsc cidi bsumr bsumi
  59                                     asumr asumi ascle argr argi aphi ang air
  60                                     aii aarg)
  61                (type (f2cl-lib:integer4) idum nw nuf nn ndai nd nai k j inu in
  62                                          iflag i))
  63       (setf nz 0)
  64       (setf nd n)
  65       (setf nlast 0)
  66       (setf cscl (/ 1.0 tol))
  67       (setf crsc tol)
  68       (setf (f2cl-lib:fref cssr (1) ((1 3))) cscl)
  69       (setf (f2cl-lib:fref cssr (2) ((1 3))) coner)
  70       (setf (f2cl-lib:fref cssr (3) ((1 3))) crsc)
  71       (setf (f2cl-lib:fref csrr (1) ((1 3))) crsc)
  72       (setf (f2cl-lib:fref csrr (2) ((1 3))) coner)
  73       (setf (f2cl-lib:fref csrr (3) ((1 3))) cscl)
  74       (setf (f2cl-lib:fref bry (1) ((1 3)))
  75               (/ (* 1000.0 (f2cl-lib:d1mach 1)) tol))
  76       (setf znr zi)
  77       (setf zni (- zr))
  78       (setf zbr zr)
  79       (setf zbi zi)
  80       (setf cidi (- coner))
  81       (setf inu (f2cl-lib:int fnu))
  82       (setf ang (* hpi (- fnu inu)))
  83       (setf c2r (cos ang))
  84       (setf c2i (sin ang))
  85       (setf car$ c2r)
  86       (setf sar c2i)
  87       (setf in (f2cl-lib:int-sub (f2cl-lib:int-add inu n) 1))
  88       (setf in (f2cl-lib:int-add (mod in 4) 1))
  89       (setf str
  90               (- (* c2r (f2cl-lib:fref cipr (in) ((1 4))))
  91                  (* c2i (f2cl-lib:fref cipi (in) ((1 4))))))
  92       (setf c2i
  93               (+ (* c2r (f2cl-lib:fref cipi (in) ((1 4))))
  94                  (* c2i (f2cl-lib:fref cipr (in) ((1 4))))))
  95       (setf c2r str)
  96       (if (> zi 0.0) (go label10))
  97       (setf znr (- znr))
  98       (setf zbi (- zbi))
  99       (setf cidi (- cidi))
 100       (setf c2i (- c2i))
 101      label10
 102       (setf fn (max fnu 1.0))
 103       (multiple-value-bind
 104             (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10
 105              var-11 var-12 var-13 var-14 var-15 var-16)
 106           (zunhj znr zni fn 1 tol phir phii argr argi zeta1r zeta1i zeta2r
 107            zeta2i asumr asumi bsumr bsumi)
 108         (declare (ignore var-0 var-1 var-2 var-3 var-4))
 109         (setf phir var-5)
 110         (setf phii var-6)
 111         (setf argr var-7)
 112         (setf argi var-8)
 113         (setf zeta1r var-9)
 114         (setf zeta1i var-10)
 115         (setf zeta2r var-11)
 116         (setf zeta2i var-12)
 117         (setf asumr var-13)
 118         (setf asumi var-14)
 119         (setf bsumr var-15)
 120         (setf bsumi var-16))
 121       (if (= kode 1) (go label20))
 122       (setf str (+ zbr zeta2r))
 123       (setf sti (+ zbi zeta2i))
 124       (setf rast (coerce (realpart (/ fn (zabs str sti))) 'double-float))
 125       (setf str (* str rast rast))
 126       (setf sti (* (- sti) rast rast))
 127       (setf s1r (- str zeta1r))
 128       (setf s1i (- sti zeta1i))
 129       (go label30)
 130      label20
 131       (setf s1r (- zeta2r zeta1r))
 132       (setf s1i (- zeta2i zeta1i))
 133      label30
 134       (setf rs1 s1r)
 135       (if (> (abs rs1) elim) (go label150))
 136      label40
 137       (setf nn (min (the f2cl-lib:integer4 2) (the f2cl-lib:integer4 nd)))
 138       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 139                     ((> i nn) nil)
 140         (tagbody
 141           (setf fn (+ fnu (f2cl-lib:int-sub nd i)))
 142           (multiple-value-bind
 143                 (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9
 144                  var-10 var-11 var-12 var-13 var-14 var-15 var-16)
 145               (zunhj znr zni fn 0 tol phir phii argr argi zeta1r zeta1i zeta2r
 146                zeta2i asumr asumi bsumr bsumi)
 147             (declare (ignore var-0 var-1 var-2 var-3 var-4))
 148             (setf phir var-5)
 149             (setf phii var-6)
 150             (setf argr var-7)
 151             (setf argi var-8)
 152             (setf zeta1r var-9)
 153             (setf zeta1i var-10)
 154             (setf zeta2r var-11)
 155             (setf zeta2i var-12)
 156             (setf asumr var-13)
 157             (setf asumi var-14)
 158             (setf bsumr var-15)
 159             (setf bsumi var-16))
 160           (if (= kode 1) (go label50))
 161           (setf str (+ zbr zeta2r))
 162           (setf sti (+ zbi zeta2i))
 163           (setf rast (coerce (realpart (/ fn (zabs str sti))) 'double-float))
 164           (setf str (* str rast rast))
 165           (setf sti (* (- sti) rast rast))
 166           (setf s1r (- str zeta1r))
 167           (setf s1i (+ (- sti zeta1i) (abs zi)))
 168           (go label60)
 169          label50
 170           (setf s1r (- zeta2r zeta1r))
 171           (setf s1i (- zeta2i zeta1i))
 172          label60
 173           (setf rs1 s1r)
 174           (if (> (abs rs1) elim) (go label120))
 175           (if (= i 1) (setf iflag 2))
 176           (if (< (abs rs1) alim) (go label70))
 177           (setf aphi (coerce (realpart (zabs phir phii)) 'double-float))
 178           (setf aarg (coerce (realpart (zabs argr argi)) 'double-float))
 179           (setf rs1
 180                   (- (+ rs1 (f2cl-lib:flog aphi))
 181                      (* 0.25 (f2cl-lib:flog aarg))
 182                      aic))
 183           (if (> (abs rs1) elim) (go label120))
 184           (if (= i 1) (setf iflag 1))
 185           (if (< rs1 0.0) (go label70))
 186           (if (= i 1) (setf iflag 3))
 187          label70
 188           (multiple-value-bind
 189                 (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7)
 190               (zairy argr argi 0 2 air aii nai idum)
 191             (declare (ignore var-0 var-1 var-2 var-3))
 192             (setf air var-4)
 193             (setf aii var-5)
 194             (setf nai var-6)
 195             (setf idum var-7))
 196           (multiple-value-bind
 197                 (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7)
 198               (zairy argr argi 1 2 dair daii ndai idum)
 199             (declare (ignore var-0 var-1 var-2 var-3))
 200             (setf dair var-4)
 201             (setf daii var-5)
 202             (setf ndai var-6)
 203             (setf idum var-7))
 204           (setf str (- (* dair bsumr) (* daii bsumi)))
 205           (setf sti (+ (* dair bsumi) (* daii bsumr)))
 206           (setf str (+ str (- (* air asumr) (* aii asumi))))
 207           (setf sti (+ sti (+ (* air asumi) (* aii asumr))))
 208           (setf s2r (- (* phir str) (* phii sti)))
 209           (setf s2i (+ (* phir sti) (* phii str)))
 210           (setf str (* (exp s1r) (f2cl-lib:fref cssr (iflag) ((1 3)))))
 211           (setf s1r (* str (cos s1i)))
 212           (setf s1i (* str (sin s1i)))
 213           (setf str (- (* s2r s1r) (* s2i s1i)))
 214           (setf s2i (+ (* s2r s1i) (* s2i s1r)))
 215           (setf s2r str)
 216           (if (/= iflag 1) (go label80))
 217           (multiple-value-bind (var-0 var-1 var-2 var-3 var-4)
 218               (zuchk s2r s2i nw (f2cl-lib:fref bry (1) ((1 3))) tol)
 219             (declare (ignore var-0 var-1 var-3 var-4))
 220             (setf nw var-2))
 221           (if (/= nw 0) (go label120))
 222          label80
 223           (if (<= zi 0.0) (setf s2i (- s2i)))
 224           (setf str (- (* s2r c2r) (* s2i c2i)))
 225           (setf s2i (+ (* s2r c2i) (* s2i c2r)))
 226           (setf s2r str)
 227           (setf (f2cl-lib:fref cyr (i) ((1 2))) s2r)
 228           (setf (f2cl-lib:fref cyi (i) ((1 2))) s2i)
 229           (setf j (f2cl-lib:int-add (f2cl-lib:int-sub nd i) 1))
 230           (setf (f2cl-lib:fref yr (j) ((1 n)))
 231                   (* s2r (f2cl-lib:fref csrr (iflag) ((1 3)))))
 232           (setf (f2cl-lib:fref yi (j) ((1 n)))
 233                   (* s2i (f2cl-lib:fref csrr (iflag) ((1 3)))))
 234           (setf str (* (- c2i) cidi))
 235           (setf c2i (* c2r cidi))
 236           (setf c2r str)
 237          label90))
 238       (if (<= nd 2) (go label110))
 239       (setf raz (coerce (realpart (/ 1.0 (zabs zr zi))) 'double-float))
 240       (setf str (* zr raz))
 241       (setf sti (* (- zi) raz))
 242       (setf rzr (* (+ str str) raz))
 243       (setf rzi (* (+ sti sti) raz))
 244       (setf (f2cl-lib:fref bry (2) ((1 3)))
 245               (/ 1.0 (f2cl-lib:fref bry (1) ((1 3)))))
 246       (setf (f2cl-lib:fref bry (3) ((1 3))) (f2cl-lib:d1mach 2))
 247       (setf s1r (f2cl-lib:fref cyr (1) ((1 2))))
 248       (setf s1i (f2cl-lib:fref cyi (1) ((1 2))))
 249       (setf s2r (f2cl-lib:fref cyr (2) ((1 2))))
 250       (setf s2i (f2cl-lib:fref cyi (2) ((1 2))))
 251       (setf c1r (f2cl-lib:fref csrr (iflag) ((1 3))))
 252       (setf ascle (f2cl-lib:fref bry (iflag) ((1 3))))
 253       (setf k (f2cl-lib:int-sub nd 2))
 254       (setf fn (coerce (the f2cl-lib:integer4 k) 'double-float))
 255       (f2cl-lib:fdo (i 3 (f2cl-lib:int-add i 1))
 256                     ((> i nd) nil)
 257         (tagbody
 258           (setf c2r s2r)
 259           (setf c2i s2i)
 260           (setf s2r (+ s1r (* (+ fnu fn) (- (* rzr c2r) (* rzi c2i)))))
 261           (setf s2i (+ s1i (* (+ fnu fn) (+ (* rzr c2i) (* rzi c2r)))))
 262           (setf s1r c2r)
 263           (setf s1i c2i)
 264           (setf c2r (* s2r c1r))
 265           (setf c2i (* s2i c1r))
 266           (setf (f2cl-lib:fref yr (k) ((1 n))) c2r)
 267           (setf (f2cl-lib:fref yi (k) ((1 n))) c2i)
 268           (setf k (f2cl-lib:int-sub k 1))
 269           (setf fn (- fn 1.0))
 270           (if (>= iflag 3) (go label100))
 271           (setf str (abs c2r))
 272           (setf sti (abs c2i))
 273           (setf c2m (max str sti))
 274           (if (<= c2m ascle) (go label100))
 275           (setf iflag (f2cl-lib:int-add iflag 1))
 276           (setf ascle (f2cl-lib:fref bry (iflag) ((1 3))))
 277           (setf s1r (* s1r c1r))
 278           (setf s1i (* s1i c1r))
 279           (setf s2r c2r)
 280           (setf s2i c2i)
 281           (setf s1r (* s1r (f2cl-lib:fref cssr (iflag) ((1 3)))))
 282           (setf s1i (* s1i (f2cl-lib:fref cssr (iflag) ((1 3)))))
 283           (setf s2r (* s2r (f2cl-lib:fref cssr (iflag) ((1 3)))))
 284           (setf s2i (* s2i (f2cl-lib:fref cssr (iflag) ((1 3)))))
 285           (setf c1r (f2cl-lib:fref csrr (iflag) ((1 3))))
 286          label100))
 287      label110
 288       (go end_label)
 289      label120
 290       (if (> rs1 0.0) (go label140))
 291       (setf (f2cl-lib:fref yr (nd) ((1 n))) zeror)
 292       (setf (f2cl-lib:fref yi (nd) ((1 n))) zeroi)
 293       (setf nz (f2cl-lib:int-add nz 1))
 294       (setf nd (f2cl-lib:int-sub nd 1))
 295       (if (= nd 0) (go label110))
 296       (multiple-value-bind
 297             (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10
 298              var-11)
 299           (zuoik zr zi fnu kode 1 nd yr yi nuf tol elim alim)
 300         (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-9
 301                          var-10 var-11))
 302         (setf nuf var-8))
 303       (if (< nuf 0) (go label140))
 304       (setf nd (f2cl-lib:int-sub nd nuf))
 305       (setf nz (f2cl-lib:int-add nz nuf))
 306       (if (= nd 0) (go label110))
 307       (setf fn (+ fnu (f2cl-lib:int-sub nd 1)))
 308       (if (< fn fnul) (go label130))
 309       (setf in (f2cl-lib:int-sub (f2cl-lib:int-add inu nd) 1))
 310       (setf in (f2cl-lib:int-add (mod in 4) 1))
 311       (setf c2r
 312               (- (* car$ (f2cl-lib:fref cipr (in) ((1 4))))
 313                  (* sar (f2cl-lib:fref cipi (in) ((1 4))))))
 314       (setf c2i
 315               (+ (* car$ (f2cl-lib:fref cipi (in) ((1 4))))
 316                  (* sar (f2cl-lib:fref cipr (in) ((1 4))))))
 317       (if (<= zi 0.0) (setf c2i (- c2i)))
 318       (go label40)
 319      label130
 320       (setf nlast nd)
 321       (go end_label)
 322      label140
 323       (setf nz -1)
 324       (go end_label)
 325      label150
 326       (if (> rs1 0.0) (go label140))
 327       (setf nz n)
 328       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 329                     ((> i n) nil)
 330         (tagbody
 331           (setf (f2cl-lib:fref yr (i) ((1 n))) zeror)
 332           (setf (f2cl-lib:fref yi (i) ((1 n))) zeroi)
 333          label160))
 334       (go end_label)
 335      end_label
 336       (return (values nil nil nil nil nil nil nil nz nlast nil nil nil nil)))))
 337
 338 (in-package #:cl-user)
 339 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 340 (eval-when (:load-toplevel :compile-toplevel :execute)
 341   (setf (gethash 'fortran-to-lisp::zuni2 fortran-to-lisp::*f2cl-function-info*)
 342           (fortran-to-lisp::make-f2cl-finfo
 343            :arg-types '((double-float) (double-float) (double-float)
 344                         (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 345                         (simple-array double-float (*))
 346                         (simple-array double-float (*))
 347                         (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 348                         (double-float) (double-float) (double-float)
 349                         (double-float))
 350            :return-values '(nil nil nil nil nil nil nil fortran-to-lisp::nz
 351                             fortran-to-lisp::nlast nil nil nil nil)
 352            :calls '(fortran-to-lisp::zuoik fortran-to-lisp::zuchk
 353                     fortran-to-lisp::zairy fortran-to-lisp::zabs
 354                     fortran-to-lisp::zunhj fortran-to-lisp::d1mach))))
 355