How would you implement a divisor function?

Question

The divisor function is the sum of divisors of a natural number.

Making a little research I found this to be a very good method if you want to find the divisor function of a given natural number N, so I tried to code it in Python:

def divisor_function(n):
    "Returns the sum of divisors of n"
    checked = [False]*100000
    factors = prime_factors(n)
    sum_of_divisors = 1 # It's = 1 because it will be the result of a product
    for x in factors:
        if checked[x]:
            continue
        else:
            count = factors.count(x)
            tmp = (x**(count+1)-1)//(x-1)
            sum_of_divisors*=tmp
            checked[x]=True
    return sum_of_divisors

It works pretty well,but I am sure that it can be improved(e.g. : I create a list with 100000 elements,but I am not using most of them).

How would you improve/implement it?

P.S. This is prime_factors:

def prime_factors(n):
    "Returns all the prime factors of a positive integer"
    factors = []
    d = 2
    while (n > 1):
        while (n%d==0):
            factors.append(d)
            n /= d
        d = d + 1
        if (d*d>n):
            if (n>1): factors.append(int(n));
            break;
    return factors

Answer 1

def sum_divisors(n):
    assert n > 0
    if n == 1:
        return 0
    sum = 1
    if n % 2 == 0:              # if n is even there is a need to go over even numbers
        i = 2
        while i < sqrt (n):
            if n % i == 0:
                sum = sum + i + (n//i)  # if i|n then n/i is an integer so we want to add it as well
            i = i + 1
        if type (sqrt (n)) == int:  # if sqrt(n)|n we would like to avoid adding it twice
            sum = sum + sqrt (n)
    else:
        i = 3
        while i < sqrt (n):     # this loop will only go over odd numbers since 2 is not a factor
            if n % i == 0:
                sum = sum + i + (n//i)  # if i|n then n/i is an integer so we want to add it as well
            i = i + 2
        if type (sqrt (n)) == int:  # if sqrt(n)|n we would like to avoid adding it twice
            sum = sum + sqrt (n)
    return sum

Answer 2

Here is what I do in my Java number utilities (extensively used for Project Euler):

• Generate the prime factorization with explicit exponents (see the answer by Gareth Rees).

• Unfold the prime factorization in the various functions based on it. I.e., use the same algorithm as for prime factorization but directly compute the desire value instead of storing the factors and exponents.

• By default test only divisors two and odd numbers. I have methods firstDivisor(n) and nextDivisor(n,d) for that.

• Optionally precompute a table of least divisors for all numbers below a bound. This is very useful if you need to factorize all or most numbers below the bound (improves speed by about sqrt(limit)). I hook the table into the firstDivisor(n) and nextDivisor(n,d) methods, so this doesn't change the factorization algorithms.

Answer 3

why use dict or set - or count() - at all, when prime_factors() is guaranteed to return the factors in ascending order? You only ever deal with a previous factor. Counting can just be a part of iteration:

def divisor_function(n):
    "Returns the sum of divisors of n"
    factors = prime_factors(n)
    sum_of_divisors = 1 
    count = 0; prev = 0;
    for x in factors:
        if x==prev:
            count += 1
        else:
            if prev: sum_of_divisors *= (prev**(count+1)-1)//(prev-1)
            count = 1; prev = x;
    if prev: sum_of_divisors *= (prev**(count+1)-1)//(prev-1)
    return sum_of_divisors

Answer 4

When computing the sum of divisors, you need the factorization of n in the form p₁^k₁ p₂^k₂ ... — that is, you need the exponent of each prime in the factorization. At the moment you are doing this by computing a flat list of prime factors, and then calling count to work out the exponent. This is a waste of time because you can easily generate the prime factorization in the format you need in the first place, like this:

def factorization(n):
    """
    Generate the prime factorization of n in the form of pairs (p, k)
    where the prime p appears k times in the factorization.

    >>> list(factorization(1))
    []
    >>> list(factorization(24))
    [(2, 3), (3, 1)]
    >>> list(factorization(1001))
    [(7, 1), (11, 1), (13, 1)]
    """
    p = 1
    while p * p < n:
        p += 1
        k = 0
        while n % p == 0:
            k += 1
            n //= p
        if k:
            yield p, k
    if n != 1:
        yield n, 1

Notes on the code above:

1. I've transformed this code so that it generates the factorization, instead of constructing a list (by repeated calls to append) and returning it. In Python, this transformation is nearly always an improvement, because it allows you to consume elements one by one as they are generated, without having to store the whole sequence in memory.

2. This is the kind of function for which doctests work well.

Now computing the sum of divisors is really simple: there's no need to store the set of checked factors or to count the number of times each factor occurs. In fact you can do it in just one line:

from operator import mul

def sum_of_divisors(n):
    """
    Return the sum of divisors of n.

    >>> sum_of_divisors(1)
    1
    >>> sum_of_divisors(33550336) // 2
    33550336
    """
    return reduce(mul, ((p**(k+1)-1) // (p-1) for p, k in factorization(n)), 1)

Answer 5

You need to change two lines only:

def divisor_function(n):
    "Returns the sum of divisors of n"
    checked = {}
    factors = prime_factors(n)
    sum_of_divisors = 1 # It's = 1 because it will be the result of a product
    for x in factors:
        if checked.get(x,False):
            continue
        else:
            count = factors.count(x)
            tmp = (x**(count+1)-1)//(x-1)
            sum_of_divisors*=tmp
            checked[x]=True
    return sum_of_divisors