R语言：素数幻方

求N阶的素数幻方。即在一个NXN 的矩阵中，每一个格填入一个数字，使每一行、每一列和两条对角线上的N 个数字所构成的N位数，均为可逆素数。

说明：对付高出3阶的幻方，对内存是一个检验。其它语言有没有更简捷的要领，递归就免了吧。

solve.cube <- function(n) {
stopifnot(n > 1)
max.n <- 10 ^ n – 1 # 较大值 99…9 (n个9)
min.n <- max.n / 9 # 最小值 11…1 (n个1)，因为幻方边上不行能存在0
prime <- rep(T, max.n – 1) # 初始化数组，以便用删除法求素数
for (i in 2:sqrt(max.n)) if (prime[i]) prime[seq(i * 2, max.n, by = i)] <- F
prime <- which(prime)
prime <- prime[prime >= min.n] # 获得 min.n 到 max.n 之间的所有素数

# 取个中的可逆素数
rev.prime <- as.integer(sapply(as.character(prime),
function(s) rawToChar(rev(charToRaw(s)))))
prime <- prime[prime %in% rev.prime]

# 拆分数位
prime.digits <- sapply(as.character(prime), function(s) as.integer(charToRaw(s)) – 48)

# 判定函数（因为列是从可逆素数中取，所以只判定行和对角线）
matrix(1:(n^2), n) -> mi
mi <- rbind(mi, diag(mi), diag(mi[n:1,]))
mf <- matrix(10^(n:1-1), n + 2, n, byrow=T)
check_matrix <- function(m) {
all(rowSums(m[mi] * mf) %in% prime)
}

# 实验开始组合，取出幻方功效
combs <- expand.grid(as.data.frame(matrix(seq_along(prime), length(prime), n)))
idx <-