[R] precision of gamma function

[Petr Savicky]

On Wed, Feb 09, 2011 at 05:23:27PM +0200, Chuse chuse wrote:
> Dear R users,
> 
> I have to calculate gamma functions for negative numbers beyond -171.4.
> e.x. gamma(-500.4)
> I got following:
> 
> > gamma(-170.4)
> [1] -5.824625e-308
> > gamma(-171.4)
> [1] 0
> Warning message:
> underflow occurred in 'gammafn'
> 
> I have tried to use a recursion getting values a little futher -180.
> How could I solve this problem? Thank you beforehand.

Hello.

As others pointed out, the function lgamma() should be used. If
the accuracy of double precision is not sufficient or in case
of doubt concerning accuracy in extreme cases, Rmpfr package may
be used. In this case, we get

  x <- c(-171.4, -170.4)
  y1 <- lgamma(x)
  y1

  [1] -712.5781 -707.4341

  library(Rmpfr)
  y2 <- lgamma(mpfr(x, precBits=200))
  y1 - y2

  2 'mpfr' numbers of precision  200   bits 
  [1] -8.0100329934347335532420397139357516023505624051256280109870069e-14
  [2] -1.0849741205998334826524274950211713854549088941953001164911718e-13

Hope this helps.

Petr Savicky.