用于生成数字所有可能因式分解的 R 算法

Question

例如，考虑数字96，它可以写成以下几种方式:

1. 96
2. 48 * 2
3. 24 * 2 * 2
4. 12 * 2 * 2 * 2
5. 6 * 2 * 2 * 2 * 2
6. 3 * 2 * 2 * 2 * 2 * 2
7. 4 * 3 * 2 * 2 * 2
8. 8 * 3 * 2 * 2
9. 6 * 4 * 2 * 2
10. 16 * 3 * 2
11. 4 * 4 * 3 * 2
12. 12 * 4 * 2
13. 8 * 6 * 2
14. 32 * 3
15. 8 * 4 * 3
16. 24 * 4
17. 6 * 4 * 4
18. 16 * 6
19. 12 * 8

我知道这与分区有关，因为任何数字写成单个素数 p 的幂 n，就是您可以写的方式的数量n。例如，要找到 2^5 的所有因式分解，我们必须找到所有写 5 的方法。它们是:

1. 1+1+1+1+1 ==>> 2^1 * 2^1 * 2^1 * 2^1 * 2^1
2. 1+1+1+2 ==>> 2^1 * 2^1 * 2^1 * 2^2
3. 1+1+3 ==>> 2^1 * 2^1 * 2^3
4. 1+2+2 ==>> 2^1 * 2^2 * 2^2
5. 1+4 ==>> 2^1 * 2^4
6. 2+3 ==>> 2^2 * 2^3
7. 5 ==>> 2^5

我发现了一篇由 Jerome Kelleher 撰写的关于分区生成算法的精彩文章 here .我已经将他的一个 python 算法改编为 R。代码如下:

library(partitions)   ## using P(n) to determine number of partitions of an integer
IntegerPartitions <- function(n) {
    a <- 0L:n
    k <- 2L
    a[2L] <- n
    MyParts <- vector("list", length=P(n))
    count <- 0L
    while (!(k==1L)) {
        x <- a[k-1L]+1L
        y <- a[k]-1L
        k <- k-1L
        while (x<=y) {a[k] <- x; y <- y-x; k <- k+1L}
        a[k] <- x+y
        count <- count+1L
        MyParts[[count]] <- a[1L:k]
    }
    MyParts
}

我试图将此方法扩展到具有不止一个质因数的数字，但我的代码变得非常笨拙。在与这个想法搏斗了一段时间之后，我决定尝试一条不同的路线。我的新算法不使用任何生成分区。它更像是一种利用已经生成的因式分解的“回溯”算法。代码如下:

FactorRepresentations <- function(n) {

MyFacts <- EfficientFactorList(n)
MyReps <- lapply(1:n, function(x) x)

    for (k in 4:n) {
        if (isprime(k)) {next}
        myset <- MyFacts[[k]]
        mylist <- vector("list")
        mylist[[1]] <- k
        count <- 1L
            for (j in 2:ceiling(length(myset)/2)) {
                count <- count+1L
                temp <- as.integer(k/myset[j])
                myvec <- sort(c(myset[j], temp), decreasing=TRUE)
                mylist[[count]] <- myvec
                MyTempRep <- MyReps[[temp]]

                if (isprime(temp) || temp==k) {next}

                if (length(MyTempRep)>1) {
                    for (i in 1:length(MyTempRep)) {
                        count <- count+1L
                        myvec <- sort(c(myset[j], MyTempRep[[i]]), decreasing=TRUE)
                        mylist[[count]] <- myvec
                    }
                }
            }
        MyReps[[k]] <- unique(mylist)
    }
    MyReps
}

上面代码中的第一个函数只是一个生成所有因子的函数。如果你好奇的话，这里是代码:

EfficientFactorList <- function(n) {
    MyFactsList <- lapply(1:n, function(x) 1)
    for (j in 2:n) {
        for (r in seq.int(j, n, j)) {MyFactsList[[r]] <- c(MyFactsList[[r]], j)}
    }
    MyFactsList
}

如果您只关心小于 10,000 的数字，我的算法还可以(它会在大约 17 秒内为每个 <= 10,000 的数字生成所有因式分解)，但它的扩展性肯定不好。我想找到一种算法，它具有为小于或等于 n 的每个数字生成所有因式分解列表的相同前提，因为我想到的一些应用程序将引用给定的因式分解倍数次，因此将它放在列表中应该比每次都动态生成它更快(我知道这里有内存成本)。

Answer 1

您的函数 EfficientFactorList 可以很好地有效地获取从 1 到 n 的每个数字的所有因子的集合，因此剩下的就是获取所有因式分解的集合。正如您所建议的，使用较小值的因式分解来计算较大值的因式分解似乎很有效。

考虑一个数字 k，其因子为 k_1、k_2、...、k_n。一种朴素的方法是组合 k/k_1、k/k_2、...、k/k_n 的分解，将 k_i 附加到 k/k_i 的每个分解以产生 k 的分解。作为一个有效的例子，考虑计算 16 的因式分解(它有非平凡的因子 2、4 和 8)。 2 有因式分解 {2}，4 有因式分解 {4, 2*2}，而 8 有因式分解 {8, 4*2, 2*2*2}，所以我们会通过首先计算 {2 *8, 4*4, 2*2*4, 8*2, 4*2*2, 2*2*2*2} 然后进行唯一分解，{8*2, 4*4, 4*2 *2, 2*2*2*2}。将 16 相加得到最终答案。

一种更有效的方法是注意我们不需要将 k_i 附加到 k/k_i 的所有因式分解。例如，我们不需要从 4 的因式分解中添加 2*2*4，因为它已经包含在 8 的因式分解中。同样，我们不需要从 2 的因式分解中添加 2*8，因为这已经包含在 8 的因式分解中。一般来说，如果因式分解中的所有值都是 k_i 或更大，我们只需要包含 k/k_i 的因式分解。

在代码中:

library(gmp)
all.fact <- function(n) {
  facts <- EfficientFactorList(n)
  facts[[1]] <- list(1)
  for (x in 2:n) {
    if (length(facts[[x]]) == 2) {
      facts[[x]] <- list(x)  # Prime number
    } else {
      x.facts <- facts[[x]][facts[[x]] != 1 & facts[[x]] <= (x^0.5+0.001)]
      allSmaller <- lapply(x.facts, function(pf) lapply(facts[[x/pf]], function(y) {
        if (all(y >= pf)) {
            return(c(pf, y))
        } else {
            return(NULL)
        }
      }))
      allSmaller <- do.call(c, allSmaller)
      facts[[x]] <- c(x, allSmaller[!sapply(allSmaller, function(y) is.null(y))])
    }
  }
  return(facts)
}

这比发布的代码快得多:

system.time(f1 <- FactorRepresentations(10000))
#    user  system elapsed 
#  13.470   0.159  13.765 
system.time(f2 <- all.fact(10000))
#    user  system elapsed 
#   1.602   0.028   1.641 

作为完整性检查，它还为每个数字返回相同数量的因式分解:

lf1 <- sapply(f1, length)
lf2 <- sapply(f2, length)
all.equal(lf1, lf2)
# [1] TRUE