Speeding up an algorithm for finding the number of distinct lines

Question

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)

• Danica: Um, that first for loop will go over 10 ** 10 / 2 values; even if you could get each call to aux down to a millisecond that'd take two months. You'll need a different algorithm.
• Wooble: You probably want to use xrange instead of range so you don't have a list of 5 billion integers. Not that it will help all that much.
• Simon Bergot: a(n) = Sum(moebius(k)*((floor(n/k)+1)^3-1), k=1..n) Don't know if it is correct though
• Joel Cornett: use sqrt(n) instead of n**0.5. It's a bit faster
• Joel Cornett: Also, if you are using python 2.x, integer division is the default. There is no need to type cast k/2 as int()

Answer 1

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 & 1 == 1):
 f = 3
 step = 2
 while(f <= 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 & 1 == 0):
 lst.append(2)
 num /= 2
 while(num & 1 == 0):
 lst.append(2)
 num /= 2
 factor = 3
 maxFactor = sqrt(num)
 while num>1 and factor <= 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 > 1):
 lst.append(num)
 return lst

For efficiency reasons you could also make generator functions.

• Osama Gamal: I tried it .. It's 20seconds slower than mine for solving N=10**6 .. I'm not sure why
• Joel Cornett: @OsamaGamal: There's an extra pow() operation at the beginning of this one. Also, the while loop is slower than iterating over the elements of range() or xrange().
• matcauthon: Thats kind of weird. The first pow operations is just for checking for perfect squares. But the while loop must be faster, due to the heavily reduced interations. I will have a look at it if i found some time
• matcauthon: Ok. Found another algorithm based on factorization. Will updating the post...
• Osama Gamal: That's better. Now it is running a little bit faster for 10**6. I just had to update the line number 5 in the aux(k) function to: divdict[k] = sorted(list(divisors(k)))