Difference in MATLAB function evaluation

According to the properties of gamma function, for any x > 0

But when I try to use it in MATLAB, it shows a very large error

gamma(150) - 149*gamma(149)

ans =

 -1.5240e+247

Can anybody explain the reason for this? I need it for an approximation algorithm I am trying to develop.

Thanks in advance.

[Star Strider]

FWIW, I get similar results for these (in R2016b):

n = 149;
q1 = gamma(n+1) - n*gamma(n);
q2 = gamma(n+1) - factorial(n);
q3 = gamma(n+1) - prod(1:n);

q1 =
   -15.2402e+246

q2 =
  -819.3656e+243

q3 =
  -819.3656e+243

I hope this is a problem with my understanding of the gamma function, and not a problem with its implementation in MATLAB.

Answer 1

[Honglei Chen]

For your use case, 150! is a very large number so the numerical precision issue starts to kick in. I'd suggest you to use gammaln instead, for example, you can do

gammaln(150)-(log(149)+gammaln(149))

and the result is indeed 0.

HTH.

Answer 2

[Jan]

While the absolute error is huge, the relative error is tiny:

(gamma(150) - 149*gamma(149)) / gamma(150)
>> -3.9868e-014

This is the expected error caused by using the IEEE754 double format. This stores about 16 valid digits.

If an approxmation algorithm is sensitive for this difference, it is instable and of limited use.