diff --git a/R/lauricella.R b/R/lauricella.R
index 8a5d3dfac1f0af5ee81562eef91449a6dab308b7..459631bfcf0174732a8007eb72031a68897e0762 100644
--- a/R/lauricella.R
+++ b/R/lauricella.R
@@ -35,29 +35,7 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
   #' \doi{10.1109/LSP.2019.2915000}
   #' @importFrom utils combn
   #' @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)
   
@@ -79,12 +57,13 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
   # # (i.e. all arrangements of n elements from {0:k})
   M <- expand.grid(rep(list(0:k), n))
 
-  Msum <- apply(M, 1, sum)
+  # Sum of the indices
+  Msum <- rowSums(M)
   
   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
+  # Product x^{m_1} * ... * x^{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 <- as.data.frame(
@@ -100,6 +79,8 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
   gridxfact <- expand.grid(xfact)
   prodxfact <- apply(gridxfact, 1, prod)
   
+  # Logarithm of the product m_1! * ... * m_n! for m_1 = 0...k, ...,  m_n = 0...k
+  # i.e. \sum_{i=0}^n{\log{m_i!}}
   lnfact <- as.data.frame(
     matrix(sapply(Munique, lnfactorial), nrow = length(Munique),
            ncol = n, byrow = FALSE, dimnames = list(Munique, 1:n)))