[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
RE: [Axiom-math] special functions
From: |
yigal |
Subject: |
RE: [Axiom-math] special functions |
Date: |
Sun, 29 Jan 2006 04:45:28 -0800 |
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,address@hidden 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,address@hidden 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