1

Incomplete gamma functions can be calculated in R with pgamma, or with gamma_inc_Q from library(gsl), or with gammainc from library(expint). However, all of these functions take only real input.

I need an implementation of the incomplete gamma function which will take complex input. Specifically, I have an integer for the first argument, and a complex number for the second argument (the limit in the integral).

This function is well-defined for complex inputs (see Wikipedia), and I've been calculating it in Mathematica. It doesn't seem to be built into R though, and I don't see it in any libraries.

So, can anyone suggest a shorter path to doing these calculations, than looking up an algorithm, implementing it in C, and writing an R interface?

(If I do have to implement it myself, here's the only algorithm for complex inputs that I've found: Kostlan & Gokhman 1987)

user54038
  • 336
  • 1
  • 9
  • 1
    R implementations of special functions are not as extensive as the Mathematica case, I guess if you didn't find what you want in gsl is pretty likely that the best thing you can do is what you have already suggested yourself. Don't be shy interacting with Rcpp ;) – Batato Dec 07 '17 at 19:22
  • I hadn't heard of Rcpp before, that looks really cool. I've had a lot of painful experiences with foreign function interfaces in R, anything to simplify it is welcome. – user54038 Dec 08 '17 at 19:29
  • If you are using RStudio, and in a unix based OS, it shouldn't be hard at all. Take a look at this [example](https://github.com/rcarcasses/schrodinger), it may be of help ;) – Batato Dec 08 '17 at 20:46

1 Answers1

1

Here is an implementation, assuming you want the lower incomplete gamma function. I've compared a couple of values with Wolfram and they match.

library(CharFun)

incgamma <- function(s,z){
  z^s * exp(-z) * hypergeom1F1(z, 1, s+1) / s
}

Perhaps the evaluation fails for a large s.

EDIT

Looks like CharFun has been removed from CRAN. You can use IncGamma in HypergeoMat:

> library(HypergeoMat)
> IncGamma(m=50, 2+2i, 5-6i)
[1] 0.3841221+0.3348439i

The result is the same on Wolfram.

Stéphane Laurent
  • 75,186
  • 15
  • 119
  • 225