diff --git a/R/lauricella.R b/R/lauricella.R
index a218f1e2354b2d9f9fa81e172d2c220f39c9cd63..5dbe91f684853f5fd4abcebd97fc23b5fb14bfa3 100644
--- a/R/lauricella.R
+++ b/R/lauricella.R
@@ -30,7 +30,7 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
   #' If this happens, the \code{attr(, "epsilon")} attribute is the precision that was really reached.
   #'
   #' @author Pierre Santagostini, Nizar Bouhlel
-  #' @references N. Bouhlel, A. Dziri, Jullback-Leibler Divergence Between Multivariate Generalized Gaussian Distributions.
+  #' @references N. Bouhlel, A. Dziri, Kullback-Leibler Divergence Between Multivariate Generalized Gaussian Distributions.
   #' IEEE Signal Processing Letters, vol. 26 no. 7, July 2019.
   #' \doi{10.1109/LSP.2019.2915000}
   #' @export
@@ -55,35 +55,37 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
   Msum <- apply(M, 1, sum)
   
   # 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 / factorial(Mi) ))
+  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
-  apoch <- sapply(Msum, function(j) pochhammer(a, j))
+  lnapoch <- sapply(Msum, function(j) lnpochhammer(a, j))
   
-  bpoch <- numeric(nrow(M))
+  lnbpoch <- numeric(nrow(M))
   for (i in 1:nrow(M)) {
     # Product pochhammer(b_1,m_1) * ...* pochhammer(b_n, m_n)
-    bpoch[i] <- prod( sapply(1:n, function(j) pochhammer(b[j], M[i, j])) )
+    lnbpoch[i] <- sum( sapply(1:n, function(j) lnpochhammer(b[j], M[i, j])) )
   }
   
   # pochhammer(c, m_1+...+m_n) for m_1 = 0...k, ...,  m_p = 0...k
-  gpoch <- sapply(Msum, function(j) pochhammer(g, j))
+  lngpoch <- sapply(Msum, function(j) lnpochhammer(g, j))
   
-  res1 <- sum( xfact * apoch * bpoch / gpoch )
+  res1 <- lnapoch + lnbpoch - lngpoch - lnfact
+  res2 <- sum(xfact * exp(res1))
   
   kstep <- 5
   
   k1 <- 1:k
   result <- 0
   
-  while (abs(res1) > eps/10 & !is.nan(res1)) {
+  while (abs(res2) > eps/10 & !is.nan(res2)) {
     
-    epsret <- signif(abs(res1), 1)*10
+    epsret <- signif(abs(res2), 1)*10
     
     k <- k1[length(k1)]
     
     k1 <- k + (1:kstep)
-    result <- result + res1
+    result <- result + res2
     
     # M: data.frame of the indices for the nested sums
     M <- as.data.frame(matrix(nrow = 0, ncol = n))
@@ -102,86 +104,27 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
     Msum <- apply(M, 1, sum)
     
     # 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 / factorial(Mi) ))
+    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
-    apoch <- sapply(Msum, function(j) pochhammer(a, j))
+    lnapoch <- sapply(Msum, function(j) lnpochhammer(a, j))
     
-    bpoch <- numeric(nrow(M))
+    lnbpoch <- numeric(nrow(M))
     for (i in 1:nrow(M)) {
       # Product pochhammer(b_1,m_1) * ...* pochhammer(b_n, m_n)
-      bpoch[i] <- prod( sapply(1:n, function(j) pochhammer(b[j], M[i, j])) )
+      lnbpoch[i] <- sum( sapply(1:n, function(j) lnpochhammer(b[j], M[i, j])) )
     }
     
     # pochhammer(c, m_1+...+m_n) for m_1 = 0...k, ...,  m_p = 0...k
-    gpoch <- sapply(Msum, function(j) pochhammer(g, j))
+    lngpoch <- sapply(Msum, function(j) lnpochhammer(g, j))
     
-    res1 <- sum( xfact * apoch * bpoch / gpoch )
-    
-  }
-  
-  if (is.nan(res1)) {
-    
-    res1 <- 0
-    k1 <- k
-    
-    while(!is.nan(res1)) {
-      
-      if (res1 > 0)
-        epsret <- abs(res1)*10
-      
-      k <- k1
-      
-      k1 <- k + 1
-      result <- result + res1
-      
-      # M: data.frame of the indices for the nested sums
-      M <- as.data.frame(matrix(nrow = 0, ncol = n))
-      if (n > 1) {
-        for (i in 1:(n-1)) {
-          Mlist <- c( rep(list(0:k), n-i), rep(list(k1), i) )
-          M <- rbind( M, expand.grid(Mlist) )
-          for (j in 1:(n-1)) {
-            Mlist <- Mlist[c(n, 1:(n-1))]
-            M <- rbind(M, expand.grid(Mlist))
-          }
-        }
-      }
-      M <- rbind( M, rep(k1, n) )
-      
-      Msum <- apply(M, 1, sum)
-      
-      # 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 / factorial(Mi) ))
-      
-      # pochhammer(a, m_1+...+m_n) for m_1 = 0...k, ...,  m_p = 0...k
-      apoch <- sapply(Msum, function(j) pochhammer(a, j))
-      
-      bpoch <- numeric(nrow(M))
-      for (i in 1:nrow(M)) {
-        # Product pochhammer(b_1,m_1) * ...* pochhammer(b_n, m_n)
-        bpoch[i] <- prod( sapply(1:n, function(j) pochhammer(b[j], M[i, j])) )
-      }
-      
-      # pochhammer(c, m_1+...+m_n) for m_1 = 0...k, ...,  m_p = 0...k
-      gpoch <- sapply(Msum, function(j) pochhammer(g, j))
-      
-      res1 <- sum( xfact * apoch * bpoch / gpoch )
-      
-    }
-  }
-  
-  if (is.nan(res1)) {
-    epsret <- signif(epsret, 1)
-    warning("\nCannot reach the precision ", eps, " due to NaN\n",
-            "Number of iteration: ", k, "\n",
-            "Precision reached:", epsret)
-    attr(result, "epsilon") <- epsret
-  } else {
-    # result <- result + res1
-    attr(result, "epsilon") <- eps
+    res1 <- lnapoch + lnbpoch - lngpoch - lnfact
+    res2 <- sum(xfact * exp(res1))
   }
   
+  result <- Re(result)
+  attr(result, "epsilon") <- eps
   attr(result, "k") <- k
   
   # Returns the result of the nested sums