求3阶的素数幻方。即在一个3X3 的矩阵中，每一个格填入一个数字，使每一行、每一列和两条对角线上的3 个数字所组成的3位数，均为可逆素数。
下面的代码比较繁琐，而且推广到4阶以上就很困难，因为计算机内存的原因。有没有改进的余地。
<br />
# 判断是否为素数<br />
is.prime <- function(m) ifelse(m==2,TRUE,sum( m%%c(2:(m-1)) ==0 ) == 0)<br />
# 求素数<br />
prime<-function(m){<br />
x<-c(2:m) ; y<-NULL<br />
repeat{<br />
z<-x[(x%%x[1])!=0] ; y<-c(y,x[1])<br />
if(x[1]>sqrt(m))break<br />
x<-z<br />
}<br />
c(y,z)<br />
}<br />
a <- prime(999)[prime(999) > 100] #计算三位素数<br />
b = data.frame(i=(a%/%100)%%10,j=(a%/%10)%%10,k=a%%10)<br />
d = transform(b,x=a,y=as.numeric(paste(b$k,b$j,b$i,sep='')))<br />
#可逆素数<br />
e <- d[sapply(d$y,is.prime),]<br />
m <- e$x<br />
n <- combn(m,3) #三位素数组合<br />
#分解数字<br />
fa <- function(x) as.integer(unlist(strsplit(as.character(x),NULL)))<br />
#合并数字<br />
fb <- function(x) x[1]*100+x[2]*10+x[3]<br />
#判断各列对角线是否为素数<br />
fc <- function(x) {<br />
ma <- matrix(fa(n[,x]),3,3,byrow=T)<br />
x1 <- fb(diag(ma[nrow(ma):1, ])) #副对角线<br />
x2 <- fb(diag(ma)) #主对角线<br />
x3 <- fb(t(ma)[1,]) #第一列<br />
x4 <- fb(t(ma)[2,]) #第二列<br />
x5 <- fb(t(ma)[3,]) #第三列<br />
xx <- c(x1,x2,x3,x4,x5)<br />
if (sum(xx %in% m) == 5)<br />
return(n[,x])<br />
}<br />
dd <- lapply(1:ncol(n),function(i) fc(i))<br />
do.call(rbind,dd)<br />
[,1] [,2] [,3]<br />
[1,] 113 151 311<br />
[2,] 113 181 311<br />
[3,] 113 751 971<br />
[4,] 191 797 919<br />
[5,] 337 353 733<br />
[6,] 337 353 739<br />
[7,] 337 383 733<br />
[8,] 337 383 739<br />
</p>
