Computing the divisor sum in bulk without division, multiplication or factorisation (SPOJ DIVSUM)

Question

SPOJ problem #0074 DIVSUM - Divisor Summation requires computing the sums of proper divisors in bulk (200000 test cases, time limit 3 seconds):

Given a natural number n (1 <= n <= 500000), please output the summation of all its proper divisors.

Definition: A proper divisor of a natural number is the divisor that is strictly less than the number.

e.g. number 20 has 5 proper divisors: 1, 2, 4, 5, 10, and the divisor summation is: 1 + 2 + 4 + 5 + 10 = 22.

The sum of proper divisors is also called the aliquot sum (OEIS A001065); it can be computed from the sum of all divisors (sigma function, OEIS A000203) by removing the number itself from the sum.

In the MathExchange topic Is there a formula to calculate the sum of all proper divisors of a number? - a treasure trove for SPOJ DIVSUM - I found an answer that mentioned an unusual recursive formula for the sigma function.

This formulation requires neither division nor multiplication or factorisation, and its structure makes it ideal for a fast cache-based implementation. And so I set about doing exactly that, because the idea of coding yet another table-based bulk factoriser seemed boring in comparison (I'm on a 'learn C#' tour of the challenge sites ATM).

The mathematical notation and explanation can be inspected in the linked answer, so here's the formula in C#:

class SPOJ_0074_DIVSUM_v1
{
    static int[] sign_sequence = { 1, 1, -1, -1 };

    static int sigma (int n)
    {
        if (n <= 1)
            return n < 1 ? 0 : 1;

        int[] gpn_index = { 1, -1 };

        for (int sum = 0, i = 0; ; ++i)
        {
            int k = gpn(gpn_index[i & 1]);

            if (k >= n)
                return sum + (k == n ? n * sign_sequence[i & 3] : 0);

            sum += sigma(n - k) * sign_sequence[i & 3];

            gpn_index[0] += i & 1;
            gpn_index[1] -= i & 1;
        }
    }

    // overflow if n outside the range [-26754, 26755] but DIVSUM needs only a fraction of that
    static long gpn (long n)
    {
        return (3 * n * n - n) >> 1;
    }
}

The function gpn() computes the generalised pentagonal numbers that are needed for the recursion, i.e. the k values.

The code above nicely shows the structure of the computation but it is exceedingly slow because of the heavy, un-memoised recursion (in fact, it is almost as good at causing stack overflows as the Fibonacci function). It is so slow that it takes several minutes just for going through the two digit numbers, and without memoisation it probably takes hours to blow the stack.

Besides memoisation, an efficient implementation needs to unroll the loop to cover the natural rythm of two sequences: the 1, -1, 2, -2, ... of the gpn() parameter, and the 1, 1, -1, -1, ... for the summation sign.

The last required bit is strength reduction on the gpn() calls. Subtracting gpn(n) from gpn(n + 1) gives a delta of 3 * n + 1, which itself requires a reduction to get rid of the multiplication which then becomes delta += 3.

The value for gpn(-n) is simply gpn(n) + n, of course.

Without further ado, here's the final code that achieves SPOJ-worthy speeds:

using System;
using System.Diagnostics;

class SPOJ_0074_DIVSUM_v3
{
    public const int LIMIT = 500000;

    static int[] sigma_cache;
    static int   sigma_limit;

    static SPOJ_0074_DIVSUM_v3 ()
    {
        sigma_cache = new int[LIMIT + 1];
        sigma_cache[0] = 0;
        sigma_cache[1] = 1;
        sigma_limit = 1;
    }

    public static int sigma (int n)
    {
        if (n > sigma_limit)
            extend_sigma_cache(n);

        return sigma_cache[n];
    }

    static void extend_sigma_cache (int new_limit)
    {
        Trace.Assert(sigma_cache != null && new_limit < sigma_cache.Length);
        Trace.Assert(sigma_limit >= 1);

        for (int n = sigma_limit + 1; n <= new_limit; ++n)
        {
            int gpn_index = 1;
            int gpn_value = 1;
            int gpn_delta = 1;
            int sum = 0, k;

            for ( ; ; )
            {
                if ((k = gpn_value) >= n) break;

                sum += sigma_cache[n - k];

                if ((k += gpn_index) >= n) break;

                sum += sigma_cache[n - k];

                ++gpn_index; gpn_value += (gpn_delta += 3);

                if ((k = gpn_value) >= n) break;

                sum -= sigma_cache[n - k];

                if ((k += gpn_index) >= n) break;

                sum -= sigma_cache[n - k];

                ++gpn_index; gpn_value += (gpn_delta += 3);
            }

            if (k == n)
                sum += (gpn_index & 1) != 0 ? n : -n;

            sigma_cache[n] = sum;
        }

        sigma_limit = new_limit;
    }
}

I chose to have the code extend the cache as needed instead of using one-shot initialisation, because that makes it easier to do experiments (see timings later). It also allows to control the initialisation from outside the class - by the simple expedient of a call like sigma(25000) - without having to expose various knobs and levers.

Here are the timings taken on my Haswell laptop. The first couple of lines show a couple of 'cold' sigma calls, which give the timings for the successive filling of the empty cache. The last call is satisfied from the cache, naturally, and I added it to avoid the misleading impression that a single call might normally take such a long time.

Then comes a line that shows the time for calling the sigma function for n ranging from 0 to the limit (500000) - done as a summation to get a handy checksum - to see what the cost of 'hot' calls is.

After that you can see the verification of the sigma results against the simplest possible (i.e. easiest to verify and most difficult to get wrong) implementation of the sigma function based on trial division. Most of those lines have been snipped because the first and last few tell all that needs telling - which is mostly that simple approaches based on trial division are way too slow for this challenge.

### SPOJ_0074_DIVSUM_v3.cs ###

sigma(100000) =  246078  //   39 ms
sigma(200000) =  496062  //   79 ms
sigma(300000) =  984312  //  108 ms
sigma(400000) =  996030  //  125 ms
sigma(500000) = 1230453  //  143 ms
sigma(500000) = 1230453  //    0 ms (total: 502 ms)

summing sigma(0)..sigma(500000) ... 205617099435 (504 ms)

verifying       1 to   25000 ...    1321 ms
verifying   25001 to   50000 ...    3659 ms
verifying   50001 to   75000 ...    7553 ms
... 
verifying  450001 to  475000 ...  280243 ms
verifying  475001 to  500000 ...  310217 ms

The initialisation times aren't bad at all (same ballpark as the C++ version) but the time for reading the cached values from the array (longer than the total initialisation time) beggars belief. That's a whole microsecond per access, whereas C++ is about 1000 times as fast.

Here are the results for compiling the code as C++ (mutatis mutandis, e.g. std::vector<unsigned> instead of int[] and so on):

*** VC++ 14.0 *** RTTI _CPPUNWIND 190023506.0

sigma(100000) =  246078  //   19.7 ms
sigma(200000) =  496062  //   42.7 ms
sigma(300000) =  984312  //   56.9 ms
sigma(400000) =  996030  //   68.6 ms
sigma(500000) = 1230453  //   78.9 ms
sigma(500000) = 1230453  //    0.0 ms (total: 267.3 ms)

summing sigma(0)..sigma(500000) ... 205617099435 (0.2 ms)

verifying       1 to   25000 ...    1025.2 ms
verifying   25001 to   50000 ...    3053.0 ms
verifying   50001 to   75000 ...    5054.7 ms

verifying  450001 to  475000 ...   37406.9 ms
verifying  475001 to  500000 ...   39903.4 ms

The combined time for full initialisation and 500000 sigma calls of the C# version amounts to about 1 second, which is sufficient for SPOJ although it does not leave a lot of reserve for the not-quite-lightning I/O of the framework. SPOJ clocked 2.02 seconds for the solution with I/O via int.Parse(Console.ReadLine()) and Console.WriteLine(); the timeout is 3 seconds.

The time for accessing the cache can only be called disappointing, though. The C# time actuall refers to naked array access instead of calls to the wrapper function but there's virtually no difference. So I guess that the array bounds check-o-matic must be slowing things down tremendously...

Here's the offending code, as per mjolka's suggestion:

var t = new Stopwatch();
t.Start();
// ...
t.Stop();
Console.WriteLine("(total: {0} ms)", t.ElapsedMilliseconds);

Console.Write("\nsumming sigma(0)..sigma({0}) ... ", sigma_limit);
long sum = 0;
t.Start();  // <--
for (int i = 0; i <= sigma_limit; ++i)
    sum += /** /sigma(i)/*/sigma_cache[i]/**/;
t.Stop();
Console.WriteLine("{0} ({1} ms)\n", sum, t.ElapsedMilliseconds);

I've compared this to mjolka's code and marked a notable point of difference with // <--.

Are there other lessons regarding this code that experienced C# knights can teach to a C# rookie/padawan?

Answer 1

I ran the code and found that summing the values sigma(0) ... sigma(500000) took 1 or 2ms, so there was a discrepancy of about 502ms compared to the timing posted.

Fortunately you included the time taken to fill the cache:

sigma(500000) = 1230453  //    0 ms (total: 502 ms)

My guess is that the timer used for summing the values is mistakenly including the time taken to fill the cache.

Answer 2

Can be faster using "sieve"

If you iterate through each number and in a sieve like fashion add itself to each of its multiples, you can generate all 500000 sigmas in a fraction of the time of any other method. It only takes 4 lines of code:

    for (int i = 1; i <= LIMIT; i++) {
        for (int j = i; j <= LIMIT;j += i) {
            sigma_cache[j] += i;
        }
    }

The theoretical running time of this solution is \$O(n \log n)\$. The theoretical time of the original post solution is approximately \$O(n^{1.5})\$ because for each sigma you compute, you sum a series of approximately \$\sqrt n\$ values.

Full implementation and timing results

I used your test code and inserted the sieve-like method, and compared it to the method in the original post. Here is the full test code:

using System;
using System.Diagnostics;

class SPOJ_0074_DIVSUM_v3
{
    public const int LIMIT = 500000;

    static void Main(string[] args)
    {
        var t = new Stopwatch();
        t.Start();

        // Use sieve-like method to compute sum of divisors:
        int[] sigma_cache = new int[LIMIT + 1];
        for (int i = 1; i <= LIMIT; i++) {
            for (int j = i; j <= LIMIT;j += i) {
                sigma_cache[j] += i;
            }
        }

        t.Stop();

        Console.WriteLine("(total: {0} ms)", t.ElapsedMilliseconds);

        Console.Write("\nsumming sigma(0)..sigma({0}) ... ", LIMIT);
        long sum = 0;
        t.Restart();
        for (int i = 0; i <= LIMIT; ++i)
            sum += sigma_cache[i];
        t.Stop();
        Console.WriteLine("{0} ({1} ms)\n", sum, t.ElapsedMilliseconds);
    }
}

And here are the results:

Sieve-like method          :  12 ms
Original post method       : 527 ms

Other methods

I also tested some other methods. I found that the order of speed from fastest to slowest was:

1. Sieve-like method (12 ms)
2. Trial division using list of primes (222 ms)
3. Original post method (527 ms)
4. Trial division using all numbers (1042 ms)

With the trial division by primes method, the theoretical running time is \$O(n * numPrimes(\sqrt n))\$.

Since the sieve is better than trial division using primes, I won't really explain the primes method other than to just give the code:

using System;
using System.Diagnostics;

class SPOJ_0074_DIVSUM_v3
{
    public const int LIMIT = 500000;

    static int [] primes = {
          2,     3,     5,     7,    11,    13,    17,    19,    23,    29,
         31,    37,    41,    43,    47,    53,    59,    61,    67,    71,
         73,    79,    83,    89,    97,   101,   103,   107,   109,   113,
        127,   131,   137,   139,   149,   151,   157,   163,   167,   173,
        179,   181,   191,   193,   197,   199,   211,   223,   227,   229,
        233,   239,   241,   251,   257,   263,   269,   271,   277,   281,
        283,   293,   307,   311,   313,   317,   331,   337,   347,   349,
        353,   359,   367,   373,   379,   383,   389,   397,   401,   409,
        419,   421,   431,   433,   439,   443,   449,   457,   461,   463,
        467,   479,   487,   491,   499,   503,   509,   521,   523,   541,
        547,   557,   563,   569,   571,   577,   587,   593,   599,   601,
        607,   613,   617,   619,   631,   641,   643,   647,   653,   659,
        661,   673,   677,   683,   691,   701,
    };

    static void Main(string[] args)
    {
        var t = new Stopwatch();
        t.Start();

        int[] sigma_cache = new int[LIMIT + 1];
        int numPrimes = primes.Length;

        // Compute sigmas using trial division by primes.
        for (int i=1; i<=LIMIT; i++) {
            int num = i;
            int result = 1;
            for (int j=0; num > 1 && j<numPrimes; j++) {
                int prime = primes[j];
                if ((num % prime) == 0) {
                    int factor = 1 + prime;
                    int curVal = prime;

                    num /= prime;
                    while ((num % prime) == 0) {
                        num /= prime;
                        curVal *= prime;
                        factor += curVal;
                    }
                    result *= factor;
                }
            }
            if (num != 1) {
                result *= num + 1;
            }
            sigma_cache[i] = result;
        }

        t.Stop();

        Console.WriteLine("(total: {0} ms)", t.ElapsedMilliseconds);

        Console.Write("\nsumming sigma(0)..sigma({0}) ... ", LIMIT);
        long sum = 0;
        t.Restart();
        for (int i = 0; i <= LIMIT; ++i)
            sum += sigma_cache[i];
        t.Stop();
        Console.WriteLine("{0} ({1} ms)\n", sum, t.ElapsedMilliseconds);
    }
}