After a couple hours of tinkering the incomplete gamma function has taken on a simpler and more useful form. Is there an easier way of making an alternating sequence, i.e. a sequence of the form: (a1,b1,a2,b2,...).
)clear all j:=20 n(i,a) == [i-1,i-a]; numg(1,a)== [1,1-a]; numg(i|i>1,a) == append(numg(i-1,a),n(i,a)); d(i,x) == [x,1]; deng(1,x)== [x,1]; deng(i|i>1,x) == append(deng(i-1,x),d(i,x)); num(a) == [numg(i,a).i for i in 1..]; den(x) == [deng(i,x).i for i in 1..]; cf(a,x) == continuedFraction(0,num(a),den(x)); ccf(a,x) == convergents cf(a,x); gam(a,x) == exp(-x)*x^a*ccf(a,x).j; -- -- a test gamma(a,x) == factorial(a-1)*exp(-x)*reduce(+,[x^i/factorial(i) for i in 0..(a-1)]); gam(7.0,20.0) gamma(7,20.0) Thanks, Yigal Weinstein On Sun, 2006-01-29 at 01:16 -0800, yigal wrote: > I made a small mistake but numerically important in my sloppy incomplete > gamma code here is a modified version that seems like its working: > > x:Float;x:=20.0 > a:Integer;a:=7 > n : (Integer) -> List Polynomial Float; > n(y) == [y-1,[EMAIL PROTECTED] Polynomial Float; > numg : (Integer) -> List Polynomial Float; > numg(1)== [1,1-a]; > numg(i|i>1) == append(numg(i-1),n(i))@List Polynomial Float; > d : (Integer) -> List Polynomial Float; > d(y) == [x,[EMAIL PROTECTED] Polynomial Float; > deng : (Integer) -> List Polynomial Float; > deng(1)== [x,1]; > deng(i|i>1) == append(deng(i-1),d(i))@List Polynomial Float; > num := [numg.i.i for i in 1..]; > den := [deng.i.i for i in 1..]; > cf := continuedFraction(0,num,den) > ccf := convergents cf; > gam(i) == exp(-x)*x^a*ccf.i; > gamma(n,x) == factorial(n-1)*exp(-x)*reduce(+,[x^i/factorial(i) for i in > 0..(n-1)]); > > where > gam.15 = 0.1836881970 165365277 > gamma(7,20.)=0.1836881970 165365277 > > thank you all for your previous help, > > Yigal Weinstein _______________________________________________ Axiom-math mailing list [email protected] http://lists.nongnu.org/mailman/listinfo/axiom-math
