The lnapoch, lnbpoch and lngpoch elements are now calculated once for each...

The lnapoch, lnbpoch and lngpoch elements are now calculated once for each value in M or Msum. This improves the execution time.

R/lauricella.R

@@ -35,6 +35,28 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
   #' \doi{10.1109/LSP.2019.2915000}
   #' @export
   
+  # # Internal functions
+  # lnfactorial <- function(n) {
+  #   # Logarithm of the factorial of a positive integer
+  #   if (n == 0) {
+  #     return(0)
+  #   }
+  #   
+  #   if (n > 0) {
+  #     return(sum(log(1:n)))
+  #   }
+  # }
+  # lnpochhammer <- function(x, n) {
+  #   # Logarithm of the Pochhammer Symbol
+  #   if (n == 0) {
+  #     return(0)
+  #   }
+  #   if (n > 0) {
+  #     y <- x + (1:n) - 1
+  #     return(sum(log(abs(y)) + (y < 0)*pi*1i))
+  #   }
+  # }
+  # 
   # Number of variables
   n <- length(x)
   
@@ -54,23 +76,41 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
   
   Msum <- apply(M, 1, sum)
   
+  Munique <- 0:max(M)
+  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)))
   
   # pochhammer(a, m_1+...+m_n) for m_1 = 0...k, ...,  m_p = 0...k
-  lnapoch <- sapply(Msum, function(j) lnpochhammer(a, j))
-  
-  lnbpoch <- numeric(nrow(M))
-  for (i in 1:nrow(M)) {
+  # lnapoch <- sapply(Msum, function(j) lnpochhammer(a, j))
+  lnapoch <- sapply(Msumunique, function(j) lnpochhammer(a, j))
+  names(lnapoch) <- Msumunique
+
+  lnbpoch <- 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)
-    lnbpoch[i] <- sum( sapply(1:n, function(j) lnpochhammer(b[j], M[i, j])) )
+    lnbpoch[i+1, j] <- lnpochhammer(b[j], i)
+  }
+  # matlnbpoch <- as.data.frame(matrix(nrow = max(Munique)+1, ncol = n))
+  # for (i in Munique) for (j in 1:n) {
+  #   # Product pochhammer(b_1,m_1) * ...* pochhammer(b_n, m_n)
+  #   matlnbpoch[i+1, j] <- lnpochhammer(b[j], Munique[i+1])
+  # }
+  # lnbpoch <- data.frame(stack(matlnbpoch), M = 0:maxM)
+  sumlnpochb <- function(ind) {
+    sum(mapply(function(i, j) lnbpoch[i, j], ind+1, 1:n))
   }
   
   # pochhammer(c, m_1+...+m_n) for m_1 = 0...k, ...,  m_p = 0...k
-  lngpoch <- sapply(Msum, function(j) lnpochhammer(g, j))
+  # lngpoch <- sapply(Msum, function(j) lnpochhammer(g, j))
+  lngpoch <- sapply(Msumunique, function(j) lnpochhammer(g, j))
+  names(lngpoch) <- Msumunique
   
-  res1 <- lnapoch + lnbpoch - lngpoch - lnfact
+  # 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))
   
   kstep <- 5
@@ -102,24 +142,44 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
     M <- rbind( M, expand.grid(rep(list(k1), n)) )
     
     Msum <- apply(M, 1, sum)
+
+    Munique <- (max(Munique)+1):max(M)
+    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)))
     
     # 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))
+    # lnapoch <- sapply(Msum, function(j) lnpochhammer(a, j))
+    lnapochsupp <- sapply(Msumunique, function(j) lnpochhammer(a, j))
+    names(lnapochsupp) <- Msumunique
+    lnapoch <- c(lnapoch, lnapochsupp)
     
-    lnbpoch <- numeric(nrow(M))
-    for (i in 1:nrow(M)) {
+    # lnbpoch <- numeric(nrow(M))
+    # for (i in 1:nrow(M)) {
+    #   # Product pochhammer(b_1,m_1) * ...* pochhammer(b_n, m_n)
+    #   lnbpoch[i] <- sum( sapply(1:n, function(j) lnpochhammer(b[j], M[i, j])) )
+    # }
+    lnbpochsupp <- 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)
-      lnbpoch[i] <- sum( sapply(1:n, function(j) lnpochhammer(b[j], M[i, j])) )
+      lnbpochsupp[i, j] <- lnpochhammer(b[j], Munique[i])
+    }
+    lnbpoch <- rbind(lnbpoch, lnbpochsupp)
+    sumlnpochb <- function(ind) {
+      sum(mapply(function(i, j) lnbpoch[i, j], ind+1, 1:n))
     }
     
     # pochhammer(c, m_1+...+m_n) for m_1 = 0...k, ...,  m_p = 0...k
-    lngpoch <- sapply(Msum, function(j) lnpochhammer(g, j))
+    # lngpoch <- sapply(Msum, function(j) lnpochhammer(g, j))
+    lngpochsupp <- sapply(Msumunique, function(j) lnpochhammer(g, j))
+    names(lngpochsupp) <- Msumunique
+    lngpoch <- c(lngpoch, lngpochsupp)
     
-    res1 <- lnapoch + lnbpoch - lngpoch - lnfact
+    # 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))
   }