correction of the calculation of factorials and of the x^m

R/lauricella.R

@@ -85,17 +85,14 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
   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]
+    xfact[i+1, j] <- x[j]^i
   }
   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)
-  }
+  lnfact <- matrix(sapply(Munique, lnfactorial), nrow = length(Munique),
+                   ncol = n, byrow = FALSE, dimnames = list(Munique, 1:n))
   sumlnfact <- function(ind) {
     sum(mapply(function(i, j) lnfact[i, j], ind+1, 1:n))
   }
@@ -165,8 +162,24 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
     Msumunique <- (max(Msumunique)+1):max(Msum)
     
     # Product x^{m_1}/m_1! * ... * x^{m_n}/m_n! for m_1, ...,  m_n given by the rows of M
-    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)))
+    xfactsupp <- matrix(nrow = length(Munique), ncol = n, dimnames = list(Munique, 1:n))
+    for (i in 1:length(Munique)) for (j in 1:n) {
+      # Product pochhammer(b_1,m_1) * ...* pochhammer(b_n, m_n)
+      xfactsupp[i, j] <- x[j]^Munique[i]
+    }
+    xfact <- rbind(xfact, xfactsupp)
+    # prodxfact <- function(ind) {
+    #   prod(mapply(function(i, j) xfact[i, j], ind+1, 1:n))
+    # }
+    
+    lnfactsupp <- matrix(sapply(Munique, lnfactorial), nrow = length(Munique),
+                         ncol = n, dimnames = list(Munique, 1:n))
+    lnfact <- rbind(lnfact, lnfactsupp)
+    # sumlnfact <- function(ind) {
+    #   sum(mapply(function(i, j) lnfact[i, j], ind+1, 1:n))
+    # }
     
     # pochhammer(a, m_1+...+m_n) for m_1, ...,  m_n given by the rows of M
     # lnapoch <- sapply(Msum, function(j) lnpochhammer(a, j))
@@ -197,8 +210,8 @@ 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] + apply(M, 1, sumlnpochb) - lngpoch[Msum] - apply(M, 1, sumlnfact)
+    res2 <- sum( apply(M, 1, prodxfact) * exp(res1) )
   }
   
   result <- Re(result)