src/numerical/slatec/dbesi.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 ':array)
  14 ;;;           (:array-slicing t) (:declare-common nil)
  15 ;;;           (:float-format double-float))
  16
  17 (in-package :slatec)
  18
  19
  20 (let ((rttpi 0.398942280401433) (inlim 80))
  21   (declare (type (double-float) rttpi) (type (f2cl-lib:integer4) inlim))
  22   (defun dbesi (x alpha kode n y nz)
  23     (declare (type (array double-float (*)) y)
  24              (type (f2cl-lib:integer4) nz n kode)
  25              (type (double-float) alpha x))
  26     (f2cl-lib:with-multi-array-data
  27         ((y double-float y-%data% y-%offset%))
  28       (prog ((temp (make-array 3 :element-type 'double-float)) (ain 0.0)
  29              (ak 0.0) (akm 0.0) (ans 0.0) (ap 0.0) (arg 0.0) (atol 0.0)
  30              (tolln 0.0) (dfn 0.0) (dtm 0.0) (dx 0.0) (earg 0.0) (elim 0.0)
  31              (etx 0.0) (flgik 0.0) (fn 0.0) (fnf 0.0) (fni 0.0) (fnp1 0.0)
  32              (fnu 0.0) (gln 0.0) (ra 0.0) (s 0.0) (sx 0.0) (sxo2 0.0) (s1 0.0)
  33              (s2 0.0) (t$ 0.0) (ta 0.0) (tb 0.0) (tfn 0.0) (tm 0.0) (tol 0.0)
  34              (trx 0.0) (t2 0.0) (xo2 0.0) (xo2l 0.0) (z 0.0) (i 0) (ialp 0)
  35              (in 0) (is 0) (i1 0) (k 0) (kk 0) (km 0) (kt 0) (nn 0) (ns 0))
  36         (declare (type (f2cl-lib:integer4) ns nn kt km kk k i1 is in ialp i)
  37                  (type (array double-float (3)) temp)
  38                  (type (double-float) z xo2l xo2 t2 trx tol tm tfn tb ta t$ s2
  39                                       s1 sxo2 sx s ra gln fnu fnp1 fni fnf fn
  40                                       flgik etx elim earg dx dtm dfn tolln atol
  41                                       arg ap ans akm ak ain))
  42         (setf nz 0)
  43         (setf kt 1)
  44         (setf ra (f2cl-lib:d1mach 3))
  45         (setf tol (max ra 1.0e-15))
  46         (setf i1 (f2cl-lib:int-sub (f2cl-lib:i1mach 15)))
  47         (setf gln (f2cl-lib:d1mach 5))
  48         (setf elim (* 2.303 (- (* i1 gln) 3.0)))
  49         (setf i1 (f2cl-lib:int-add (f2cl-lib:i1mach 14) 1))
  50         (setf tolln (* 2.303 gln i1))
  51         (setf tolln (min tolln 34.5388))
  52         (f2cl-lib:arithmetic-if (f2cl-lib:int-sub n 1)
  53                                 (go label590)
  54                                 (go label10)
  55                                 (go label20))
  56        label10
  57         (setf kt 2)
  58        label20
  59         (setf nn n)
  60         (if (or (< kode 1) (> kode 2)) (go label570))
  61         (f2cl-lib:arithmetic-if x (go label600) (go label30) (go label80))
  62        label30
  63         (f2cl-lib:arithmetic-if alpha (go label580) (go label40) (go label50))
  64        label40
  65         (setf (f2cl-lib:fref y-%data% (1) ((1 *)) y-%offset%) 1.0)
  66         (if (= n 1) (go end_label))
  67         (setf i1 2)
  68         (go label60)
  69        label50
  70         (setf i1 1)
  71        label60
  72         (f2cl-lib:fdo (i i1 (f2cl-lib:int-add i 1))
  73                       ((> i n) nil)
  74           (tagbody
  75             (setf (f2cl-lib:fref y-%data% (i) ((1 *)) y-%offset%) 0.0)
  76            label70))
  77         (go end_label)
  78        label80
  79         (if (< alpha 0.0) (go label580))
  80         (setf ialp (f2cl-lib:int alpha))
  81         (setf fni
  82                 (coerce
  83                  (the f2cl-lib:integer4
  84                       (f2cl-lib:int-sub (f2cl-lib:int-add ialp n) 1))
  85                  'double-float))
  86         (setf fnf (- alpha ialp))
  87         (setf dfn (+ fni fnf))
  88         (setf fnu dfn)
  89         (setf in 0)
  90         (setf xo2 (* x 0.5))
  91         (setf sxo2 (* xo2 xo2))
  92         (setf etx
  93                 (coerce (the f2cl-lib:integer4 (f2cl-lib:int-sub kode 1))
  94                         'double-float))
  95         (setf sx (* etx x))
  96         (if (<= sxo2 (+ fnu 1.0)) (go label90))
  97         (if (<= x 12.0) (go label110))
  98         (setf fn (* 0.55 fnu fnu))
  99         (setf fn (max 17.0 fn))
 100         (if (>= x fn) (go label430))
 101         (setf ans (max (- 36.0 fnu) 0.0))
 102         (setf ns (f2cl-lib:int ans))
 103         (setf fni (+ fni ns))
 104         (setf dfn (+ fni fnf))
 105         (setf fn dfn)
 106         (setf is kt)
 107         (setf km (f2cl-lib:int-add (f2cl-lib:int-sub n 1) ns))
 108         (if (> km 0) (setf is 3))
 109         (go label120)
 110        label90
 111         (setf fn fnu)
 112         (setf fnp1 (+ fn 1.0))
 113         (setf xo2l (f2cl-lib:flog xo2))
 114         (setf is kt)
 115         (if (<= x 0.5) (go label230))
 116         (setf ns 0)
 117        label100
 118         (setf fni (+ fni ns))
 119         (setf dfn (+ fni fnf))
 120         (setf fn dfn)
 121         (setf fnp1 (+ fn 1.0))
 122         (setf is kt)
 123         (if (> (f2cl-lib:int-add (f2cl-lib:int-sub n 1) ns) 0) (setf is 3))
 124         (go label230)
 125        label110
 126         (setf xo2l (f2cl-lib:flog xo2))
 127         (setf ns (f2cl-lib:int (- sxo2 fnu)))
 128         (go label100)
 129        label120
 130         (if (= kode 2) (go label130))
 131         (if (< alpha 1.0) (go label150))
 132         (setf z (/ x alpha))
 133         (setf ra (f2cl-lib:fsqrt (+ 1.0 (* z z))))
 134         (setf gln (f2cl-lib:flog (/ (+ 1.0 ra) z)))
 135         (setf t$ (+ (* ra (- 1.0 etx)) (/ etx (+ z ra))))
 136         (setf arg (* alpha (- t$ gln)))
 137         (if (> arg elim) (go label610))
 138         (if (= km 0) (go label140))
 139        label130
 140         (setf z (/ x fn))
 141         (setf ra (f2cl-lib:fsqrt (+ 1.0 (* z z))))
 142         (setf gln (f2cl-lib:flog (/ (+ 1.0 ra) z)))
 143         (setf t$ (+ (* ra (- 1.0 etx)) (/ etx (+ z ra))))
 144         (setf arg (* fn (- t$ gln)))
 145        label140
 146         (if (< arg (- elim)) (go label280))
 147         (go label190)
 148        label150
 149         (if (> x elim) (go label610))
 150         (go label130)
 151        label160
 152         (if (/= km 0) (go label170))
 153         (setf (f2cl-lib:fref y-%data% (1) ((1 *)) y-%offset%)
 154                 (f2cl-lib:fref temp (3) ((1 3))))
 155         (go end_label)
 156        label170
 157         (setf (f2cl-lib:fref temp (1) ((1 3)))
 158                 (f2cl-lib:fref temp (3) ((1 3))))
 159         (setf in ns)
 160         (setf kt 1)
 161         (setf i1 0)
 162        label180
 163         (setf is 2)
 164         (setf fni (- fni 1.0))
 165         (setf dfn (+ fni fnf))
 166         (setf fn dfn)
 167         (if (= i1 2) (go label350))
 168         (setf z (/ x fn))
 169         (setf ra (f2cl-lib:fsqrt (+ 1.0 (* z z))))
 170         (setf gln (f2cl-lib:flog (/ (+ 1.0 ra) z)))
 171         (setf t$ (+ (* ra (- 1.0 etx)) (/ etx (+ z ra))))
 172         (setf arg (* fn (- t$ gln)))
 173        label190
 174         (setf i1 (abs (f2cl-lib:int-sub 3 is)))
 175         (setf i1 (max (the f2cl-lib:integer4 i1) (the f2cl-lib:integer4 1)))
 176         (setf flgik 1.0)
 177         (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7)
 178             (dasyik x fn kode flgik ra arg i1
 179              (f2cl-lib:array-slice temp double-float (is) ((1 3))))
 180           (declare (ignore var-0 var-1 var-2 var-3 var-6 var-7))
 181           (setf ra var-4)
 182           (setf arg var-5))
 183         (f2cl-lib:computed-goto (label180 label350 label510) is)
 184        label230
 185         (setf gln (dlngam fnp1))
 186         (setf arg (- (* fn xo2l) gln sx))
 187         (if (< arg (- elim)) (go label300))
 188         (setf earg (exp arg))
 189        label240
 190         (setf s 1.0)
 191         (if (< x tol) (go label260))
 192         (setf ak 3.0)
 193         (setf t2 1.0)
 194         (setf t$ 1.0)
 195         (setf s1 fn)
 196         (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
 197                       ((> k 17) nil)
 198           (tagbody
 199             (setf s2 (+ t2 s1))
 200             (setf t$ (/ (* t$ sxo2) s2))
 201             (setf s (+ s t$))
 202             (if (< (abs t$) tol) (go label260))
 203             (setf t2 (+ t2 ak))
 204             (setf ak (+ ak 2.0))
 205             (setf s1 (+ s1 fn))
 206            label250))
 207        label260
 208         (setf (f2cl-lib:fref temp (is) ((1 3))) (* s earg))
 209         (f2cl-lib:computed-goto (label270 label350 label500) is)
 210        label270
 211         (setf earg (/ (* earg fn) xo2))
 212         (setf fni (- fni 1.0))
 213         (setf dfn (+ fni fnf))
 214         (setf fn dfn)
 215         (setf is 2)
 216         (go label240)
 217        label280
 218         (setf (f2cl-lib:fref y-%data% (nn) ((1 *)) y-%offset%) 0.0)
 219         (setf nn (f2cl-lib:int-sub nn 1))
 220         (setf fni (- fni 1.0))
 221         (setf dfn (+ fni fnf))
 222         (setf fn dfn)
 223         (f2cl-lib:arithmetic-if (f2cl-lib:int-sub nn 1)
 224                                 (go label340)
 225                                 (go label290)
 226                                 (go label130))
 227        label290
 228         (setf kt 2)
 229         (setf is 2)
 230         (go label130)
 231        label300
 232         (setf (f2cl-lib:fref y-%data% (nn) ((1 *)) y-%offset%) 0.0)
 233         (setf nn (f2cl-lib:int-sub nn 1))
 234         (setf fnp1 fn)
 235         (setf fni (- fni 1.0))
 236         (setf dfn (+ fni fnf))
 237         (setf fn dfn)
 238         (f2cl-lib:arithmetic-if (f2cl-lib:int-sub nn 1)
 239                                 (go label340)
 240                                 (go label310)
 241                                 (go label320))
 242        label310
 243         (setf kt 2)
 244         (setf is 2)
 245        label320
 246         (if (<= sxo2 fnp1) (go label330))
 247         (go label130)
 248        label330
 249         (setf arg (+ (- arg xo2l) (f2cl-lib:flog fnp1)))
 250         (if (< arg (- elim)) (go label300))
 251         (go label230)
 252        label340
 253         (setf nz (f2cl-lib:int-sub n nn))
 254         (go end_label)
 255        label350
 256         (setf nz (f2cl-lib:int-sub n nn))
 257        label360
 258         (if (= kt 2) (go label420))
 259         (setf s1 (f2cl-lib:fref temp (1) ((1 3))))
 260         (setf s2 (f2cl-lib:fref temp (2) ((1 3))))
 261         (setf trx (/ 2.0 x))
 262         (setf dtm fni)
 263         (setf tm (* (+ dtm fnf) trx))
 264         (if (= in 0) (go label390))
 265         (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 266                       ((> i in) nil)
 267           (tagbody
 268             (setf s s2)
 269             (setf s2 (+ (* tm s2) s1))
 270             (setf s1 s)
 271             (setf dtm (- dtm 1.0))
 272             (setf tm (* (+ dtm fnf) trx))
 273            label380))
 274         (setf (f2cl-lib:fref y-%data% (nn) ((1 *)) y-%offset%) s1)
 275         (if (= nn 1) (go end_label))
 276         (setf (f2cl-lib:fref y-%data%
 277                              ((f2cl-lib:int-sub nn 1))
 278                              ((1 *))
 279                              y-%offset%)
 280                 s2)
 281         (if (= nn 2) (go end_label))
 282         (go label400)
 283        label390
 284         (setf (f2cl-lib:fref y-%data% (nn) ((1 *)) y-%offset%) s1)
 285         (setf (f2cl-lib:fref y-%data%
 286                              ((f2cl-lib:int-sub nn 1))
 287                              ((1 *))
 288                              y-%offset%)
 289                 s2)
 290         (if (= nn 2) (go end_label))
 291        label400
 292         (setf k (f2cl-lib:int-add nn 1))
 293         (f2cl-lib:fdo (i 3 (f2cl-lib:int-add i 1))
 294                       ((> i nn) nil)
 295           (tagbody
 296             (setf k (f2cl-lib:int-sub k 1))
 297             (setf (f2cl-lib:fref y-%data%
 298                                  ((f2cl-lib:int-sub k 2))
 299                                  ((1 *))
 300                                  y-%offset%)
 301                     (+
 302                      (* tm
 303                         (f2cl-lib:fref y-%data%
 304                                        ((f2cl-lib:int-sub k 1))
 305                                        ((1 *))
 306                                        y-%offset%))
 307                      (f2cl-lib:fref y-%data% (k) ((1 *)) y-%offset%)))
 308             (setf dtm (- dtm 1.0))
 309             (setf tm (* (+ dtm fnf) trx))
 310            label410))
 311         (go end_label)
 312        label420
 313         (setf (f2cl-lib:fref y-%data% (1) ((1 *)) y-%offset%)
 314                 (f2cl-lib:fref temp (2) ((1 3))))
 315         (go end_label)
 316        label430
 317         (setf earg (/ rttpi (f2cl-lib:fsqrt x)))
 318         (if (= kode 2) (go label440))
 319         (if (> x elim) (go label610))
 320         (setf earg (* earg (exp x)))
 321        label440
 322         (setf etx (* 8.0 x))
 323         (setf is kt)
 324         (setf in 0)
 325         (setf fn fnu)
 326        label450
 327         (setf dx (+ fni fni))
 328         (setf tm 0.0)
 329         (if (and (= fni 0.0) (< (abs fnf) tol)) (go label460))
 330         (setf tm (* 4.0 fnf (+ fni fni fnf)))
 331        label460
 332         (setf dtm (* dx dx))
 333         (setf s1 etx)
 334         (setf trx (- dtm 1.0))
 335         (setf dx (/ (- (+ trx tm)) etx))
 336         (setf t$ dx)
 337         (setf s (+ 1.0 dx))
 338         (setf atol (* tol (abs s)))
 339         (setf s2 1.0)
 340         (setf ak 8.0)
 341         (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
 342                       ((> k 25) nil)
 343           (tagbody
 344             (setf s1 (+ s1 etx))
 345             (setf s2 (+ s2 ak))
 346             (setf dx (- dtm s2))
 347             (setf ap (+ dx tm))
 348             (setf t$ (/ (* (- t$) ap) s1))
 349             (setf s (+ s t$))
 350             (if (<= (abs t$) atol) (go label480))
 351             (setf ak (+ ak 8.0))
 352            label470))
 353        label480
 354         (setf (f2cl-lib:fref temp (is) ((1 3))) (* s earg))
 355         (if (= is 2) (go label360))
 356         (setf is 2)
 357         (setf fni (- fni 1.0))
 358         (setf dfn (+ fni fnf))
 359         (setf fn dfn)
 360         (go label450)
 361        label500
 362         (setf akm (max (- 3.0 fn) 0.0))
 363         (setf km (f2cl-lib:int akm))
 364         (setf tfn (+ fn km))
 365         (setf ta
 366                 (/ (+ (- (+ gln tfn) 0.9189385332) (/ -0.0833333333 tfn))
 367                    (+ tfn 0.5)))
 368         (setf ta (- xo2l ta))
 369         (setf tb (/ (- (+ 1.0 (/ (* -1 1.0) tfn))) tfn))
 370         (setf ain
 371                 (+ (/ tolln (- (f2cl-lib:fsqrt (- (* ta ta) (* tolln tb))) ta))
 372                    1.5))
 373         (setf in (f2cl-lib:int ain))
 374         (setf in (f2cl-lib:int-add in km))
 375         (go label520)
 376        label510
 377         (setf t$ (/ 1.0 (* fn ra)))
 378         (setf ain
 379                 (+
 380                  (/ tolln
 381                     (+ gln (f2cl-lib:fsqrt (+ (* gln gln) (* t$ tolln)))))
 382                  1.5))
 383         (setf in (f2cl-lib:int ain))
 384         (if (> in inlim) (go label160))
 385        label520
 386         (setf trx (/ 2.0 x))
 387         (setf dtm (+ fni in))
 388         (setf tm (* (+ dtm fnf) trx))
 389         (setf ta 0.0)
 390         (setf tb tol)
 391         (setf kk 1)
 392        label530
 393         (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 394                       ((> i in) nil)
 395           (tagbody
 396             (setf s tb)
 397             (setf tb (+ (* tm tb) ta))
 398             (setf ta s)
 399             (setf dtm (- dtm 1.0))
 400             (setf tm (* (+ dtm fnf) trx))
 401            label540))
 402         (if (/= kk 1) (go label550))
 403         (setf ta (* (/ ta tb) (f2cl-lib:fref temp (3) ((1 3)))))
 404         (setf tb (f2cl-lib:fref temp (3) ((1 3))))
 405         (setf kk 2)
 406         (setf in ns)
 407         (if (/= ns 0) (go label530))
 408        label550
 409         (setf (f2cl-lib:fref y-%data% (nn) ((1 *)) y-%offset%) tb)
 410         (setf nz (f2cl-lib:int-sub n nn))
 411         (if (= nn 1) (go end_label))
 412         (setf tb (+ (* tm tb) ta))
 413         (setf k (f2cl-lib:int-sub nn 1))
 414         (setf (f2cl-lib:fref y-%data% (k) ((1 *)) y-%offset%) tb)
 415         (if (= nn 2) (go end_label))
 416         (setf dtm (- dtm 1.0))
 417         (setf tm (* (+ dtm fnf) trx))
 418         (setf km (f2cl-lib:int-sub k 1))
 419         (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 420                       ((> i km) nil)
 421           (tagbody
 422             (setf (f2cl-lib:fref y-%data%
 423                                  ((f2cl-lib:int-sub k 1))
 424                                  ((1 *))
 425                                  y-%offset%)
 426                     (+ (* tm (f2cl-lib:fref y-%data% (k) ((1 *)) y-%offset%))
 427                        (f2cl-lib:fref y-%data%
 428                                       ((f2cl-lib:int-add k 1))
 429                                       ((1 *))
 430                                       y-%offset%)))
 431             (setf dtm (- dtm 1.0))
 432             (setf tm (* (+ dtm fnf) trx))
 433             (setf k (f2cl-lib:int-sub k 1))
 434            label560))
 435         (go end_label)
 436        label570
 437         (xermsg "SLATEC" "DBESI" "SCALING OPTION, KODE, NOT 1 OR 2." 2 1)
 438         (go end_label)
 439        label580
 440         (xermsg "SLATEC" "DBESI" "ORDER, ALPHA, LESS THAN ZERO." 2 1)
 441         (go end_label)
 442        label590
 443         (xermsg "SLATEC" "DBESI" "N LESS THAN ONE." 2 1)
 444         (go end_label)
 445        label600
 446         (xermsg "SLATEC" "DBESI" "X LESS THAN ZERO." 2 1)
 447         (go end_label)
 448        label610
 449         (xermsg "SLATEC" "DBESI" "OVERFLOW, X TOO LARGE FOR KODE = 1." 6 1)
 450         (go end_label)
 451        end_label
 452         (return (values nil nil nil nil nil nz))))))
 453
 454 (in-package #:cl-user)
 455 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 456 (eval-when (:load-toplevel :compile-toplevel :execute)
 457   (setf (gethash 'fortran-to-lisp::dbesi fortran-to-lisp::*f2cl-function-info*)
 458           (fortran-to-lisp::make-f2cl-finfo
 459            :arg-types '((double-float) (double-float)
 460                         (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 461                         (array double-float (*)) (fortran-to-lisp::integer4))
 462            :return-values '(nil nil nil nil nil fortran-to-lisp::nz)
 463            :calls '(fortran-to-lisp::xermsg fortran-to-lisp::dlngam
 464                     fortran-to-lisp::dasyik fortran-to-lisp::i1mach
 465                     fortran-to-lisp::d1mach))))
 466