埃及分数问题

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. //
   
3. 
   
4. #include "stdafx.h"
   
5. #define OUTPUT
   
6. 
   
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.     }
    
46. 
    
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.     }
    
79. 
    
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

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