用于生成数字所有可能因式分解的 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 ==>> 2^1 * 2^1 * 2^1 * 2^1 * 2^1
1+1+1+2 ==>> 2^1 * 2^1 * 2^1 * 2^2
1+1+3 ==>> 2^1 * 2^1 * 2^3
1+2+2 ==>> 2^1 * 2^2 * 2^2
1+4 ==>> 2^1 * 2^4
2+3 ==>> 2^2 * 2^3
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
}

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

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

作为完整性检查，它还 returns 每个数字的分解次数相同：

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

Answer 2

万一有人对生成一个数字 n 的乘法分区感兴趣，下面是两个算法可以做到这一点（函数 IntegerPartition 来自上面的问题):

library(gmp)
library(partitions)
get_Factorizations1 <- function(MyN) {
    pfs <- function (x1) {
        n1 <- length(x1)
        y1 <- x1[-1L] != x1[-n1]
        i <- c(which(y1), n1)
        list(lengths = diff(c(0L, i)), values = x1[i], uni = sum(y1)+1L)
    }

    if (MyN==1L) return(MyN)
    else {
        pfacs <- pfs(as.integer(factorize(MyN)))
        unip <- pfacs$values
        pv <- pfacs$lengths
        n <- pfacs$uni
        mySort <- order(pv, decreasing = TRUE)
        pv <- pv[mySort]
        unip <- unip[mySort]
        myReps <- lapply(IntegerPartitions(pv[1L]), function(y) unip[1L]^y)
        if (n > 1L) {
            mySet <- unlist(lapply(2L:n, function(x) rep(unip[x],pv[x])))
            for (p in mySet) {
                myReps <- unique(do.call(c,
                    lapply(myReps, function(j) { 
                        dupJ <- duplicated(j)
                        nDupJ <- !dupJ
                        SetJ <- j[which(nDupJ)]
                        lenJ <- sum(nDupJ)
                        if (any(dupJ)) {v1 <- j[which(dupJ)]} else {v1 <- vector(mode="integer")}
                        tList <- vector("list", length=lenJ+1L)
                        tList[[1L]] <- sort(c(j,p))

                        if (lenJ > 1L) {c2 <- 1L
                            for (a in 1:lenJ) {tList[[c2 <- c2+1L]] <- sort(c(v1,SetJ[-a],SetJ[a]*p))}
                        } else {
                            tList[[2L]] <- sort(c(v1,p*SetJ))
                        }
                        tList
                    }
                )))
            }
        }
    }
    myReps
}

下面是来自上面的 josliber 的代码，用于处理单个案例。函数 MyFactors 来自这个 post（它 returns 给定数字的所有因数）。

library(gmp)
get_Factorizations2 <- function(n) {
    myFacts <- as.integer(MyFactors(n))
    facts <- lapply(myFacts, function(x) 1L)
    numFacs <- length(myFacts)
    facts[[numFacs]] <- myFacts
    names(facts) <- facts[[numFacs]]
    for (j in 2L:numFacs) {
        x <- myFacts[j]
        if (isprime(x)>0L) {
            facts[[j]] <- list(x)
        } else {
            facts[[j]] <- myFacts[which(x%%myFacts[myFacts <= x]==0L)]
            x.facts <- facts[[j]][facts[[j]] != 1 & facts[[j]] <= (x^0.5+0.001)]
            allSmaller <- lapply(x.facts, function(pf) lapply(facts[[which(names(facts)==(x/pf))]], function(y) {
                if (all(y >= pf)) {
                    return(c(pf, y))
                } else {
                    return(NULL)
                }
            }))
            allSmaller <- do.call(c, allSmaller)
            facts[[j]] <- c(x, allSmaller[!sapply(allSmaller, function(y) is.null(y))])
        }
    }
    facts[[numFacs]]
}

以下是一些基准：

set.seed(101)
samp <- sample(10^7, 10^4)
library(rbenchmark)
benchmark(getFacs1=sapply(samp, get_Factorizations),
            getFacs2=sapply(samp, get_Factorizations2),
            replications=5,
            columns = c("test", "replications", "elapsed", "relative"),
            order = "relative")
test replications elapsed relative
1 getFacs1            5  117.68    1.000
2 getFacs2            5  216.39    1.839


system.time(t2 <- get_Factorizations(25401600))
 user  system elapsed 
10.89    0.03   10.97
system.time(t2 <- get_Factorizations2(25401600))
 user  system elapsed 
21.08    0.00   21.12 

length(t1)==length(t2)
[1] TRUE

object.size(t1)
28552768 bytes
object.size(t2)
20908768 bytes

尽管 get_Factorizations1 更快，但第二种方法更直观（请参阅上面 josliber 的出色解释）并且它生成的对象更小。对于感兴趣的 reader，here 是一篇关于该主题的非常好的论文。

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

R Algorithm for generating all possible factorizations of a number

algorithm

r

factorization