src/numerical/slatec/zrati.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 ((czeror 0.0)
  21       (czeroi 0.0)
  22       (coner 1.0)
  23       (conei 0.0)
  24       (rt2 1.4142135623730951))
  25   (declare (type (double-float) czeror czeroi coner conei rt2))
  26   (defun zrati (zr zi fnu n cyr cyi tol)
  27     (declare (type (simple-array double-float (*)) cyi cyr)
  28              (type (f2cl-lib:integer4) n)
  29              (type (double-float) tol fnu zi zr))
  30     (prog ((i 0) (id 0) (idnu 0) (inu 0) (itime 0) (k 0) (kk 0) (magz 0)
  31            (ak 0.0) (amagz 0.0) (ap1 0.0) (ap2 0.0) (arg 0.0) (az 0.0)
  32            (cdfnui 0.0) (cdfnur 0.0) (dfnu 0.0) (fdnu 0.0) (flam 0.0)
  33            (fnup 0.0) (pti 0.0) (ptr 0.0) (p1i 0.0) (p1r 0.0) (p2i 0.0)
  34            (p2r 0.0) (rak 0.0) (rap1 0.0) (rho 0.0) (rzi 0.0) (rzr 0.0)
  35            (test 0.0) (test1 0.0) (tti 0.0) (ttr 0.0) (t1i 0.0) (t1r 0.0))
  36       (declare (type (double-float) t1r t1i ttr tti test1 test rzr rzi rho rap1
  37                                     rak p2r p2i p1r p1i ptr pti fnup flam fdnu
  38                                     dfnu cdfnur cdfnui az arg ap2 ap1 amagz ak)
  39                (type (f2cl-lib:integer4) magz kk k itime inu idnu id i))
  40       (setf az (coerce (realpart (zabs zr zi)) 'double-float))
  41       (setf inu (f2cl-lib:int fnu))
  42       (setf idnu (f2cl-lib:int-sub (f2cl-lib:int-add inu n) 1))
  43       (setf magz (f2cl-lib:int az))
  44       (setf amagz
  45               (coerce (the f2cl-lib:integer4 (f2cl-lib:int-add magz 1))
  46                       'double-float))
  47       (setf fdnu (coerce (the f2cl-lib:integer4 idnu) 'double-float))
  48       (setf fnup (max amagz fdnu))
  49       (setf id (f2cl-lib:int-sub idnu magz 1))
  50       (setf itime 1)
  51       (setf k 1)
  52       (setf ptr (/ 1.0 az))
  53       (setf rzr (* ptr (+ zr zr) ptr))
  54       (setf rzi (* (- ptr) (+ zi zi) ptr))
  55       (setf t1r (* rzr fnup))
  56       (setf t1i (* rzi fnup))
  57       (setf p2r (- t1r))
  58       (setf p2i (- t1i))
  59       (setf p1r coner)
  60       (setf p1i conei)
  61       (setf t1r (+ t1r rzr))
  62       (setf t1i (+ t1i rzi))
  63       (if (> id 0) (setf id 0))
  64       (setf ap2 (coerce (realpart (zabs p2r p2i)) 'double-float))
  65       (setf ap1 (coerce (realpart (zabs p1r p1i)) 'double-float))
  66       (setf arg (/ (+ ap2 ap2) (* ap1 tol)))
  67       (setf test1 (f2cl-lib:fsqrt arg))
  68       (setf test test1)
  69       (setf rap1 (/ 1.0 ap1))
  70       (setf p1r (* p1r rap1))
  71       (setf p1i (* p1i rap1))
  72       (setf p2r (* p2r rap1))
  73       (setf p2i (* p2i rap1))
  74       (setf ap2 (* ap2 rap1))
  75      label10
  76       (setf k (f2cl-lib:int-add k 1))
  77       (setf ap1 ap2)
  78       (setf ptr p2r)
  79       (setf pti p2i)
  80       (setf p2r (- p1r (- (* t1r ptr) (* t1i pti))))
  81       (setf p2i (- p1i (+ (* t1r pti) (* t1i ptr))))
  82       (setf p1r ptr)
  83       (setf p1i pti)
  84       (setf t1r (+ t1r rzr))
  85       (setf t1i (+ t1i rzi))
  86       (setf ap2 (coerce (realpart (zabs p2r p2i)) 'double-float))
  87       (if (<= ap1 test) (go label10))
  88       (if (= itime 2) (go label20))
  89       (setf ak (coerce (realpart (* (zabs t1r t1i) 0.5)) 'double-float))
  90       (setf flam (+ ak (f2cl-lib:fsqrt (- (* ak ak) 1.0))))
  91       (setf rho (min (/ ap2 ap1) flam))
  92       (setf test (* test1 (f2cl-lib:fsqrt (/ rho (- (* rho rho) 1.0)))))
  93       (setf itime 2)
  94       (go label10)
  95      label20
  96       (setf kk (f2cl-lib:int-sub (f2cl-lib:int-add k 1) id))
  97       (setf ak (coerce (the f2cl-lib:integer4 kk) 'double-float))
  98       (setf t1r ak)
  99       (setf t1i czeroi)
 100       (setf dfnu (+ fnu (f2cl-lib:int-sub n 1)))
 101       (setf p1r (/ 1.0 ap2))
 102       (setf p1i czeroi)
 103       (setf p2r czeror)
 104       (setf p2i czeroi)
 105       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 106                     ((> i kk) nil)
 107         (tagbody
 108           (setf ptr p1r)
 109           (setf pti p1i)
 110           (setf rap1 (+ dfnu t1r))
 111           (setf ttr (* rzr rap1))
 112           (setf tti (* rzi rap1))
 113           (setf p1r (+ (- (* ptr ttr) (* pti tti)) p2r))
 114           (setf p1i (+ (* ptr tti) (* pti ttr) p2i))
 115           (setf p2r ptr)
 116           (setf p2i pti)
 117           (setf t1r (- t1r coner))
 118          label30))
 119       (if (or (/= p1r czeror) (/= p1i czeroi)) (go label40))
 120       (setf p1r tol)
 121       (setf p1i tol)
 122      label40
 123       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 124           (zdiv p2r p2i p1r p1i (f2cl-lib:fref cyr (n) ((1 n)))
 125            (f2cl-lib:fref cyi (n) ((1 n))))
 126         (declare (ignore var-0 var-1 var-2 var-3))
 127         (setf (f2cl-lib:fref cyr (n) ((1 n))) var-4)
 128         (setf (f2cl-lib:fref cyi (n) ((1 n))) var-5))
 129       (if (= n 1) (go end_label))
 130       (setf k (f2cl-lib:int-sub n 1))
 131       (setf ak (coerce (the f2cl-lib:integer4 k) 'double-float))
 132       (setf t1r ak)
 133       (setf t1i czeroi)
 134       (setf cdfnur (* fnu rzr))
 135       (setf cdfnui (* fnu rzi))
 136       (f2cl-lib:fdo (i 2 (f2cl-lib:int-add i 1))
 137                     ((> i n) nil)
 138         (tagbody
 139           (setf ptr
 140                   (+ cdfnur
 141                      (- (* t1r rzr) (* t1i rzi))
 142                      (f2cl-lib:fref cyr ((f2cl-lib:int-add k 1)) ((1 n)))))
 143           (setf pti
 144                   (+ cdfnui
 145                      (+ (* t1r rzi) (* t1i rzr))
 146                      (f2cl-lib:fref cyi ((f2cl-lib:int-add k 1)) ((1 n)))))
 147           (setf ak (coerce (realpart (zabs ptr pti)) 'double-float))
 148           (if (/= ak czeror) (go label50))
 149           (setf ptr tol)
 150           (setf pti tol)
 151           (setf ak (* tol rt2))
 152          label50
 153           (setf rak (/ coner ak))
 154           (setf (f2cl-lib:fref cyr (k) ((1 n))) (* rak ptr rak))
 155           (setf (f2cl-lib:fref cyi (k) ((1 n))) (* (- rak) pti rak))
 156           (setf t1r (- t1r coner))
 157           (setf k (f2cl-lib:int-sub k 1))
 158          label60))
 159       (go end_label)
 160      end_label
 161       (return (values nil nil nil nil nil nil nil)))))
 162
 163 (in-package #:cl-user)
 164 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 165 (eval-when (:load-toplevel :compile-toplevel :execute)
 166   (setf (gethash 'fortran-to-lisp::zrati fortran-to-lisp::*f2cl-function-info*)
 167           (fortran-to-lisp::make-f2cl-finfo
 168            :arg-types '((double-float) (double-float) (double-float)
 169                         (fortran-to-lisp::integer4)
 170                         (simple-array double-float (*))
 171                         (simple-array double-float (*)) (double-float))
 172            :return-values '(nil nil nil nil nil nil nil)
 173            :calls '(fortran-to-lisp::zdiv fortran-to-lisp::zabs))))
 174