R编程效率——数字的阶乘分解

Question

我正在尝试从 R 中称为阶乘分解的代码战中解决一个卡塔。卡塔的目的是分解 n! (阶乘 n) 转化为它的质因数。该函数应返回一个类似

的字符串

decomp(12) -> "2^10 * 3^5 * 5^2 * 7 * 11"

我能够解决它，但达到了服务器超时(通过 74 个分配)。我尝试对其进行一些优化(lapply of pointwise())，但我无法更改基本核心(for-while-loop)。

任何帮助将不胜感激，因为我已经投入了比我应该拥有的更多的时间。

##' A function for a factorial decomposition of a number
##' @title decomp
##' @param n integer
##' @return a String with the factorial decomposition
##' @author krisselack

decomp <- function(n) {

  # https://stackoverflow.com/questions/19767408/prime-number-function-in-r
  is.prime <- function(n) n == 2L || all(n %% 2L:ceiling(sqrt(n)) != 0)

  p <- 2:n
  primes <- p[as.logical(vapply(p, is.prime, 1))]
  erg <- NULL

  pointwise <- function(x) {

    primloop <- primes[primes<=x]

    for(j in primloop){

      while(x %% j == 0){
        x <- x/j
        erg <- c(erg, j)
      }
    }

    if(length(erg)>0)
      return(erg)
  }

  erg2 <- unlist(lapply(p, pointwise))

  ergfin <- table(erg2)

  namen <- paste(ifelse(ergfin>1, paste0(names(ergfin), "^", ergfin),
                        paste(names(ergfin))),
                 collapse = " * ")

  return(namen)
}

decomp(5) # -> "2^3 * 3 * 5"
decomp(12) # -> "2^10 * 3^5 * 5^2 * 7 * 11"
decomp(17) # -> "2^15 * 3^6 * 5^3 * 7^2 * 11 * 13 * 17"
decomp(25) # -> "2^22 * 3^10 * 5^6 * 7^3 * 11^2 * 13 * 17 * 19 * 23"

Answer 1

library("purrr")

# https://stackoverflow.com/questions/19767408/prime-number-function-in-r
is.prime <- function(n) n == 2L || all(n %% 2L:ceiling(sqrt(n)) != 0)

#' Multiplicity of prime p in the decomp of n!
#' @param p A prime
#' @param n An integer
multiplicity_in_factorial <- function(p, n) {
  # Adding epsilon to avoid rounding errors.
  # For example at p = 3, n = 243
  max_mul <- floor(log(n) / log(p) + 0.0001)
  prime_mul <- p ^ (1:max_mul)
  how_many_of_each <- map_dbl(prime_mul, ~ floor(n / .))
  sum(how_many_of_each)
}



decomp2 <- function(n) {
  p <- 2:n
  primes <- p[as.logical(vapply(p, is.prime, 1))]

  primes_mul <- map_dbl(primes, multiplicity_in_factorial, n)

  namen <- paste(ifelse(primes_mul > 1,
                        paste0(primes, "^", primes_mul),
                        primes),
                 collapse = " * ")

  return(namen)
}

check <- function(n) {
  decomp(n) == decomp2(n)
}

想法是循环n 下的素数，并计算出它们出现在阶乘中的频率。

The key is that the multiplicity of p in n! is the sum of the multiplicities of p in k for k = 1..n.

为了说明，n = 100 且 p = 2。在 1 到 100 之间有 50 个 2 的倍数。但这并没有考虑多重性 > 1 的因素。

我们还必须考虑 4(有 25)、8(有 12)、16(有 6)、32(有 3)和 64(有 1)的倍数。

这就是 multiplicity in factorial 中发生的情况。剩下的就很简单了。

高值的瓶颈是素数的计算，这可以通过使用 Eratosthenes 筛法得到改进。

# https://gist.github.com/seankross/5946396
microbenchmark::microbenchmark(
   sieve = sieveOfEratosthenes(N),
   naive_filter = {
     p <- 2:N
     primes <- p[as.logical(vapply(p, is.prime, 1))]
   }
 )

Unit: microseconds
          expr      min        lq      mean    median       uq      max neval
         sieve  395.010  405.4015  423.2184  410.8445  439.629   584.71   100
  naive_filter 2875.782 2936.5195 3268.4576 2979.4925 3016.060 16875.81   100

但我不会打扰:真正的瓶颈将是众所周知的缓慢的字符串粘贴。

在我的笔记本电脑上，decomp2(10e5) 需要几秒钟，decomp2(10e6) 需要大约 2 分钟。我 99% 确定字符串粘贴实际上是这种情况下的瓶颈。