素数幻方

mathe

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

1. #include <stdio.h>
   
2. 
   
3. #define BUF_SIZE ((1<<20)>>5)
   
4. #define HALF_TABLE (1<<10)
   
5. #define GET_WORD(x) ((x)>>5)
   
6. #define GET_INDEX(x) ((x)&31)
   
7. #define SET_BIT(x)  pbuf[GET_WORD(x)]|=1<<GET_INDEX(x)
   
8. #define TEST_BIT(x) (pbuf[GET_WORD(x)]&(1<<GET_INDEX(x)))
   
9. int pbuf[BUF_SIZE];
   
10. #define IS_PRIME(x) (!TEST_BIT(x))
    
11. #define SET_NOT_PRIME(x) SET_BIT(x)
    
12. #define MCOUNT 68906
    
13. int p1[MCOUNT];//用于保存所有6位素数(要求逆向读也是素数)
    
14. int p2[MCOUNT];//用于保存所有可以用于4条边上的素数
    
15. int index1[100];///根据6位素数的最高位和最低位取值，分组
    
16. int index2[10];///根据最高位的值分组
    
17. int c1,c2;
    
18. int is_contain_zero(int x);
    
19. int last_column(int x);
    
20. int is_reverse_prime(int x);
    
21. int is_short_reverse(int x, int k);
    
22. int is_edge_prime(int x);
    
23. #define LOWER 100000
    
24. #define UPPER 999999
    
25. 
    
26. int cmp_p(const void *p, const void *q){
    
27.     const int *cp=(const int *)p;
    
28.     const int *cq=(const int *)q;
    
29.     int sp=*cp%10;
    
30.     int ep=*cp/LOWER;
    
31.     int sq=*cq%10;
    
32.     int eq=*cq/LOWER;
    
33.     if(ep<eq)return -1;
    
34.     if(ep>eq)return 1;
    
35.     if(sp<sq)return -1;
    
36.     if(sp>sq)return 1;
    
37.     return *cp-*cq;
    
38. }
    
39. 
    
40. void init_prime()
    
41. {
    
42.     int i,j,s,u;
    
43.     SET_NOT_PRIME(0);
    
44.     SET_NOT_PRIME(1);
    
45.     for(i=2;i<HALF_TABLE;i++){
    
46.         if(IS_PRIME(i)){
    
47.            for(j=i*i;j<=UPPER;j+=i){
    
48.                SET_NOT_PRIME(j);
    
49.            }
    
50.         }
    
51.     }
    
52.     s=1;u=10;
    
53.     for(i=1;i<=5;i++){
    
54.         for(j=s;j<u;j++){
    
55.             if(IS_PRIME(j)&&!is_short_reverse(j,i)){
    
56.                 SET_NOT_PRIME(j);
    
57.             }
    
58.         }
    
59.         s=u;u*=10;
    
60.     }
    
61. 
    
62.     c1=c2=0;
    
63.     for(i=LOWER;i<=UPPER;i++){
    
64.         if(IS_PRIME(i)&&is_reverse_prime(i)){
    
65.             p1[c1++]=i;
    
66.             if(is_edge_prime(i)){
    
67.                 p2[c2++]=i;
    
68.             }
    
69.         }else{
    
70.             SET_NOT_PRIME(i);
    
71.         }
    
72.     }
    
73. 
    
74.     qsort(p1,c1,sizeof(int),cmp_p);
    
75.     for(i=0,j=1,index2[0]=0;i<c2;i++){
    
76.         if(p2[i]>=(j+1)*LOWER){
    
77.             index2[j]=i;
    
78.             j++;
    
79.         }
    
80.         while(p2[i]>=(j+1)*LOWER){
    
81.             index2[j]=i;
    
82.             j++;
    
83.         }
    
84.     }
    
85.     for(;j<10;j++)index2[j]=c2;
    
86.     for(i=0,j=1,index1[0]=0;i<c1;i++){
    
87.         int high=j/10,lower=j%10;
    
88.         int p=p1[i];
    
89.         int ph=p/LOWER,pl=p%10;
    
90.         if(ph!=high||pl!=lower){
    
91.             int next_j=ph*10+pl;
    
92.             for(;j<next_j;j++){
    
93.                 index1[j]=i;
    
94.             }
    
95.         }
    
96.     }
    
97.     for(;j<100;j++)index1[j]=c1;
    
98.     fprintf(stderr,"c1=%d,c2=%d\n",c1,c2);
    
99. }
    
100. 
     
101. int is_short_reverse(int x, int k)
     
102. {
     
103.     int i,y;
     
104.     y=0;
     
105.     for(i=0;i<k;i++){
     
106.         y*=10;
     
107.         y+=x%10;
     
108.         x/=10;
     
109.     }
     
110.     return IS_PRIME(y);
     
111. }
     
112. 
     
113. int is_reverse_prime(int x)
     
114. {
     
115.     int i,y;
     
116.     y=0;
     
117.     for(i=0;i<6;i++){
     
118.         y*=10;
     
119.         y+=x%10;
     
120.         x/=10;
     
121.     }
     
122.     return IS_PRIME(y);
     
123. }
     
124. 
     
125. int is_contain_zero(int x)
     
126. {
     
127.     int i;
     
128.     for(i=0;i<6;i++){
     
129.         if(x%10==0)return 1;
     
130.         x/=10;
     
131.     }
     
132.     return 0;
     
133. }
     
134. 
     
135. int is_edge_prime(int x)
     
136. {
     
137.     int last=x%10;
     
138.     if(!IS_PRIME(last))
     
139.         return 0;
     
140.     if(!IS_PRIME(x/100000))
     
141.         return 0;
     
142.     if(!last_column(x))
     
143.         return 0;
     
144.     return !is_contain_zero(x);
     
145. }
     
146. 
     
147. int last_column(int x)
     
148. {
     
149.     int valid[10]={0,1,0,1,0,0,0,1,0,1};
     
150.     int i;
     
151.     for(i=0;i<6;i++){
     
152.         int d=x%10;
     
153.         if(!valid[d])return 0;
     
154.         x/=10;
     
155.     }
     
156.     return 1;
     
157. }
     
158. 
     
159. char board[6][6];
     
160. 
     
161. void check_row()
     
162. {
     
163.     int i,j,s;
     
164.     for(i=1;i<6;i++){
     
165.         int x=0;
     
166.         for(j=0;j<6;j++){
     
167.             x*=10;
     
168.             x+=board[i][5-j];
     
169.         }
     
170.         if(!IS_PRIME(x))
     
171.             return;
     
172.     }
     
173.     s=0;
     
174.     for(i=0;i<6;i++){
     
175.         s*=10;
     
176.         s+=board[i][i];
     
177.     }
     
178.     if(!IS_PRIME(s))
     
179.         return;
     
180.     s=0;
     
181.     for(i=0;i<6;i++){
     
182.         s*=10;
     
183.         s+=board[i][5-i];
     
184.     }
     
185.     if(!IS_PRIME(s))
     
186.         return;
     
187.     for(i=0;i<6;i++){///first check s+t=i
     
188.         int sp=0;
     
189.         for(j=0;j<=i;j++){
     
190.             sp*=10;
     
191.             sp+=board[j][i-j];
     
192.         }
     
193.         if(!IS_PRIME(sp))
     
194.             return;
     
195.         sp=0;
     
196.         for(j=i;j<=5;j++){///check s+t = 5+i
     
197.             sp*=10;
     
198.             sp+=board[j][5+i-j];
     
199.         }
     
200.         if(!IS_PRIME(sp))
     
201.             return;
     
202.         sp=0;
     
203.         for(j=0;j<=5-i;j++){///check s-t=i;
     
204.             sp*=10;
     
205.             sp+=board[i+j][j];
     
206.         }
     
207.         if(!IS_PRIME(sp))
     
208.             return;
     
209.         sp=0;
     
210.         for(j=0;j<=5-i;j++){///Check s-t=-i;
     
211.             sp*=10;
     
212.             sp+=board[j][i+j];
     
213.         }
     
214.         if(!IS_PRIME(sp))
     
215.             return;
     
216.     }
     
217.     for(i=0;i<6;i++){
     
218.         for(j=0;j<6;j++){
     
219.             printf("%d",board[i][j]);
     
220.         }
     
221.         printf("\n");
     
222.     }
     
223.     printf("\n");
     
224.     fflush(stdout);
     
225. }
     
226. void search()
     
227. {
     
228.     int p,q;
     
229.     int pp,qq,ss,tt;
     
230.     int i,j;
     
231.     for(pp=0;pp<c2;pp++){
     
232.         int lead;
     
233.         int s;
     
234.         p=p2[pp];
     
235.         lead=p/LOWER;
     
236.         s=p;
     
237.         for(i=0;i<6;i++){
     
238.             board[0][5-i]=s%10;
     
239.             s/=10;
     
240.         }
     
241.         for(qq=index2[lead-1];qq<index2[lead];qq++){
     
242.             int j1,j2,j3,j4;
     
243.             int h1,h2,h3,h4;
     
244.             int sp;
     
245.             q=p2[qq];
     
246.             s=q;
     
247.             for(j=0;j<6;j++){
     
248.                 board[5-j][0]=s%10;
     
249.                 s/=10;
     
250.             }
     
251.             sp=10*board[1][0]+board[0][1];
     
252.             if(!IS_PRIME(sp))
     
253.                 continue;
     
254.             fprintf(stderr,"begin with edge <%d,%d,...>\n",p2[pp],p2[qq]);
     
255.             for(ss=index2[board[0][5]-1];ss<index2[board[0][5]];ss++){
     
256.                 s=p2[ss];
     
257.                 for(i=0;i<6;i++){
     
258.                     board[5-i][5]=s%10;
     
259.                     s/=10;
     
260.                 }
     
261.                 sp=10*board[0][4]+board[1][5];
     
262.                 if(!IS_PRIME(sp))
     
263.                     continue;
     
264.                 for(tt=index2[board[5][0]-1];tt<index2[board[5][0]];tt++){
     
265.                     s=p2[tt];
     
266.                     if(s%10!=board[5][5])
     
267.                         continue;
     
268.                     for(i=0;i<6;i++){
     
269.                         board[5][5-i]=s%10;
     
270.                         s/=10;
     
271.                     }
     
272.                     sp =10*board[4][0]+board[5][1];
     
273.                     if(!IS_PRIME(sp))
     
274.                         continue;
     
275.                     sp = 10*board[4][5]+board[5][4];
     
276.                     if(!IS_PRIME(sp))
     
277.                         continue;
     
278.                     h1=board[0][1]*10+board[5][1];
     
279.                     h2=board[0][2]*10+board[5][2];
     
280.                     h3=board[0][3]*10+board[5][3];
     
281.                     h4=board[0][4]*10+board[5][4];
     
282.                     for(j1=index1[h1-1];j1<index1[h1];j1++){
     
283.                         s=p1[j1];
     
284.                         for(j=0;j<5;j++){
     
285.                             board[5-j][1]=s%10;
     
286.                             s/=10;
     
287.                         }
     
288.                         sp=board[2][0]*100+board[1][1]*10+board[0][2];
     
289.                         if(!IS_PRIME(sp))
     
290.                             continue;
     
291.                         sp=board[3][0]*100+board[4][1]*10+board[5][2];
     
292.                         if(!IS_PRIME(sp))
     
293.                             continue;
     
294.                         for(j2=index1[h2-1];j2<index1[h2];j2++){
     
295.                             s=p1[j2];
     
296.                             for(j=0;j<5;j++){
     
297.                                 board[5-j][2]=s%10;
     
298.                                 s/=10;
     
299.                             }
     
300.                             sp=board[3][0]*1000+board[2][1]*100+board[1][2]*10+board[0][3];
     
301.                             if(!IS_PRIME(sp))
     
302.                                 continue;
     
303.                             sp=board[2][0]*1000+board[3][1]*100+board[4][2]*10+board[5][3];
     
304.                             if(!IS_PRIME(sp))
     
305.                                 continue;
     
306.                             for(j3=index1[h3-1];j3<index1[h3];j3++){
     
307.                                 s=p1[j3];
     
308.                                 for(j=0;j<5;j++){
     
309.                                     board[5-j][3]=s%10;
     
310.                                     s/=10;
     
311.                                 }
     
312.                                 sp=board[4][0]*10000+board[3][1]*1000+board[2][2]*100+board[1][3]*10+board[0][4];
     
313.                                 if(!IS_PRIME(sp))
     
314.                                     continue;
     
315.                                 sp=board[1][0]*10000+board[2][1]*1000+board[3][2]*100+board[4][3]*10+board[5][4];
     
316.                                 if(!IS_PRIME(sp))
     
317.                                     continue;
     
318.                                 for(j4=index1[h4-1];j4<index1[h4];j4++){
     
319.                                     s=p1[j4];
     
320.                                     for(j=0;j<5;j++){
     
321.                                         board[5-j][4]=s%10;
     
322.                                         s/=10;
     
323.                                     }
     
324.                                     check_row();
     
325.                                 }
     
326.                             }
     
327.                         }
     
328.                     }
     
329.                 }
     
330.             }
     
331.         }
     
332.     }
     
333. }
     
334. 
     
335. int main()
     
336. {
     
337.     init_prime();
     
338.     search();
     
339. }

复制代码

mathe

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

无心人

那就去掉吧

要不来个低阶的看下？

mathe

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

mathe

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

无心人

那就满世界发帖子求一个超完美的
限定在16阶内

估计只要边界数字限定在1， 3， 7， 9就可以

葡萄糖

找到一张图片
有趣的是，其定和为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}