[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
RE: [Axiommath] special functions
From: 
yigal 
Subject: 
RE: [Axiommath] 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;
n(y) == [1,address@hidden Polynomial Float;

 alternating generator for numerator

numg : (Integer) > List Polynomial Float;
numg(1)== [1,1a];
numg(ii>1) == append(numg(i1),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(ii>1) == append(deng(i1),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