src/numerical/slatec/zasyi.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 ((pi$ 3.141592653589793)
  21       (rtpi 0.15915494309189535)
  22       (zeror 0.0)
  23       (zeroi 0.0)
  24       (coner 1.0)
  25       (conei 0.0))
  26   (declare (type (double-float) pi$ rtpi zeror zeroi coner conei))
  27   (defun zasyi (zr zi fnu kode n yr yi nz rl tol elim alim)
  28     (declare (type (simple-array double-float (*)) yi yr)
  29              (type (f2cl-lib:integer4) nz n kode)
  30              (type (double-float) alim elim tol rl fnu zi zr))
  31     (prog ((i 0) (ib 0) (il 0) (inu 0) (j 0) (jl 0) (k 0) (koded 0) (m 0)
  32            (nn 0) (aa 0.0) (aez 0.0) (ak 0.0) (ak1i 0.0) (ak1r 0.0) (arg 0.0)
  33            (arm 0.0) (atol 0.0) (az 0.0) (bb 0.0) (bk 0.0) (cki 0.0) (ckr 0.0)
  34            (cs1i 0.0) (cs1r 0.0) (cs2i 0.0) (cs2r 0.0) (czi 0.0) (czr 0.0)
  35            (dfnu 0.0) (dki 0.0) (dkr 0.0) (dnu2 0.0) (ezi 0.0) (ezr 0.0)
  36            (fdn 0.0) (p1i 0.0) (p1r 0.0) (raz 0.0) (rtr1 0.0) (rzi 0.0)
  37            (rzr 0.0) (s 0.0) (sgn 0.0) (sqk 0.0) (sti 0.0) (str 0.0) (s2i 0.0)
  38            (s2r 0.0) (tzi 0.0) (tzr 0.0))
  39       (declare (type (double-float) tzr tzi s2r s2i str sti sqk sgn s rzr rzi
  40                                     rtr1 raz p1r p1i fdn ezr ezi dnu2 dkr dki
  41                                     dfnu czr czi cs2r cs2i cs1r cs1i ckr cki bk
  42                                     bb az atol arm arg ak1r ak1i ak aez aa)
  43                (type (f2cl-lib:integer4) nn m koded k jl j inu il ib i))
  44       (setf nz 0)
  45       (setf az (coerce (realpart (zabs zr zi)) 'double-float))
  46       (setf arm (* 1000.0 (f2cl-lib:d1mach 1)))
  47       (setf rtr1 (f2cl-lib:fsqrt arm))
  48       (setf il (min (the f2cl-lib:integer4 2) (the f2cl-lib:integer4 n)))
  49       (setf dfnu (+ fnu (f2cl-lib:int-sub n il)))
  50       (setf raz (/ 1.0 az))
  51       (setf str (* zr raz))
  52       (setf sti (* (- zi) raz))
  53       (setf ak1r (* rtpi str raz))
  54       (setf ak1i (* rtpi sti raz))
  55       (multiple-value-bind (var-0 var-1 var-2 var-3)
  56           (zsqrt$ ak1r ak1i ak1r ak1i)
  57         (declare (ignore var-0 var-1))
  58         (setf ak1r var-2)
  59         (setf ak1i var-3))
  60       (setf czr zr)
  61       (setf czi zi)
  62       (if (/= kode 2) (go label10))
  63       (setf czr zeror)
  64       (setf czi zi)
  65      label10
  66       (if (> (abs czr) elim) (go label100))
  67       (setf dnu2 (+ dfnu dfnu))
  68       (setf koded 1)
  69       (if (and (> (abs czr) alim) (> n 2)) (go label20))
  70       (setf koded 0)
  71       (multiple-value-bind (var-0 var-1 var-2 var-3)
  72           (zexp czr czi str sti)
  73         (declare (ignore var-0 var-1))
  74         (setf str var-2)
  75         (setf sti var-3))
  76       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
  77           (zmlt ak1r ak1i str sti ak1r ak1i)
  78         (declare (ignore var-0 var-1 var-2 var-3))
  79         (setf ak1r var-4)
  80         (setf ak1i var-5))
  81      label20
  82       (setf fdn 0.0)
  83       (if (> dnu2 rtr1) (setf fdn (* dnu2 dnu2)))
  84       (setf ezr (* zr 8.0))
  85       (setf ezi (* zi 8.0))
  86       (setf aez (* 8.0 az))
  87       (setf s (/ tol aez))
  88       (setf jl (f2cl-lib:int (+ rl rl 2)))
  89       (setf p1r zeror)
  90       (setf p1i zeroi)
  91       (if (= zi 0.0) (go label30))
  92       (setf inu (f2cl-lib:int fnu))
  93       (setf arg (* (- fnu inu) pi$))
  94       (setf inu (f2cl-lib:int-sub (f2cl-lib:int-add inu n) il))
  95       (setf ak (- (sin arg)))
  96       (setf bk (cos arg))
  97       (if (< zi 0.0) (setf bk (- bk)))
  98       (setf p1r ak)
  99       (setf p1i bk)
 100       (if (= (mod inu 2) 0) (go label30))
 101       (setf p1r (- p1r))
 102       (setf p1i (- p1i))
 103      label30
 104       (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
 105                     ((> k il) nil)
 106         (tagbody
 107           (setf sqk (- fdn 1.0))
 108           (setf atol (* s (abs sqk)))
 109           (setf sgn 1.0)
 110           (setf cs1r coner)
 111           (setf cs1i conei)
 112           (setf cs2r coner)
 113           (setf cs2i conei)
 114           (setf ckr coner)
 115           (setf cki conei)
 116           (setf ak 0.0)
 117           (setf aa 1.0)
 118           (setf bb aez)
 119           (setf dkr ezr)
 120           (setf dki ezi)
 121           (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 122                         ((> j jl) nil)
 123             (tagbody
 124               (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 125                   (zdiv ckr cki dkr dki str sti)
 126                 (declare (ignore var-0 var-1 var-2 var-3))
 127                 (setf str var-4)
 128                 (setf sti var-5))
 129               (setf ckr (* str sqk))
 130               (setf cki (* sti sqk))
 131               (setf cs2r (+ cs2r ckr))
 132               (setf cs2i (+ cs2i cki))
 133               (setf sgn (- sgn))
 134               (setf cs1r (+ cs1r (* ckr sgn)))
 135               (setf cs1i (+ cs1i (* cki sgn)))
 136               (setf dkr (+ dkr ezr))
 137               (setf dki (+ dki ezi))
 138               (setf aa (/ (* aa (abs sqk)) bb))
 139               (setf bb (+ bb aez))
 140               (setf ak (+ ak 8.0))
 141               (setf sqk (- sqk ak))
 142               (if (<= aa atol) (go label50))
 143              label40))
 144           (go label110)
 145          label50
 146           (setf s2r cs1r)
 147           (setf s2i cs1i)
 148           (if (>= (+ zr zr) elim) (go label60))
 149           (setf tzr (+ zr zr))
 150           (setf tzi (+ zi zi))
 151           (multiple-value-bind (var-0 var-1 var-2 var-3)
 152               (zexp (- tzr) (- tzi) str sti)
 153             (declare (ignore var-0 var-1))
 154             (setf str var-2)
 155             (setf sti var-3))
 156           (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 157               (zmlt str sti p1r p1i str sti)
 158             (declare (ignore var-0 var-1 var-2 var-3))
 159             (setf str var-4)
 160             (setf sti var-5))
 161           (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 162               (zmlt str sti cs2r cs2i str sti)
 163             (declare (ignore var-0 var-1 var-2 var-3))
 164             (setf str var-4)
 165             (setf sti var-5))
 166           (setf s2r (+ s2r str))
 167           (setf s2i (+ s2i sti))
 168          label60
 169           (setf fdn (+ fdn (* 8.0 dfnu) 4.0))
 170           (setf p1r (- p1r))
 171           (setf p1i (- p1i))
 172           (setf m (f2cl-lib:int-add (f2cl-lib:int-sub n il) k))
 173           (setf (f2cl-lib:fref yr (m) ((1 n))) (- (* s2r ak1r) (* s2i ak1i)))
 174           (setf (f2cl-lib:fref yi (m) ((1 n))) (+ (* s2r ak1i) (* s2i ak1r)))
 175          label70))
 176       (if (<= n 2) (go end_label))
 177       (setf nn n)
 178       (setf k (f2cl-lib:int-sub nn 2))
 179       (setf ak (coerce (the f2cl-lib:integer4 k) 'double-float))
 180       (setf str (* zr raz))
 181       (setf sti (* (- zi) raz))
 182       (setf rzr (* (+ str str) raz))
 183       (setf rzi (* (+ sti sti) raz))
 184       (setf ib 3)
 185       (f2cl-lib:fdo (i ib (f2cl-lib:int-add i 1))
 186                     ((> i nn) nil)
 187         (tagbody
 188           (setf (f2cl-lib:fref yr (k) ((1 n)))
 189                   (+
 190                    (* (+ ak fnu)
 191                       (-
 192                        (* rzr
 193                           (f2cl-lib:fref yr ((f2cl-lib:int-add k 1)) ((1 n))))
 194                        (* rzi
 195                           (f2cl-lib:fref yi
 196                                          ((f2cl-lib:int-add k 1))
 197                                          ((1 n))))))
 198                    (f2cl-lib:fref yr ((f2cl-lib:int-add k 2)) ((1 n)))))
 199           (setf (f2cl-lib:fref yi (k) ((1 n)))
 200                   (+
 201                    (* (+ ak fnu)
 202                       (+
 203                        (* rzr
 204                           (f2cl-lib:fref yi ((f2cl-lib:int-add k 1)) ((1 n))))
 205                        (* rzi
 206                           (f2cl-lib:fref yr
 207                                          ((f2cl-lib:int-add k 1))
 208                                          ((1 n))))))
 209                    (f2cl-lib:fref yi ((f2cl-lib:int-add k 2)) ((1 n)))))
 210           (setf ak (- ak 1.0))
 211           (setf k (f2cl-lib:int-sub k 1))
 212          label80))
 213       (if (= koded 0) (go end_label))
 214       (multiple-value-bind (var-0 var-1 var-2 var-3)
 215           (zexp czr czi ckr cki)
 216         (declare (ignore var-0 var-1))
 217         (setf ckr var-2)
 218         (setf cki var-3))
 219       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 220                     ((> i nn) nil)
 221         (tagbody
 222           (setf str
 223                   (- (* (f2cl-lib:fref yr (i) ((1 n))) ckr)
 224                      (* (f2cl-lib:fref yi (i) ((1 n))) cki)))
 225           (setf (f2cl-lib:fref yi (i) ((1 n)))
 226                   (+ (* (f2cl-lib:fref yr (i) ((1 n))) cki)
 227                      (* (f2cl-lib:fref yi (i) ((1 n))) ckr)))
 228           (setf (f2cl-lib:fref yr (i) ((1 n))) str)
 229          label90))
 230       (go end_label)
 231      label100
 232       (setf nz -1)
 233       (go end_label)
 234      label110
 235       (setf nz -2)
 236       (go end_label)
 237      end_label
 238       (return (values nil nil nil nil nil nil nil nz nil nil nil nil)))))
 239
 240 (in-package #:cl-user)
 241 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 242 (eval-when (:load-toplevel :compile-toplevel :execute)
 243   (setf (gethash 'fortran-to-lisp::zasyi fortran-to-lisp::*f2cl-function-info*)
 244           (fortran-to-lisp::make-f2cl-finfo
 245            :arg-types '((double-float) (double-float) (double-float)
 246                         (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 247                         (simple-array double-float (*))
 248                         (simple-array double-float (*))
 249                         (fortran-to-lisp::integer4) (double-float)
 250                         (double-float) (double-float) (double-float))
 251            :return-values '(nil nil nil nil nil nil nil fortran-to-lisp::nz nil
 252                             nil nil nil)
 253            :calls '(fortran-to-lisp::zdiv fortran-to-lisp::zmlt
 254                     fortran-to-lisp::zexp fortran-to-lisp::zsqrt$
 255                     fortran-to-lisp::d1mach fortran-to-lisp::zabs))))
 256