素数幻方

mathe

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

1. #include <stdio.h>
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
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. }
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. }
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. }
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. }
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. }
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. }
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. }
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. }
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. }
159. char board[6][6];
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. }
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}