c++algorithmc++11mathperfect-numbers

Best way to find 'quite good' numbers up to 1 million?


I am working on an assignment involving 'quite good' numbers. The task describes them as:

"A "quite good" number is an integer whose "badness" – the size of the difference between the sum of its divisors and the number itself – is not greater than a specified value. For example, if the maximum badness is set at 3, there are 12 "quite good" numbers less than 100: 2, 3, 4, 6, 8, 10, 16, 18, 20, 28, 32, and 64; Your task is to write a C++ program, quitegood, that determines numbers of a specified maximum badness that are less than a specified value. The limiting value and maximum badness are specified as command-line arguments when the program is executed."

The task asks me to write a program that prints perfect numbers with a specified badness limit up to a million. So, the command line argument of quitegood 1000000 1 should print 2 4 6 8 16 28 32 64 128 256 496 512 1024 2048 4096 8128 8192 16384 32768 65536 131072 262144 524288.

I have gotten this to work with the following code

#include <iostream>

using namespace std;

int main(int argc, char *argv[]) {

    const int limit = argc > 1 ? atoi(argv[1]) : 1000000;
    const int badness = argc > 2 ? atoi(argv[2]) : 10;

    for(int number = 2; number < limit; number++) {
        int sum = 1;
        for (int factor = 2; factor < number; factor++){
            if (number % factor == 0) {
                sum += factor;
            }
        }

        if (number >= (sum - badness) && number <= (sum + badness)) {
            cout << number << " ";
        }
    }

    return 0;
}

The only issue is that this code is far too slow finding the 'quite good' numbers up to 1 million. Is there any way of optimising this?

Thank you


Solution

  • If f is a factor of n then so is n/f (although when f is the square-root of n, f and n/f are the same factor). So you can make the code a lot faster by counting factors only up to sqrt(number), and then when you find one also include the matching factor number/factor (except for the square-root case).

    for (int factor = 2; factor * factor <= number; factor++){
        if (number % factor == 0) {
            sum += factor;
            if (factor * factor != number) {
                sum += number / factor;
            }
        }
    }
    

    This code runs in 1.554s on my machine in the case of limit being 1 million, and badness 1. I got bored after several minutes waiting for the original code to complete.

    To make the code even faster, you can find the prime factorization of the number, and use the formula for the sum of the divisors based on the prime factorization.

    Even without pre-computing the primes, using this method runs in 0.713s on my machine. Here's my code to compute sum from number:

    int n = number;
    int i = 2;
    while (n > 1) {
        if (i * i > n) {
            sum *= (n + 1);
            break;
        }
        int pp = i;
        while (n % i == 0) {
            pp *= i;
            n /= i;
        }
        sum *= (pp - 1) / (i - 1);
        i += 1;
    }
    sum -= number;
    

    It finds all prime powers that divide number, and for each p^m multiplies sum by (p^(m+1) - 1) / (p - 1). Like the first solution, it stops early, when i*i > n, which at that point means n is a prime.

    It's a lot faster than the first solution in the average case, because although we're still doing trial division, n gets smaller as prime factors are found.

    If you have precomputed a large enough list of primes (that is, it includes at least one larger than the square root of limit), you can be a little more efficient again in computing sum:

    int n = number;
    for (int i = 0; primes[i] * primes[i] <= n; ++i) {
        int pp = primes[i];
        while (n % primes[i] == 0) {
            pp *= primes[i];
            n /= primes[i];
        }
        sum *= (pp - 1) / (primes[i] - 1);
    }
    if (n > 1) sum *= (n + 1);
    sum -= number;
    

    The code with this way of computing sum runs in 0.189s on my machine.