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.

    &gt;&gt;&gt; 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) &amp;= (1 + 3 + 5 + 15) + (2 + 6 + 10 + 30) + (4 + 12 + 20 + 60) \\ &amp;= (1 + 3 + 5 + 15) + 2(1 + 3 + 5 + 15) + 2^2(1 + 3 + 5 + 15) \\ &amp;= (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) &amp;= (1 + 5) + (3 + 15) \\ &amp;= (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.

    &gt;&gt;&gt; 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 &lt; 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.