Sum of divisors summatory function with Erathosthenes' sieve

Question

I ran across the following problem from an online problem bank: there are up to $~10^5~$ queries each of which asks to compute the sum $$\sum_{k = L}^{R} \sigma(k)$$ where $\sigma(k)$ is the sum of divisors of $k$. It is given that $1 \leq L \leq R \leq 5\cdot 10^6$.

My solution (described below) is based on Erathosthenes's sieve. I've implemented it in C++ and it works in about $0.9$ seconds on average which is too slow. I know that this problem can be solved at least twice faster but don't know how.

So here is my solution (arrays are 0-based):

M = 5 * 1e6
M = array of zeroes of size M + 1
A[1] = 1
for (k = 2; k <= M; k += 1)
    for (j = k; j <= M; j += k)
        A[j] += k

I precalculate $\sigma(k)$ via Erathosthenes' sieve for each $k$ below max possible value. When the main loop reaches $k$, $A[k]$ keeps the value of $\sigma(k)$. Then I reassign $A[k]$ to be $\sum_{i=1}^{k}\sigma(i)$. After such preprocessing all queries can be computed in $O(1)$ time by computing $A[R] - A[L-1]$.

How can I make it faster? I know two formulas: $$ (a) ~~~~~ \sigma(p_{1}^{a_1} \cdots p_{s}^{a_s}) = \prod_{i=1}^{s} \frac{p_{i}^{a_i + 1} - 1}{p_{i} - 1}$$ $$ (b) ~~~~~ \sum_{k=1}^{n} \sigma(k) = \sum_{k=1}^{n} k \left \lfloor \frac{n}{k} \right \rfloor$$

The problem with (a) is that computing it (at least in my implementation) is slower than given above. The problem with (b) is that I don't understand how to compute prefix sum with such approach faster than in $O(n^2)$ time.

Is there a more efficient algorithm for this problem?

(The problem bank credits the original source of the problem as 2012 Kharkiv, Winter School, Day of Sergey Kopelovich, Problem H.)

Answer 1

This isn't really computer science...

You create a table d where you store the sum of the divisors of k, for k = 1 to M, where M = $5 · 10^6$. That's the part that is time critical. Then you create a table s where you store the sum of divisors for all 1 ≤ j ≤ k, for k = 1 to M. That's easy, $s_0 = 0$, $s_{k+1} = s_k + d_{k+1}$. And then f (L, R) = $s_R - s_{L-1}$.

The first table is the problem. You handle this in $O (n \log n)$. And you only need a factor two, you say...

You will have an array d with 5 million entries, likely 4 byte per entry = 20 Megabyte. On a typical processor that you would have in your home computer, 20 Megabyte doesn't fit into any cache. And your code does lots of accesses to elements of that array in quasi random order. For each potential divisor k, you visit all numbers that are divisible by k, and increase the sum of divisors by k.

Let's do that with fewer visits: When you visit j which is divisible by k, add the two divisors k and j/k. But when you do that, start with $j = k^2$, adding only k (because k = j / k, and you don't want to count the divisor twice), and then add k and j / k for further j. You don't need to divide, because j / k will equal k+1, k+2, k+3 etc. We initialise the array for the case k = 1, that is setting A [j] = 1 + j / 1 for j ≥ 2.

A [1] = 1
for (j = 2; j ≤ M; j += 1)
    A [j] = 1 + j

for (k = 2; k*k ≤ M; k += 1)
    j = k*k
    A [j] += k
    j += k
    s = k + (k + 1)
    while j ≤ M
        A [j] += s
        j += k
        s += 1 // s equals k + j / k

You don't save operations. However, you are now accessing the array A in a much more regular pattern, so this will save you time because access to the items will be faster. j will be smaller, making the number of iterations for each j larger, which will make branch prediction work better.

For more improvement, you would find out how many array items fit into the processor cache in your computer, then perform the whole code for subranges of the array only (for example changing only A [0] to A [99999], then changing A [100000] to A [199999], and so on). That way, most memory accesses will only access cache memory, which might be substantially faster.

You are doing N lookups in a table of size M. If M is substantially larger than N, then you should probably think about approaches that don't build this table, and that may be a lot slower per lookup, but faster overall due to the small number of lookups. Even in the case here where N ≤ 100,000 and M = 5,000,000, you might for example not count divisors 1, 2, 3, 4, j/1, j/2, j/3, j/4 in the table (which makes it a bit faster to build), and handle that during the lookup.

Or you could add the sum of the divisors for odd numbers only, then calculate the sum of divisors for even numbers (if the sum of the divisors of an odd k is s, then the sum for 2k is 3s, for 4k it is 7s, for 8k it is 15s etc), which would save almost a factor 2.

PS. I measured it... making the algorithm for counting all sums of divisors more cache friendly by adding both j and k/j doubled the speed. Calculating the sum of divisors for odd k first, then calculating even k from the odd values, makes it a total of 7 times faster. Obviously all just constant factors.

Answer 2

So let me rearrange your problem a bit: using prime sieve should be helpful, but normal Erathostenes sieve is not good enough.

What you need is prime sieve working in linear time, hitting every number just once.
One description of linear time prime sieve shows how to cross every number just once.
What are benefits? Well, if instead of crossing numbers we insert the sum of divisors there, we have fast algorithm placing divisors (please remember about $1$ as a divisor).

Also there is one additional step, prime numbers are not calculated, so encountering one we should write it's divisor as this number + 1.

Next there should be cummulative pass (going through array adding last item to make it sum of all previous divisors).

This way every number should be written exactly once, so this is for sure better than original attempt.

What else could be done?
Since there are less queries than numbers, I thought maybe we can ommit the calculation of the whole array?

This can be done in at least two ways: obvious one is to make partial (or even whole) array offline (not during time measurement), making program bigger, but there was no limit of size.

Another one is to calculate whole array of cumulative divisors, and then fit some functions retreiving results from indices.

The functions itself might be a bit complicated or to make thinks easier we can divide it into ranges - making them shorter and easier to find.
The huge complexity behind that is done offline, and during runtime only queries time matters, since there is no sieve at all.

Answer 3

You may store precalculated results for intervals {L=1, R = k*10^4} and brute-force only about 2*10^4 numbers