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