Computing hypergeometric function of two variables

[Quantox]

Hi there,

I'm about to implement an option pricing model based on a variance gamma process
published in a paper by Madan, Carr and Chang.
There needs to be computed the degenerative hypergeometric function of two variables,
defined by Humbert, denoted as H(a,b,c;x,y).
I couldn't find any numerical solution or straight formula for doing this.
Does anybody know how to compute it or which
Python/R-libraries to use?
Thanks in advance!