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