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

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 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 //

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:

• Sieve-like method (12 ms)

• Trial division using list of primes (222 ms)

• Original post method (527 ms)

• 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 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);
    }
}

StackExchange Code Review Q#120572, answer score: 3