埃及数问题

这个问题是从数学研发网（bbs.emath.ac.cn)看来的，

定义：

设有N个自然数n1, n2, nN, 如果所有这些数的倒数相加，结果等于1, 则称这一组数为一组埃及数。例如：（2,3,6）是一组埃及数，因为 1/2+1/3+1/6=1. 又如， （2, 3, 12, 13, 156）也是一组埃及数， 因为 1/2+1/3+1/12+1/13+1/156==1。 埃及数很少吗？ 非也非也，实际上很多。 比如， 搜索在128以内的， 8个数一组的埃及数， 总共有9855个， 但如果把它们全找出来， 计算量是很大的。 下面这个程序就是查找埃及数的。

 

//

// (c) DGU版权所有， 自由软件， 无限制。

//全文转载请注明出处。

//如有意见或建议,请email gdyxgzy@hotmail.com

//

 

#include <iostream>
#include <vector>
#include <algorithm>
#include <functional>

using namespace std;

static vector<int> primes;

//get all prime numbers less than n, these numbers are used to find GCD
vector<int> get_prime_nums(int n)
{
    vector<int> primes;

    for(int i=2; i<=n; ++i)
        primes.push_back(i);

    int len=primes.size();
    int idx = 0;
    while(idx<len)
    {
        int x=primes[idx];
        for(int i=x+x; i<=n; i+=x)
            primes[i-2] = 0;
        idx++;
        while(idx<len && primes[idx]<=0)
            ++idx;
    }

    primes.erase(remove_if(primes.begin(), primes.end(), bind2nd(less_equal<int>(), 0)), primes.end());
    return primes;
}

int find_gcd(int& up, int& low)
{
    int gcd=1;
    for(int i=0; i<primes.size(); ++i)
    {
        int factor=primes[i];
        if(up<factor)
            break;

        while(up%factor==0 && low%factor==0)
        {
            gcd*=factor;
            up /= factor;
            low /= factor;
        }
    }
    return gcd;
}

vector<vector<int> > egyptian(int Max, int n, int start, int up, int low)
{
    vector<vector<int> >result;
    if(n==1)
    {
        if(low%up==0 && low/up>=start)
            result.push_back(vector<int>(1, low/up));

        return result;
    }

    int end_point = (unsigned long)low*n/up+1;
    if(end_point>Max)
        end_point = Max;

    for(int i=start; i<end_point; ++i)
    {
        int candidate=i;
        if(candidate*up <= low)
            continue;
        int new_up = candidate*up - low;
        int new_low = candidate*low;

        //must > 1/Max
        if(Max*new_up <= new_low)
            continue;

        int gcd = find_gcd(new_up, new_low);

        vector<vector<int> > temp = egyptian(Max, n-1, candidate+1, new_up, new_low);

        if(temp.size()>0)
        {
            for(int i=0; i<temp.size(); ++i)
            {
                vector<int> t = temp[i];
                t.push_back(candidate);
                result.push_back(t);
            }
        }
    }
    return result;
}

int main()
{
    //The total results increases terribly as max_number, and, especially num_parts, increase.
    //Setting max_number to 128 and num_parts to 8 seems a good one, this will give out 9855 results.
    //If max_number=200, and num_parts=5, there will be 57 results.

 

    int max_number = 200;
    int num_parts = 5;

    primes = get_prime_nums(max_number);

    vector<vector<int> > result = egyptian(max_number, num_parts, 2, 1, 1);

    for(int i=0; i<result.size(); ++i)
    {
        vector<int> temp = result[i];
        for(int j=temp.size()-1; j>=0; --j)
            cout<<temp[j]<<", ";
        cout<<endl;
    }
    cout<<"---------------------------------------------------"<<endl;
    cout<<"There're altogether "<<result.size()<<" solutions found."<<endl;
    cout<<"---------------------------------------------------"<<endl;
    cout<<endl;

    return 0;
}