2
\$\begingroup\$

I have written a solution to problem 23 of Project Euler, and after doing some testing and editing, I've managed to get the run time down from 1.7443 seconds to 1.56616 seconds. The problem is

Find the sum of all the positive integers which cannot be written as the sum of two abundant numbers.

What are some other ways that I could improve this program? I'm pasting it entirely, and it should be noted that I'm compiling with g++ -o3.

#include <chrono>
#include <cmath>
#include <iostream>
#include <ctime>
#include <ratio>

int main()
{
//Value given in the question
const int LIMIT = 28124;
//Used to measure the runtime
std::chrono::high_resolution_clock::time_point t1 = std::chrono::high_resolution_clock::now();

//Sieve of Erasthotenes to remove a significant number of iterations in the next step
int rootMax = (int)sqrt((double)LIMIT);

bool set[LIMIT+1] = {0};

for(int i = 2; i <= rootMax; ++i)
{
    if(!set[i])
    {
        for(int j = (i*i); j <= LIMIT; j += i)
        {
            set[j] = true;
        }
    }
}

bool isSumOfAb[LIMIT] = {0};

int abundant[6965];
int count = 0;
for(int i = 12; i < LIMIT; ++i) //This can start at 12, as we know that 12 is the 1st abundant number
{
    if(!set[i]) // Prime numbers won't be abundant, so next number
    {
        continue;
    }

    int sum = 0;
    for(int j = 1; j < i; ++j)
    {
        if(i % j == 0)
        {
            sum += j;   
        }
    } 

    if(sum > i)
    {
        abundant[count] = i;
        ++count;
    }
}   

//Lots of work to do here.
//Establish all of the sums and compare
int a, c;
for(int i = 0; i < 6965; ++i)
{
    a = abundant[i];

    for(int j = i; j < 6965; ++j)
    {
        c = a + abundant[j];
        //Big speed increase when made like this as opposed to `if(c <= LIMIT) isSumOfAb[c] = true;`
        if(c > LIMIT)
        {
            break;  
        }else
        {
            isSumOfAb[c] = true;
        }
    }
}

int sum = 0;
for(int i = 24; i < LIMIT; ++i) //can start at 24, as that is the 1st sum of 2 abundant numbers
{
    if(!isSumOfAb[i])
    {
        sum += i;
    }
}

std::cout << "Sum: " <<  sum;
std::chrono::high_resolution_clock::time_point t2 = std::chrono::high_resolution_clock::now();
std::chrono::duration<double> time_span = std::chrono::duration_cast<std::chrono::duration<double>>(t2 - t1);
std::cout << "\nIt took me " << time_span.count() << " seconds.";
}
\$\endgroup\$
  • \$\begingroup\$ The best solution to this problem is to pre-compute the result (because it does not change with input). Then the application you submit just prints the number. \$\endgroup\$ – Martin York Sep 8 '14 at 20:50
  • \$\begingroup\$ @LokiAstari of course the best solution is to pre-compute it, but that's missing the point. I'm iffy about knowing the number of numbers so that a static array can be used. The aim of the program is to compute it, not print it out, otherwise a lot of project euler's program solutions would print out wolfram-alpha.com \$\endgroup\$ – Yann Sep 8 '14 at 20:56
2
\$\begingroup\$

After some tests (adding more time-points), we can find that the double for-loop used to find all abundant numbers is the slowest. Now we can think about speeding it up, by shortening the inner loop. We start with sum=1 and try every number we find can divide i upto sqrt(i). We immediatelly know that the product will of course divide the number as well, but need to take special care for numbers of the form i=j*j.

//  stop is the biggest integer such that stop*stop <= i
//  the reasoning behind using it is similar to the reasoning for using it in first phase
//  when marking all primes, but we need to take care for div=i/j when i%j == 0
    int sum = 1, stop = sqrt(i);
    for(int j = 2;;)
    {
        int div = i/j;         // divison and modulo will probably be optimized
        int mod = i%j;         // as one instruction (we could use some divmod built-in)

        if(mod == 0)           // when this condition is met we know that:
        {                      // i = j*div
            if(j == div)
            {                  // but can be i = j*j when j=div
                sum += j;      // we need to sum it only once in such case (e.g. 16=4*4)
                break;         // and can exit the loop, because j=stop here
            }
            sum += j+div;      // this is the tricky part when we sum both j and div
        }                      // e.g. for 15: 1 + 3+5 (j=3, div=5, stop=3 here)
        if(++j > stop)
            break;             // moving the condition here is another optimization
    }                          // ...and could be rewritten as do-while loop

The original code needed more than 8.8s on my old PC, but the current version needs only 0.11s. Full code:

#include <chrono>
#include <iostream>
#include <algorithm>

using the_clock = std::chrono::high_resolution_clock;
using time_point = the_clock::time_point;
using time_span = std::chrono::duration<double>;

// all integers greater than 28123 can be written as the sum of two abundant numbers
static const int LIMIT = 28124;
// smallest abundant number is 12
static const int minAbundant = 12;

// all prime numbers will have false in this array
static bool not_prime[LIMIT+1]; // zero initialized
static const int maxPrimeRoot = (int)sqrt(LIMIT);
static inline bool is_prime(int x) {
    return !not_prime[x];
}

// all abundant numbers
static const int abundantSz = 6965; // we know that, see final output
static int abundant[abundantSz];

// all numbers that are sum of abundant numbers will be marked true here
static bool isSumOfAb[LIMIT];

int main()
{
//  1. Mark all primes:
    time_point t1 = the_clock::now();

    for(int i = 2; i <= maxPrimeRoot; i++)
    {
        if(is_prime(i))
        {
            for(int j = (i*i); j <= LIMIT; j += i)
                not_prime[j] = true;
        }
    }

//  2. Find all abundant numbers
    time_point t2 = the_clock::now();

    abundant[0] = minAbundant;
    int count = 1;
    for(int i = minAbundant+1; i < LIMIT; ++i)
    {
        if(is_prime(i)) // Prime numbers won't be abundant, so next number
            continue;

    //  BIG CHANGE HERE
    //  we sum both the number we find can divide i
    //  and the product + stop after sqrt(i)
    //  we need to take special care for numbers of the form i*i (e.g. 16=4*4)
        int sum = 1, stop = sqrt(i);
        for(int j = 2;;)
        {
            int div = i/j;
            int mod = i%j;
            if(mod == 0) {
                if(j == div) {
                    sum += j;
                    break; }
                sum += j+div;
            }
            if(++j > stop)
                break;
        }

        if(sum > i)
        {
            abundant[count] = i;
            ++count;
        }
    }

//  3. Mark are numbers that are sum of abundant numbers
    time_point t3 = the_clock::now();

    for(int i = 0; i < count; ++i)
    {
        int a = abundant[i];

        for(int j = i; j < count; j++)
        {
            int c = a + abundant[j];
            if(c > LIMIT)
                break;
            isSumOfAb[c] = true;
        }
    }

//  4. Sum them all
    time_point t4 = the_clock::now();

    int sum = 0;
    for(int i = 24; i < LIMIT; ++i) //can start at 24, as that is the 1st sum of 2 abundant numbers
    {
        if(!isSumOfAb[i])
            sum += i;
    }

//  Output
    time_point t5 = the_clock::now();
    using std::cout;
    cout << "Sum: " <<  sum;
    cout << "\nIt took me " << time_span(t5-t1).count() << " seconds.";
    cout << "\nTimes: " << time_span(t2 - t1).count()
        << ", " << time_span(t3 - t2).count()
        << ", " << time_span(t4 - t3).count()
        << ", " << time_span(t5 - t4).count();
    cout << "\nAbundant count: " << count; // just a check

    if(sum != 4179595)
        cout << "\nERROR - bad sum\n";
}

You can see I have moved few variables outside main() and added few using to reduce the code a bit. But I am really not the one to talk about code-styling :D

\$\endgroup\$
  • \$\begingroup\$ I'll definitely have a proper look at this tomorrow and likely accept as an answer. I'd assumed that the double for loop was the slowest bit, but I couldn't figure out how to fix it \$\endgroup\$ – Yann Sep 8 '14 at 20:01
  • \$\begingroup\$ After a while, and asking a friend, I understand why the magical inner loop works and you only need to go up to sqrt(i). I think that it'd be helpful for others if you could edit in an explanation, if possible? \$\endgroup\$ – Yann Sep 8 '14 at 20:58
  • \$\begingroup\$ I hope to find some time tomorrow (already in bed using smartphone). The sqrt is like with the prime numbers, we know it is the point where all bigger values have already been processed, because we used sum += j+div. I will try to describe it better, tomorrow. 15&16 should be good examples to show how it works: 1 + 3+5, 1 + 2+8 + 4, if I am correct. \$\endgroup\$ – user52292 Sep 8 '14 at 21:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.