Project Euler Problem 521 Smallest prime factor - Solution



Smallest prime factor

   

Problem 521

Let smpf(n) be the smallest prime factor of n.
smpf(91)=7 because 91=7×13 and smpf(45)=3 because 45=3×3×5.
Let S(n) be the sum of smpf(i) for 2 ≤ i ≤ n.
E.g. S(100)=1257.
Find S(1012) mod 109.

If you didn't already, read Hint 1 and Hint 2 first. 
Below is an implementation based on ideas from Hint 1 and Hint 2, with some performance tuning in place.

Below c++ code reaches correct answer in ~1 hour.


#include <iostream>
#include <vector>

using namespace std;
int main()
{
    const long long maxp = 1000000;
    const long long maxbillion = 100;
    const long long maxp2 = 10000000000;
    const long long maxarrsize = maxp2 + 1;
    const long long mod = 1000000000;
    std::vector<bool> primemap( maxp + 1, true );
    std::vector<long long> arrp = {};

    long long sum1 = 0;
    for (long long i = 2; i < maxp + 1; i++) {        
        if (primemap[i] == true) {            
            arrp.push_back(i);
            for (long long j = 2; j < maxp / i + 1; j++) {
                if (primemap[i * j] == true) {
                    primemap[i * j] = false;                    
                }
            }
        }
    }
#pragma omp parallel for num_threads(4)
    for (long long i = 0; i < maxbillion; i++) {
        long long localsum = 0;
#pragma omp critical 
        {
            cout << "Iteration: " << i << endl;
        }
        std::vector<bool> mprimemap(maxarrsize, true);

        for (long long j = 1; j < long long(arrp.size()); j++) {
            long long prime = arrp[j];
            long long start = ((i * maxp2) / prime) + 1;
            long long end = ((i + 1) * maxp2) / prime;
            long long count = 0;
            
            for (long long k = start; k <= end; k++) {
                if (k % 2 == 1) {
                    long long sub = prime * k - i * maxp2;
                    if (mprimemap[sub] == true) {
                        mprimemap[sub] = false;
                        count++;
                    }
                }
            }
            localsum = (localsum + prime * count) % mod;            
        }
        for (long long j = 1; j < maxp2 + 1; j+=2) {
            if (mprimemap[j] == true) {
                localsum = (j + localsum) % mod;             
            }
        }
#pragma omp critical 
        {
            sum1 += localsum;
            sum1 %= mod;
        }
    }    
    cout << "Answer: " << sum1 << endl;
    return 0;
}

Comments

Popular posts from this blog

Project Euler Problem 684 Inverse Digit Sum - Solution

Project Euler Problem 686 Powers of Two - Solution

Project Euler Problem 717 Summation of a Modular Formula - Solution