素数幻方

mathe · 发表于 2008-4-16 17:45:20

下面程序在我的计算机上运行了将近一个小时，没有找到任何结果。如果程序没有错误（很有可能哪里不小心弄错了一点），那么超完美的就是无解了。

#include <stdio.h>
#define BUF_SIZE ((1<<20)>>5)
#define HALF_TABLE (1<<10)
#define GET_WORD(x) ((x)>>5)
#define GET_INDEX(x) ((x)&31)
#define SET_BIT(x) pbuf[GET_WORD(x)]|=1<<GET_INDEX(x)
#define TEST_BIT(x) (pbuf[GET_WORD(x)]&(1<<GET_INDEX(x)))
int pbuf[BUF_SIZE];
#define IS_PRIME(x) (!TEST_BIT(x))
#define SET_NOT_PRIME(x) SET_BIT(x)
#define MCOUNT 68906
int p1[MCOUNT];//用于保存所有6位素数(要求逆向读也是素数)
int p2[MCOUNT];//用于保存所有可以用于4条边上的素数
int index1[100];///根据6位素数的最高位和最低位取值，分组
int index2[10];///根据最高位的值分组
int c1,c2;
int is_contain_zero(int x);
int last_column(int x);
int is_reverse_prime(int x);
int is_short_reverse(int x, int k);
int is_edge_prime(int x);
#define LOWER 100000
#define UPPER 999999
int cmp_p(const void *p, const void *q){
const int *cp=(const int *)p;
const int *cq=(const int *)q;
int sp=*cp%10;
int ep=*cp/LOWER;
int sq=*cq%10;
int eq=*cq/LOWER;
if(ep<eq)return -1;
if(ep>eq)return 1;
if(sp<sq)return -1;
if(sp>sq)return 1;
return *cp-*cq;
}
void init_prime()
{
int i,j,s,u;
SET_NOT_PRIME(0);
SET_NOT_PRIME(1);
for(i=2;i<HALF_TABLE;i++){
if(IS_PRIME(i)){
for(j=i*i;j<=UPPER;j+=i){
SET_NOT_PRIME(j);
}
}
}
s=1;u=10;
for(i=1;i<=5;i++){
for(j=s;j<u;j++){
if(IS_PRIME(j)&&!is_short_reverse(j,i)){
SET_NOT_PRIME(j);
}
}
s=u;u*=10;
}
c1=c2=0;
for(i=LOWER;i<=UPPER;i++){
if(IS_PRIME(i)&&is_reverse_prime(i)){
p1[c1++]=i;
if(is_edge_prime(i)){
p2[c2++]=i;
}
}else{
SET_NOT_PRIME(i);
}
}
qsort(p1,c1,sizeof(int),cmp_p);
for(i=0,j=1,index2[0]=0;i<c2;i++){
if(p2[i]>=(j+1)*LOWER){
index2[j]=i;
j++;
}
while(p2[i]>=(j+1)*LOWER){
index2[j]=i;
j++;
}
}
for(;j<10;j++)index2[j]=c2;
for(i=0,j=1,index1[0]=0;i<c1;i++){
int high=j/10,lower=j%10;
int p=p1[i];
int ph=p/LOWER,pl=p%10;
if(ph!=high||pl!=lower){
int next_j=ph*10+pl;
for(;j<next_j;j++){
index1[j]=i;
}
}
}
for(;j<100;j++)index1[j]=c1;
fprintf(stderr,"c1=%d,c2=%d\n",c1,c2);
}
int is_short_reverse(int x, int k)
{
int i,y;
y=0;
for(i=0;i<k;i++){
y*=10;
y+=x%10;
x/=10;
}
return IS_PRIME(y);
}
int is_reverse_prime(int x)
{
int i,y;
y=0;
for(i=0;i<6;i++){
y*=10;
y+=x%10;
x/=10;
}
return IS_PRIME(y);
}
int is_contain_zero(int x)
{
int i;
for(i=0;i<6;i++){
if(x%10==0)return 1;
x/=10;
}
return 0;
}
int is_edge_prime(int x)
{
int last=x%10;
if(!IS_PRIME(last))
return 0;
if(!IS_PRIME(x/100000))
return 0;
if(!last_column(x))
return 0;
return !is_contain_zero(x);
}
int last_column(int x)
{
int valid[10]={0,1,0,1,0,0,0,1,0,1};
int i;
for(i=0;i<6;i++){
int d=x%10;
if(!valid[d])return 0;
x/=10;
}
return 1;
}
char board[6][6];
void check_row()
{
int i,j,s;
for(i=1;i<6;i++){
int x=0;
for(j=0;j<6;j++){
x*=10;
x+=board[i][5-j];
}
if(!IS_PRIME(x))
return;
}
s=0;
for(i=0;i<6;i++){
s*=10;
s+=board[i][i];
}
if(!IS_PRIME(s))
return;
s=0;
for(i=0;i<6;i++){
s*=10;
s+=board[i][5-i];
}
if(!IS_PRIME(s))
return;
for(i=0;i<6;i++){///first check s+t=i
int sp=0;
for(j=0;j<=i;j++){
sp*=10;
sp+=board[j][i-j];
}
if(!IS_PRIME(sp))
return;
sp=0;
for(j=i;j<=5;j++){///check s+t = 5+i
sp*=10;
sp+=board[j][5+i-j];
}
if(!IS_PRIME(sp))
return;
sp=0;
for(j=0;j<=5-i;j++){///check s-t=i;
sp*=10;
sp+=board[i+j][j];
}
if(!IS_PRIME(sp))
return;
sp=0;
for(j=0;j<=5-i;j++){///Check s-t=-i;
sp*=10;
sp+=board[j][i+j];
}
if(!IS_PRIME(sp))
return;
}
for(i=0;i<6;i++){
for(j=0;j<6;j++){
printf("%d",board[i][j]);
}
printf("\n");
}
printf("\n");
fflush(stdout);
}
void search()
{
int p,q;
int pp,qq,ss,tt;
int i,j;
for(pp=0;pp<c2;pp++){
int lead;
int s;
p=p2[pp];
lead=p/LOWER;
s=p;
for(i=0;i<6;i++){
board[0][5-i]=s%10;
s/=10;
}
for(qq=index2[lead-1];qq<index2[lead];qq++){
int j1,j2,j3,j4;
int h1,h2,h3,h4;
int sp;
q=p2[qq];
s=q;
for(j=0;j<6;j++){
board[5-j][0]=s%10;
s/=10;
}
sp=10*board[1][0]+board[0][1];
if(!IS_PRIME(sp))
continue;
fprintf(stderr,"begin with edge <%d,%d,...>\n",p2[pp],p2[qq]);
for(ss=index2[board[0][5]-1];ss<index2[board[0][5]];ss++){
s=p2[ss];
for(i=0;i<6;i++){
board[5-i][5]=s%10;
s/=10;
}
sp=10*board[0][4]+board[1][5];
if(!IS_PRIME(sp))
continue;
for(tt=index2[board[5][0]-1];tt<index2[board[5][0]];tt++){
s=p2[tt];
if(s%10!=board[5][5])
continue;
for(i=0;i<6;i++){
board[5][5-i]=s%10;
s/=10;
}
sp =10*board[4][0]+board[5][1];
if(!IS_PRIME(sp))
continue;
sp = 10*board[4][5]+board[5][4];
if(!IS_PRIME(sp))
continue;
h1=board[0][1]*10+board[5][1];
h2=board[0][2]*10+board[5][2];
h3=board[0][3]*10+board[5][3];
h4=board[0][4]*10+board[5][4];
for(j1=index1[h1-1];j1<index1[h1];j1++){
s=p1[j1];
for(j=0;j<5;j++){
board[5-j][1]=s%10;
s/=10;
}
sp=board[2][0]*100+board[1][1]*10+board[0][2];
if(!IS_PRIME(sp))
continue;
sp=board[3][0]*100+board[4][1]*10+board[5][2];
if(!IS_PRIME(sp))
continue;
for(j2=index1[h2-1];j2<index1[h2];j2++){
s=p1[j2];
for(j=0;j<5;j++){
board[5-j][2]=s%10;
s/=10;
}
sp=board[3][0]*1000+board[2][1]*100+board[1][2]*10+board[0][3];
if(!IS_PRIME(sp))
continue;
sp=board[2][0]*1000+board[3][1]*100+board[4][2]*10+board[5][3];
if(!IS_PRIME(sp))
continue;
for(j3=index1[h3-1];j3<index1[h3];j3++){
s=p1[j3];
for(j=0;j<5;j++){
board[5-j][3]=s%10;
s/=10;
}
sp=board[4][0]*10000+board[3][1]*1000+board[2][2]*100+board[1][3]*10+board[0][4];
if(!IS_PRIME(sp))
continue;
sp=board[1][0]*10000+board[2][1]*1000+board[3][2]*100+board[4][3]*10+board[5][4];
if(!IS_PRIME(sp))
continue;
for(j4=index1[h4-1];j4<index1[h4];j4++){
s=p1[j4];
for(j=0;j<5;j++){
board[5-j][4]=s%10;
s/=10;
}
check_row();
}
}
}
}
}
}
}
}
}
int main()
{
init_prime();
search();
}

复制代码

mathe · 发表于 2008-4-16 17:49:33

如果将最外层的4个循环结果事先处理一下，去掉对称情况，应该可以提高将近8倍的速度

无心人 · 发表于 2008-4-16 17:54:31

那就去掉吧要不来个低阶的看下？

mathe · 发表于 2008-4-16 18:06:06

呵呵，最好两阶，这个显然有解: $[(3,7),(7,3)]$

mathe · 发表于 2008-4-16 18:06:40

哦，两阶都不行呀对角线不行

无心人 · 发表于 2008-4-16 20:23:44

那就满世界发帖子求一个超完美的限定在16阶内估计只要边界数字限定在1， 3， 7， 9就可以

葡萄糖 · 发表于 2019-3-17 22:27:07

找到一张图片
有趣的是，其定和为666

\begin{array}{|r|r|r|r|}
\hline 3 & 107 & 5 & 131 &109 &311 \\
\hline 7 & 331 & 193 & 11 &83 &41 \\
\hline 103 & 53 & 71 & 89 &151 &199 \\
\hline 113 & 61 & 97 & 197 &167 &31 \\
\hline 367 & 13 & 173 & 59 &17 &37 \\
\hline 73 & 101 & 127 & 179 &139 &47 \\
\hline \end{array}

账号		自动登录	找回密码
密码			欢迎注册

[擂台] 素数幻方