src/numerical/slatec/dbsknu.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 ((x1 2.0)
  21       (x2 17.0)
  22       (pi$ 3.14159265358979)
  23       (rthpi 1.2533141373155)
  24       (cc
  25        (make-array 8
  26                    :element-type 'double-float
  27                    :initial-contents '(0.577215664901533 -0.0420026350340952
  28                                        -0.0421977345555443 0.007218943246663
  29                                        -2.152416741149e-4 -2.01348547807e-5
  30                                        1.133027232e-6 6.116095e-9))))
  31   (declare (type (double-float) x1 x2 pi$ rthpi)
  32            (type (simple-array double-float (8)) cc))
  33   (defun dbsknu (x fnu kode n y nz)
  34     (declare (type (simple-array double-float (*)) y)
  35              (type (f2cl-lib:integer4) nz n kode)
  36              (type (double-float) fnu x))
  37     (prog ((a (make-array 160 :element-type 'double-float))
  38            (b (make-array 160 :element-type 'double-float)) (ak 0.0) (a1 0.0)
  39            (a2 0.0) (bk 0.0) (ck 0.0) (coef 0.0) (cx 0.0) (dk 0.0) (dnu 0.0)
  40            (dnu2 0.0) (elim 0.0) (etest 0.0) (ex 0.0) (f 0.0) (fc 0.0)
  41            (fhs 0.0) (fk 0.0) (fks 0.0) (flrx 0.0) (fmu 0.0) (g1 0.0) (g2 0.0)
  42            (p 0.0) (pt 0.0) (p1 0.0) (p2 0.0) (q 0.0) (rx 0.0) (s 0.0)
  43            (smu 0.0) (sqk 0.0) (st 0.0) (s1 0.0) (s2 0.0) (tm 0.0) (tol 0.0)
  44            (t1 0.0) (t2 0.0) (i 0) (iflag 0) (inu 0) (j 0) (k 0) (kk 0)
  45            (koded 0) (nn 0))
  46       (declare (type (f2cl-lib:integer4) nn koded kk k j inu iflag i)
  47                (type (double-float) t2 t1 tol tm s2 s1 st sqk smu s rx q p2 p1
  48                                     pt p g2 g1 fmu flrx fks fk fhs fc f ex
  49                                     etest elim dnu2 dnu dk cx coef ck bk a2 a1
  50                                     ak)
  51                (type (simple-array double-float (160)) b a))
  52       (setf kk (f2cl-lib:int-sub (f2cl-lib:i1mach 15)))
  53       (setf elim (* 2.303 (- (* kk (f2cl-lib:d1mach 5)) 3.0)))
  54       (setf ak (f2cl-lib:d1mach 3))
  55       (setf tol (max ak 1.0e-15))
  56       (if (<= x 0.0) (go label350))
  57       (if (< fnu 0.0) (go label360))
  58       (if (or (< kode 1) (> kode 2)) (go label370))
  59       (if (< n 1) (go label380))
  60       (setf nz 0)
  61       (setf iflag 0)
  62       (setf koded kode)
  63       (setf rx (/ 2.0 x))
  64       (setf inu (f2cl-lib:int (+ fnu 0.5)))
  65       (setf dnu (- fnu inu))
  66       (if (= (abs dnu) 0.5) (go label120))
  67       (setf dnu2 0.0)
  68       (if (< (abs dnu) tol) (go label10))
  69       (setf dnu2 (* dnu dnu))
  70      label10
  71       (if (> x x1) (go label120))
  72       (setf a1 (- 1.0 dnu))
  73       (setf a2 (+ 1.0 dnu))
  74       (setf t1 (/ 1.0 (dgamma a1)))
  75       (setf t2 (/ 1.0 (dgamma a2)))
  76       (if (> (abs dnu) 0.1) (go label40))
  77       (setf s (f2cl-lib:fref cc (1) ((1 8))))
  78       (setf ak 1.0)
  79       (f2cl-lib:fdo (k 2 (f2cl-lib:int-add k 1))
  80                     ((> k 8) nil)
  81         (tagbody
  82           (setf ak (* ak dnu2))
  83           (setf tm (* (f2cl-lib:fref cc (k) ((1 8))) ak))
  84           (setf s (+ s tm))
  85           (if (< (abs tm) tol) (go label30))
  86          label20))
  87      label30
  88       (setf g1 (- s))
  89       (go label50)
  90      label40
  91       (setf g1 (/ (- t1 t2) (+ dnu dnu)))
  92      label50
  93       (setf g2 (* (+ t1 t2) 0.5))
  94       (setf smu 1.0)
  95       (setf fc 1.0)
  96       (setf flrx (f2cl-lib:flog rx))
  97       (setf fmu (* dnu flrx))
  98       (if (= dnu 0.0) (go label60))
  99       (setf fc (* dnu pi$))
 100       (setf fc (/ fc (sin fc)))
 101       (if (/= fmu 0.0) (setf smu (/ (sinh fmu) fmu)))
 102      label60
 103       (setf f (* fc (+ (* g1 (cosh fmu)) (* g2 flrx smu))))
 104       (setf fc (exp fmu))
 105       (setf p (/ (* 0.5 fc) t2))
 106       (setf q (/ 0.5 (* fc t1)))
 107       (setf ak 1.0)
 108       (setf ck 1.0)
 109       (setf bk 1.0)
 110       (setf s1 f)
 111       (setf s2 p)
 112       (if (or (> inu 0) (> n 1)) (go label90))
 113       (if (< x tol) (go label80))
 114       (setf cx (* x x 0.25))
 115      label70
 116       (setf f (/ (+ (* ak f) p q) (- bk dnu2)))
 117       (setf p (/ p (- ak dnu)))
 118       (setf q (/ q (+ ak dnu)))
 119       (setf ck (/ (* ck cx) ak))
 120       (setf t1 (* ck f))
 121       (setf s1 (+ s1 t1))
 122       (setf bk (+ bk ak ak 1.0))
 123       (setf ak (+ ak 1.0))
 124       (setf s (/ (abs t1) (+ 1.0 (abs s1))))
 125       (if (> s tol) (go label70))
 126      label80
 127       (setf (f2cl-lib:fref y (1) ((1 *))) s1)
 128       (if (= koded 1) (go end_label))
 129       (setf (f2cl-lib:fref y (1) ((1 *))) (* s1 (exp x)))
 130       (go end_label)
 131      label90
 132       (if (< x tol) (go label110))
 133       (setf cx (* x x 0.25))
 134      label100
 135       (setf f (/ (+ (* ak f) p q) (- bk dnu2)))
 136       (setf p (/ p (- ak dnu)))
 137       (setf q (/ q (+ ak dnu)))
 138       (setf ck (/ (* ck cx) ak))
 139       (setf t1 (* ck f))
 140       (setf s1 (+ s1 t1))
 141       (setf t2 (* ck (- p (* ak f))))
 142       (setf s2 (+ s2 t2))
 143       (setf bk (+ bk ak ak 1.0))
 144       (setf ak (+ ak 1.0))
 145       (setf s (+ (/ (abs t1) (+ 1.0 (abs s1))) (/ (abs t2) (+ 1.0 (abs s2)))))
 146       (if (> s tol) (go label100))
 147      label110
 148       (setf s2 (* s2 rx))
 149       (if (= koded 1) (go label170))
 150       (setf f (exp x))
 151       (setf s1 (* s1 f))
 152       (setf s2 (* s2 f))
 153       (go label170)
 154      label120
 155       (setf coef (/ rthpi (f2cl-lib:fsqrt x)))
 156       (if (= koded 2) (go label130))
 157       (if (> x elim) (go label330))
 158       (setf coef (* coef (exp (- x))))
 159      label130
 160       (if (= (abs dnu) 0.5) (go label340))
 161       (if (> x x2) (go label280))
 162       (setf etest (/ (cos (* pi$ dnu)) (* pi$ x tol)))
 163       (setf fks 1.0)
 164       (setf fhs 0.25)
 165       (setf fk 0.0)
 166       (setf ck (+ x x 2.0))
 167       (setf p1 0.0)
 168       (setf p2 1.0)
 169       (setf k 0)
 170      label140
 171       (setf k (f2cl-lib:int-add k 1))
 172       (setf fk (+ fk 1.0))
 173       (setf ak (/ (- fhs dnu2) (+ fks fk)))
 174       (setf bk (/ ck (+ fk 1.0)))
 175       (setf pt p2)
 176       (setf p2 (- (* bk p2) (* ak p1)))
 177       (setf p1 pt)
 178       (setf (f2cl-lib:fref a (k) ((1 160))) ak)
 179       (setf (f2cl-lib:fref b (k) ((1 160))) bk)
 180       (setf ck (+ ck 2.0))
 181       (setf fks (+ fks fk fk 1.0))
 182       (setf fhs (+ fhs fk fk))
 183       (if (> etest (* fk p1)) (go label140))
 184       (setf kk k)
 185       (setf s 1.0)
 186       (setf p1 0.0)
 187       (setf p2 1.0)
 188       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 189                     ((> i k) nil)
 190         (tagbody
 191           (setf pt p2)
 192           (setf p2
 193                   (/ (- (* (f2cl-lib:fref b (kk) ((1 160))) p2) p1)
 194                      (f2cl-lib:fref a (kk) ((1 160)))))
 195           (setf p1 pt)
 196           (setf s (+ s p2))
 197           (setf kk (f2cl-lib:int-sub kk 1))
 198          label150))
 199       (setf s1 (* coef (/ p2 s)))
 200       (if (or (> inu 0) (> n 1)) (go label160))
 201       (go label200)
 202      label160
 203       (setf s2 (/ (* s1 (+ x dnu 0.5 (/ (- p1) p2))) x))
 204      label170
 205       (setf ck (/ (+ dnu dnu 2.0) x))
 206       (if (= n 1) (setf inu (f2cl-lib:int-sub inu 1)))
 207       (if (> inu 0) (go label180))
 208       (if (> n 1) (go label200))
 209       (setf s1 s2)
 210       (go label200)
 211      label180
 212       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 213                     ((> i inu) nil)
 214         (tagbody
 215           (setf st s2)
 216           (setf s2 (+ (* ck s2) s1))
 217           (setf s1 st)
 218           (setf ck (+ ck rx))
 219          label190))
 220       (if (= n 1) (setf s1 s2))
 221      label200
 222       (if (= iflag 1) (go label220))
 223       (setf (f2cl-lib:fref y (1) ((1 *))) s1)
 224       (if (= n 1) (go end_label))
 225       (setf (f2cl-lib:fref y (2) ((1 *))) s2)
 226       (if (= n 2) (go end_label))
 227       (f2cl-lib:fdo (i 3 (f2cl-lib:int-add i 1))
 228                     ((> i n) nil)
 229         (tagbody
 230           (setf (f2cl-lib:fref y (i) ((1 *)))
 231                   (+ (* ck (f2cl-lib:fref y ((f2cl-lib:int-sub i 1)) ((1 *))))
 232                      (f2cl-lib:fref y ((f2cl-lib:int-sub i 2)) ((1 *)))))
 233           (setf ck (+ ck rx))
 234          label210))
 235       (go end_label)
 236      label220
 237       (setf s (- (f2cl-lib:flog s1) x))
 238       (setf (f2cl-lib:fref y (1) ((1 *))) 0.0)
 239       (setf nz 1)
 240       (if (< s (- elim)) (go label230))
 241       (setf (f2cl-lib:fref y (1) ((1 *))) (exp s))
 242       (setf nz 0)
 243      label230
 244       (if (= n 1) (go end_label))
 245       (setf s (- (f2cl-lib:flog s2) x))
 246       (setf (f2cl-lib:fref y (2) ((1 *))) 0.0)
 247       (setf nz (f2cl-lib:int-add nz 1))
 248       (if (< s (- elim)) (go label240))
 249       (setf nz (f2cl-lib:int-sub nz 1))
 250       (setf (f2cl-lib:fref y (2) ((1 *))) (exp s))
 251      label240
 252       (if (= n 2) (go end_label))
 253       (setf kk 2)
 254       (if (< nz 2) (go label260))
 255       (f2cl-lib:fdo (i 3 (f2cl-lib:int-add i 1))
 256                     ((> i n) nil)
 257         (tagbody
 258           (setf kk i)
 259           (setf st s2)
 260           (setf s2 (+ (* ck s2) s1))
 261           (setf s1 st)
 262           (setf ck (+ ck rx))
 263           (setf s (- (f2cl-lib:flog s2) x))
 264           (setf nz (f2cl-lib:int-add nz 1))
 265           (setf (f2cl-lib:fref y (i) ((1 *))) 0.0)
 266           (if (< s (- elim)) (go label250))
 267           (setf (f2cl-lib:fref y (i) ((1 *))) (exp s))
 268           (setf nz (f2cl-lib:int-sub nz 1))
 269           (go label260)
 270          label250))
 271       (go end_label)
 272      label260
 273       (if (= kk n) (go end_label))
 274       (setf s2 (+ (* s2 ck) s1))
 275       (setf ck (+ ck rx))
 276       (setf kk (f2cl-lib:int-add kk 1))
 277       (setf (f2cl-lib:fref y (kk) ((1 *))) (exp (- (f2cl-lib:flog s2) x)))
 278       (if (= kk n) (go end_label))
 279       (setf kk (f2cl-lib:int-add kk 1))
 280       (f2cl-lib:fdo (i kk (f2cl-lib:int-add i 1))
 281                     ((> i n) nil)
 282         (tagbody
 283           (setf (f2cl-lib:fref y (i) ((1 *)))
 284                   (+ (* ck (f2cl-lib:fref y ((f2cl-lib:int-sub i 1)) ((1 *))))
 285                      (f2cl-lib:fref y ((f2cl-lib:int-sub i 2)) ((1 *)))))
 286           (setf ck (+ ck rx))
 287          label270))
 288       (go end_label)
 289      label280
 290       (setf nn 2)
 291       (if (and (= inu 0) (= n 1)) (setf nn 1))
 292       (setf dnu2 (+ dnu dnu))
 293       (setf fmu 0.0)
 294       (if (< (abs dnu2) tol) (go label290))
 295       (setf fmu (* dnu2 dnu2))
 296      label290
 297       (setf ex (* x 8.0))
 298       (setf s2 0.0)
 299       (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
 300                     ((> k nn) nil)
 301         (tagbody
 302           (setf s1 s2)
 303           (setf s 1.0)
 304           (setf ak 0.0)
 305           (setf ck 1.0)
 306           (setf sqk 1.0)
 307           (setf dk ex)
 308           (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 309                         ((> j 30) nil)
 310             (tagbody
 311               (setf ck (/ (* ck (- fmu sqk)) dk))
 312               (setf s (+ s ck))
 313               (setf dk (+ dk ex))
 314               (setf ak (+ ak 8.0))
 315               (setf sqk (+ sqk ak))
 316               (if (< (abs ck) tol) (go label310))
 317              label300))
 318          label310
 319           (setf s2 (* s coef))
 320           (setf fmu (+ fmu (* 8.0 dnu) 4.0))
 321          label320))
 322       (if (> nn 1) (go label170))
 323       (setf s1 s2)
 324       (go label200)
 325      label330
 326       (setf koded 2)
 327       (setf iflag 1)
 328       (go label120)
 329      label340
 330       (setf s1 coef)
 331       (setf s2 coef)
 332       (go label170)
 333      label350
 334       (xermsg "SLATEC" "DBSKNU" "X NOT GREATER THAN ZERO" 2 1)
 335       (go end_label)
 336      label360
 337       (xermsg "SLATEC" "DBSKNU" "FNU NOT ZERO OR POSITIVE" 2 1)
 338       (go end_label)
 339      label370
 340       (xermsg "SLATEC" "DBSKNU" "KODE NOT 1 OR 2" 2 1)
 341       (go end_label)
 342      label380
 343       (xermsg "SLATEC" "DBSKNU" "N NOT GREATER THAN 0" 2 1)
 344       (go end_label)
 345      end_label
 346       (return (values nil nil nil nil nil nz)))))
 347
 348 (in-package #:cl-user)
 349 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 350 (eval-when (:load-toplevel :compile-toplevel :execute)
 351   (setf (gethash 'fortran-to-lisp::dbsknu
 352                  fortran-to-lisp::*f2cl-function-info*)
 353           (fortran-to-lisp::make-f2cl-finfo
 354            :arg-types '((double-float) (double-float)
 355                         (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 356                         (simple-array double-float (*))
 357                         (fortran-to-lisp::integer4))
 358            :return-values '(nil nil nil nil nil fortran-to-lisp::nz)
 359            :calls '(fortran-to-lisp::xermsg fortran-to-lisp::dgamma
 360                     fortran-to-lisp::d1mach fortran-to-lisp::i1mach))))
 361