看到一个月前版主提到素数幻方矩阵。
参考资料
http://cos.name/cn/topic/107777#post-236019
突发奇想在之前版主的算法上加上一些幻方性质来加速幻方计算。
本人计算机配置如下图
单线程能体现算法性能，故所有运算均以单线程计算
目前：
单线程4阶幻方计算内存消耗200M以内，计算时间2秒左右。
单线程5阶幻方计算内存消耗1300M以内，计算时间50小时左右,根据我的计算方法，5阶运算被分为10000份，跑完1份的时间20秒左右，5阶幻方计算需要约3000分钟左右。理论上使用我的小破本50小时，可算出所有5阶幻方，还可继续优化。
为方便说明算法：
定义:分族,将可逆素数根据边缘进行分类，取首数字与末数字作为数字特征。例如：3319，取边缘值39，作为特征。根据首尾数字只有可能是1，3，7，9。故共有16类族群，全体可逆素数一定属于16族之一。
根据素数幻方矩阵边缘的矩阵的性质，当确定上下两行元素的组合后，4列数据分别属于的族就确定了。
此时我们只需要在对应族群里选取变量组合即可。
同理，左右两列确定后，对行的族选择也有限制作用。
由于我现在只使用上下两行+第一列制作分族运算(3边缘分族)。
对于更高阶幻方可以使用上下两列+左右两列进行分族运算，以减少内存消耗。
本文只使用3个边缘进行分族，计算效率已经飞跃性提升。
因时间有限4边缘分族，给出计算复杂度估计：
空间占用：
4阶幻方内存占用到50M左右。
5阶幻方内存占用到130M左右。
6阶幻方内存占用到300M左右。
原理：
1、通过第一行、最后一行分组，限制候选数集。
2、由于素数性质，第一行、最后一行、第一列、最后一列一定由1、3、7、9组成
3、利用数组查询代替in
4、使用矩阵运算代替apply运算
5、去除一些无意义的计时运算
3阶幻方计算时间如下：很慢(由于算法为计算高阶幻方，加入了大量逻辑运算，降低了低阶幻方运算速度)
4阶幻方计算时间如下，很快吧，哈哈哈
图片看不到这里也有
http://blog.sina.com.cn/s/blog_3f16beff0102v3bv.html
<br />
solve.cube <- function(n) {<br />
stopifnot(n > 1)<br />
max.n <- 10 ^ n - 1 # 最大值 99...9 (n个9)<br />
min.n <- max.n / 9 # 最小值 11...1 (n个1)，因为幻方边上不可能存在0<br />
prime <- rep(T, max.n - 1) # 初始化数组，以便用删除法求素数<br />
for (i in 2:sqrt(max.n)) {<br />
if (prime[i]) {<br />
prime[seq(i * 2, max.n, by = i)] <- F<br />
}<br />
}<br />
prime <- which(prime)<br />
prime <- prime[prime >= min.n] # 得到 min.n 到 max.n 之间的所有素数</p>
<p> # 取其中的可逆素数<br />
rev.prime <- as.integer(sapply(as.character(prime),<br />
function(s) rawToChar(rev(charToRaw(s)))))<br />
prime <- prime[prime %in% rev.prime]<br />
# 拆分数位<br />
prime.digits <- sapply(as.character(prime), function(s) as.integer(charToRaw(s)) - 48)<br />
# 标记特征头尾<br />
prime.r.r <- c(1:ncol(prime.digits))[apply(matrix(apply(matrix(prime.digits[2:(n-1),],n-2),2,function(v) v%in%c(1,3,7,9)),n-2),2,all)]<br />
#变量分族<br />
prime.class<-apply(prime.digits,2,function(v) v[1]*10+v[n])<br />
#prime.r<-as.matrix(t(prime.digits[,prime.digits.r]))%*%mf<br />
a<-matrix(0,10^n)<br />
a[prime]=1<br />
# 判断函数（因为列是从可逆素数中取，所以只判断行和对角线）<br />
matrix(1:(n^2), n) -> mi<br />
mi <- diag(mi[n:1,])<br />
mf <- matrix(10^(n:1-1))<br />
check_matrix <- function(m) {<br />
#通过数组索引识别一般来说效率高于in<br />
all(a[m%*%mf]==1)<br />
#a[m[mi[i,]]%*%mf]==1<br />
}<br />
check_matrix_d <- function(m) {<br />
all(a[matrix(m[mi],ncol=n)%*%mf]==1,a[diag(m)%*%mf]==1)<br />
#a[m[mi[i,]]%*%mf]==1<br />
}<br />
val_class <- function(m){<br />
c(1:ncol(prime.digits))[prime.class==lim.range[m]]<br />
}<br />
#构建首尾约束<br />
prime.lim<-expand.grid(prime.r.r,prime.r.r)<br />
#构建候选集<br />
prime.list<-vector("list",n)<br />
prime.list[[n]]<-prime.r.r<br />
for(i in 2:n-1){<br />
prime.list[[i]]<-seq_along(prime)<br />
}<br />
lim.range<-matrix(0,1,n)<br />
flag2=1<br />
for(k in 1:length(prime.r.r)^2){<br />
prime.lim[k,]<br />
lim.up<-prime.digits[,prime.lim[k,][[1]]]<br />
lim.down<-prime.digits[,prime.lim[k,][[2]]]<br />
#构建候选集计算族<br />
lim.range<-lim.up*10+lim.down<br />
#以下表达式计算n族位置<br />
#c(1:ncol(prime.digits))[prime.class==lim.range[n]]<br />
prime.list<-vector("list",n)<br />
prime.list[[n]]<-prime.r.r[prime.r.r%in%val_class(n)]<br />
flag=1<br />
if(length(prime.list[[n]])==0){<br />
flag=0<br />
}<br />
for(i in 2:n-1){<br />
prime.list[[i]]<-val_class(i)<br />
if(length(prime.list[[i]])==0){<br />
flag=0<br />
}<br />
}<br />
if(flag==0||length(prime.r.r[prime.r.r%in%val_class(1)])==0){<br />
next<br />
}</p>
<p> for(i in 1:length(prime.r.r[prime.r.r%in%val_class(1)])){<br />
#合成前可以先过滤函数<br />
prime.list[[1]] <- prime.r.r[prime.r.r%in%val_class(1)][i]<br />
combs <- do.call("cbind",expand.grid(prime.list))<br />
#4层就70+万个候选矩阵了,理论上生成过程可以剪枝，不过那就要自己写函数了。<br />
#行检验可以矩阵运算得到计算速度较快，消耗时间比例（全行检验,安位置抽取对角线，diag函数--1:2:10）<br />
#经测试对角线、副对角线，diag函数耗时极高只能按位置抽取位置抽取耗时较高，故先通过行检验减少数据<br />
#apply函数效率太低估替换为矩阵运算<br />
#idx <- apply(combs, 1, function(v) check_matrix(prime.digits[,v]))<br />
tmp<-combs<br />
if(nrow(tmp)<1){<br />
next<br />
}<br />
for(j in 2:n-1){<br />
m<-matrix(prime.digits[j,tmp],ncol=n)<br />
idx<-a[m%*%mf]==1<br />
if(all(!idx)){<br />
break<br />
}<br />
tmp<-matrix(tmp[idx,],ncol=n)<br />
}<br />
if(all(!idx)){<br />
next<br />
}<br />
if(length(dim(tmp))<=0){<br />
next<br />
}<br />
#nrow(tmp)#筛选后只有500+个<br />
#开始对角线检验<br />
#prime.digits[,tmp[1,]]<br />
#idxt <- all(a[matrix(m[mi],ncol=n)%*%mf]==1,a[diag(m)%*%mf]==1)<br />
idxt <- apply(tmp, 1, function(v) check_matrix_d(prime.digits[,v]))<br />
if(all(!idxt)){<br />
next<br />
}<br />
#tmp<-tmp[idxt,]<br />
#if(length(tmp)==0){<br />
# next<br />
#}<br />
if(flag2==1){<br />
res<-tmp[idxt,]<br />
flag2=0<br />
}else{<br />
res=rbind(res,tmp[idxt,])<br />
}<br />
}<br />
}<br />
results<-sapply(apply(res, 1, function(v) list(prime.digits[,v])), function(x) x)<br />
}<br />
</p>
复制过来代码格式都乱了，为了编辑代码格式，摁空格手指都点抽筋。。。
