I need to compute a sum in R that involves the gamma function for each element. When the arguments of the gamma increase I get NaN as a result, and I suspect that the problem is numerical with the evaluation of the gamma. I already read the documentation of Rmpfr and gmp but I only found factiorial for integers. I also post the code here, maybe you have a better idea about the source of the error.
V1<-w1*(b1^a1)/gamma(a1)
VV1<-outer(V1,V1)
a1mat<-outer(a1, a1-1, FUN="+")
b1mat<-outer(b1, b1, FUN="+")
A<-sum(VV1*gamma(a1mat)/(b1mat^a1mat))
w1 is an array of positive real numbers that sums up to 1, a1 and b1 are vectors of positive values. A becames NaN when a1 (and a1mat) becomes long and with high values (~150).
Try working in log-space:
w1 <- c(0.2, 0.3, 0.2, 0.1, 0.1, 0.1)
a1 <- 3:8
b1 <- rep(4, 6)
a1mat <- outer(a1, a1 - 1, "+")
b1mat <- outer(b1, b1, "+")
# working in log-space
logV1 <- log(w1) + a1*log(b1) - lgamma(a1)
logVV1 <- outer(logV1, logV1, "+")
sum(exp(logVV1 + lgamma(a1mat) - a1mat*log(b1mat)))
#> [1] 0.4614941
# compared to original
V1 <- w1*(b1^a1)/gamma(a1)
VV1 <- outer(V1, V1)
sum(VV1*gamma(a1mat)/(b1mat^a1mat))
#> [1] 0.4614941