用你的方法，仅仅为了计算精度问题，在LCM大于64位整数的情况，可以采用多个互素的整数，
比如取三个互素的32位整数p1,p2,p3,使得p1*p2*p2大于LCM.
然后所有的求和过程改成求和并分别对p1,p2,p3求模（实际上如果只有加法，那么除法运算都不需要）。
当然，这样，搜索过程中和太大和太小本来可以提前剪枝的，这个方法不行。所以可以通过移位方法来做删减。
比如假设LCM比64位整数大8个比特，那么我们可以将所有的数除以256（分别求上取整和下取整）。在求和过程中同时计算这些上取整和下取整的和，如果发现它们越界，也可以提前裁减。
不过我上面介绍的方法，将会采用不同的方法。
原帖由 无心人 于 2008-3-7 13:55 发表 http://bbs.emath.ac.cn/images/common/back.gif
这个题目要求技巧高些而已
还有输出多点
应该超过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 1179414567 秒
按此估计，即使以后每增加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稍后有更好的算法
原帖由 liangbch 于 2008-3-10 14:39 发表 http://bbs.emath.ac.cn/images/common/back.gif
mathe 和 楼主， 能否告诉我， 你们预测的解的个数是多少个，如果解的个数太多，建议降低搜索层数，比如从16减低10或者更低。
我采用简单的剪枝方法（9#中的规则 “i)我们可以用部分和进行剪枝:”），用自定义 ...
我没有真正计算过，不过也感觉解的数目会太多。所以我首先的目标是就算解的数目，而不是将所有解输出。
现在介绍一下我的思路，主要数学依据还是在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};
27应该对应模式应该是63个，而不是64个，上面写错了，下面代码可以打印出所有这些模式:// flter.cpp : Defines the entry point for the console application.
//
#include "stdafx.h"
#define OUTPUT
int mask[]={1,1,2,6,12,60,60,420,840,2520,2520,27720,27720,360360,360360,360360,720720};
int cands[]={17,19,23,29,31,37,41,43,47,53,59,61};
int c2[]={5,7};
int c3[]={1,2,4,5,6,7,2};
#define NUMOF(x)(sizeof(x)/sizeof(x))
int _tmain(int argc, _TCHAR* argv[])
{
int i,j;
for(i=0;i<NUMOF(cands);i++){
int p=cands;
int mt=256/p;
int count=0;
int lm=mask;
for(j=1;j<(1<<mt);j++){
int sum=0;
int k;
for(k=0;k<mt;k++){
if(j&(1<<k)){
sum+=lm/(k+1);
sum%=p;
}
}
if(sum==0){
int start=1;
count++;
#ifdef OUTPUT
for(k=0;k<mt;k++){
if(j&(1<<k)){
if(!start)fprintf(stderr,"+");
fprintf(stderr,"1/%d",(k+1)*p);
start=0;
}
}
fprintf(stderr,"\n");
#endif
}
}
printf("Count %d for prime %d\n",count,p);
}
for(i=0;i<NUMOF(c2);i++){
int p=c2;
int p2=p*p;
int count=0;
int mt=256/p2;
int lm=mask;
for(j=1;j<(1<<mt);j++){
int sum=0;
int k;
for(k=0;k<mt;k++){
if(j&(1<<k)){
sum+=lm/(k+1);
sum%=p;
}
}
if(sum==0){
int start=1;
count++;
#ifdef OUTPUT
for(k=0;k<mt;k++){
if(j&(1<<k)){
if(!start)fprintf(stderr,"+");
fprintf(stderr,"1/%d",(k+1)*p2);
start=0;
}
}
fprintf(stderr,"\n");
#endif
}
}
printf("Count %d for %d\n",count,p2);
}
{
int p=3;
int p3=27;
int mt=7;
int count=0;
int lm=mask;
for(j=1;j<(1<<mt);j++){
int sum=0;
int k;
for(k=0;k<mt;k++){
if(j&(1<<k)){
sum+=lm/c3;
sum%=p;
}
}
if(sum==0){
int start=1;
count++;
#ifdef OUTPUT
for(k=0;k<mt;k++){
if(j&(1<<k)){
if(!start)fprintf(stderr,"+");
if(k==mt-1){
fprintf(stderr,"1/81+1/162");
}else{
fprintf(stderr,"1/%d",c3*p3);
}
start=0;
}
}
fprintf(stderr,"\n");
#endif
}
}
printf("Count %d for 27\n",count);
lm = 1*3*5*7;
count=0;
for(j=1;j<16;j++){
int sum=0;
int k;
for(k=0;k<4;k++){
if(j&(1<<k)){
sum+=lm/(2*k+1);
sum%=2;
}
}
if(sum==0){
int start=1;
count++;
#ifdef OUTPUT
for(k=0;k<4;k++){
if(j&(1<<k)){
if(!start)fprintf(stderr,"+");
fprintf(stderr,"1/%d",(2*k+1)*32);
start=0;
}
}
fprintf(stderr,"\n");
#endif
}
}
printf("Count %d for 32\n",count);
}
return 0;
}我们首先需要将它们转化为我上面说明的结构。
上面程序还有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;memset(tc,0,sizeof(tc));
for(i1=0;i1<1928;i1++)for(i2=0;i2<431;i2++)for(i3=0;i3<260;i3++){
int count=L17.count+L19.count+L25.count;
int sum=L17.value+L19.value+L25.value;
if(0<count&&count<=16&&sum<=720720){
tc++;
}
}
然后对于余下的87个模式列表，我们按任意顺序事先排列好，假设为
ELEMENT *L;
其中每项模式数目为int c;
然后我们同样创建另外一个数组tc2;
后面的程序中，需要频繁交换两个数组的位置，我们可以采用来个指针来轮换它们地位，具体做法这里就不写出了。
现在假设tc数组中已经保存了使用前3个模式列表在加上87个模式列表中t个模式列表(0<=t<87)中的模式使用了k个数达到和x的数目保存在tc中
那么我们可以使用动态规划计算再加上的t个模式后的结果为memcpy(tc2,tc,sizeof(tc2));
for(i=0;i<c;i++){
int c=L.count;
int v=L.value;
for(j=c;j<=16;j++)for(s=v;s<=720720;s++){
tc2+=tc;
}
}这个过程对于所有的t处理完的复杂度可以估计为16*720720*Sum(c[.])=16*720720*(88+8+6+1+1+1+4+3+7+1+77)=2.27G,现在的计算机处理起来应该很快。
全部计算完以后，那么tc2就保存了分别使用1,2,...,16个埃及分数达到1的组合确切数目。
当然上面的过程，如果我们使用了更大的内存，显然不难可以计算出16个以上埃及分数和为1的确切数目（主要问题在于内存使用量）。
现在我们开始考虑另外一个问题，如果我们需要将所有组合都输出，那么使用上面的方法，是否可以提供一种非常有效的方法。
我们想象，如果我们能够将上面计算过程中的所有tc2数组的内容都保存下来（共88个这么大的数组），然后我们再使用动态规划的方法，反向再做一次（这次的目的是将那些不能达到1的中间状态全部找到）
比如已经使用了87个模式列表，得到一个tc[.][.],再使用最后一个模式列表的时候，我们应该可以发现，对于很多tc[][]中的数据，根本对数据tc2没有贡献，对于那些项，我们完全可以将它们的计数置0.这个方向过程可以将大量计数置0。有了这个过程以后，那么接下去穷举过程中，我们可以直接用这些tc2[.][.]数组作为指导，如果穷举到某层，对应的状态达到的计数tc2[.][.]是0，那么对于这个分支，我们就没有必要穷举下去。反之，如果对应计数非0，说明通过这个分支继续下去必然有解，所以这种剪枝方法是非常有效的（只要没有被剪掉的分支，必然继续穷举下去还有解），几乎没有任何浪费。
可惜的是数组这88个tc2[][]数组都很大，几乎不可能保存下来。
不过根据前面的分析，上面计算过程中，实际上我们关心的只是每个tc2[.][.]是否是0，而对于它具体的值并不关心。所以实际上，我们只需要保存88个长度为720720*16的比特向量就可以了。而88个长度为720720*16的比特向量占用的内存空间为88*720720*16/32=31.7M,这个内存空间一点也不大。
所以我们可以发现，这种方法完全可行。
而且实际计算过程中，我们更本不需要计算tc2[.][.],而是直接计算这些比特向量就行了，只要将上面所有用到的代码中+改成|，那么整个动态规划的正向过程一点也不需要改变。同样，对于动态规划的逆向过程，代码没有任何变化，只是操作的对象由整数数组改成比特数组了。
至于这个计算过程的速度如何，我没有实际实现过，无法评估（这个跟实际数据量有关系），但是应该已经是非常有效的方法。
有谁有兴趣实现一下上面的算法（分别实现计数和搜索的代码）？