euler0347.cc

// https://projecteuler.net/problem=347

#include <vector>
#include <cmath>
#include <iostream>

std::vector getSieveOfEratosthenes(const int& n)
{
    if (n <= 0) return {};

    std::vector sieve(n, false);
    int i, j, inc;
    double end = sqrt(n);

    if (n >= 2) sieve[2] = true;

    for (i = 3; i <= n; i += 2)
        sieve[i] = true;

    for (i = 3; i <= end; i += 2)
        if (sieve[i])
        {
            inc = i * 2;
            for (j = i * i; j <= n; j += inc)
                sieve[j] = false;
        }

    return sieve;
}

std::vector getPrimes(const int& n)
{
    std::vector isPrime = getSieveOfEratosthenes(n);

    std::vector primes;
    int s = isPrime.size();

    for (int i = 0; i < s; i++)
        if (isPrime[i])
            primes.push_back(i);

    return primes;
}

unsigned long maxPrimeDividend(const unsigned long& limit)
{
    if (limit <= 0) return 0;

    std::vector primes = getPrimes(limit/2);

    unsigned long sum = 0, maxProduct, product, current;
    int n = primes.size(), i, j;
    double end = sqrt(limit);

    for (i = 0; i < n; i++) // first prime
    {
        if (primes[i] > end) break;

        for (j = i + 1; j < n; j++) // second prime
        {
            product = primes[i] * primes[j]; // initialization: i^1 * j^1

            if (product > limit) break;

            maxProduct = 0;  // largest (k, h) pair such that i^k * j^h <= limit
            do
            {
                current = product;

                // find maximum exponent for j
                while (current * primes[j] <= limit) current *= primes[j];

                if (maxProduct < current) maxProduct = current;

                product *= primes[i]; // increment i's exponent by one
            } while (product <= limit);

            // sum all maximum products
            sum += maxProduct;
        }
    }

    return sum;
}

int main()
{
    std::cout << "S(10e6) = " << maxPrimeDividend(10e6) << std::endl;
    return 0;
}