(defun besj (x alpha n y nz) (declare (type float x)) (declare (type float alpha)) (declare (type fixnum n)) (declare (type (simple-array double-float (*)) y)) (declare (type fixnum nz)) (prog ((gamln 0.0d0) (spmpar 0.0d0) (ak 0.0d0) (akm 0.0d0) (ans 0.0d0) (ap 0.0d0) (arg 0.0d0) (coef 0.0d0) (dalpha 0.0d0) (dfn 0.0d0) (dtm 0.0d0) (earg 0.0d0) (elim 0.0d0) (etx 0.0d0) (fidal 0.0d0) (flgjy 0.0d0) (fn 0.0d0) (fnf 0.0d0) (fni 0.0d0) (fnp1 0.0d0) (fnu 0.0d0) (fnulim (make-array '(2) :element-type 'double-float)) (gln 0.0d0) (pdf 0.0d0) (pidt 0.0d0) (pp (make-array '(4) :element-type 'double-float)) (rden 0.0d0) (relb 0.0d0) (rttp 0.0d0) (rtwo 0.0d0) (rtx 0.0d0) (rzden 0.0d0) (s 0.0d0) (sa 0.0d0) (sb 0.0d0) (sxo2 0.0d0) (s1 0.0d0) (s2 0.0d0) (t_ 0.0d0) (ta 0.0d0) (tau 0.0d0) (tb 0.0d0) (temp (make-array '(3) :element-type 'double-float)) (tfn 0.0d0) (tm 0.0d0) (tol 0.0d0) (tolln 0.0d0) (trx 0.0d0) (tx 0.0d0) (t1 0.0d0) (t2 0.0d0) (wk (make-array '(7) :element-type 'double-float)) (xo2 0.0d0) (xo2l 0.0d0) (ipmpar 0) (i 0) (ialp 0) (idalp 0) (iflw 0) (in 0) (inlim 0) (is 0) (i1 0) (i2 0) (k 0) (kk 0) (km 0) (kt 0) (nn 0) (ns 0) ) (declare (type float gamln)) (declare (type float spmpar)) (declare (type float ak)) (declare (type float akm)) (declare (type float ans)) (declare (type float ap)) (declare (type float arg)) (declare (type float coef)) (declare (type float dalpha)) (declare (type float dfn)) (declare (type float dtm)) (declare (type float earg)) (declare (type float elim)) (declare (type float etx)) (declare (type float fidal)) (declare (type float flgjy)) (declare (type float fn)) (declare (type float fnf)) (declare (type float fni)) (declare (type float fnp1)) (declare (type float fnu)) (declare (type (simple-array double-float (*)) fnulim)) (declare (type float gln)) (declare (type float pdf)) (declare (type float pidt)) (declare (type (simple-array double-float (*)) pp)) (declare (type float rden)) (declare (type float relb)) (declare (type float rttp)) (declare (type float rtwo)) (declare (type float rtx)) (declare (type float rzden)) (declare (type float s)) (declare (type float sa)) (declare (type float sb)) (declare (type float sxo2)) (declare (type float s1)) (declare (type float s2)) (declare (type float t_)) (declare (type float ta)) (declare (type float tau)) (declare (type float tb)) (declare (type (simple-array double-float (*)) temp)) (declare (type float tfn)) (declare (type float tm)) (declare (type float tol)) (declare (type float tolln)) (declare (type float trx)) (declare (type float tx)) (declare (type float t1)) (declare (type float t2)) (declare (type (simple-array double-float (*)) wk)) (declare (type float xo2)) (declare (type float xo2l)) (declare (type fixnum ipmpar)) (declare (type fixnum i)) (declare (type fixnum ialp)) (declare (type fixnum idalp)) (declare (type fixnum iflw)) (declare (type fixnum in)) (declare (type fixnum inlim)) (declare (type fixnum is)) (declare (type fixnum i1)) (declare (type fixnum i2)) (declare (type fixnum k)) (declare (type fixnum kk)) (declare (type fixnum km)) (declare (type fixnum kt)) (declare (type fixnum nn)) (declare (type fixnum ns)) (setq rtwo 1.3483998) (replace pp '(8.729092 0.26569393 0.12457858 7.7013375E-4) :end 3) (setq inlim 150) (replace fnulim '(100.0 60.0) :end 1) (setf tb (ipmpar 4)) (setf ta (/ (spmpar 1) tb)) (if (= tb 2.0) (go label1)) (if (= tb 8.0) (go label2)) (if (= tb 16.0) (go label3)) (setf tb (alog tb)) (go label5) label1 (setf tb 0.69315) (go label5) label2 (setf tb 2.07944) (go label5) label3 (setf tb 2.77259) label5 (setf tol (amax1 ta 1.0E-15)) (setf i1 (ipmpar 5)) (setf i2 (ipmpar 6)) (setf elim (+ (* (float (- i2)) tb) (- 6.90776))) (setf tolln (* (float i1) tb)) (setf tolln (amin1 tolln 34.5388)) (setf nz 0) (setf kt 1) (arithmetic-if (+ n (- 1)) (go label720) (go label10) (go label20)) label10 (setf kt 2) label20 (setf nn n) (arithmetic-if x (go label730) (go label30) (go label80)) label30 (arithmetic-if alpha (go label710) (go label40) (go label50)) label40 (fset (fref y 1) 1.0) (if (= n 1) (go end_label)) (setf i1 2) (go label60) label50 (setf i1 1) label60 (fdo ((i i1 (+ i 1))) ((> i n) nil) (tagbody (fset (fref y i) 0.0))) (go end_label) label80 (if (< alpha 0.0) (go label710)) (setf ialp (int alpha)) (setf fni (float (+ (+ ialp n) (- 1)))) (setf fnf (+ alpha (- (float ialp)))) (setf dfn (+ fni fnf)) (setf fnu dfn) (setf xo2 (* x 0.5)) (setf sxo2 (* xo2 xo2)) (if (<= sxo2 (+ fnu 1.0)) (go label90)) (setf ta (amax1 20.0 fnu)) (if (> x ta) (go label120)) (if (> x 12.0) (go label110)) (setf xo2l (alog xo2)) (setf ns (+ (int (+ sxo2 (- fnu))) 1)) (go label100) label90 (setf fn fnu) (setf fnp1 (+ fn 1.0)) (setf xo2l (alog xo2)) (setf is kt) (if (<= x 0.5) (go label330)) (setf ns 0) label100 (setf fni (+ fni (float ns))) (setf dfn (+ fni fnf)) (setf fn dfn) (setf fnp1 (+ fn 1.0)) (setf is kt) (if (> (+ (+ n (- 1)) ns) 0) (setf is 3)) (go label330) label110 (setf ans (amax1 (+ 36.0 (- fnu)) 0.0)) (setf ns (int ans)) (setf fni (+ fni (float ns))) (setf dfn (+ fni fnf)) (setf fn dfn) (setf is kt) (if (> (+ (+ n (- 1)) ns) 0) (setf is 3)) (go label130) label120 (setf rtx (sqrt x)) (setf tau (* rtwo rtx)) (setf ta (+ tau (fref fnulim kt))) (if (<= fnu ta) (go label480)) (setf fn fnu) (setf is kt) label130 (setf i1 (iabs (+ 3 (- is)))) (setf i1 (max0 i1 1)) (setf flgjy 1.0) (multiple-value-setq (jairy x fn flgjy i1 tol elim dummy_var wk iflw) (asjy jairy x fn flgjy i1 tol elim (fref temp is) wk iflw) ) (if (/= iflw 0) (go label380)) (case is (1) (2 label450) (3)) label310 (fset (fref temp 1) (fref temp 3)) (setf kt 1) (setf is 2) (setf fni (+ fni (- 1.0))) (setf dfn (+ fni fnf)) (setf fn dfn) (if (= i1 2) (go label450)) (go label130) label330 (setf gln (gamln fnp1)) (setf arg (+ (* fn xo2l) (- gln))) (if (< arg (- elim)) (go label400)) (setf earg (exp arg)) label340 (setf s 1.0) (if (< x tol) (go label360)) (setf ak 3.0) (setf t2 1.0) (setf t_ 1.0) (setf s1 fn) (fdo ((k 1 (+ k 1))) ((> k 17) nil) (tagbody (setf s2 (+ t2 s1)) (setf t_ (/ (* (* -1 t_) sxo2) s2)) (setf s (+ s t_)) (if (< (abs t_) tol) (go label360)) (setf t2 (+ t2 ak)) (setf ak (+ ak 2.0)) (setf s1 (+ s1 fn)) )) label360 (fset (fref temp is) (* s earg)) (case is (1) (2 label450) (3)) (setf earg (/ (* earg fn) xo2)) (setf fni (+ fni (- 1.0))) (setf dfn (+ fni fnf)) (setf fn dfn) (setf is 2) (go label340) label380 (fset (fref y nn) 0.0) (setf nn (+ nn (- 1))) (setf fni (+ fni (- 1.0))) (setf dfn (+ fni fnf)) (setf fn dfn) (arithmetic-if (+ nn (- 1)) (go label440) (go label390) (go label130)) label390 (setf kt 2) (setf is 2) (go label130) label400 (fset (fref y nn) 0.0) (setf nn (+ nn (- 1))) (setf fnp1 fn) (setf fni (+ fni (- 1.0))) (setf dfn (+ fni fnf)) (setf fn dfn) (arithmetic-if (+ nn (- 1)) (go label440) (go label410) (go label420)) label410 (setf kt 2) (setf is 2) label420 (if (<= sxo2 fnp1) (go label430)) (go label130) label430 (setf arg (+ (+ arg (- xo2l)) (alog fnp1))) (if (< arg (- elim)) (go label400)) (go label330) label440 (setf nz (+ n (- nn))) (go end_label) label450 (setf nz (+ n (- nn))) (if (= kt 2) (go label470)) (fset (fref y nn) (fref temp 1)) (fset (fref y (+ nn (- 1))) (fref temp 2)) (if (= nn 2) (go end_label)) (setf trx (/ 2.0 x)) (setf dtm fni) (setf tm (* (+ dtm fnf) trx)) (setf k (+ nn 1)) (fdo ((i 3 (+ i 1))) ((> i nn) nil) (tagbody (setf k (+ k (- 1))) (fset (fref y (+ k (- 2))) (+ (* tm (fref y (+ k (- 1)))) (- (fref y k)))) (setf dtm (+ dtm (- 1.0))) (setf tm (* (+ dtm fnf) trx)) )) (go end_label) label470 (fset (fref y 1) (fref temp 2)) (go end_label) label480 (setf in (int (+ (+ alpha (- tau)) 2.0))) (if (<= in 0) (go label490)) (setf idalp (+ (+ ialp (- in)) (- 1))) (setf kt 1) (go label500) label490 (setf idalp ialp) (setf in 0) label500 (setf is kt) (setf fidal (float idalp)) (setf dalpha (+ fidal fnf)) (setf arg (+ (+ x (* (* -1 pidt) dalpha)) (- pdf))) (setf sa (sin arg)) (setf sb (cos arg)) (setf coef (/ rttp rtx)) (setf etx (* 8.0 x)) label510 (setf dtm (+ fidal fidal)) (setf dtm (* dtm dtm)) (setf tm 0.0) (if (and (= fidal 0.0) (< (abs fnf) tol)) (go label520)) (setf tm (* (* 4.0 fnf) (+ (+ fidal fidal) fnf))) label520 (setf trx (+ dtm (- 1.0))) (setf t2 (/ (+ trx tm) etx)) (setf s2 t2) (setf relb (* tol (abs t2))) (setf t1 etx) (setf s1 1.0) (setf fn 1.0) (setf ak 8.0) (fdo ((k 1 (+ k 1))) ((> k 13) nil) (tagbody (setf t1 (+ t1 etx)) (setf fn (+ fn ak)) (setf trx (+ dtm (- fn))) (setf ap (+ trx tm)) (setf t2 (/ (* (* -1 t2) ap) t1)) (setf s1 (+ s1 t2)) (setf t1 (+ t1 etx)) (setf ak (+ ak 8.0)) (setf fn (+ fn ak)) (setf trx (+ dtm (- fn))) (setf ap (+ trx tm)) (setf t2 (/ (* t2 ap) t1)) (setf s2 (+ s2 t2)) (if (<= (abs t2) relb) (go label540)) (setf ak (+ ak 8.0)) )) label540 (fset (fref temp is) (* coef (+ (* s1 sb) (* (* -1 s2) sa)))) (if (= is 2) (go label560)) (setf fidal (+ fidal 1.0)) (setf dalpha (+ fidal fnf)) (setf is 2) (setf tb sa) (setf sa (- sb)) (setf sb tb) (go label510) label560 (if (= kt 2) (go label470)) (setf s1 (fref temp 1)) (setf s2 (fref temp 2)) (setf tx (/ 2.0 x)) (setf tm (* dalpha tx)) (if (= in 0) (go label580)) (fdo ((i 1 (+ i 1))) ((> i in) nil) (tagbody (setf s s2) (setf s2 (+ (* tm s2) (- s1))) (setf tm (+ tm tx)) (setf s1 s) )) (if (= nn 1) (go label600)) (setf s s2) (setf s2 (+ (* tm s2) (- s1))) (setf tm (+ tm tx)) (setf s1 s) label580 (fset (fref y 1) s1) (fset (fref y 2) s2) (if (= nn 2) (go end_label)) (fdo ((i 3 (+ i 1))) ((> i nn) nil) (tagbody (fset (fref y i) (+ (* tm (fref y (+ i (- 1)))) (- (fref y (+ i (- 2)))))) (setf tm (+ tm tx)) )) (go end_label) label600 (fset (fref y 1) s2) (go end_label) (setf akm (amax1 (+ 3.0 (- fn)) 0.0)) (setf km (int akm)) (setf tfn (+ fn (float km))) (setf ta (/ (+ (+ (+ gln tfn) (- 0.9189385)) (/ (* -1 0.083333336) tfn)) (+ tfn 0.5)) ) (setf ta (+ xo2l (- ta))) (setf tb (/ (* -1 (+ 1.0 (/ (* -1 1.5) tfn))) tfn)) (setf akm (+ (/ tolln (+ (- ta) (sqrt (+ (* ta ta) (* (* -1 tolln) tb))))) 1.5) ) (setf in (+ km (int akm))) (go label660) (setf gln (+ (fref wk 3) (fref wk 2))) (if (> (fref wk 6) 30.0) (go label640)) (setf rden (+ (* (+ (* (fref pp 4) (fref wk 6)) (fref pp 3)) (fref wk 6)) 1.0) ) (setf rzden (+ (fref pp 1) (* (fref pp 2) (fref wk 6)))) (setf ta (/ rzden rden)) (if (< (fref wk 1) 0.1) (go label630)) (setf tb (/ gln (fref wk 5))) (go label650) label630 (setf tb (/ (+ 1.2599211 (* (+ 0.16798948 (* 0.08879443 (fref wk 1))) (fref wk 1))) (fref wk 7) )) (go label650) label640 (setf ta (/ (* 0.5 tolln) (fref wk 4))) (setf ta (* (* (+ (* (+ (* 0.049382716 ta) (- 0.11111111)) ta) 0.6666667) ta) (fref wk 6) )) (if (< (fref wk 1) 0.1) (go label630)) (setf tb (/ gln (fref wk 5))) label650 (setf in (int (+ (/ ta tb) 1.5))) (if (> in inlim) (go label310)) label660 (setf dtm (+ fni (float in))) (setf trx (/ 2.0 x)) (setf tm (* (+ dtm fnf) trx)) (setf ta 0.0) (setf tb tol) (setf kk 1) label670 (fdo ((i 1 (+ i 1))) ((> i in) nil) (tagbody (setf s tb) (setf tb (+ (* tm tb) (- ta))) (setf ta s) (setf dtm (+ dtm (- 1.0))) (setf tm (* (+ dtm fnf) trx)) )) (if (/= kk 1) (go label690)) (setf ta (* (/ ta tb) (fref temp 3))) (setf tb (fref temp 3)) (setf kk 2) (setf in ns) (if (/= ns 0) (go label670)) label690 (fset (fref y nn) tb) (setf nz (+ n (- nn))) (if (= nn 1) (go end_label)) (setf k (+ nn (- 1))) (fset (fref y k) (+ (* tm tb) (- ta))) (if (= nn 2) (go end_label)) (setf dtm (+ dtm (- 1.0))) (setf tm (* (+ dtm fnf) trx)) (setf km (+ k (- 1))) (fdo ((i 1 (+ i 1))) ((> i km) nil) (tagbody (fset (fref y (+ k (- 1))) (+ (* tm (fref y k)) (- (fref y (+ k 1))))) (setf dtm (+ dtm (- 1.0))) (setf tm (* (+ dtm fnf) trx)) (setf k (+ k (- 1))) )) (go end_label) label710 (setf nz (- 2)) (go end_label) label720 (setf nz (- 3)) (go end_label) label730 (setf nz (- 1)) (go end_label) end_label (return (values x alpha n y nz)) ))