src/numerical/slatec/zbknu.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 ((kmax 30)
  21       (czeror 0.0)
  22       (czeroi 0.0)
  23       (coner 1.0)
  24       (conei 0.0)
  25       (ctwor 2.0)
  26       (r1 2.0)
  27       (dpi 3.141592653589793)
  28       (rthpi 1.2533141373155003)
  29       (spi 1.909859317102744)
  30       (hpi 1.5707963267948966)
  31       (fpi 1.8976999933151775)
  32       (tth 0.6666666666666666)
  33       (cc
  34        (make-array 8
  35                    :element-type 'double-float
  36                    :initial-contents '(0.5772156649015329 -0.04200263503409524
  37                                        -0.04219773455554433 0.0072189432466631
  38                                        -2.1524167411495098e-4
  39                                        -2.013485478078824e-5
  40                                        1.133027231981696e-6
  41                                        6.116095104481416e-9))))
  42   (declare (type (f2cl-lib:integer4) kmax)
  43            (type (double-float) czeror czeroi coner conei ctwor r1 dpi rthpi
  44                                 spi hpi fpi tth)
  45            (type (simple-array double-float (8)) cc))
  46   (defun zbknu (zr zi fnu kode n yr yi nz tol elim alim)
  47     (declare (type (simple-array double-float (*)) yi yr)
  48              (type (f2cl-lib:integer4) nz n kode)
  49              (type (double-float) alim elim tol fnu zi zr))
  50     (prog ((cssr (make-array 3 :element-type 'double-float))
  51            (csrr (make-array 3 :element-type 'double-float))
  52            (bry (make-array 3 :element-type 'double-float))
  53            (cyr (make-array 2 :element-type 'double-float))
  54            (cyi (make-array 2 :element-type 'double-float)) (i 0) (iflag 0)
  55            (inu 0) (k 0) (kflag 0) (kk 0) (koded 0) (idum 0) (j 0) (ic 0)
  56            (inub 0) (nw 0) (aa 0.0) (ak 0.0) (ascle 0.0) (a1 0.0) (a2 0.0)
  57            (bb 0.0) (bk 0.0) (caz 0.0) (cbi 0.0) (cbr 0.0) (cchi 0.0)
  58            (cchr 0.0) (cki 0.0) (ckr 0.0) (coefi 0.0) (coefr 0.0) (crscr 0.0)
  59            (csclr 0.0) (cshi 0.0) (cshr 0.0) (csi 0.0) (csr 0.0) (czi 0.0)
  60            (czr 0.0) (dnu 0.0) (dnu2 0.0) (etest 0.0) (fc 0.0) (fhs 0.0)
  61            (fi 0.0) (fk 0.0) (fks 0.0) (fmui 0.0) (fmur 0.0) (fr 0.0) (g1 0.0)
  62            (g2 0.0) (pi$ 0.0) (pr 0.0) (pti 0.0) (ptr 0.0) (p1i 0.0) (p1r 0.0)
  63            (p2i 0.0) (p2m 0.0) (p2r 0.0) (qi 0.0) (qr 0.0) (rak 0.0) (rcaz 0.0)
  64            (rzi 0.0) (rzr 0.0) (s 0.0) (smui 0.0) (smur 0.0) (sti 0.0)
  65            (str 0.0) (s1i 0.0) (s1r 0.0) (s2i 0.0) (s2r 0.0) (tm 0.0) (t1 0.0)
  66            (t2 0.0) (elm 0.0) (celmr 0.0) (zdr 0.0) (zdi 0.0) (as 0.0)
  67            (alas 0.0) (helim 0.0))
  68       (declare (type (simple-array double-float (2)) cyi cyr)
  69                (type (simple-array double-float (3)) cssr csrr bry)
  70                (type (double-float) helim alas as zdi zdr celmr elm t2 t1 tm
  71                                     s2r s2i s1r s1i str sti smur smui s rzr rzi
  72                                     rcaz rak qr qi p2r p2m p2i p1r p1i ptr pti
  73                                     pr pi$ g2 g1 fr fmur fmui fks fk fi fhs fc
  74                                     etest dnu2 dnu czr czi csr csi cshr cshi
  75                                     csclr crscr coefr coefi ckr cki cchr cchi
  76                                     cbr cbi caz bk bb a2 a1 ascle ak aa)
  77                (type (f2cl-lib:integer4) nw inub ic j idum koded kk kflag k inu
  78                                          iflag i))
  79       (setf caz (coerce (realpart (zabs zr zi)) 'double-float))
  80       (setf csclr (/ 1.0 tol))
  81       (setf crscr tol)
  82       (setf (f2cl-lib:fref cssr (1) ((1 3))) csclr)
  83       (setf (f2cl-lib:fref cssr (2) ((1 3))) 1.0)
  84       (setf (f2cl-lib:fref cssr (3) ((1 3))) crscr)
  85       (setf (f2cl-lib:fref csrr (1) ((1 3))) crscr)
  86       (setf (f2cl-lib:fref csrr (2) ((1 3))) 1.0)
  87       (setf (f2cl-lib:fref csrr (3) ((1 3))) csclr)
  88       (setf (f2cl-lib:fref bry (1) ((1 3)))
  89               (/ (* 1000.0 (f2cl-lib:d1mach 1)) tol))
  90       (setf (f2cl-lib:fref bry (2) ((1 3)))
  91               (/ 1.0 (f2cl-lib:fref bry (1) ((1 3)))))
  92       (setf (f2cl-lib:fref bry (3) ((1 3))) (f2cl-lib:d1mach 2))
  93       (setf nz 0)
  94       (setf iflag 0)
  95       (setf koded kode)
  96       (setf rcaz (/ 1.0 caz))
  97       (setf str (* zr rcaz))
  98       (setf sti (* (- zi) rcaz))
  99       (setf rzr (* (+ str str) rcaz))
 100       (setf rzi (* (+ sti sti) rcaz))
 101       (setf inu (f2cl-lib:int (+ fnu 0.5)))
 102       (setf dnu (- fnu inu))
 103       (if (= (abs dnu) 0.5) (go label110))
 104       (setf dnu2 0.0)
 105       (if (> (abs dnu) tol) (setf dnu2 (* dnu dnu)))
 106       (if (> caz r1) (go label110))
 107       (setf fc 1.0)
 108       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4)
 109           (zlog rzr rzi smur smui idum)
 110         (declare (ignore var-0 var-1))
 111         (setf smur var-2)
 112         (setf smui var-3)
 113         (setf idum var-4))
 114       (setf fmur (* smur dnu))
 115       (setf fmui (* smui dnu))
 116       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 117           (zshch fmur fmui cshr cshi cchr cchi)
 118         (declare (ignore var-0 var-1))
 119         (setf cshr var-2)
 120         (setf cshi var-3)
 121         (setf cchr var-4)
 122         (setf cchi var-5))
 123       (if (= dnu 0.0) (go label10))
 124       (setf fc (* dnu dpi))
 125       (setf fc (/ fc (sin fc)))
 126       (setf smur (/ cshr dnu))
 127       (setf smui (/ cshi dnu))
 128      label10
 129       (setf a2 (+ 1.0 dnu))
 130       (setf t2
 131               (exp
 132                (-
 133                 (multiple-value-bind (ret-val var-0 var-1)
 134                     (dgamln a2 idum)
 135                   (declare (ignore var-0))
 136                   (setf idum var-1)
 137                   ret-val))))
 138       (setf t1 (/ 1.0 (* t2 fc)))
 139       (if (> (abs dnu) 0.1) (go label40))
 140       (setf ak 1.0)
 141       (setf s (f2cl-lib:fref cc (1) ((1 8))))
 142       (f2cl-lib:fdo (k 2 (f2cl-lib:int-add k 1))
 143                     ((> k 8) nil)
 144         (tagbody
 145           (setf ak (* ak dnu2))
 146           (setf tm (* (f2cl-lib:fref cc (k) ((1 8))) ak))
 147           (setf s (+ s tm))
 148           (if (< (abs tm) tol) (go label30))
 149          label20))
 150      label30
 151       (setf g1 (- s))
 152       (go label50)
 153      label40
 154       (setf g1 (/ (- t1 t2) (+ dnu dnu)))
 155      label50
 156       (setf g2 (* (+ t1 t2) 0.5))
 157       (setf fr (* fc (+ (* cchr g1) (* smur g2))))
 158       (setf fi (* fc (+ (* cchi g1) (* smui g2))))
 159       (multiple-value-bind (var-0 var-1 var-2 var-3)
 160           (zexp fmur fmui str sti)
 161         (declare (ignore var-0 var-1))
 162         (setf str var-2)
 163         (setf sti var-3))
 164       (setf pr (/ (* 0.5 str) t2))
 165       (setf pi$ (/ (* 0.5 sti) t2))
 166       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 167           (zdiv 0.5 0.0 str sti ptr pti)
 168         (declare (ignore var-0 var-1 var-2 var-3))
 169         (setf ptr var-4)
 170         (setf pti var-5))
 171       (setf qr (/ ptr t1))
 172       (setf qi (/ pti t1))
 173       (setf s1r fr)
 174       (setf s1i fi)
 175       (setf s2r pr)
 176       (setf s2i pi$)
 177       (setf ak 1.0)
 178       (setf a1 1.0)
 179       (setf ckr coner)
 180       (setf cki conei)
 181       (setf bk (- 1.0 dnu2))
 182       (if (or (> inu 0) (> n 1)) (go label80))
 183       (if (< caz tol) (go label70))
 184       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 185           (zmlt zr zi zr zi czr czi)
 186         (declare (ignore var-0 var-1 var-2 var-3))
 187         (setf czr var-4)
 188         (setf czi var-5))
 189       (setf czr (* 0.25 czr))
 190       (setf czi (* 0.25 czi))
 191       (setf t1 (* 0.25 caz caz))
 192      label60
 193       (setf fr (/ (+ (* fr ak) pr qr) bk))
 194       (setf fi (/ (+ (* fi ak) pi$ qi) bk))
 195       (setf str (/ 1.0 (- ak dnu)))
 196       (setf pr (* pr str))
 197       (setf pi$ (* pi$ str))
 198       (setf str (/ 1.0 (+ ak dnu)))
 199       (setf qr (* qr str))
 200       (setf qi (* qi str))
 201       (setf str (- (* ckr czr) (* cki czi)))
 202       (setf rak (/ 1.0 ak))
 203       (setf cki (* (+ (* ckr czi) (* cki czr)) rak))
 204       (setf ckr (* str rak))
 205       (setf s1r (+ (- (* ckr fr) (* cki fi)) s1r))
 206       (setf s1i (+ (* ckr fi) (* cki fr) s1i))
 207       (setf a1 (* a1 t1 rak))
 208       (setf bk (+ bk ak ak 1.0))
 209       (setf ak (+ ak 1.0))
 210       (if (> a1 tol) (go label60))
 211      label70
 212       (setf (f2cl-lib:fref yr (1) ((1 n))) s1r)
 213       (setf (f2cl-lib:fref yi (1) ((1 n))) s1i)
 214       (if (= koded 1) (go end_label))
 215       (multiple-value-bind (var-0 var-1 var-2 var-3)
 216           (zexp zr zi str sti)
 217         (declare (ignore var-0 var-1))
 218         (setf str var-2)
 219         (setf sti var-3))
 220       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 221           (zmlt s1r s1i str sti (f2cl-lib:fref yr (1) ((1 n)))
 222            (f2cl-lib:fref yi (1) ((1 n))))
 223         (declare (ignore var-0 var-1 var-2 var-3))
 224         (setf (f2cl-lib:fref yr (1) ((1 n))) var-4)
 225         (setf (f2cl-lib:fref yi (1) ((1 n))) var-5))
 226       (go end_label)
 227      label80
 228       (if (< caz tol) (go label100))
 229       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 230           (zmlt zr zi zr zi czr czi)
 231         (declare (ignore var-0 var-1 var-2 var-3))
 232         (setf czr var-4)
 233         (setf czi var-5))
 234       (setf czr (* 0.25 czr))
 235       (setf czi (* 0.25 czi))
 236       (setf t1 (* 0.25 caz caz))
 237      label90
 238       (setf fr (/ (+ (* fr ak) pr qr) bk))
 239       (setf fi (/ (+ (* fi ak) pi$ qi) bk))
 240       (setf str (/ 1.0 (- ak dnu)))
 241       (setf pr (* pr str))
 242       (setf pi$ (* pi$ str))
 243       (setf str (/ 1.0 (+ ak dnu)))
 244       (setf qr (* qr str))
 245       (setf qi (* qi str))
 246       (setf str (- (* ckr czr) (* cki czi)))
 247       (setf rak (/ 1.0 ak))
 248       (setf cki (* (+ (* ckr czi) (* cki czr)) rak))
 249       (setf ckr (* str rak))
 250       (setf s1r (+ (- (* ckr fr) (* cki fi)) s1r))
 251       (setf s1i (+ (* ckr fi) (* cki fr) s1i))
 252       (setf str (- pr (* fr ak)))
 253       (setf sti (- pi$ (* fi ak)))
 254       (setf s2r (+ (- (* ckr str) (* cki sti)) s2r))
 255       (setf s2i (+ (* ckr sti) (* cki str) s2i))
 256       (setf a1 (* a1 t1 rak))
 257       (setf bk (+ bk ak ak 1.0))
 258       (setf ak (+ ak 1.0))
 259       (if (> a1 tol) (go label90))
 260      label100
 261       (setf kflag 2)
 262       (setf a1 (+ fnu 1.0))
 263       (setf ak (* a1 (abs smur)))
 264       (if (> ak alim) (setf kflag 3))
 265       (setf str (f2cl-lib:fref cssr (kflag) ((1 3))))
 266       (setf p2r (* s2r str))
 267       (setf p2i (* s2i str))
 268       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 269           (zmlt p2r p2i rzr rzi s2r s2i)
 270         (declare (ignore var-0 var-1 var-2 var-3))
 271         (setf s2r var-4)
 272         (setf s2i var-5))
 273       (setf s1r (* s1r str))
 274       (setf s1i (* s1i str))
 275       (if (= koded 1) (go label210))
 276       (multiple-value-bind (var-0 var-1 var-2 var-3)
 277           (zexp zr zi fr fi)
 278         (declare (ignore var-0 var-1))
 279         (setf fr var-2)
 280         (setf fi var-3))
 281       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 282           (zmlt s1r s1i fr fi s1r s1i)
 283         (declare (ignore var-0 var-1 var-2 var-3))
 284         (setf s1r var-4)
 285         (setf s1i var-5))
 286       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 287           (zmlt s2r s2i fr fi s2r s2i)
 288         (declare (ignore var-0 var-1 var-2 var-3))
 289         (setf s2r var-4)
 290         (setf s2i var-5))
 291       (go label210)
 292      label110
 293       (multiple-value-bind (var-0 var-1 var-2 var-3)
 294           (zsqrt$ zr zi str sti)
 295         (declare (ignore var-0 var-1))
 296         (setf str var-2)
 297         (setf sti var-3))
 298       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 299           (zdiv rthpi czeroi str sti coefr coefi)
 300         (declare (ignore var-0 var-1 var-2 var-3))
 301         (setf coefr var-4)
 302         (setf coefi var-5))
 303       (setf kflag 2)
 304       (if (= koded 2) (go label120))
 305       (if (> zr alim) (go label290))
 306       (setf str (* (exp (- zr)) (f2cl-lib:fref cssr (kflag) ((1 3)))))
 307       (setf sti (* (- str) (sin zi)))
 308       (setf str (* str (cos zi)))
 309       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 310           (zmlt coefr coefi str sti coefr coefi)
 311         (declare (ignore var-0 var-1 var-2 var-3))
 312         (setf coefr var-4)
 313         (setf coefi var-5))
 314      label120
 315       (if (= (abs dnu) 0.5) (go label300))
 316       (setf ak (cos (* dpi dnu)))
 317       (setf ak (abs ak))
 318       (if (= ak czeror) (go label300))
 319       (setf fhs (abs (- 0.25 dnu2)))
 320       (if (= fhs czeror) (go label300))
 321       (setf t1
 322               (coerce
 323                (the f2cl-lib:integer4
 324                     (f2cl-lib:int-sub (f2cl-lib:i1mach 14) 1))
 325                'double-float))
 326       (setf t1 (* t1 (f2cl-lib:d1mach 5) 3.321928094))
 327       (setf t1 (max t1 12.0))
 328       (setf t1 (min t1 60.0))
 329       (setf t2 (- (* tth t1) 6.0))
 330       (if (/= zr 0.0) (go label130))
 331       (setf t1 hpi)
 332       (go label140)
 333      label130
 334       (setf t1 (f2cl-lib:datan (/ zi zr)))
 335       (setf t1 (abs t1))
 336      label140
 337       (if (> t2 caz) (go label170))
 338       (setf etest (/ ak (* dpi caz tol)))
 339       (setf fk coner)
 340       (if (< etest coner) (go label180))
 341       (setf fks ctwor)
 342       (setf ckr (+ caz caz ctwor))
 343       (setf p1r czeror)
 344       (setf p2r coner)
 345       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 346                     ((> i kmax) nil)
 347         (tagbody
 348           (setf ak (/ fhs fks))
 349           (setf cbr (/ ckr (+ fk coner)))
 350           (setf ptr p2r)
 351           (setf p2r (- (* cbr p2r) (* p1r ak)))
 352           (setf p1r ptr)
 353           (setf ckr (+ ckr ctwor))
 354           (setf fks (+ fks fk fk ctwor))
 355           (setf fhs (+ fhs fk fk))
 356           (setf fk (+ fk coner))
 357           (setf str (* (abs p2r) fk))
 358           (if (< etest str) (go label160))
 359          label150))
 360       (go label310)
 361      label160
 362       (setf fk (+ fk (* spi t1 (f2cl-lib:fsqrt (/ t2 caz)))))
 363       (setf fhs (abs (- 0.25 dnu2)))
 364       (go label180)
 365      label170
 366       (setf a2 (f2cl-lib:fsqrt caz))
 367       (setf ak (/ (* fpi ak) (* tol (f2cl-lib:fsqrt a2))))
 368       (setf aa (/ (* 3.0 t1) (+ 1.0 caz)))
 369       (setf bb (/ (* 14.7 t1) (+ 28.0 caz)))
 370       (setf ak
 371               (/
 372                (+ (f2cl-lib:flog ak)
 373                   (/ (* caz (cos aa)) (+ 1.0 (* 0.008 caz))))
 374                (cos bb)))
 375       (setf fk (+ (/ (* 0.12125 ak ak) caz) 1.5))
 376      label180
 377       (setf k (f2cl-lib:int fk))
 378       (setf fk (coerce (the f2cl-lib:integer4 k) 'double-float))
 379       (setf fks (* fk fk))
 380       (setf p1r czeror)
 381       (setf p1i czeroi)
 382       (setf p2r tol)
 383       (setf p2i czeroi)
 384       (setf csr p2r)
 385       (setf csi p2i)
 386       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 387                     ((> i k) nil)
 388         (tagbody
 389           (setf a1 (- fks fk))
 390           (setf ak (/ (+ fks fk) (+ a1 fhs)))
 391           (setf rak (/ 2.0 (+ fk coner)))
 392           (setf cbr (* (+ fk zr) rak))
 393           (setf cbi (* zi rak))
 394           (setf ptr p2r)
 395           (setf pti p2i)
 396           (setf p2r (* (- (* ptr cbr) (* pti cbi) p1r) ak))
 397           (setf p2i (* (- (+ (* pti cbr) (* ptr cbi)) p1i) ak))
 398           (setf p1r ptr)
 399           (setf p1i pti)
 400           (setf csr (+ csr p2r))
 401           (setf csi (+ csi p2i))
 402           (setf fks (+ (- a1 fk) coner))
 403           (setf fk (- fk coner))
 404          label190))
 405       (setf tm (coerce (realpart (zabs csr csi)) 'double-float))
 406       (setf ptr (/ 1.0 tm))
 407       (setf s1r (* p2r ptr))
 408       (setf s1i (* p2i ptr))
 409       (setf csr (* csr ptr))
 410       (setf csi (* (- csi) ptr))
 411       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 412           (zmlt coefr coefi s1r s1i str sti)
 413         (declare (ignore var-0 var-1 var-2 var-3))
 414         (setf str var-4)
 415         (setf sti var-5))
 416       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 417           (zmlt str sti csr csi s1r s1i)
 418         (declare (ignore var-0 var-1 var-2 var-3))
 419         (setf s1r var-4)
 420         (setf s1i var-5))
 421       (if (or (> inu 0) (> n 1)) (go label200))
 422       (setf zdr zr)
 423       (setf zdi zi)
 424       (if (= iflag 1) (go label270))
 425       (go label240)
 426      label200
 427       (setf tm (coerce (realpart (zabs p2r p2i)) 'double-float))
 428       (setf ptr (/ 1.0 tm))
 429       (setf p1r (* p1r ptr))
 430       (setf p1i (* p1i ptr))
 431       (setf p2r (* p2r ptr))
 432       (setf p2i (* (- p2i) ptr))
 433       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 434           (zmlt p1r p1i p2r p2i ptr pti)
 435         (declare (ignore var-0 var-1 var-2 var-3))
 436         (setf ptr var-4)
 437         (setf pti var-5))
 438       (setf str (- (+ dnu 0.5) ptr))
 439       (setf sti (- pti))
 440       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 441           (zdiv str sti zr zi str sti)
 442         (declare (ignore var-0 var-1 var-2 var-3))
 443         (setf str var-4)
 444         (setf sti var-5))
 445       (setf str (+ str 1.0))
 446       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 447           (zmlt str sti s1r s1i s2r s2i)
 448         (declare (ignore var-0 var-1 var-2 var-3))
 449         (setf s2r var-4)
 450         (setf s2i var-5))
 451      label210
 452       (setf str (+ dnu 1.0))
 453       (setf ckr (* str rzr))
 454       (setf cki (* str rzi))
 455       (if (= n 1) (setf inu (f2cl-lib:int-sub inu 1)))
 456       (if (> inu 0) (go label220))
 457       (if (> n 1) (go label215))
 458       (setf s1r s2r)
 459       (setf s1i s2i)
 460      label215
 461       (setf zdr zr)
 462       (setf zdi zi)
 463       (if (= iflag 1) (go label270))
 464       (go label240)
 465      label220
 466       (setf inub 1)
 467       (if (= iflag 1) (go label261))
 468      label225
 469       (setf p1r (f2cl-lib:fref csrr (kflag) ((1 3))))
 470       (setf ascle (f2cl-lib:fref bry (kflag) ((1 3))))
 471       (f2cl-lib:fdo (i inub (f2cl-lib:int-add i 1))
 472                     ((> i inu) nil)
 473         (tagbody
 474           (setf str s2r)
 475           (setf sti s2i)
 476           (setf s2r (+ (- (* ckr str) (* cki sti)) s1r))
 477           (setf s2i (+ (* ckr sti) (* cki str) s1i))
 478           (setf s1r str)
 479           (setf s1i sti)
 480           (setf ckr (+ ckr rzr))
 481           (setf cki (+ cki rzi))
 482           (if (>= kflag 3) (go label230))
 483           (setf p2r (* s2r p1r))
 484           (setf p2i (* s2i p1r))
 485           (setf str (abs p2r))
 486           (setf sti (abs p2i))
 487           (setf p2m (max str sti))
 488           (if (<= p2m ascle) (go label230))
 489           (setf kflag (f2cl-lib:int-add kflag 1))
 490           (setf ascle (f2cl-lib:fref bry (kflag) ((1 3))))
 491           (setf s1r (* s1r p1r))
 492           (setf s1i (* s1i p1r))
 493           (setf s2r p2r)
 494           (setf s2i p2i)
 495           (setf str (f2cl-lib:fref cssr (kflag) ((1 3))))
 496           (setf s1r (* s1r str))
 497           (setf s1i (* s1i str))
 498           (setf s2r (* s2r str))
 499           (setf s2i (* s2i str))
 500           (setf p1r (f2cl-lib:fref csrr (kflag) ((1 3))))
 501          label230))
 502       (if (/= n 1) (go label240))
 503       (setf s1r s2r)
 504       (setf s1i s2i)
 505      label240
 506       (setf str (f2cl-lib:fref csrr (kflag) ((1 3))))
 507       (setf (f2cl-lib:fref yr (1) ((1 n))) (* s1r str))
 508       (setf (f2cl-lib:fref yi (1) ((1 n))) (* s1i str))
 509       (if (= n 1) (go end_label))
 510       (setf (f2cl-lib:fref yr (2) ((1 n))) (* s2r str))
 511       (setf (f2cl-lib:fref yi (2) ((1 n))) (* s2i str))
 512       (if (= n 2) (go end_label))
 513       (setf kk 2)
 514      label250
 515       (setf kk (f2cl-lib:int-add kk 1))
 516       (if (> kk n) (go end_label))
 517       (setf p1r (f2cl-lib:fref csrr (kflag) ((1 3))))
 518       (setf ascle (f2cl-lib:fref bry (kflag) ((1 3))))
 519       (f2cl-lib:fdo (i kk (f2cl-lib:int-add i 1))
 520                     ((> i n) nil)
 521         (tagbody
 522           (setf p2r s2r)
 523           (setf p2i s2i)
 524           (setf s2r (+ (- (* ckr p2r) (* cki p2i)) s1r))
 525           (setf s2i (+ (* cki p2r) (* ckr p2i) s1i))
 526           (setf s1r p2r)
 527           (setf s1i p2i)
 528           (setf ckr (+ ckr rzr))
 529           (setf cki (+ cki rzi))
 530           (setf p2r (* s2r p1r))
 531           (setf p2i (* s2i p1r))
 532           (setf (f2cl-lib:fref yr (i) ((1 n))) p2r)
 533           (setf (f2cl-lib:fref yi (i) ((1 n))) p2i)
 534           (if (>= kflag 3) (go label260))
 535           (setf str (abs p2r))
 536           (setf sti (abs p2i))
 537           (setf p2m (max str sti))
 538           (if (<= p2m ascle) (go label260))
 539           (setf kflag (f2cl-lib:int-add kflag 1))
 540           (setf ascle (f2cl-lib:fref bry (kflag) ((1 3))))
 541           (setf s1r (* s1r p1r))
 542           (setf s1i (* s1i p1r))
 543           (setf s2r p2r)
 544           (setf s2i p2i)
 545           (setf str (f2cl-lib:fref cssr (kflag) ((1 3))))
 546           (setf s1r (* s1r str))
 547           (setf s1i (* s1i str))
 548           (setf s2r (* s2r str))
 549           (setf s2i (* s2i str))
 550           (setf p1r (f2cl-lib:fref csrr (kflag) ((1 3))))
 551          label260))
 552       (go end_label)
 553      label261
 554       (setf helim (* 0.5 elim))
 555       (setf elm (exp (- elim)))
 556       (setf celmr elm)
 557       (setf ascle (f2cl-lib:fref bry (1) ((1 3))))
 558       (setf zdr zr)
 559       (setf zdi zi)
 560       (setf ic -1)
 561       (setf j 2)
 562       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 563                     ((> i inu) nil)
 564         (tagbody
 565           (setf str s2r)
 566           (setf sti s2i)
 567           (setf s2r (+ (- (* str ckr) (* sti cki)) s1r))
 568           (setf s2i (+ (* sti ckr) (* str cki) s1i))
 569           (setf s1r str)
 570           (setf s1i sti)
 571           (setf ckr (+ ckr rzr))
 572           (setf cki (+ cki rzi))
 573           (setf as (coerce (realpart (zabs s2r s2i)) 'double-float))
 574           (setf alas (f2cl-lib:flog as))
 575           (setf p2r (- alas zdr))
 576           (if (< p2r (- elim)) (go label263))
 577           (multiple-value-bind (var-0 var-1 var-2 var-3 var-4)
 578               (zlog s2r s2i str sti idum)
 579             (declare (ignore var-0 var-1))
 580             (setf str var-2)
 581             (setf sti var-3)
 582             (setf idum var-4))
 583           (setf p2r (- str zdr))
 584           (setf p2i (- sti zdi))
 585           (setf p2m (/ (exp p2r) tol))
 586           (setf p1r (* p2m (cos p2i)))
 587           (setf p1i (* p2m (sin p2i)))
 588           (multiple-value-bind (var-0 var-1 var-2 var-3 var-4)
 589               (zuchk p1r p1i nw ascle tol)
 590             (declare (ignore var-0 var-1 var-3 var-4))
 591             (setf nw var-2))
 592           (if (/= nw 0) (go label263))
 593           (setf j (f2cl-lib:int-sub 3 j))
 594           (setf (f2cl-lib:fref cyr (j) ((1 2))) p1r)
 595           (setf (f2cl-lib:fref cyi (j) ((1 2))) p1i)
 596           (if (= ic (f2cl-lib:int-sub i 1)) (go label264))
 597           (setf ic i)
 598           (go label262)
 599          label263
 600           (if (< alas helim) (go label262))
 601           (setf zdr (- zdr elim))
 602           (setf s1r (* s1r celmr))
 603           (setf s1i (* s1i celmr))
 604           (setf s2r (* s2r celmr))
 605           (setf s2i (* s2i celmr))
 606          label262))
 607       (if (/= n 1) (go label270))
 608       (setf s1r s2r)
 609       (setf s1i s2i)
 610       (go label270)
 611      label264
 612       (setf kflag 1)
 613       (setf inub (f2cl-lib:int-add i 1))
 614       (setf s2r (f2cl-lib:fref cyr (j) ((1 2))))
 615       (setf s2i (f2cl-lib:fref cyi (j) ((1 2))))
 616       (setf j (f2cl-lib:int-sub 3 j))
 617       (setf s1r (f2cl-lib:fref cyr (j) ((1 2))))
 618       (setf s1i (f2cl-lib:fref cyi (j) ((1 2))))
 619       (if (<= inub inu) (go label225))
 620       (if (/= n 1) (go label240))
 621       (setf s1r s2r)
 622       (setf s1i s2i)
 623       (go label240)
 624      label270
 625       (setf (f2cl-lib:fref yr (1) ((1 n))) s1r)
 626       (setf (f2cl-lib:fref yi (1) ((1 n))) s1i)
 627       (if (= n 1) (go label280))
 628       (setf (f2cl-lib:fref yr (2) ((1 n))) s2r)
 629       (setf (f2cl-lib:fref yi (2) ((1 n))) s2i)
 630      label280
 631       (setf ascle (f2cl-lib:fref bry (1) ((1 3))))
 632       (multiple-value-bind
 633             (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10
 634              var-11)
 635           (zkscl zdr zdi fnu n yr yi nz rzr rzi ascle tol elim)
 636         (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-7 var-8 var-9
 637                          var-10 var-11))
 638         (setf nz var-6))
 639       (setf inu (f2cl-lib:int-sub n nz))
 640       (if (<= inu 0) (go end_label))
 641       (setf kk (f2cl-lib:int-add nz 1))
 642       (setf s1r (f2cl-lib:fref yr (kk) ((1 n))))
 643       (setf s1i (f2cl-lib:fref yi (kk) ((1 n))))
 644       (setf (f2cl-lib:fref yr (kk) ((1 n)))
 645               (* s1r (f2cl-lib:fref csrr (1) ((1 3)))))
 646       (setf (f2cl-lib:fref yi (kk) ((1 n)))
 647               (* s1i (f2cl-lib:fref csrr (1) ((1 3)))))
 648       (if (= inu 1) (go end_label))
 649       (setf kk (f2cl-lib:int-add nz 2))
 650       (setf s2r (f2cl-lib:fref yr (kk) ((1 n))))
 651       (setf s2i (f2cl-lib:fref yi (kk) ((1 n))))
 652       (setf (f2cl-lib:fref yr (kk) ((1 n)))
 653               (* s2r (f2cl-lib:fref csrr (1) ((1 3)))))
 654       (setf (f2cl-lib:fref yi (kk) ((1 n)))
 655               (* s2i (f2cl-lib:fref csrr (1) ((1 3)))))
 656       (if (= inu 2) (go end_label))
 657       (setf t2 (+ fnu (f2cl-lib:int-sub kk 1)))
 658       (setf ckr (* t2 rzr))
 659       (setf cki (* t2 rzi))
 660       (setf kflag 1)
 661       (go label250)
 662      label290
 663       (setf koded 2)
 664       (setf iflag 1)
 665       (setf kflag 2)
 666       (go label120)
 667      label300
 668       (setf s1r coefr)
 669       (setf s1i coefi)
 670       (setf s2r coefr)
 671       (setf s2i coefi)
 672       (go label210)
 673      label310
 674       (setf nz -2)
 675       (go end_label)
 676      end_label
 677       (return (values nil nil nil nil nil nil nil nz nil nil nil)))))
 678
 679 (in-package #:cl-user)
 680 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 681 (eval-when (:load-toplevel :compile-toplevel :execute)
 682   (setf (gethash 'fortran-to-lisp::zbknu fortran-to-lisp::*f2cl-function-info*)
 683           (fortran-to-lisp::make-f2cl-finfo
 684            :arg-types '((double-float) (double-float) (double-float)
 685                         (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 686                         (simple-array double-float (*))
 687                         (simple-array double-float (*))
 688                         (fortran-to-lisp::integer4) (double-float)
 689                         (double-float) (double-float))
 690            :return-values '(nil nil nil nil nil nil nil fortran-to-lisp::nz nil
 691                             nil nil)
 692            :calls '(fortran-to-lisp::zkscl fortran-to-lisp::zuchk
 693                     fortran-to-lisp::i1mach fortran-to-lisp::zsqrt$
 694                     fortran-to-lisp::zmlt fortran-to-lisp::zdiv
 695                     fortran-to-lisp::zexp fortran-to-lisp::dgamln
 696                     fortran-to-lisp::zshch fortran-to-lisp::zlog
 697                     fortran-to-lisp::d1mach fortran-to-lisp::zabs))))
 698