Sum of proper divisors of every number up to N

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.

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.

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.