Sum of proper divisors of every number up to N

Question

I wrote a function which provides the sum of proper divisors of every number up-to N (N is a natural number) and I want to improve its performance.

For example, if N = 10 then the output is:

[0, 0, 1, 1, 3, 1, 6, 1, 7, 4, 8]

This is my proper divisor sum function:

def proper_divisor_sum(N):
    N += 1
    divSum = [1] * N
    divSum[0], divSum[1] = 0, 0
    for i in range(2, N//2+1):
        for j in range(i+i, N, i):
            divSum[j] += i
    return divSum

I am certain that there is an algorithmic/mathematical optimization for this problem but I am having trouble thinking about it because I am not good at math.

Answer 1

1. Analysis

The code in the post loops over \$2 ≤ i ≤ {n\over 2}\$ and then over the multiples of \$i\$ below \$n\$:

for i in range(2, N//2+1):
    for j in range(i+i, N, i):
        div_sum[j] += i + j//i if i != j//i else i  # (A)

This means that the line (A) is executed at most $$ {n\over 2} + {n\over 3} + {n\over 4} \dotsb + {n \over n/2} $$ times, that is $$ n\left({1\over 2} + {1\over 3} + {1\over 4} \dotsb {1\over n/2}\right) $$ which is \$n(H_{n/2} - 1)\$ (where \$H_n\$ is the \$n\$th harmonic number), and this is \$Θ(n \log n)\$.

2. Improvement

The first thing to note is that if we had a function that computed the sum of all the divisors, then we could compute the sum of proper divisors by subtracting the number itself:

def sum_proper_divisors(n):
    """Return list of the sums of proper divisors of the numbers below n.

    >>> sum_proper_divisors(10) # https://oeis.org/A001065
    [0, 0, 1, 1, 3, 1, 6, 1, 7, 4]

    """
    result = sum_divisors(n)
    for i in range(n):
        result[i] -= i
    return result

Let's call the sum-of-all-divisors function \$σ\$, and consider for example \$σ(60)\$: $$ σ(60) = 1 + 2 + 3 + 4 + 5 + 6 + 10 + 12 + 15 + 20 + 30 + 60. $$ Collect together the divisors by the largest power of 2 in each divisor: $$ \eqalign{σ(60) &= (1 + 3 + 5 + 15) + (2 + 6 + 10 + 30) + (4 + 12 + 20 + 60) \\ &= (1 + 3 + 5 + 15) + 2(1 + 3 + 5 + 15) + 2^2(1 + 3 + 5 + 15) \\ &= (1 + 2 + 2^2)(1 + 3 + 5 + 15).} $$ Now, the remaining factor \$1 + 3 + 5 + 15\$ is \$σ(15)\$ so we can repeat the process, collecting together the divisors by the largest power of 3: $$ \eqalign{σ(15) &= (1 + 5) + (3 + 15) \\ &= (1 + 3)(1 + 5).}$$ And so $$ σ(60) = (1 + 2 + 2^2)(1 + 3)(1 + 5). $$ This is obviously connected to the fact that \$ 60 = 2^2·3·5 \$. In general, if we can factorize \$n\$ as: $$ n = 2^a·3^b·5^c\dotsb $$ then $$ σ(n) = (1 + 2 + \dotsb + 2^a)(1 + 3 + \dotsb + 3^b)(1 + 5 + \dotsb + 5^c)\dotsm $$ These multipliers occur many times, for example \$(1 + 2 + 2^2)\$ occurs in the sum of divisors of every number divisible by 4 but not by 8, so it's most efficient to sieve for the sums of many divisors at once, like this:

def sum_divisors(n):
    """Return a list of the sums of divisors for the numbers below n.

    >>> sum_divisors(10) # https://oeis.org/A000203
    [0, 1, 3, 4, 7, 6, 12, 8, 15, 13]

    """
    result = [1] * n
    result[0] = 0
    for p in range(2, n):
        if result[p] == 1: # p is prime
            p_power, last_m = p, 1
            while p_power < n:
                m = last_m + p_power
                for i in range(p_power, n, p_power):
                    result[i] //= last_m    # (B)
                    result[i] *= m          # (B)
                last_m = m
                p_power *= p
    return result

3. Analysis

The lines marked (B) are executed at most $${n\over 2} + {n\over 2^2} + \dotsb + {n\over 3} + {n\over 3^2} + \dotsb $$ times, that is $$n\left({1 \over 2} + {1 \over 2^2} + \dotsb + {1\over 2^{\lfloor \log_2 n\rfloor}} + {1 \over 3} + {1 \over 3^2} + \dotsb + {1\over 3^{\lfloor \log_3 n\rfloor}} + \ldots\right)$$ which is less than $$n\left({1\over 2-1} + {1\over3-1} + {1\over5-1} + \dotsb + {1\over p-1}\right)$$ where \$p\$ is the largest prime less than or equal to \$n\$, and this is \$Θ(n \log \log n)\$ (see divergence of the sum of the reciprocals of the primes) and so asymptotically better than the code in the post.

Answer 2

Non-Performance Changes

Note that I haven't really made that many style or other recommendations as you didn't ask for any. However, I did still change two things. First, I avoided modifying an inputted variable, as this could have unintended consequences, and second I changed the name of divSum to div_sum.

Performance Changes

I noticed a couple of optimizations that you could make to your code to improve performance. The first change I made after noticing that there was a bit of unnecessary repetition in certain areas (Please forgive me for being slightly vague as it's currently 4 in the morning). Given any divisor of a number and that number itself, a second divisor can be calculated. I used this information to optimize the code to this form.

def proper_divisor_sum(N):
     div_sum = [1] * (N+1)
     div_sum[0], div_sum[1] = 0, 0

     for i in range(2, int(math.sqrt(N))+1):
         for j in range(i*i, N+1, i):
             div_sum[j] += i + j//i if i != j//i else i
     return div_sum

This code is roughly 33% faster than the original function. I originally thought that I needed to have an if-statement in there to handle perfect squares and the double-counting of those squares.

I then used this code to analyze the states of each individual iteration. This is what I found:

[[0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  1, 1,  1, 1],
 [0, 0, 0, 0, 2, 0, 5, 0, 6, 0, 7, 0, 8, 0, 9, 0, 10, 0, 11, 0, 12],
 [0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 7, 0, 0, 8, 0,  0, 9,  0, 0],
 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4,  0, 0,  0, 9]]

What I noticed was that the leading coefficients of each row progressed linearly, but the coefficients that followed obeyed a pattern. I'm sure you can work out the other patterns I saw in these numbers for yourself, but the primary thing this told me was that my previous assumption (about having to handle perfect squares) was actually a slip-up, as I had at some point forgotten that a number can't have more than one square, or I had at least accounted for the possibility in my code (Don't write code at 4 AM). Realizing this I was able to simplify the code down into its current form.

def proper_divisor_sum(N):
     div_sum = [1] * (N+1)
     div_sum[0], div_sum[1] = 0, 0

     for i in range(2, int(math.sqrt(N))+1):
         div_sum[i*i] += i
         for j in range(i*i+i, N+1, i):
             div_sum[j] += i + j//i
     return div_sum

This final function is roughly 40% faster than your original function.