I am currently working on a model which requires the use of an upper incomplete gamma function. For some parameter combinations the argument of this incomplete gamma can become negative. Unfortunately, the gamma function in Stan/boost only takes non-negative values as an input. A colleague uses Matlab and there the implemented incomplete gamma function takes negatives inputs(the documentation however makes no mention of how exactly that is computed).
I have come across this paper about an algorithm of computing the incomplete gamma with negative position and argument(Algorithm 969). It has been written in Fortran and uses series expansion if I understand it correctly(this slightly outside of my expertise and comfort zone).
I was wondering if it is straight forward to incorporate these Fortran function into my Stan code via its capability to link to functions in other languages or are there particular things to consider.
Otherwise, does anyone have any other suggestion of going about this problem? Amy help is appreciated :)