incomplete gamma function: evaluation and taylor series

[suvra...@gmail.com]

I am new to Sage; trying to explore open source alternatives to Mathematica. 

However, I seem to be having trouble with the incomplete gamma function. Here are two difficulties. First, in trying to evaluate the incomplete gamma function at a point where the result should be very small, I just get zero even if I increase the precision arbitrarily. In particular consider 

numerical_approx(gamma(9, 10^(-3))-gamma(9), digits=40)

the value of this number is approximate -1.1 \times 10^(-28), but I just get 0.00000 ...0 for the input above.

Second, I was trying to do a series expansion of the incomplete gamma function. (In fact this series expansion is the simplest way of seeing the correct value above.)

So, I wanted to consider

taylor((gamma(p,z)-gamma(p))/z^p, z, 0, 1)

and Sage gives me 0.

Actually, the right answer is: -1/p + z/(1 + p)

I can get this answer by putting a specific value of p above in Sage i.e.:

taylor((gamma(4,z)-gamma(4))/z^p, z, 0, 1)

but not otherwise. In fact, for symbolic p everything I tried failed:

for example,

taylor(gamma(p,z), z, 0, 1)

gives me gamma(p) and refuses to give any further terms.

Could you please help? What is the right way to series expand the incomplete gamma function, and evaluate it numerically at points where it gives a small value?

Suvrat

[Dima Pasechnik]

On 2014-05-03, suvra...@gmail.com <suvra...@gmail.com> wrote:
> I am new to Sage; trying to explore open source alternatives to
> Mathematica.
>
> However, I seem to be having trouble with the incomplete gamma function.
> Here are two difficulties. First, in trying to evaluate the incomplete
> gamma function at a point where the result should be very small, I just get
> zero even if I increase the precision arbitrarily. In particular consider
>
> numerical_approx(gamma(9, 10^(-3))-gamma(9), digits=40)
>
> the value of this number is approximate -1.1 \times 10^(-28), but I just
> get 0.00000 ...0 for the input above.

I think this is a bug in the Sage's interface to PARI's incomplete
gamma-function implementation.
The latter has a parameter 'precision' that actually allows you do
increase the precision of the computation here as much as you like.
E.g.

sage: pari(9).incgam(1/1000,precision=100).sage() - gamma(9)
-1.1101098718046915201731717980618372763e-28

(the part pari(9).incgam(1/1000,precision=100) calls the PARI's
incomplete gamma-function code, and .sage() pulls it back into Sage.)

and you can increase the precision, e.g.:
sage: pari(9).incgam(1/1000,precision=200).sage()-gamma(9)
-1.110111565517708813007363808370583494179636586117301130595952112736737061589e-28
sage: pari(9).incgam(1/1000,precision=500).sage()-gamma(9)
-1.11011156551770881300736380837058349417963658951616739908084575093982182228943211704794146886994896523236690487484369923202470676572472807219205439592102e-28

etc.

Any takes to fix this?
(meanwhile, you can use the workaround as above)

HTH,
Dmitrii

[Peter Bruin]

Hello,

> However, I seem to be having trouble with the incomplete gamma function.
> Here are two difficulties. First, in trying to evaluate the incomplete
> gamma function at a point where the result should be very small, I just get
> zero even if I increase the precision arbitrarily. In particular consider
>
> numerical_approx(gamma(9, 10^(-3))-gamma(9), digits=40)
>
> the value of this number is approximate -1.1 \times 10^(-28), but I just
> get 0.00000 ...0 for the input above.

I think this is a bug in the Sage's interface to PARI's incomplete
gamma-function implementation.
The latter has a parameter 'precision' that actually allows you do
increase the precision of the computation here as much as you like.

 

From looking at the source code, I suspect that the problem is in sage.functions.other.Function_gamma_inc, where the method _evalf_() does not use its argument "parent" (and hence does not know about the precision of that parent).

If you don't want to use PARI explicitly, you can also avoid the bug by typing

sage: C = ComplexField(400)
sage: C(9).gamma_inc(1/1000) - gamma(9)
-1.11011156551770881154239109830530085698810642670976295680548580836149189576782725366794514487711609262987622059881687164e-28

Peter

[Peter Bruin]

This is in fact a long-standing bug, reported here:

http://trac.sagemath.org/ticket/7099  (serious incomplete gamma function precision bugs)

[Peter Bruin]

Correction: it is not exactly the same bug, but the two are certainly related.

[kcrisman]

On Saturday, May 3, 2014 7:17:04 AM UTC-4, suvra...@gmail.com wrote:

I am new to Sage; trying to explore open source alternatives to Mathematica. 

However, I seem to be having trouble with the incomplete gamma function. Here are two difficulties. First, in trying to evaluate the incomplete gamma function at a point where the result should be very small, I just get zero even if I increase the precision arbitrarily. In particular consider 

numerical_approx(gamma(9, 10^(-3))-gamma(9), digits=40)

the value of this number is approximate -1.1 \times 10^(-28), but I just get 0.00000 ...0 for the input above.

I believe mpmath (in Sage) should also be able to do this fine, though I haven't tested your particular case:

http://mpmath.googlecode.com/svn/trunk/doc/build/functions/expintegrals.html#gammainc