Project Euler Problem 87 (Python和C++代码实现和解析)

Problem 87 : Prime power triples

The smallest number expressible as the sum of a prime square, prime cube, and prime fourth power is 28. In fact, there are exactly four numbers below fifty that can be expressed in such a way:

28 = 2² + 2³ + 2⁴
33 = 3² + 2³ + 2⁴
49 = 5² + 2³ + 2⁴
47 = 2² + 3³ + 2⁴

How many numbers below fifty million can be expressed as the sum of a prime square, prime cube, and prime fourth power?

1. 欧拉项目第87道题 质数幂三元组

最小的可以表示为一个质数的平方，加上一个质数的立方，再加上一个质数的四次方的数是28。实际上，在小于50的数中，一共有4个数满足这一性质：

28 = 2² + 2³ + 2⁴
33 = 3² + 2³ + 2⁴
49 = 5² + 2³ + 2⁴
47 = 2² + 3³ + 2⁴

在小于五千万的数中，有多少个数可以表示为一个质数的平方，加上一个质数的立方，再加上一个质数的四次方？

2. 求解分析

求解这道题有几个点：第一，使用质数筛子获取质数列表，第二, 一个质数的平方，立方和四次方的幂之和要小于5千万，第三，符合条件的幂之和有可能重复，需要剔除重复的，但如果使用集合，会比较耗时。

第一步：我们找到最大的质数p的质数列表 [2,3, …, p]（p² + 2³ + 2⁴ < 50000000)
第二步：枚举判断质数三元组(p₁, p₂, p₃) 且 p₁² + p₂² + p₃² < 50000000
第三步：剔除重复的幂之和， 统计不重复的幂之和的个数。

3. Python 代码实现

Python 代码

# PE0087.py
import math

def createPrimeSieve(n):
    """ create prime sieve and return primes list """
    sieve = [True] * n
    sieve[0] = sieve[1] = False

    for i in range(2, int(n**0.5)+1):
        if True == sieve[i]:
            for j in range(i*i, n, i):
                sieve[j] = False

    return [ i for i in range(1, n) if True==sieve[i] ]
    
def getPrimesList(max_number):
    '''
    find primes list which sum of a prime square, prime cube, 
    and prime fourth power is below fifty million
    '''
    max_prime = int(math.sqrt(max_number - 2**3 - 2**4))
    primes_list = createPrimeSieve(max_prime+1)
    return primes_list

def getNumOfPrimePowerTriples(max_number):
    '''
    max_number is fifty million
    '''
    primes_powers_sum_list = [ 0 ] * max_number
    primes_list = getPrimesList(max_number)
    for prime3 in primes_list:
        prime_fourth_power = prime3**4
        for prime2 in primes_list:
            prime_cube = prime2**3
            sum32 = prime_fourth_power + prime_cube
            if sum32 >= max_number:
                break
            for prime1 in primes_list:
                prime_square = prime1**2
                sum321 = sum32 + prime_square
                if sum321 >= max_number:
                    break
                else:
                    primes_powers_sum_list[sum321] = 1
   
    return sum(primes_powers_sum_list)
                            
def main():
    total = getNumOfPrimePowerTriples(50*10**6)
    print(total, "numbers below fifty million can be expressed as the")
    print("sum of a prime square, prime cube, and prime fourth power.")
    
if  __name__ == '__main__':
    main()

４. C++ 代码实现

C++采用了求解分析的方法，首先通过质数筛子找到所有小于 50万的平方根的质数列表，保存质数列表的时候使用STL vector或者申请内存来保存，然后枚举判断满足条件的质数三元组，当然统计的时候，要剔除重复的幂之和。

实际结果表明：使用STL vector或者申请内存来保存质数，运行时间相差不大。

Ｃ++代码　(编译开关为　USE_STL_VECTOR)

#include <iostream>
#include <vector>
#include <cmath>    // pow()
#include <ctime>    // clock()
#include <cstring>  // memset()

using namespace std;

//#define USE_STL_VECTOR

class PE0087
{
private:
    static const int max_number_int = 50000000; // 50 million
    bool *m_primesSieve;

    vector<int> primes_vec;
    int  *m_primesArray = nullptr;
    int  m_numOfPrimes;

    void createPrimeSieve(int max_n);

    int getPrimesList();

public:
    PE0087() { getPrimesList(); };
    ~PE0087() { delete[] m_primesArray; };

    int getNumOfPrimePowerTriples();
    int getNumOfPrimePowerTriples1();
};

// create a prime Sieve of Eratosthenes below max value
void PE0087::createPrimeSieve(int max_n)
{
    m_primesSieve[0] = m_primesSieve[1] = false;

    for (int i = 2; i < (int)sqrt((double)max_n); i++)
    {
        if (true == m_primesSieve[i])
        {
            for (int j = i * i; j < max_n; j += i)
            {
                m_primesSieve[j] = false;
            }
        }
    }
}

int PE0087::getPrimesList()
{
    int max_square_root = (int)sqrt((double)max_number_int);

    m_primesSieve = new bool[max_square_root + 1];
    memset(m_primesSieve, true, (max_square_root + 1) * sizeof(bool));

    createPrimeSieve(max_square_root);

#ifndef USE_STL_VECTOR
    m_primesArray = new int[max_square_root + 1];
    m_numOfPrimes = 0;
#endif

    for (int n = 2; n < max_square_root; n++)
    {
        if (true == m_primesSieve[n])
        {
#ifdef USE_STL_VECTOR
            primes_vec.push_back(n);
#else
            m_primesArray[m_numOfPrimes++] = n;
#endif
        }
    }

    delete[] m_primesSieve;

    return 0;
}

int PE0087::getNumOfPrimePowerTriples1()
{
    double max_number = (double)max_number_int;

    int *PrimePowerTriples = new int[max_number_int + 1];
    memset(PrimePowerTriples, 0, sizeof(int)*(max_number_int + 1));

    for (int i = 0; i < m_numOfPrimes; i++)  // primes vector
    {
        double prime_square = pow((double)m_primesArray[i], 2);
        for (int j = 0; j < m_numOfPrimes; j++)  // prime cube
        {
            double prime_cube = pow((double)m_primesArray[j], 3);

            if (prime_square + prime_cube >= max_number)
            {
                break;
            }

            for (int k = 0; k < m_numOfPrimes; k++) //prime fourth power
            {
                double number = prime_square + prime_cube + pow((double)m_primesArray[k], 4);
                if (number > 0 && number < max_number)
                {
                    PrimePowerTriples[(int)number] = 1;
                }
                else
                {
                    break;
                }
            }
        }
    }

    int total = 0;
    for (int n = 1; n < max_number_int; n++)
    {
        if (1 == PrimePowerTriples[n])
        {
            total++;
        }
    }

    delete[] PrimePowerTriples;

    return total;
}

int PE0087::getNumOfPrimePowerTriples()
{
    double max_number = (double)max_number_int;

    int *PrimePowerTriples = new int[max_number_int];
    memset(PrimePowerTriples, 0, sizeof(int)*max_number_int);

    for (auto prime1 : primes_vec)  // primes vector
    {
        double prime_square = prime1 * prime1;
        for (auto prime2 : primes_vec)  // prime cube
        {
            double prime_cube = pow((double)prime2, 3);

            if (prime_square + prime_cube >= max_number)
            {
                break;
            }

            for (auto prime3 : primes_vec) //prime fourth power
            {
                double number = prime_square + prime_cube + pow((double)prime3, 4);
                if (number < max_number)
                {
                    PrimePowerTriples[(int)number] = 1;
                }
                else
                {
                    break;
                }
            }
        }
    }

    int total = 0;
    for (int n = 0; n < max_number_int; n++)
    {
        if (1 == PrimePowerTriples[n])
        {
            total++;
        }
    }

    delete [] PrimePowerTriples;

    return total;
}

int main()
{
    clock_t start = clock();

    PE0087 pe0087;

#ifdef USE_STL_VECTOR
    int total = pe0087.getNumOfPrimePowerTriples();
#else
    int total = pe0087.getNumOfPrimePowerTriples1();
#endif

    cout << total << " below fifty million can be expressed as the sum" << endl;
    cout << "of a prime square, prime cube, and prime fourth power" << endl;

    clock_t finish = clock();
    double duration = (double)(finish - start) / CLOCKS_PER_SEC;
    cout << "C/C++ running time: " << duration << " seconds" << endl;

    return 0;
}