4
\$\begingroup\$

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.

\$\endgroup\$
7
\$\begingroup\$

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.

\$\endgroup\$
5
\$\begingroup\$

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.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.