R编程效率 - 一个数字的阶乘分解

Question

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

的字符串

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

我能够解决它，但达到了服务器超时（通过 74 次分配）。我尝试对其进行一些优化（pointwise() 的 lapply），但我无法更改基本核心（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 以下的素数上循环，并找出它们在阶乘中出现的频率。

关键是p在n中的多重性！是 k = 1..n 时 p 在 k 中的重数之和。

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

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

这就是multiplicity in factorial 中发生的事情。其余的很简单。

高值的瓶颈是primes 的计算，可以通过使用 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% 确定字符串粘贴实际上是瓶颈。

Answer 2

您可以使用包 {primefactr} 来做到这一点：

> primefactr::ReducePrime(c(rep(0, 999), 1), out.summary = TRUE)
       [,1] [,2]
primes    2    5
power     3    3
> primefactr::ReducePrime(c(rep(0, 9999), 1), out.summary = TRUE)
       [,1] [,2]
primes    2    5
power     4    4
> primefactr::ReducePrime(c(rep(0, 999999), 1), out.summary = TRUE)
       [,1] [,2]
primes    2    5
power     6    6