Speeding up an algorithm for finding the number of distinct lines Is there a way to speed up the following python code:

auxdict = {0:0}
divdict = {}

def divisors(n):
    divisors = []
    for i in range(1, int(n**0.5) + 1):
        if (n % i == 0):
            divisors.extend([i, n/i])
    divisors = list(set(divisors))
    divisors.sort()
    return divisors

def aux(k):
    if k in auxdict:
        return auxdict[k]
    else:
        if not k in divdict:
            divdict[k] = divisors(k)
        r = sum(map(aux, divdict[k][:len(divdict[k])-1]))
        auxdict[k] = (k+1)**3 - k**3 - r
        return auxdict[k]

def lines(k):
    result = 0
    for i in range(int(k/2),0,-1):
        result += int(k/i-1) * aux(i)
    return (k+1)**3 - result - 1

for i in [10**6, 10**10]:
    print i, lines(i)


It is a modified version of a Mathematica code I found. The code finds the maximum number of distinct lines in a 3D cube of lattices of size n: http://oeis.org/A090025

The problem is the code is very slow and doesn't work for larger values. I want to evaluate lines(10**10)
You can improve your divisors function:

Version 1: (not includes n)

def divisors(n):
    """Return the divisors of n (including 1)"""
    if(n == 1):
        return (0)

    divs = set([1])
    r = floor(sqrt(n))
    if(r**2 == n):
        divs.add(r)

    f = 2
    step = 1
    if(n &amp; 1 == 1):
        f = 3
        step = 2

    while(f &lt;= r):
        if(n % f == 0):
            divs.add(n/f)
            divs.add(f)
        f = f + step

    return sorted(divs)


Version 2: (includes n) from here

def divisors2(n) :
    """
    Generates all divisors, unordered, from the prime factorization.
    """
    factors = factorI(n)
    ps = sorted(set(factors))
    omega = len(ps)

    def rec_gen(n = 0) :
        if n == omega :
            yield 1
        else :
            pows = [1]
            for j in xrange(factors.count(ps[n])) :
                pows += [pows[-1] * ps[n]]
            for q in rec_gen(n + 1) :
                for p in pows :
                    yield p * q

    for p in rec_gen() :
        yield p


where factorI is:

def factorI(n):
    """Returns the factors of n."""
    lst = []
    num = int(n)
    if(num &amp; 1 == 0):
        lst.append(2)
        num /= 2
        while(num &amp; 1 == 0):
            lst.append(2)
            num /= 2

    factor = 3
    maxFactor = sqrt(num)
    while num&gt;1 and factor &lt;= maxFactor:
        if(num % factor == 0):
            num /= factor
            lst.append(factor)
            while(num % factor == 0):
                lst.append(factor)
                num /= factor
            maxFactor = sqrt(num)
        factor += 2

    if(num &gt; 1):
        lst.append(num)
    return lst


For efficiency reasons you could also make generator functions.