看到一个月前版主提到素数幻方矩阵。

参考资料

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>

复制过来代码格式都乱了，为了编辑代码格式，摁空格手指都点抽筋。。。