axiom-math
[Top][All Lists]

## RE: [Axiom-math] special functions

 From: yigal Subject: RE: [Axiom-math] special functions Date: Sat, 28 Jan 2006 20:15:50 -0800

```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;
--
-- 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;
--
-- 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;