Change for the computation of Mi! and x^{Mi}

R/lauricella.R

@@ -80,8 +80,25 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
   Msumunique <- 0:max(Msum)
   
   # Product x^{m_1}/m_1! * ... * x^{m_n}/m_n! for m_1 = 0...k, ...,  m_n = 0...k
-  xfact <- apply(M, 1, function(Mi) prod( x^Mi ))
-  lnfact <- apply(M, 1, function(Mi) sum(sapply(Mi, lnfactorial)))
+  # xfact <- apply(M, 1, function(Mi) prod( x^Mi ))
+  # lnfact <- apply(M, 1, function(Mi) sum(sapply(Mi, lnfactorial)))
+  xfact <- matrix(nrow = max(M) + 1, ncol = n, dimnames = list(Munique, 1:n))
+  for (i in Munique) for (j in 1:n) {
+    # Product pochhammer(b_1,m_1) * ...* pochhammer(b_n, m_n)
+    xfact[i+1, j] <- x[j]^Munique[i+1]
+  }
+  prodxfact <- function(ind) {
+    prod(mapply(function(i, j) xfact[i, j], ind+1, 1:n))
+  }
+  
+  lnfact <- matrix(nrow = max(M) + 1, ncol = n, dimnames = list(Munique, 1:n))
+  for (i in Munique) {
+    # Product pochhammer(b_1,m_1) * ...* pochhammer(b_n, m_n)
+    lnfact[i+1, ] <- lnfactorial(i)
+  }
+  sumlnfact <- function(ind) {
+    sum(mapply(function(i, j) lnfact[i, j], ind+1, 1:n))
+  }
   
   # pochhammer(a, m_1+...+m_n) for m_1 = 0...k, ...,  m_p = 0...k
   # lnapoch <- sapply(Msum, function(j) lnpochhammer(a, j))
@@ -110,8 +127,9 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
   
   # res1 <- lnapoch + lnbpoch - lngpoch - lnfact
   # res2 <- sum(xfact * exp(res1))
-  res1 <- lnapoch[Msum+1] + apply(M, 1, sumlnpochb) - lngpoch[Msum+1] - lnfact
-  res2 <- sum(xfact * exp(res1))
+  res1 <- lnapoch[Msum+1] + apply(M, 1, sumlnpochb) - lngpoch[Msum+1] -
+    apply(M, 1, sumlnfact)
+  res2 <- sum( apply(M, 1, prodxfact) * exp(res1) )
   
   kstep <- 5