求能计算所有方阵的程序

libinary

改进。
因为N行和N列可以用SUM减前N-1项求得，所以我们只要穷举(N - 1) * (N - 1)的方阵就可以了，
代码增加了一个gen函数，改变了func函数，
因为行和列的和已经保证等于SUM，所以testsum函数中去掉了行和列的检查。
另：因为我的程序是用define定义的N，所以改方阵阶数的时候要同时改
#define N 3
#define N2 9
#define SUM 15
三行。

1. 
   
2. #include <stdio.h>
   
3. 
   
4. #define N 3
   
5. #define N2 9
   
6. #define SUM 15
   
7. #define FALSE 0
   
8. #define TRUE 1
   
9. #define WIDTH 2
   
10. 
    
11. int testsum(void);
    
12. void print_a(void);
    
13. void gen(void);
    
14. void func(int);
    
15. 
    
16. int a[N2];
    
17. int num[N2];
    
18. char line[(WIDTH + 1) * N + 1];
    
19. int count;
    
20. 
    
21. int main(void)
    
22. {
    
23.   int i;
    
24. 
    
25.   memset(line, '-', (WIDTH + 1) * N);
    
26.   line[(WIDTH + 1) * N] = '\0';
    
27.   for(i = 0; i < N2; ++i)
    
28.     num[i] = 0;
    
29. 
    
30.   func(0);
    
31.   printf("count: %d\n", count);
    
32. 
    
33.   exit(0);
    
34. }
    
35. 
    
36. int
    
37. testsum(void)
    
38. {
    
39.   int x, y, sum;
    
40. 
    
41.   sum = 0;
    
42.   for(x = 0, y = 0; x < N; ++x, ++y)
    
43.     sum += *(a + y * N + x);
    
44.   if(sum != SUM)
    
45.     return(FALSE);
    
46.   sum = 0;
    
47.   for(x = N - 1, y = 0; y < N; --x, ++y)
    
48.     sum += *(a + y * N + x);
    
49.   if(sum != SUM)
    
50.     return(FALSE);
    
51. 
    
52.   ++count;
    
53.   return(TRUE);
    
54. }
    
55. 
    
56. void
    
57. print_a(void)
    
58. {
    
59.   int x, y, i;
    
60. 
    
61.   for(y = 0; y < N; ++y){
    
62.     for(x = 0; x < N; ++x)
    
63.       printf("%*d ", WIDTH, *(a + y * N + x));
    
64.     printf("\n");
    
65.   }
    
66.   printf("%s\n", line);
    
67. }
    
68. 
    
69. void
    
70. gen(void)
    
71. {
    
72.   int al[N2];
    
73.   int i, j, sum, val;
    
74. 
    
75.   memcpy(al, num, N2 * sizeof(num[0]));
    
76.   for(i = N2 - N; i < N2 - 1; ++i){
    
77.     sum = 0;
    
78.     for(j = i % N; j < i; j += N)
    
79.       sum += a[j];
    
80.     val = SUM - sum;
    
81.     if(val < 1 || val > N2 || al[val - 1] == 1)
    
82.       return;
    
83.     else{
    
84.       al[val - 1] = 1;
    
85.       a[i] = val;
    
86.     }
    
87.   }
    
88.   for(i = N - 1; i < N2; i += N){
    
89.     sum = 0;
    
90.     for(j = i - N + 1; j < i; ++j)
    
91.       sum += a[j];
    
92.     val = SUM - sum;
    
93.     if(val < 1 || val > N2 || al[val - 1] == 1)
    
94.       return;
    
95.     else{
    
96.       al[val - 1] = 1;
    
97.       a[i] = val;
    
98.     }
    
99.   }
    
100.   if(testsum())
     
101.     print_a();
     
102. }
     
103. 
     
104. void
     
105. func(int n)
     
106. {
     
107.   int i;
     
108. 
     
109.   for(i = 0; i < N2; ++i){
     
110.     if(num[i] == 0){
     
111.       num[i] = 1;
     
112.       a[n] = i + 1;
     
113.       if(n == N2 - N - 2)
     
114.         gen();
     
115.       else if(n % N == N - 2)
     
116.         func(n + 2);
     
117.       else
     
118.         func(n + 1);
     
119.       num[i] = 0;
     
120.     }
     
121.   }
     
122. }
     

复制代码

time：
N=3：0.010 (比完全穷举快了26倍)
N=4：39m25s 7040个结果

to：lordbyorn
我运行过你的程序，N=4的时候好象没有30s就出来吧，是不是30m？
我运行程序的时候是在X下，开了十几个窗口，还有人在我的ftp上传一个电影，所以有点慢，呵呵

再改进，想法是：
对N-1行和N-1列的取值是有范围的，最小值是SUM - N2 - (前N-2项的和)
因为SUM - 第N项的值 == 前N-1项的和
又因为第N项的值最大是N2，所以SUM - N2 == 前N-2项的和 + 第N-1项的和
回头再改改。

libinary

按上面的想法改了func函数，其他代码没变：

1. 
   
2. void
   
3. func(int n)
   
4. {
   
5.   int i;
   
6.   int min = 0, max = N2, sum = 0;
   
7. 
   
8.   if(n % N == N - 2){
   
9.     for(i = n - N + 2; i < n; ++i)
   
10.       sum += a[i];
    
11.     min = SUM - N2 - sum - 1 > 0 ? SUM - N2 - sum - 1 : 0;
    
12.     max = SUM - sum - 1 < N2 ? SUM - sum - 1 : N2;
    
13.   }else if(n / N == N - 2){
    
14.     for(i = n % N; i < n; i += N)
    
15.       sum += a[i];
    
16.     min = SUM - N2 - sum - 1 > 0 ? SUM - N2 - sum - 1 : 0;
    
17.     max = SUM - sum - 1 < N2 ? SUM - sum - 1 : N2;
    
18.   }
    
19. 
    
20.   for(i = min; i < max; ++i){
    
21.     if(num[i] == 0){
    
22.       num[i] = 1;
    
23.       a[n] = i + 1;
    
24.       if(n == N2 - N - 2)
    
25.         gen();
    
26.       else if(n % N == N - 2)
    
27.         func(n + 2);
    
28.       else
    
29.         func(n + 1);
    
30.       num[i] = 0;
    
31.     }
    
32.   }
    
33. }
    

复制代码

time：
N=3：0.006
N=4：7m39s (呵呵，快了5倍) 7040个结果

lordbyorn

我的想法是这样的。对每一个数都要检查。我的range函数就是这样的作用。
为什么要这样做呢？？在(N×N)！中，每减小一个乘数，程序的速度就会快很多。
我的是这样的：
用sum减去同列(行，对角线)中已经填了的数，若再减去可能填的最大数的和，可以求出可能填的最小数(若再减去可能填的最小数的和，可以求出可能填的最大数)
如
1 2 - - -
要填第三个数，
sum=65
minnum=65-1-2-24-25=13
maxnum=65-1-2-3-4=55
因为55>25;所以maxnum=25
所以第三个数的范围是 13<=num<=25
很好吧？？？

我现在想的是怎样优化range函数。
libinary，你从这个起点出发吧，除非你有更好的办法。
我想到一点办法了。开工！！！！

lordbyorn

最初由 libinary 发表

我运行过你的程序，N=4的时候好象没有30s就出来吧，是不是30m？

什么意思？？30s出不来？？什么30M?
我的最高记录是23s(N=4).
用O2,在console下，没有开其它的程序。

kj501

考虑到方阵的对称性，可以用旋转和对称操作来减少运算量。一个满足条件的幻方分别旋转90度、180度和270度就变成了另外一个满足条件的幻方。再看看对称，軕对称方式有两种，分别是上下对称和左右对称。由于上下对称和旋转180度等价。实际上只有左右对称。这样组合起来，一个满足条件的幻方就可以通过旋转和对称操作来变成8个满足条件的幻方。事实上算出八分之一，然后通过旋转和对称操作就可以得到全部的幻方。以N＝4为例，只要算出880个，然后通过旋转和对称操作就可以得到880X4=7040个幻方了。

lordbyorn

问题是怎样判断递归结束?

libinary

程序正在继续改进，今天看了一下你的程序，我的程序简直快变成和你的一模一样了，呼呼，晕。

第一次看你的程序的时候感觉看不懂，就没好好看，后来写的求等式数目的程序还附了一堆说明，原来这个不加注释你也懂呀，又晕。

lordbyorn

呵呵，我还没看就懂了。

第一次看你的程序的时候感觉看不懂

我第二次看我的程序时都没用看懂（range函数不懂）

lordbyorn

我准备求完美魔方.如果一个n阶幻方任一行任一列及任一泛对角线上各数的和相等， 则称此幻方为n阶完美幻方。

lordbyorn

exit 函数在哪里？？
我的n阶完美幻方程序写好了。

1. 
   
2. #include<time.h>
   
3. #include<stdio.h>
   
4. 
   
5. #define Max 20
   
6. #define MAX(x,y) (x)>(y)?(x) : (y)
   
7. 
   
8. int count=0;
   
9. int END;
   
10. int M,MM;
    
11. int sum;
    
12. int maxnum,minnum;
    
13. int data[Max*Max];
    
14. int flags[Max*Max];
    
15. 
    
16. FILE *outfile;
    
17. time_t start,end;
    
18. int min_not_used[Max+1]={0};
    
19. int max_not_used[Max+1]={0};
    
20. 
    
21. class STACKDATA{
    
22. private:
    
23.         int buff[Max+1];
    
24.         int cursor;
    
25. public:
    
26.         STACKDATA(){
    
27.                 cursor=-1;
    
28.                 sum=0;
    
29.         }
    
30.         ~STACKDATA(){
    
31.         }
    
32.         int push(int newdata){
    
33.                 buff[++cursor]=newdata;
    
34.                 sum+=newdata;
    
35.                 return 1;
    
36.         }
    
37.         int pop(){
    
38.                 if(cursor<0)return -1;
    
39.                 sum-=buff[cursor];
    
40.                 return buff[cursor--];
    
41.         }
    
42.         int sum;
    
43. };
    
44. 
    
45. void init();
    
46. void bye();
    
47. void range(int,int,int,int );
    
48. void magic(int);
    
49. void output(void);
    
50. 
    
51. 
    
52. STACKDATA rowbuff[Max+1];
    
53. STACKDATA colbuff[Max+1];
    
54. STACKDATA ldbuff[Max+1];
    
55. STACKDATA rdbuff[Max+1];
    
56. 
    
57. /***********  bigin main   *************/
    
58. 
    
59. int main(int argc,char *argv[]) {
    
60. if(argc!=2){
    
61.          fprintf(stderr,"usage: %s num \n",argv[0]);
    
62.          return 1;
    
63. }
    
64. sscanf(argv[1],"%d",&M);
    
65. init();
    
66. magic(0);
    
67. bye();
    
68. return 0;
    
69. }
    
70. 
    
71. /******    end of main  *****************/
    
72. 
    
73. void init()
    
74. {
    
75. 
    
76. int i;
    
77. outfile=fopen("result.txt","w");
    
78. for(i=0;i<Max*Max;i++)flags[i]=0;
    
79. if(M>Max){puts("Max M exceed!!\n");}
    
80. MM=M*M;
    
81. sum=(MM+1)*M/2;
    
82. END=MM-1;
    
83. start=time(&start);
    
84. fprintf(outfile,"start at %s\n//////////////////////////////////\n",ctime(&start));
    
85. 
    
86. }
    
87. 
    
88. /************************************************/
    
89. void bye()
    
90. {
    
91. end=time(&end);
    
92. fprintf(outfile,"\n//////////////////////////////////\nend at %s\n",ctime(&end));
    
93. fprintf(outfile,"%d matchs found!",count);
    
94. fclose(outfile);
    
95. }
    
96. /***************************************************/
    
97. void magic(int n)
    
98. {
    
99. int i,ii;
    
100. int max;
     
101. int col,row,diag;
     
102. 
     
103. row=n/M;
     
104. col=n%M;
     
105. diag=(col-row+M)%M;
     
106. range(n,row,col,diag);
     
107. 
     
108.         for(i=minnum,max=maxnum;i<=max;i++) {
     
109.                 if(flags[i]==0) {
     
110.                         ii=i+1;
     
111.                         flags[i]=1;
     
112.                         data[n]=ii;
     
113.                         rowbuff[row].push(ii);
     
114.                         colbuff[col].push(ii);
     
115.                         ldbuff[diag].push(ii);
     
116.                         rdbuff[diag].push(ii);
     
117.                         if(n==(END)){
     
118.                                 output();
     
119.                         }
     
120.                         else magic(n+1);
     
121.                         flags[i]=0;
     
122.                         rowbuff[row].pop();
     
123.                         colbuff[col].pop();
     
124.                         ldbuff[diag].pop();
     
125.                         rdbuff[diag].pop();
     
126.                 }
     
127.         }
     
128.         return;
     
129. }
     
130. 
     
131. /**********************************************************/
     
132. void output()
     
133. {
     
134.         int i,j;
     
135. 
     
136.         count++;
     
137.         fprintf(outfile,"%d:\n",count);
     
138.         for(i=0;i<M;i++){
     
139.                 for(j=0;j<M;j++){
     
140.                         fprintf(outfile,"%5d",data[i*M+j]);
     
141.                 }
     
142.                 fputc('\n',outfile);
     
143.         }
     
144.         fflush(outfile);
     
145. }
     
146. 
     
147. 
     
148. 
     
149. /*******************  This function to large!!!        ************/
     
150. void range(int n,int row,int col,int diag) {
     
151. static int used;
     
152. static int bottom,top,left;
     
153. static int c;
     
154. static int tmp;
     
155. static int i,j,k;
     
156. 
     
157. k=MAX(M-col,M-row);
     
158. 
     
159. if(k>M)k=M;
     
160. for(i=0,c=1,tmp=0;c<=k;i++){
     
161.          if(flags[i]==0){
     
162.                 tmp+=i+1;
     
163.                 c++;
     
164.                 min_not_used[c]=tmp;
     
165.          }
     
166. }
     
167. for(i=END,c=1,tmp=0;c<=k;i--){
     
168.          if(flags[i]==0){
     
169.                 tmp+=i+1;
     
170.                 c++;
     
171.                 max_not_used[c]=tmp;
     
172.          }
     
173. }
     
174. 
     
175. used=rowbuff[row].sum;
     
176. left=sum-used-1;
     
177. bottom=max_not_used[M-col];
     
178. top   =min_not_used[M-col];
     
179.         minnum=left-bottom;
     
180.         maxnum=left-top;
     
181.          
     
182. used=colbuff[col].sum;
     
183. left=sum-used-1;
     
184. bottom=max_not_used[M-row];
     
185. top   =min_not_used[M-row];
     
186. minnum=minnum<left-bottom?left-bottom:minnum;
     
187. maxnum=maxnum>left-top?left-top:maxnum;
     
188. 
     
189. used=ldbuff[diag].sum;
     
190. left=sum-used-1;
     
191. bottom=max_not_used[M-row];
     
192. top   =min_not_used[M-row];
     
193. minnum=minnum<left-bottom?left-bottom:minnum;
     
194. maxnum=maxnum>left-top?left-top:maxnum;
     
195. 
     
196. used=rdbuff[diag].sum;
     
197. left=sum-used-1;
     
198. bottom=max_not_used[M-row];
     
199. top   =min_not_used[M-row];
     
200. minnum=minnum<left-bottom?left-bottom:minnum;
     
201. maxnum=maxnum>left-top?left-top:maxnum;
     
202. 
     
203. if(minnum<0)minnum=0;
     
204. if(maxnum>=MM)maxnum=MM-1;
     
205. }
     
206. 
     
207. /*****************  end of program   *****************************************/
     

复制代码