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 mssigma(200000) =  496062  //   79 mssigma(300000) =  984312  //  108 mssigma(400000) =  996030  //  125 mssigma(500000) = 1230453  //  143 mssigma(500000) = 1230453  //    0 ms (total: 502 ms)summing sigma(0)..sigma(500000) ... 205617099435 (504 ms)verifying       1 to   25000 ...    1321 msverifying   25001 to   50000 ...    3659 msverifying   50001 to   75000 ...    7553 ms... verifying  450001 to  475000 ...  280243 msverifying  475001 to  500000 ...  310217 msThe 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.0sigma(100000) =  246078  //   19.7 mssigma(200000) =  496062  //   42.7 mssigma(300000) =  984312  //   56.9 mssigma(400000) =  996030  //   68.6 mssigma(500000) = 1230453  //   78.9 mssigma(500000) = 1230453  //    0.0 ms (total: 267.3 ms)summing sigma(0)..sigma(500000) ... 205617099435 (0.2 ms)verifying       1 to   25000 ...    1025.2 msverifying   25001 to   50000 ...    3053.0 msverifying   50001 to   75000 ...    5054.7 msverifying  450001 to  475000 ...   37406.9 msverifying  475001 to  500000 ...   39903.4 msThe 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?
Computing the divisor sum in bulk without division, multiplication or factorisation (SPOJ DIVSUM)
c#;c++;performance;programming challenge;primes
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.