RE: [Axiom-math] special functions

[yigal]

I am trying to use a continued fraction representation for the
incomplete gamma function gamma(x,y) with mixed success.  Does anyone
know why this error is obtained:

  >> Error detected within library code:
   Denominators must be greater than 0.

For this code:

num:=cons(1,[i^2 for i in 1..])
den:=cons(2,[-2*i for i in 2..])
continuedFraction(0,num,den)

(This is known as Gompertz constant and is equal to gamma(0,1))

It is true that the denominator has the negative sequence [-2*i for i in
2..] in it and that changing the - to a plus i.e. 

num:=cons(1,[i^2 for i in 1..])
den:=cons(2,[2*i for i in 2..])
continuedFraction(0,num,den)

results in 
   (26)
       1 |     1 |     4 |     9 |     16 |     25 |     36 |      |
|
     +---+ + +---+ + +---+ + +---+ + +----+ + +----+ + +----+ + ...
     | 2     | 4     | 6     | 8     | 10     | 12     | 14       
                                      Type: ContinuedFraction Integer

A perfectly good continued fraction although nothing I can actually use.
Is there a way for Axiom to use the first code.  

I created today a code that uses both positive numerators and
denominators, it works for gamma(0,x), x>0 but converges slowly, it is
also very poorly coded as I am a true newbie to Axiom.  I don't know how
to efficiently create an alternating sequence.

)clear all 
x:Float;x:=1.0
a:Integer;a:=0
n : (Integer) -> List Polynomial Float;
n(y) == [1,address@hidden Polynomial Float;
--
-- alternating generator for numerator
--
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,address@hidden Polynomial Float; 
--
-- alternating generater for denominator
--
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
-- bellow is the approximation to gamma(a,x) 
gam(i) == exp(-x)*x^a*ccf.i; 

Thank you for reading,

Yigal Weinstein