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

Question

我正在尝试解决代码战中称为 R 中的阶乘分解的套路。套路的目的是分解 n！ (阶乘 n) 转化为它的质因数。该函数应该 return 一个像

这样的字符串

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) {

  # 
  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")

# 
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 中发生的事情。剩下的就简单了。

高值的瓶颈是 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

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

R Programming Efficiency - factorial decomposition of a number

performance

r

factorial