指定的向量大小对于 32*32 矩阵来说太大
vector size specified is too large for 32*32 matrix
我在 R 中尝试了以下代码来计算 matrix 的永久性。但是我得到这个错误：
Error in vector("list", gamma(n + 1)) : vector size specified is too large.
我尝试了很多解决方案，但都无法解决问题。
library(combinat)
perm <- function(A)
{
stopifnot(is.matrix(A))
n <- nrow(A)
if(n != ncol(A)) stop("Matrix is not square.")
if(n < 1) stop("Matrix has a dimension of size 0.")
sum(sapply(combinat::permn(n), function(sigma) prod(sapply(1:n, function(i) A[i, sigma[i]]))))
}
#We copy our test cases from the Python example.
testData <- list("A" = rbind(c(0.5,0.3,0.4,0.4,0.4,0.3,0.3,0.4,0.3,0.4,0.5,0.6,0.5,0,0,0.3,0,0,0,0.6,0.6,0.4,0.1,0.2,0.4,0.7,0.4,0.7,0.8,0.9,0.4,0.4),
c(0.7,0.7,0.7,0.7,0.7,0.5,0.6,0.5,0.6,0.7,0.8,0.8,0.8,0,0,0.7,0,0.6,0.7,0.9,0.8,0.7,0.5,0.5,0.7,0,0.7,0.8,0,0.9,0.6,0.7),
c(0.6,0.3,0.5,0.7,0.7,0.4,0.6,0.6,0.6,0.6,0.7,0.7,0.8,0,0,0.6,0,0.6,0.6,0.8,0.8,0.6,0.4,0.5,0.8,0.9,0.6,0.8,0,1,0.6,0.7),
c(0.6,0.3,0.3,0.5,0.7,0.4,0.5,0.7,0.4,0.6,0.7,0.7,0.8,0.4,0.4,0.5,0,0.5,0.6,0.8,0.7,0.5,0.4,0.4,0.7,0.9,0.6,0.8,0,0.9,0.6,0.7),
c(0.6,0.3,0.3,0.3,0.5,0.3,0.4,0.4,0.4,0.4,0.5,0.6,0.6,0,0,0.3,0,0,0,0.7,0.6,0.4,0.3,0.3,0.4,0,0.4,0.6,0,0.7,0.5,0.4),
c(0.7,0.5,0.6,0.6,0.7,0.6,0.6,0.6,0.6,0.6,0.8,0.8,0.8,0,0,0.6,0,0.6,0.6,0.9,0.7,0.6,0.5,0.6,0.7,0.9,0.6,0.8,0,0.9,0.6,0.7),
c(0.7,0.4,0.4,0.5,0.6,0.4,0.4,0.5,0.4,0.5,0.7,0.7,0.7,0,0,0.4,0,0.4,0.7,0.8,0.8,0.6,0.3,0.5,0.7,0.9,0.6,0.7,0,0.9,0.6,0.7),
c(0.6,0.5,0.4,0.3,0.6,0.4,0.5,0.6,0.4,0.5,0.7,0.7,0.7,0,0,0.6,0,0,0,0.8,0.7,0.6,0.4,0.5,0.7,0,0.6,0.7,0,0,0.6,0.6),
c(0.7,0.4,0.4,0.6,0.6,0.4,0.6,0.6,0.5,0.7,0.7,0.7,0.7,0,0,0.6,0,0.5,0.6,0.8,0.7,0.6,0.4,0.5,0.7,0,0.6,0.7,0.9,0.9,0.7,0.8),
c(0.6,0.3,0.4,0.4,0.6,0.4,0.5,0.5,0.3,0.6,0.7,0.7,0.7,0,0,0.6,0,0,0,0.8,0.6,0.6,0.5,0.5,0.6,0.8,0.6,0.6,0,0,0.6,0.7),
c(0.5,0.2,0.3,0.3,0.5,0.2,0.3,0.3,0.3,0.3,0.4,0.6,0.5,0,0,0.4,0,0.3,0.4,0.7,0.5,0.4,0.2,0.2,0.4,0,0.4,0.5,0.7,0.7,0.4,0.4),
c(0.4,0.2,0.3,0.3,0.4,0.2,0.3,0.3,0.3,0.3,0.4,0.4,0.5,0,0,0.3,0,0.4,0.4,0.7,0.5,0.3,0.3,0.2,0.4,0,0.4,0.5,0.8,0.9,0.4,0.4),
c(0.5,0.2,0.2,0.2,0.4,0.2,0.3,0.3,0.3,0.3,0.5,0.5,0.3,0,0,0.3,0,0,0,0.8,0.4,0.4,0.1,0.1,0.4,0,0.4,0.4,0.8,0.8,0.4,0.4),
c(1,1,1,0.6,1,1,1,1,1,1,1,1,1,0.5,0.7,0.7,0,0.7,0.7,0.9,0.7,0.7,0.5,0.6,0.7,0.9,0.7,0.8,0.9,1,0.7,0.8),
c(1,1,1,0.6,1,1,1,1,0,1,1,0,1,0,0.4,0.5,0,0.5,0.5,0.9,0,0.6,0.4,0.5,0.7,0.9,0.6,0,0,0.9,0.6,0.7),
c(0.7,0.3,0.4,0.5,0.7,0.4,0.6,0.4,0,0.4,0.6,0,0.7,0,0,0.6,0,0.6,0.7,0.8,0.8,0.6,0.4,0.5,0.7,0,0.6,0,0,0.9,0.7,0.7),
c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0.3,0.7,0.7,1,0.7,0.7,0.4,0.6,0.7,1,0.7,0.9,1,1,0.7,0.7),
c(1,0.4,0.4,0.5,1,0.4,0.6,1,0.5,1,0.7,0.6,1,0,0,0.4,0,0.2,0.6,0.8,0.7,0.7,0.4,0.5,0.7,0,0.7,0.8,0,1,0.8,0.7),
c(1,0.3,0.4,0.4,1,0.4,0.3,1,0.4,1,0.6,0,1,0,0,0.3,0,0.4,0.1,0.8,0,0.6,0.4,0.6,0.7,0,0.7,0,0,0.9,0.7,0.7),
c(0.4,0.1,0.2,0.2,0.3,0.1,0.2,0.2,0.2,0.2,0.3,0,0.2,0,0,0.2,0,0.2,0.2,0.5,0,0.3,0.1,0.1,0.3,0.7,0.4,0.4,0.7,0.8,0.4,0.4),
c(0.4,0.2,0.2,0.3,0.4,0.3,0.2,0.3,0.3,0.4,0.5,0.5,0.6,0.3,1,0.2,0,0.3,1,1,0.4,0.4,0.3,0.3,0.4,0.8,0.5,0.7,0.9,0.9,0.4,0.4),
c(0.6,0.3,0.4,0.5,0.6,0.4,0.4,0.4,0,0.4,0.6,0,0.6,0,0,0.4,0,0.3,0.4,0.7,0,0.5,0.3,0.4,0.5,0,0.5,0,0,1,0.6,0.6),
c(0.9,0.5,0.6,0.6,0.7,0.5,0.7,0.6,0.6,0.5,0.8,0.7,0.9,0,0,0.6,0,0.6,0.6,0.9,0.7,0.7,0.5,0.7,0.8,0,0.6,0.8,0,0,0.6,0.7),
c(0.8,0.5,0.5,0.6,0.7,0.4,0.5,0.5,0.5,0.5,0.8,0.8,0.9,0,0,0.5,0,0,0,0.9,0.7,0.6,0.3,0.6,0.8,0,0.7,0.8,0.9,0.9,0.7,0.8),
c(0.6,0.3,0.2,0.3,0.6,0.3,0.3,0.3,0.3,0.4,0.6,0.6,0.6,0.3,0.3,0.3,0,0.3,0.3,0.7,0.6,0.5,0.2,0.2,0.4,0.8,0.5,0.7,0.9,0.9,0.5,0.6),
c(0.3,1,0.1,0,0,0.1,0.1,0,0,0,0,0,0,0,0,0,0,0,0,0.3,0,0,0,0,0.2,0.2,0.3,0.3,0.5,0.5,0,0),
c(0.6,0.3,0.4,0,0,0.4,0.4,0.4,0.4,0.4,0.6,0,0.6,0,0,0.4,0,0.3,0.3,0.6,0,0.5,0.4,0.3,0.5,0,0.2,0,0,0.8,0.4,0.5),
c(0.3,0.2,0.2,0.2,0.4,0.2,0.3,0.3,0.3,0.4,0.5,0.5,0.6,0,0,1,0,0,0,0.6,0.3,1,0.2,0.2,0.3,0.7,1,0.4,0.9,0.9,0.5,0.5),
c(0.2,1,0,0,0,1,0,1,0.1,1,0.3,0.2,0.2,0,0,0,0,0,0,0,0,0,1,0.1,0.1,0.5,0,0.1,0.5,0.5,0,0),
c(0.1,0.1,0,0,0,0.1,0,1,0.1,1,0.3,0.1,0.2,0,0,0,0,0,0,0,0,0,1,0.1,0.1,0.5,0,0.1,0.5,0.4,0,0),
c(0.6,0.4,0.4,0.4,0.5,0.4,0.4,0.4,0.3,0.4,0.6,0.6,0.6,0,0,0.3,0,0,0,0.6,0.6,0.4,0.4,0.3,0.5,0,0.6,0.5,1,1,0.4,0.7),
c(0.6,0.3,0.3,0,0,0.3,0.3,0.4,0.2,0.3,0.6,0,0.6,0,0,0.3,0,0,0.3,0.6,0,0.4,0.3,0.2,0.4,0,0.5,0.5,0,1,0.3,0.3)))
print(sapply(testData, function(x) list( Permanent = perm(x))))
计算矩阵永久物是一个计算难题（例如，“计算永久物”有其专用的 Wikipedia article 与操作的数学定义分开......）
您实施的原始方法（已讨论 here）将需要计算长度为 32 的向量！ = 2.6e35（这将占用超过 10^24 TB 的内存，更不用说完成所有计算所花费的时间......）
幸运的是，自从上次在 SO 上讨论这个问题（2014 年）以来，有人已经将 BosonSampling 包发布到 CRAN，它使用编译代码（和一种高效的算法！）来完成计算。 tl;dr BosonSampling::rePerm(testData$A) 在我的笔记本电脑上不到 2 分钟就完成了。
设置
## repeated definition from above, for convenience
library(BosonSampling)
library(combinat)
perm <- function(A) {
stopifnot(is.matrix(A))
n <- nrow(A)
if(n != ncol(A)) stop("Matrix is not square.")
if(n < 1) stop("Matrix has a dimension of size 0.")
sum(sapply(combinat::permn(n), function(sigma) prod(sapply(1:n, function(i) A[i, sigma[i]]))))
}
检查小数据上简单和复杂解决方案的等价性
set.seed(101)
t1 <- matrix(rnorm(16), 4, 4)
t2 <- matrix(rnorm(36), 6, 6)
stopifnot(all.equal(perm(t1), rePerm(t1)))
stopifnot(all.equal(perm(t2), rePerm(t2)))
评估时间缩放
（在我的笔记本电脑上这需要几分钟）
f <- function(n, type = "naive") {
set.seed(101)
m <- matrix(rnorm(n^2), n, n)
t <- system.time(if (type == "naive") perm(m) else rePerm(m))
return(t[["elapsed"]])
}
n0 <- 3; n1 <- 10; n2 <- 32
ntimes <- sapply(n0:n1, f)
ptimes <- sapply(n0:n2, f, type = "rePerm")
dd <- data.frame(n = c(n0:n1, n0:n2),
time = c(ntimes, ptimes),
type = rep(c("naive", "rePerm"), c(n1-n0+1, n2-n0+1)))
saveRDS(dd, file = "SO72004136.rds")
library(ggplot2); theme_set(theme_bw())
(ggplot(dd) +
aes(n, time, colour = type) +
geom_point() + geom_line() +
scale_y_log10() +
scale_x_log10() +
expand_limits(x=32) +
labs(y = "elapsed time (sec)") +
geom_smooth(method = "lm",
lty = 2,
colour = "black",
data = subset(dd, type == "rePerm" & n>16))
)
lm(log(time) ~ log(n),
data = subset(dd, type == "rePerm" & n>16))
运行 n=10 的天真方法需要 190 秒； rePerm 对于 n=29 需要 18 秒。计算时间尺度大致为 time \propto n^17（即 log-log 回归的斜率约为 17），这非常糟糕...
一旦我发现这个问题可能需要几分钟而不是几个小时，我就在 115 秒内计算出 rePerm(testData$A) = 1.271369e+20。 （我真诚地希望你不必做比这更大的计算！）
我在 R 中尝试了以下代码来计算 matrix 的永久性。但是我得到这个错误：
Error in vector("list", gamma(n + 1)) : vector size specified is too large.
我尝试了很多解决方案，但都无法解决问题。
library(combinat)
perm <- function(A)
{
stopifnot(is.matrix(A))
n <- nrow(A)
if(n != ncol(A)) stop("Matrix is not square.")
if(n < 1) stop("Matrix has a dimension of size 0.")
sum(sapply(combinat::permn(n), function(sigma) prod(sapply(1:n, function(i) A[i, sigma[i]]))))
}
#We copy our test cases from the Python example.
testData <- list("A" = rbind(c(0.5,0.3,0.4,0.4,0.4,0.3,0.3,0.4,0.3,0.4,0.5,0.6,0.5,0,0,0.3,0,0,0,0.6,0.6,0.4,0.1,0.2,0.4,0.7,0.4,0.7,0.8,0.9,0.4,0.4),
c(0.7,0.7,0.7,0.7,0.7,0.5,0.6,0.5,0.6,0.7,0.8,0.8,0.8,0,0,0.7,0,0.6,0.7,0.9,0.8,0.7,0.5,0.5,0.7,0,0.7,0.8,0,0.9,0.6,0.7),
c(0.6,0.3,0.5,0.7,0.7,0.4,0.6,0.6,0.6,0.6,0.7,0.7,0.8,0,0,0.6,0,0.6,0.6,0.8,0.8,0.6,0.4,0.5,0.8,0.9,0.6,0.8,0,1,0.6,0.7),
c(0.6,0.3,0.3,0.5,0.7,0.4,0.5,0.7,0.4,0.6,0.7,0.7,0.8,0.4,0.4,0.5,0,0.5,0.6,0.8,0.7,0.5,0.4,0.4,0.7,0.9,0.6,0.8,0,0.9,0.6,0.7),
c(0.6,0.3,0.3,0.3,0.5,0.3,0.4,0.4,0.4,0.4,0.5,0.6,0.6,0,0,0.3,0,0,0,0.7,0.6,0.4,0.3,0.3,0.4,0,0.4,0.6,0,0.7,0.5,0.4),
c(0.7,0.5,0.6,0.6,0.7,0.6,0.6,0.6,0.6,0.6,0.8,0.8,0.8,0,0,0.6,0,0.6,0.6,0.9,0.7,0.6,0.5,0.6,0.7,0.9,0.6,0.8,0,0.9,0.6,0.7),
c(0.7,0.4,0.4,0.5,0.6,0.4,0.4,0.5,0.4,0.5,0.7,0.7,0.7,0,0,0.4,0,0.4,0.7,0.8,0.8,0.6,0.3,0.5,0.7,0.9,0.6,0.7,0,0.9,0.6,0.7),
c(0.6,0.5,0.4,0.3,0.6,0.4,0.5,0.6,0.4,0.5,0.7,0.7,0.7,0,0,0.6,0,0,0,0.8,0.7,0.6,0.4,0.5,0.7,0,0.6,0.7,0,0,0.6,0.6),
c(0.7,0.4,0.4,0.6,0.6,0.4,0.6,0.6,0.5,0.7,0.7,0.7,0.7,0,0,0.6,0,0.5,0.6,0.8,0.7,0.6,0.4,0.5,0.7,0,0.6,0.7,0.9,0.9,0.7,0.8),
c(0.6,0.3,0.4,0.4,0.6,0.4,0.5,0.5,0.3,0.6,0.7,0.7,0.7,0,0,0.6,0,0,0,0.8,0.6,0.6,0.5,0.5,0.6,0.8,0.6,0.6,0,0,0.6,0.7),
c(0.5,0.2,0.3,0.3,0.5,0.2,0.3,0.3,0.3,0.3,0.4,0.6,0.5,0,0,0.4,0,0.3,0.4,0.7,0.5,0.4,0.2,0.2,0.4,0,0.4,0.5,0.7,0.7,0.4,0.4),
c(0.4,0.2,0.3,0.3,0.4,0.2,0.3,0.3,0.3,0.3,0.4,0.4,0.5,0,0,0.3,0,0.4,0.4,0.7,0.5,0.3,0.3,0.2,0.4,0,0.4,0.5,0.8,0.9,0.4,0.4),
c(0.5,0.2,0.2,0.2,0.4,0.2,0.3,0.3,0.3,0.3,0.5,0.5,0.3,0,0,0.3,0,0,0,0.8,0.4,0.4,0.1,0.1,0.4,0,0.4,0.4,0.8,0.8,0.4,0.4),
c(1,1,1,0.6,1,1,1,1,1,1,1,1,1,0.5,0.7,0.7,0,0.7,0.7,0.9,0.7,0.7,0.5,0.6,0.7,0.9,0.7,0.8,0.9,1,0.7,0.8),
c(1,1,1,0.6,1,1,1,1,0,1,1,0,1,0,0.4,0.5,0,0.5,0.5,0.9,0,0.6,0.4,0.5,0.7,0.9,0.6,0,0,0.9,0.6,0.7),
c(0.7,0.3,0.4,0.5,0.7,0.4,0.6,0.4,0,0.4,0.6,0,0.7,0,0,0.6,0,0.6,0.7,0.8,0.8,0.6,0.4,0.5,0.7,0,0.6,0,0,0.9,0.7,0.7),
c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0.3,0.7,0.7,1,0.7,0.7,0.4,0.6,0.7,1,0.7,0.9,1,1,0.7,0.7),
c(1,0.4,0.4,0.5,1,0.4,0.6,1,0.5,1,0.7,0.6,1,0,0,0.4,0,0.2,0.6,0.8,0.7,0.7,0.4,0.5,0.7,0,0.7,0.8,0,1,0.8,0.7),
c(1,0.3,0.4,0.4,1,0.4,0.3,1,0.4,1,0.6,0,1,0,0,0.3,0,0.4,0.1,0.8,0,0.6,0.4,0.6,0.7,0,0.7,0,0,0.9,0.7,0.7),
c(0.4,0.1,0.2,0.2,0.3,0.1,0.2,0.2,0.2,0.2,0.3,0,0.2,0,0,0.2,0,0.2,0.2,0.5,0,0.3,0.1,0.1,0.3,0.7,0.4,0.4,0.7,0.8,0.4,0.4),
c(0.4,0.2,0.2,0.3,0.4,0.3,0.2,0.3,0.3,0.4,0.5,0.5,0.6,0.3,1,0.2,0,0.3,1,1,0.4,0.4,0.3,0.3,0.4,0.8,0.5,0.7,0.9,0.9,0.4,0.4),
c(0.6,0.3,0.4,0.5,0.6,0.4,0.4,0.4,0,0.4,0.6,0,0.6,0,0,0.4,0,0.3,0.4,0.7,0,0.5,0.3,0.4,0.5,0,0.5,0,0,1,0.6,0.6),
c(0.9,0.5,0.6,0.6,0.7,0.5,0.7,0.6,0.6,0.5,0.8,0.7,0.9,0,0,0.6,0,0.6,0.6,0.9,0.7,0.7,0.5,0.7,0.8,0,0.6,0.8,0,0,0.6,0.7),
c(0.8,0.5,0.5,0.6,0.7,0.4,0.5,0.5,0.5,0.5,0.8,0.8,0.9,0,0,0.5,0,0,0,0.9,0.7,0.6,0.3,0.6,0.8,0,0.7,0.8,0.9,0.9,0.7,0.8),
c(0.6,0.3,0.2,0.3,0.6,0.3,0.3,0.3,0.3,0.4,0.6,0.6,0.6,0.3,0.3,0.3,0,0.3,0.3,0.7,0.6,0.5,0.2,0.2,0.4,0.8,0.5,0.7,0.9,0.9,0.5,0.6),
c(0.3,1,0.1,0,0,0.1,0.1,0,0,0,0,0,0,0,0,0,0,0,0,0.3,0,0,0,0,0.2,0.2,0.3,0.3,0.5,0.5,0,0),
c(0.6,0.3,0.4,0,0,0.4,0.4,0.4,0.4,0.4,0.6,0,0.6,0,0,0.4,0,0.3,0.3,0.6,0,0.5,0.4,0.3,0.5,0,0.2,0,0,0.8,0.4,0.5),
c(0.3,0.2,0.2,0.2,0.4,0.2,0.3,0.3,0.3,0.4,0.5,0.5,0.6,0,0,1,0,0,0,0.6,0.3,1,0.2,0.2,0.3,0.7,1,0.4,0.9,0.9,0.5,0.5),
c(0.2,1,0,0,0,1,0,1,0.1,1,0.3,0.2,0.2,0,0,0,0,0,0,0,0,0,1,0.1,0.1,0.5,0,0.1,0.5,0.5,0,0),
c(0.1,0.1,0,0,0,0.1,0,1,0.1,1,0.3,0.1,0.2,0,0,0,0,0,0,0,0,0,1,0.1,0.1,0.5,0,0.1,0.5,0.4,0,0),
c(0.6,0.4,0.4,0.4,0.5,0.4,0.4,0.4,0.3,0.4,0.6,0.6,0.6,0,0,0.3,0,0,0,0.6,0.6,0.4,0.4,0.3,0.5,0,0.6,0.5,1,1,0.4,0.7),
c(0.6,0.3,0.3,0,0,0.3,0.3,0.4,0.2,0.3,0.6,0,0.6,0,0,0.3,0,0,0.3,0.6,0,0.4,0.3,0.2,0.4,0,0.5,0.5,0,1,0.3,0.3)))
print(sapply(testData, function(x) list( Permanent = perm(x))))
计算矩阵永久物是一个计算难题（例如，“计算永久物”有其专用的 Wikipedia article 与操作的数学定义分开......）
您实施的原始方法（已讨论 here）将需要计算长度为 32 的向量！ = 2.6e35（这将占用超过 10^24 TB 的内存，更不用说完成所有计算所花费的时间......）
幸运的是，自从上次在 SO 上讨论这个问题（2014 年）以来，有人已经将 BosonSampling 包发布到 CRAN，它使用编译代码（和一种高效的算法！）来完成计算。 tl;dr BosonSampling::rePerm(testData$A) 在我的笔记本电脑上不到 2 分钟就完成了。
设置
## repeated definition from above, for convenience
library(BosonSampling)
library(combinat)
perm <- function(A) {
stopifnot(is.matrix(A))
n <- nrow(A)
if(n != ncol(A)) stop("Matrix is not square.")
if(n < 1) stop("Matrix has a dimension of size 0.")
sum(sapply(combinat::permn(n), function(sigma) prod(sapply(1:n, function(i) A[i, sigma[i]]))))
}
检查小数据上简单和复杂解决方案的等价性
set.seed(101)
t1 <- matrix(rnorm(16), 4, 4)
t2 <- matrix(rnorm(36), 6, 6)
stopifnot(all.equal(perm(t1), rePerm(t1)))
stopifnot(all.equal(perm(t2), rePerm(t2)))
评估时间缩放
（在我的笔记本电脑上这需要几分钟）
f <- function(n, type = "naive") {
set.seed(101)
m <- matrix(rnorm(n^2), n, n)
t <- system.time(if (type == "naive") perm(m) else rePerm(m))
return(t[["elapsed"]])
}
n0 <- 3; n1 <- 10; n2 <- 32
ntimes <- sapply(n0:n1, f)
ptimes <- sapply(n0:n2, f, type = "rePerm")
dd <- data.frame(n = c(n0:n1, n0:n2),
time = c(ntimes, ptimes),
type = rep(c("naive", "rePerm"), c(n1-n0+1, n2-n0+1)))
saveRDS(dd, file = "SO72004136.rds")
library(ggplot2); theme_set(theme_bw())
(ggplot(dd) +
aes(n, time, colour = type) +
geom_point() + geom_line() +
scale_y_log10() +
scale_x_log10() +
expand_limits(x=32) +
labs(y = "elapsed time (sec)") +
geom_smooth(method = "lm",
lty = 2,
colour = "black",
data = subset(dd, type == "rePerm" & n>16))
)
lm(log(time) ~ log(n),
data = subset(dd, type == "rePerm" & n>16))
运行 n=10 的天真方法需要 190 秒； rePerm 对于 n=29 需要 18 秒。计算时间尺度大致为 time \propto n^17（即 log-log 回归的斜率约为 17），这非常糟糕...
一旦我发现这个问题可能需要几分钟而不是几个小时，我就在 115 秒内计算出 rePerm(testData$A) = 1.271369e+20。 （我真诚地希望你不必做比这更大的计算！）