一道数字操作题

mathe

上面的代码是如此的低效,以致于对于n=5的时候(这时空间复杂度还不是大问题,4维的rch数组现在的计算机还能够忍受),计算速度已经让人无法忍受了. 可以知道,实际计算过程计算机需要对于每个状态都调用一次上面的pass函数,而上面的pass函数实现过程使用了一个递归过程,显然里面会出现很多冗余的状态(很多状态会被多次遍历到,入对于"加一"操作两次分别添加到第1和第2个元素的操作,上面代码会两次遇到,分别是先加到第一个元素和先加到第二个元素. 为此,我对这个部分做了简单的优化,结果计算n=5就不是问题了: 思路很简单,我们先将所有的操作计算出来,保存在一个数组中,将重复的操作都剔除掉,然后实际使用过程分别试验使用这些预先保存的操作就可以提高很多效率了. 对于n=5的情况,我们可以将所有可能的操作进行编码,其中每个变换后的元素(4个元素)每一个都可以是变换前某个元素(5个选择,可以用3比特表示)通过0次(只能用于自身),1次,2次或3次"加一"操作得到(不同加法可以用2比特表示) 所以所有操作可以用最多4*(3+2)=20比特来表示 这部分代码同前面代码很像,它只是创建一个code数组

2. int ocount;
3. #define MAX_OC 10000
4. int code[MAX_OC];
5. #define GET_OFFSET(x) ((x)&3)
6. #define GET_TARGET(x) (((x)>>2)&3)
7. #define CREATE_VALUE(offset,target) (((target)<<2)|(offset))
8. void add(int x, int y, int z, int w)
9. {
10. int i;
11. int u=TARGET(x,y,z,w);
12. for(i=0;i<ocount;i++){
13. if(code[i]==u)
14. break;
15. }
16. if(i==ocount){
17. code[ocount++]=u;
18. }
19. }
21. #define enum_code_0(a,b,c,d) add(a,b,c,d)
22. void enum_code_1(int x, int y, int z, int w)
23. {
24. int xu,yu,zu,wu;
25. int xa,ya,za,wa;
26. int xn,yn,zn,wn;
27. add(x,y,z,w);
28. xu=GET_TARGET(x);
29. yu=GET_TARGET(y);
30. zu=GET_TARGET(z);
31. wu=GET_TARGET(w);
32. xa=GET_OFFSET(x);
33. ya=GET_OFFSET(y);
34. za=GET_OFFSET(z);
35. wa=GET_OFFSET(w);
36. xn=CREATE_VALUE(xa+1,xu);
37. yn=CREATE_VALUE(ya+1,yu);
38. zn=CREATE_VALUE(za+1,zu);
39. wn=CREATE_VALUE(wa+1,wu);
40. enum_code_0(xn,y,z,w);
41. enum_code_0(x,yn,z,w);
42. enum_code_0(x,y,zn,w);
43. enum_code_0(x,y,z,wn);
44. enum_code_0(x,xn,z,w);
45. enum_code_0(x,y,xn,w);
46. enum_code_0(x,y,z,xn);
47. enum_code_0(yn,y,z,w);
48. enum_code_0(x,y,yn,w);
49. enum_code_0(x,y,z,yn);
50. enum_code_0(zn,y,z,w);
51. enum_code_0(x,zn,z,w);
52. enum_code_0(x,y,z,zn);
53. enum_code_0(wn,y,z,w);
54. enum_code_0(x,wn,z,w);
55. enum_code_0(x,y,wn,w);
56. }
58. void enum_code_2(int x, int y, int z, int w)
59. {
60. int xu,yu,zu,wu;
61. int xa,ya,za,wa;
62. int xn,yn,zn,wn;
63. add(x,y,z,w);
64. xu=GET_TARGET(x);
65. yu=GET_TARGET(y);
66. zu=GET_TARGET(z);
67. wu=GET_TARGET(w);
68. xa=GET_OFFSET(x);
69. ya=GET_OFFSET(y);
70. za=GET_OFFSET(z);
71. wa=GET_OFFSET(w);
72. xn=CREATE_VALUE(xa+1,xu);
73. yn=CREATE_VALUE(ya+1,yu);
74. zn=CREATE_VALUE(za+1,zu);
75. wn=CREATE_VALUE(wa+1,wu);
76. enum_code_1(xn,y,z,w);
77. enum_code_1(x,yn,z,w);
78. enum_code_1(x,y,zn,w);
79. enum_code_1(x,y,z,wn);
80. enum_code_1(x,xn,z,w);
81. enum_code_1(x,y,xn,w);
82. enum_code_1(x,y,z,xn);
83. enum_code_1(yn,y,z,w);
84. enum_code_1(x,y,yn,w);
85. enum_code_1(x,y,z,yn);
86. enum_code_1(zn,y,z,w);
87. enum_code_1(x,zn,z,w);
88. enum_code_1(x,y,z,zn);
89. enum_code_1(wn,y,z,w);
90. enum_code_1(x,wn,z,w);
91. enum_code_1(x,y,wn,w);
92. }
94. void enum_code_3(int x, int y, int z, int w)
95. {
96. int xu,yu,zu,wu;
97. int xa,ya,za,wa;
98. int xn,yn,zn,wn;
99. add(x,y,z,w);
100. xu=GET_TARGET(x);
101. yu=GET_TARGET(y);
102. zu=GET_TARGET(z);
103. wu=GET_TARGET(w);
104. xa=GET_OFFSET(x);
105. ya=GET_OFFSET(y);
106. za=GET_OFFSET(z);
107. wa=GET_OFFSET(w);
108. xn=CREATE_VALUE(xa+1,xu);
109. yn=CREATE_VALUE(ya+1,yu);
110. zn=CREATE_VALUE(za+1,zu);
111. wn=CREATE_VALUE(wa+1,wu);
112. enum_code_2(xn,y,z,w);
113. enum_code_2(x,yn,z,w);
114. enum_code_2(x,y,zn,w);
115. enum_code_2(x,y,z,wn);
116. enum_code_2(x,xn,z,w);
117. enum_code_2(x,y,xn,w);
118. enum_code_2(x,y,z,xn);
119. enum_code_2(yn,y,z,w);
120. enum_code_2(x,y,yn,w);
121. enum_code_2(x,y,z,yn);
122. enum_code_2(zn,y,z,w);
123. enum_code_2(x,zn,z,w);
124. enum_code_2(x,y,z,zn);
125. enum_code_2(wn,y,z,w);
126. enum_code_2(x,wn,z,w);
127. enum_code_2(x,y,wn,w);
128. }
130. void init()
131. {
132. memset(rch,-1,sizeof(rch));
133. enum_code_3(CREATE_VALUE(0,0),CREATE_VALUE(0,1),CREATE_VALUE(0,2),CREATE_VALUE(0,3));
134. printf("total ocount=%d\n",ocount);
135. rch[0][0][0][0]=TARGET(0,0,0,0);
136. }

复制代码

此后,n=5对应的pass函数可以修改为:

2. int test(int x[4], int t)
3. {
4. int i;
5. int y[4];
6. for(i=0;i<ocount;i++){
7. int u=code[i];
8. int xu,yu,zu,wu;
9. int xa,ya,za,wa;
10. int x1,x2,x3,x4;
11. x1=(u>>24)&0xFF;
12. x2=(u>>16)&0xFF;
13. x3=(u>>8)&0xFF;
14. x4=u&0xFF;
15. xu=GET_TARGET(x1);
16. xa=GET_OFFSET(x1);
17. yu=GET_TARGET(x2);
18. ya=GET_OFFSET(x2);
19. zu=GET_TARGET(x3);
20. za=GET_OFFSET(x3);
21. wu=GET_TARGET(x4);
22. wa=GET_OFFSET(x4);
23. y[0]=x[xu]+xa;
24. y[1]=x[yu]+ya;
25. y[2]=x[zu]+za;
26. y[3]=x[wu]+wa;
27. if(y[0]>t||y[1]>t||y[2]>t||y[3]>t)
28. continue;
29. if(y[0]<t&&y[1]<t&&y[2]<t&&y[3]<t)
30. continue;
31. if(rch[y[0]][y[1]][y[2]][y[3]]>=0)
32. break;
33. }
34. if(i<ocount)return i;
35. else return -1;
36. }
38. int pass(int a,int b, int c, int d)
39. {
40. int x[4];
41. int i;
42. x[0]=a,x[1]=b,x[2]=c,x[3]=d;
43. for(i=0;i<4;i++){
44. if(x[i]==d)x[i]=0;
45. }
46. if(d==0)return 0;
47. return test(x,d-1);
48. }

复制代码

mathe

程序开头还有定义:

2. #include <stdio.h>
3. #include <string.h>
5. #define LIM 100
6. #define N 5
7. int rch[LIM][LIM][LIM][LIM];
9. #define SET(x,y,z,w,v) rch[x][y][z][w]=rch[x][z][y][w]=rch[y][x][z][w]=\
10. rch[y][z][x][w]=rch[z][x][y][w]=rch[z][y][x][w]=\
11. rch[x][y][w][z]=rch[x][w][y][z]=rch[y][x][w][z]=\
12. rch[y][w][x][z]=rch[w][x][y][z]=rch[w][y][x][z]=\
13. rch[x][z][w][y]=rch[x][w][z][y]=rch[z][x][w][y]=\
14. rch[z][w][x][y]=rch[w][x][z][y]=rch[w][z][x][y]=\
15. rch[z][y][w][x]=rch[z][w][y][x]=rch[y][z][w][x]=\
16. rch[y][w][z][x]=rch[w][z][y][x]=rch[w][y][z][x]=\
17. (v)
19. #define TARGET(x,y,z,w) (((x)<<24)|((y)<<16)|((z)<<8)|(w))

复制代码

再把main函数贴出,这就是完整的n=5时候的代码了

2. int main()
3. {
4. int i,j,k,h,r;
5. int x[4];
6. init();
7. for(i=0;i<LIM;i++){
8. for(j=0;j<=i;j++){
9. for(k=0;k<=j;k++){
10. for(h=0;h<=k;h++){
11. if((r=pass(h, k,j,i))>=0){
12. SET(h,k,j,i,r);
13. }
14. }
15. }
16. }
17. }
18. for(i=LIM-1;i>=0;i--){
19. for(j=i;j>=0;j--){
20. for(k=j;k>=0;k--){
21. for(h=k;h>=0;h--){
22. if(rch[h][k][j][i]>=0){
23. printf("Find 0,%d,%d,%d,%d\n",h,k,j,i);
24. goto end;
25. }
26. }
27. }
28. }
29. }
30. end:
31. x[0]=h,x[1]=k,x[2]=j,x[3]=i;
32. do{
33. int u=code[rch[x[0]][x[1]][x[2]][x[3]]];
34. int y[4];
35. int xu,yu,zu,wu,xa,ya,za,wa;
36. int x1,x2,x3,x4;
37. x1=(u>>24)&0xFF;
38. x2=(u>>16)&0xFF;
39. x3=(u>>8)&0xFF;
40. x4=(u)&0xFF;
41. xu=GET_TARGET(x1);
42. xa=GET_OFFSET(x1);
43. yu=GET_TARGET(x2);
44. ya=GET_OFFSET(x2);
45. zu=GET_TARGET(x3);
46. za=GET_OFFSET(x3);
47. wu=GET_TARGET(x4);
48. wa=GET_OFFSET(x4);
49. x[3]=0;
50. y[0]=x[xu]+xa;
51. y[1]=x[yu]+ya;
52. y[2]=x[zu]+za;
53. y[3]=x[wu]+wa;
54. for(i=0;i<4;i++)x[i]=y[i];
55. for(i=0;i<4;i++){
56. for(x1=0;x1<i;x1++){
57. if(x[x1]>x[i]){
58. x2=x[x1];
59. x[x1]=x[i];
60. x[i]=x2;
61. }
62. }
63. }
64. printf("From 0 %d %d %d %d\n",x[0],x[1],x[2],x[3]);
65. if(x[3]==0)break;
66. if(rch[x[0]][x[1]][x[2]][x[3]]<0){
67. printf("ERROR\n");
68. return -1;
69. }
70. }while(1);
71. return 0;
72. }

复制代码

mathe

对于$n>=6$上面的方法显然已经不能再使用了,我们没有办法分配一个5维的而且每一个维数都要到达100左右的数组了. 考虑到实际的能够到达0的数列中,数列中最大和第二大数的差值只能是1,2,3,4之一. 而最大同第三大数的差值通常也不会太大,我们可以预先设置几个宏 L1,L2,L3,L4 如果n==6时数列中5个数分别时x[0],x[1],x[2],x[3],x[4](递增数列) 我们只搜索x[4]-x[3]

2. #define LIM 100
3. #define N 6
4. #define MAX_RULE N
5. #define L1 5
6. #define L2 25
7. #define L3 40
8. #define L4 70
9. int rch[L1*L2*L3*L4*LIM];
10. int muls[]={LIM,L1,L2,L3,L4};
12. int encode(int x[N])
13. {
14. int i;
15. int s=x[N-2];
16. for(i=N-3;i>=0;i--){
17. s*=muls[N-2-i];
18. s+=x[N-2]-x[i];
19. }
20. return s;
21. }
23. void decode(int c, int x[N])
24. {
25. int i;
26. for(i=0;i<N-2;i++){
27. x[i]=c%muls[N-2-i];
28. c/=muls[N-2-i];
29. }
30. x[i]=c;
31. for(i=N-3;i>=0;i--){
32. x[i]=x[N-2]-x[i];
33. }
34. }

复制代码这个方法可以降低了内存的使用量,代价是增加了编码转化过程,从而计算速度变慢,而且有可能会漏掉一些解(除非我们能够证明使用的L1,L2,L3,L4等都的的确确足够大了) 此外,对于前面提到的状态转化数组,考虑到有些变换只是位置的交换,它们应该是完全等价的,我们还可以预先淘汰更加多的等价状态转化过程,比如 将第一个元素替换为第三个元素加1,并且将第二个元素加2 和将第一个元素替换为第二个元素加2,并且将第二个元素替换为将第三个元素加1 是完全等价的结果. 所以程序里面添加了先对这些转化数据排序的过程,代码见附件,其中函数push_code(添加一个代码)中调用了个sort_x函数,用来对状态转化过程标准化. ens66.tar.gz (1.57 KB, 下载次数: 1)

2008-6-6 15:58 上传

点击文件名下载附件

mathe

上面的方法在n=7时候的内存开销已经有点让人无法忍受了 为此在n=7时,我将里面的rch数组改成用比特位数组来表示, 通过这种修改,这个版本的程序能够勉勉强强运行n=7的情况了,只是速度非常慢. 见附件中ens77.cpp文件 另外,我注意到状态转化数组codes含有比较多的状态转化,但是如果我们能够将需要转化的数组更具它前面几个数的情况进行分类,那么就可以对于每个类别只尝试较少的状态转化过程 比如次大元素和最大元素只差1的情况,我们只需要选择那些保留第二大元素的状态转化就可以了. 这个思路可以将n=7情况的搜索过程提高至少10倍以上的速度. 附件中的ens68.cpp文件是我想到一个猜想来用来验证猜想的代码(结果没有什么意义),里面首次用了这个想法 可以看到里面有个trans_codes函数用来分类这个数组到几个小数组. 其次,我还有一个想法,其实我们完全没有必要查看那些将非零元素替换成(越界的元素总是看成0)其它元素加多少的变换,也就是我们只要看: i)将某个0加若干 ii)将某个元素加若干 而为了保险,还允许另外一种情况,将连续若干比较小的数完全看成0(也就是可以任意被其它数加若干覆盖) 这两种情况就可以了, 所以就添加了一个well_form函数,只允许使用含这三种情况的变换. 不过实际计算表明这个优化作用用处可能不大. ens77.tar.gz (2.55 KB, 下载次数: 0)

2008-6-6 16:14 上传

点击文件名下载附件

mathe

上面的算法就是计算n=7的情况也需要数小时（而对于n=6只需要很短时间)，而且找出来的结果无法保证是最优的， 这也是为什么第一次对于n=6先找出了84，后来重新调整参数才找到了85。 而对于n=7,第一次通过冗长运算找到了142，但是我一直固执的认为可能最优结果是144(这样就可以符合一个4次多项式了） 而对于更大的n,上面程序已经基本无能为力了。 其实在写上面代码过程中，我已经逐渐现成后面将要介绍的方法，只是算法需要用到树结构，比较复杂，所以一直没有动手。 对于一个最大数字为M的序列，如果我们能够按照题目定义的方法一次释放M,M-1,...,0;我们就称这个序列是合格序列，那么 我们先可以考虑一些定理: 定理一:如果一个逆序序列$a_1,a_2,...,a_n=0$是长度为n的合格序列，其中含有的数字0不唯一，那么我们将其中一个0替换成任意其它不超过a_1的序列还是合格的 定理二:如果一个合格序列中某个数字重复出现，我们将其中一个重复的数字替换成0，结果依旧是合格序列 定理三:如果一个逆序序列$a_1,a_2,...,a_n=0$是长度为n,最大元为$a_1$的合格序列，我们将其中每个元素减去k($k

mathe

我们现在开始讨论这个新算法 对于一个合格逆序序列，我们先查看最大数字是4的序列 由于我们知道 4,0,0,...,0 是合格的，我们可以得出任意以4开头的逆序序列都是合格的。 我们记这样的序列为 4e,其中后缀e表示后面可以添加一些比4小的逆序数（最后还有若干0） 同样对于5开头的逆序序列，只要第二个元素是1到4的数，那么第三个以后的元素无论是多少都是合格序列，我们可以记为 5 1-4e 同样6开头的逆序序列可以记为 6 2-5e 也就是6开头的逆序序列，如果第二个元素是0或1，那么就不是合格的，不然都是合格的。 然后还有 7 3-6e 8 4-7e 9 8 1-7e 9 5-7e 10 9 2-8e 10 6-8e 11 10 3-9e 11 7-9e 12 11 4-10e 12 10 1-9e 12 8-9e 13 12 11 1-10e 13 12 5-10e 13 11 2-10e 13 9-10e 可以看出，对于给定开头数字的所有这些序列，我们都可以把它们表示成树的形式，其中有后缀e的数为叶节点（可以看出叶节点有上界和下界，而中间节点都只含一个数字） 而我们的任务就是依次对于每个开头数字构造对应的树。 而定理三给出了使用开头数字为m的树构造开头数字为m+1的树的一种比较好的方案。 我们只要将开头数字为m+1的树上每个数字加1，然后检查所有对应序列是否是合格的， 而对于每个叶节点，由于它包含的是一个范围，根据猜想二，我们只需要从这个范围的两边开始检查，只要遇到一个合格的数据就不需要继续检查下去了。 然后这个范围就可以分裂成最多三个部分，中间部分对应合格序列，两边对应不合格序列 对于中间的合格序列，我们保留成一个叶节点就可以了（只是范围缩小了） 对于两边那些不合格的序列，每个数字都需要分裂成一个中间节点，然后对它们添加一个叶节点1 -- (k-1),其中k是对应分裂出来节点含的数字。 比如5 1-4e构造6开头的序列，我们可以马上从6 2-5e开始，由于6 2 0...0和6 5 0..0都是合格序列，我们得到这个不需要继续分裂了。 而从4e构造5开头序列，我们可以从5e开始，由于5 0 .. 0不是合格序列，我们需要添加叶节点1-4e成为 5 1-4e.然后我们需要再次判断这个是不是已经完全是合格序列了，这时我们需要判断5 1 0 ... 0 和5 4 0...0都是合格序列就可以了。（而我的第一个程序由于我猜测这时构造序列总是合格的，没有继续检测，结果得到了偏大的范围，所以对于n=8得到不超过222的结论） 而对于m+1开头序列，判断序列是否合格，我们通过一次变换去掉最大数据后，将得到结果用m开头的树来匹配就可以了，只要能够部分匹配，就是合格的，不然就是不合格的。

mathe

最后给出对应的源程序 ense.tar.gz (4.21 KB, 下载次数: 0)

2008-6-6 18:03 上传

点击文件名下载附件

其中包含3个程序 ense.cpp程序会将所有的树创建出来,以文件形式保存下来. 程序最后结束时创建编号最大数目的文件就是能够达到的序列开始数字的最大值. rt.c程序可以将其中任意一棵树(文件形式保存)用文本形式显示出来 ensr.cpp程序可以将一个结果给构造出来(前提时所有的树对应的文件已经创建) 比如ense在n=8时打印出来最大编号为221 我们运行./ensr 221就可以打印出一个221开头的序列,然后将它一次变换到5 1 0 ... 0(后面几步就很显然了,没有打印)

mathe

现在还想到一个问题，我前面一个程序中定义了L1,L2,L3,L4分别为合格数列最大元素和次大元素差值的上界，最大元素和第三大元素差值的上界等等， 很显然，L1=4就可以遍历所有的解了。那么请问L2,L3,L4分别应该是多少，才不会漏掉合格序列（也就是可以保证找到最优解）

mathe

原帖由 无心人 于 2008-6-6 11:25 发表 2G的完全可以利用到2000M 只要光用linux的字符界面

在n=9的时候,果然代码运行到第198层(花费4天多时间)时就内存不够了(预划分1.6G内存),昨天重新修改了一下代码,将双缓冲改成单缓冲,而结果数据直接写入文件.现在重新从197层开始计算了,也许再过一周左右时间能够给出n=9的最终结果?(还要希望不要停电)

无心人

你还真懒得写状态保存部分阿