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