Optimizing Divisor Sieve

Question

Update: The optimizations so far have been excellent. And I've used Python 64 bit to get around the memory issue, but I came to the problem that 50,000,000+ causes my computer to use up more than it's 8 gigs of memory and freeze. I've update the times for 20,00,000 to 40,000,000 for divisorSieve. 40,000,000 capped out at just over 7 gigs of memory. So something will have to be done, perhaps by writing to file or serializing part of the array, to free up memory.

I've never run into this issue before so suggestions are welcome. I'll also be looking into this on my own so if I come up with something, I will update again. I really thank people for their help so far!

---

Original Post (just with results and code updated): I have two sieves that I wrote in python and would like help optimizing them if at all possible. The divisorSieve calculates the divisors of all numbers up to n. Each index of the list contains a list of its divisors. The numDivisorSieve just counts the number of divisors each index has but doesn't store the divisors themselves. These sieves work in a similar way as you would do a Sieve of Eratosthenes to calculate all prime numbers up to n.

Note: divs[i * j].append(i) changed from divs[i * j] += [i] with speed increase thanks to a member over at stackoverflow. I updated the table below with the new times for divisorSieve. It was suggested to use this board instead so I look forward to your input.

def divisorSieve(n):
    divs = [[1] for x in xrange(0, n + 1)]
    divs[0] = [0]
    for i in xrange(2, n + 1):
        for j in xrange(1, n / i + 1, i):
            divs[j].append(i)
    return divs

    def numDivisorSieve(n):
        divs = [1] * (n + 1)
        divs[0] = 0
        for i in xrange(2, n + 1):
            for j in xrange(1, n / i + 1, i):
                divs[j] += 1
        return divs

#Timer test for function
if __name__=='__main__':
    from timeit import Timer
    n = ...
    t1 = Timer(lambda: divisorSieve(n))
    print n, t1.timeit(number=1)

Results:

 -----n-----|--time(divSieve)--|--time(numDivSieve)--
    100,000 | 0.118300055981   | 0.0687864651365
    200,000 | 0.196426572122   | 0.133795795235
    300,000 | 0.279482801518   | 0.20909715743
    400,000 | 0.377094918206   | 0.273495055048
    500,000 | 0.52753871991    | 0.33962928407
    600,000 | 0.566675063756   | 0.406102506687
    700,000 | 0.685204018163   | 0.471622193195
    800,000 | 0.755363562854   | 0.538575604901
    900,000 | 0.85135473036    | 0.620326242458
  1,000,000 | 0.962831374013   | 0.675673633141
  2,000,000 | 1.90850283013    | 1.34772096784
  3,000,000 | 2.83754773901    | 2.03887603409
  4,000,000 | 3.81808421969    | 2.71636126143
  5,000,000 | 4.77563537974    | 3.3678397711
  6,000,000 | 5.71348354793    | 4.0420467127
  7,000,000 | 6.64301299994    | 4.71918143932
  8,000,000 | 7.62517291783    | 5.38929091248
  9,000,000 | 8.60579119239    | 6.06912203769
 10,000,000 | 9.57554985383    | 6.76640791192
 20,000,000 | 16.231225975     |  13.4803200224
 30,000,000 | 24.2593779934    |  20.3057849723
 40,000,000 | 32.8265861168    |  27.0092311999
 50,000,000 | Out of Memory    |  33.7464085339
 60,000,000 | Out of Memory    |  40.5103277975
 70,000,000 | Out of Memory    |  47.2521350376
 80,000,000 | Out of Memory    |  54.1790253157
 90,000,000 | Out of Memory    |  60.716789824
100,000,000 | Out of Memory    |  67.4218930771

Results are pretty good and I'm happy I was able to get it this far, but I'm looking to get it even faster. If at all possible, I'd like to get 100,000,000 at a reasonable speed with the divisorSieve. Although this also brings into the issue that anything over 19,000,000+ throws a MemoryError at divs = [[1] for x in xrange(0, n + 1)]) in divisorSieve. numDivisorSieve does allow the full 100,000,000 to run. If you could also help get past the memory error, that would be great.

I've tried replacing numDivisorSieve's divs = [1] * (n + 1) with both divs = array.array('i', [1] * (n + 1)) and divs = numpy.ones((n + 1), dtype='int') but both resulted in a loss of speed (slight difference for array, much larger difference for numpy). I expect that since numDivisorSieve had a loss in efficiency, then so would divisorSieve. Of course there's always the chance I'm using one or both of these incorrectly since I'm not used to either of them.

I would appreciate any help you can give me. I hope I have provided enough details. Thank you.

Answer 1

You can use xrange's third param to do the stepping for you to shave off a little bit of time (not huge).

Changing:

for j in xrange(1, n / i + 1):
    divs[i * j].append(i)

To:

for j in xrange(i, n + 1, i):
    divs[j].append(i)

For n=100000, I go from 0.522774934769 to 0.47496509552. This difference is bigger when made to numDivisorSieve, but as I understand, you're looking for speedups in divisorSieve

Answer 2

EDIT: map(lambda s: s.append(i) , [divs[ind] for ind in xrange(i, n + 1, i)]) Seems to be ~0.2% faster ~2 times slower than Adam Wagner's (for n=1000000)

The infamous 'test the unit test' problem.

Answer 3

The following offers a very very small improvement to divisorSieve and a good improvement to numdivisorSieve. But the factors will not be sorted inside each list. For example the factors list of of 16 will be [4, 2, 8, 1, 16].

def divisorSieve(n):
    divs = [[] for j in xrange(n + 1)]
    nsqrt = int(sqrt(n))
    for i in xrange(1, nsqrt + 1):
        divs[i*i].append(i)
        for j in xrange(i, i*i, i):
            divs[j].append(j/i) #If j/i is replaced by i, a good improvement is seen. Of course, that would be wrong.
            divs[j].append(i)
    for i in xrange(i+1, n+1):
        for j in xrange(i, n+1, i):
            divs[j].append(j/i)
            divs[j].append(i)
    return divs

def numdivisorSieve(n):
    divs = [1] * (n + 1)
    divs[0] = 0
    nsqrt = int(sqrt(n))
    for i in xrange(2, nsqrt + 1):
        divs[i*i] += 1
        for j in xrange(i, i*i, i):
            divs[j] += 2
    for i in xrange(i+1, n+1):
        for j in xrange(i, n+1, i):
            divs[j] += 2
    return divs

Unfortunately, modifying this definition to create two lists divsmall ([4,2,1]) and divlarge ([8,16]) and in the end doing divsmall[j].reverse(); divsmall[j].extend(divlarge[j]); return divsmall makes it slightly slower than the original.

Also, I think it makes more sense for divs[0] to be [] instead of [0]