埃及分数问题

mathe

用你的方法，仅仅为了计算精度问题，在LCM大于64位整数的情况，可以采用多个互素的整数， 比如取三个互素的32位整数p1,p2,p3,使得p1*p2*p2大于LCM. 然后所有的求和过程改成求和并分别对p1,p2,p3求模（实际上如果只有加法，那么除法运算都不需要）。 当然，这样，搜索过程中和太大和太小本来可以提前剪枝的，这个方法不行。所以可以通过移位方法来做删减。 比如假设LCM比64位整数大8个比特，那么我们可以将所有的数除以256（分别求上取整和下取整）。在求和过程中同时计算这些上取整和下取整的和，如果发现它们越界，也可以提前裁减。 不过我上面介绍的方法，将会采用不同的方法。

liangbch

原帖由 无心人 于 2008-3-7 13:55 发表 这个题目要求技巧高些而已 还有输出多点 应该超过10万吧 ：）

mathe 和 楼主， 能否告诉我， 你们预测的解的个数是多少个，如果解的个数太多，建议降低搜索层数，比如从16减低10或者更低。 我采用简单的剪枝方法（9#中的规则 “i)我们可以用部分和进行剪枝:”），用自定义的大整数验证，当搜索层数为 2-9是，解的个数为 level 解个数 时间 2 1 3 3 4 12 5 117 6 1250 7 12986 8 126498 14 秒 9 1179414 567 秒 按此估计，即使以后每增加1层，解的个数增加5倍，最终解的个数也将超过90G, 而如果每增加1层，解的个数增加8倍，最终解的个数也将超过2000G。由于解的数目太多，将导致计算时间太长，验证也很困难。 顺便说一下，我在14楼和 20楼给出的程序 有bug,统计出的个数不对。

无心人

你先剔除23, 29, 31, 37, 41, 47等大数字 2^6 * 3^4 * 5^2 * 7^2 * 11 * 13 * 17 * 19 = 2933 1862 5600 用2933 1862 5600 / n 代替 1 / n 则整个算法成为整数加的算法 然后将2-256中23, 29, 31, 37, 41, 43的倍数剔除 再剔除128, 256, 243, 125, 250几个 得到一个表 对每个数字n , 使用2933 1862 5600 / n 代替 则使用穷举就能得到最后答案 只要注意当前面数字大到剩余数字用最小的1/255添充都太大 或者小到使用当前数字填充都太小 的时候就要终止 时间应该远小于你原来时间 至于添加23, 29, 31, 37, 41, 47的倍数的情况 mathe稍后有更好的算法

mathe

原帖由 liangbch 于 2008-3-10 14:39 发表 mathe 和 楼主， 能否告诉我， 你们预测的解的个数是多少个，如果解的个数太多，建议降低搜索层数，比如从16减低10或者更低。 我采用简单的剪枝方法（9#中的规则 “i)我们可以用部分和进行剪枝:”），用自定义 ...

我没有真正计算过，不过也感觉解的数目会太多。所以我首先的目标是就算解的数目，而不是将所有解输出。

mathe

现在介绍一下我的思路，主要数学依据还是在36楼。 根据36楼，我们可以知道，比如对于所有分母为17的倍数的埃及分数： ${1/17, 1/{2*17}, ...,1/{15*17} }$ 它们出现的模式（不考虑其它埃及分数是否出现）总共有1927种，比如: $1/17(1/1+1/2+1/5)$ $1/17(1/1+1/4+1/6)$ $1/17(1/1+1/3+1/5+1/6)$ $1/17(1/1+1/2+1/3+1/4+1/7)$ $1/17(1/2+1/4+1/5+1/7)$ $1/17(1/2+1/6+1/7)$ 等等 搜索这些模式的程序很简单 我们现在先考虑如何比较有效的保存它们。 对于每个模式，我们可以考虑用64比特信息保存 其中32比特用来保存模式中埃及分数和乘上720720的值（其实20比特就够了） 比如对于上面列出那些数据，结果分别为 72072 60060 72072 94380 46332 34320 另外我们需要知道每个数据使用了多少个埃及分数，这个数值不会超过15，所以4比特足够了，也采用8比特吧。 另外，我们希望在最后每个搜索结果出来后，可以很容易从模式还原会这些埃及序列，这个使用两部分数据表示 i) 模式是在那个基数的倍数，比如这里是17（这个数不会超过255，所以8比特足够） ii)模式使用了基数的那些倍数，最多使用了15个倍数。所以我们可以采用15位比特位来表示 比如这里对于第一个模式，分别是17的1倍，2倍，5倍被使用了，所以我们可以采用二进制数0010011来表示这个部分 而我们还知道，对于27的倍数，还有一种特殊情况，就是还可以同时出现$1/81+1/162$,我们在添加一个比特位，比特位为真表示 还要添加$1/81+1/162$,所以这部分也使用16比特。 struct ELEMENT{ unsigned value; unsigned char count; unsigned char base; unsigned short mask; }; 所以这里，第一个模式可以表示成 struct ELEMENT x={72072, 3, 17, 0010011}; 而对于27的倍数的模式: $1/27(1/1+1/2+1/4+1/5+1/7)+1/81+1/162$ 可以表示为 struct ELEMENT x2={69212,7,27,0x805B};

mathe

27应该对应模式应该是63个，而不是64个，上面写错了，下面代码可以打印出所有这些模式:

1. // flter.cpp : Defines the entry point for the console application.
2. //
4. #include "stdafx.h"
5. #define OUTPUT
7. int mask[]={1,1,2,6,12,60,60,420,840,2520,2520,27720,27720,360360,360360,360360,720720};
8. int cands[]={17,19,23,29,31,37,41,43,47,53,59,61};
9. int c2[]={5,7};
10. int c3[]={1,2,4,5,6,7,2};
11. #define NUMOF(x) (sizeof(x)/sizeof(x[0]))
12. int _tmain(int argc, _TCHAR* argv[])
13. {
14. int i,j;
15. for(i=0;i<NUMOF(cands);i++){
16. int p=cands[i];
17. int mt=256/p;
18. int count=0;
19. int lm=mask[mt];
20. for(j=1;j<(1<<mt);j++){
21. int sum=0;
22. int k;
23. for(k=0;k<mt;k++){
24. if(j&(1<<k)){
25. sum+=lm/(k+1);
26. sum%=p;
27. }
28. }
29. if(sum==0){
30. int start=1;
31. count++;
32. #ifdef OUTPUT
33. for(k=0;k<mt;k++){
34. if(j&(1<<k)){
35. if(!start)fprintf(stderr,"+");
36. fprintf(stderr,"1/%d",(k+1)*p);
37. start=0;
38. }
39. }
40. fprintf(stderr,"\n");
41. #endif
42. }
43. }
44. printf("Count %d for prime %d\n",count,p);
45. }
47. for(i=0;i<NUMOF(c2);i++){
48. int p=c2[i];
49. int p2=p*p;
50. int count=0;
51. int mt=256/p2;
52. int lm=mask[mt];
53. for(j=1;j<(1<<mt);j++){
54. int sum=0;
55. int k;
56. for(k=0;k<mt;k++){
57. if(j&(1<<k)){
58. sum+=lm/(k+1);
59. sum%=p;
60. }
61. }
62. if(sum==0){
63. int start=1;
64. count++;
65. #ifdef OUTPUT
66. for(k=0;k<mt;k++){
67. if(j&(1<<k)){
68. if(!start)fprintf(stderr,"+");
69. fprintf(stderr,"1/%d",(k+1)*p2);
70. start=0;
71. }
72. }
73. fprintf(stderr,"\n");
74. #endif
75. }
76. }
77. printf("Count %d for %d\n",count,p2);
78. }
80. {
81. int p=3;
82. int p3=27;
83. int mt=7;
84. int count=0;
85. int lm=mask[mt];
86. for(j=1;j<(1<<mt);j++){
87. int sum=0;
88. int k;
89. for(k=0;k<mt;k++){
90. if(j&(1<<k)){
91. sum+=lm/c3[k];
92. sum%=p;
93. }
94. }
95. if(sum==0){
96. int start=1;
97. count++;
98. #ifdef OUTPUT
99. for(k=0;k<mt;k++){
100. if(j&(1<<k)){
101. if(!start)fprintf(stderr,"+");
102. if(k==mt-1){
103. fprintf(stderr,"1/81+1/162");
104. }else{
105. fprintf(stderr,"1/%d",c3[k]*p3);
106. }
107. start=0;
108. }
109. }
110. fprintf(stderr,"\n");
111. #endif
112. }
113. }
114. printf("Count %d for 27\n",count);
115. lm = 1*3*5*7;
116. count=0;
117. for(j=1;j<16;j++){
118. int sum=0;
119. int k;
120. for(k=0;k<4;k++){
121. if(j&(1<<k)){
122. sum+=lm/(2*k+1);
123. sum%=2;
124. }
125. }
126. if(sum==0){
127. int start=1;
128. count++;
129. #ifdef OUTPUT
130. for(k=0;k<4;k++){
131. if(j&(1<<k)){
132. if(!start)fprintf(stderr,"+");
133. fprintf(stderr,"1/%d",(2*k+1)*32);
134. start=0;
135. }
136. }
137. fprintf(stderr,"\n");
138. #endif
139. }
140. }
141. printf("Count %d for 32\n",count);
142. }
143. return 0;
144. }

复制代码

我们首先需要将它们转化为我上面说明的结构。

mathe

上面程序还有64的倍数$1/64+1/{3*64}$没有输出 这样我们总共得到13个长度非0的模式列表，此外还有77个720720的因子，每个单独构成一个模式列表。 所以总共90个模式列表。 其中有三个模式列表最长，分别为 17的倍数1927个模式L17 19的倍数430个模式L19 25的倍数255个模式L25 对于上面三个模式列表，我们在列表末尾添上一个特殊数据ELEMENT{0,0,*,0}; 让后穷举它们的组合,并且定义一个长度为16*720720的数组tc[720720][16];

1. memset(tc,0,sizeof(tc));
2. for(i1=0;i1<1928;i1++)for(i2=0;i2<431;i2++)for(i3=0;i3<260;i3++){
3. int count=L17[i1].count+L19[i2].count+L25[i3].count;
4. int sum=L17[i1].value+L19[i2].value+L25[i3].value;
5. if(0<count&&count<=16&&sum<=720720){
6. tc[sum-1][count-1]++;
7. }
8. }

复制代码

mathe

然后对于余下的87个模式列表，我们按任意顺序事先排列好，假设为 ELEMENT *L[87]; 其中每项模式数目为int c[87]; 然后我们同样创建另外一个数组tc2[720720][16]; 后面的程序中，需要频繁交换两个数组的位置，我们可以采用来个指针来轮换它们地位，具体做法这里就不写出了。 现在假设tc数组中已经保存了使用前3个模式列表在加上87个模式列表中t个模式列表(0<=t<87)中的模式使用了k个数达到和x的数目保存在tc[x-1][k-1]中 那么我们可以使用动态规划计算再加上的t个模式后的结果为

1. memcpy(tc2,tc,sizeof(tc2));
2. for(i=0;i<c[t];i++){
3. int c=L[t][i].count;
4. int v=L[t][i].value;
5. for(j=c;j<=16;j++)for(s=v;s<=720720;s++){
6. tc2[s][j-1]+=tc[s-v][j-c];
7. }
8. }

复制代码

这个过程对于所有的t处理完的复杂度可以估计为16*720720*Sum(c[.])=16*720720*(88+8+6+1+1+1+4+3+7+1+77)=2.27G,现在的计算机处理起来应该很快。 全部计算完以后，那么tc2[720719][0..15]就保存了分别使用1,2,...,16个埃及分数达到1的组合确切数目。

mathe

当然上面的过程，如果我们使用了更大的内存，显然不难可以计算出16个以上埃及分数和为1的确切数目（主要问题在于内存使用量）。 现在我们开始考虑另外一个问题，如果我们需要将所有组合都输出，那么使用上面的方法，是否可以提供一种非常有效的方法。 我们想象，如果我们能够将上面计算过程中的所有tc2数组的内容都保存下来（共88个这么大的数组），然后我们再使用动态规划的方法，反向再做一次（这次的目的是将那些不能达到1的中间状态全部找到） 比如已经使用了87个模式列表，得到一个tc[.][.],再使用最后一个模式列表的时候，我们应该可以发现，对于很多tc[][]中的数据，根本对数据tc2[720720-1][16-1]没有贡献，对于那些项，我们完全可以将它们的计数置0.这个方向过程可以将大量计数置0。有了这个过程以后，那么接下去穷举过程中，我们可以直接用这些tc2[.][.]数组作为指导，如果穷举到某层，对应的状态达到的计数tc2[.][.]是0，那么对于这个分支，我们就没有必要穷举下去。反之，如果对应计数非0，说明通过这个分支继续下去必然有解，所以这种剪枝方法是非常有效的（只要没有被剪掉的分支，必然继续穷举下去还有解），几乎没有任何浪费。 可惜的是数组这88个tc2[][]数组都很大，几乎不可能保存下来。 不过根据前面的分析，上面计算过程中，实际上我们关心的只是每个tc2[.][.]是否是0，而对于它具体的值并不关心。所以实际上，我们只需要保存88个长度为720720*16的比特向量就可以了。而88个长度为720720*16的比特向量占用的内存空间为88*720720*16/32=31.7M,这个内存空间一点也不大。 所以我们可以发现，这种方法完全可行。 而且实际计算过程中，我们更本不需要计算tc2[.][.],而是直接计算这些比特向量就行了，只要将上面所有用到的代码中+改成|，那么整个动态规划的正向过程一点也不需要改变。同样，对于动态规划的逆向过程，代码没有任何变化，只是操作的对象由整数数组改成比特数组了。 至于这个计算过程的速度如何，我没有实际实现过，无法评估（这个跟实际数据量有关系），但是应该已经是非常有效的方法。

mathe

有谁有兴趣实现一下上面的算法（分别实现计数和搜索的代码）？