Divisors

When I wrote the first post of the FFT series I confessed it was out of boredom from number theory. That statement got me into a coffee break discussion with a friend and mathematician, who claimed that the FFT is as number theory as it gets, and went on to explain why. Since I usually don’t understand a word of what he is saying when he gets started on one of this rants, which invariably lead to Riemann’s zeta function, I just pretended to acknowledge his doubts, and then forgot about the conversation altogether.

Well, a few weeks down the road it seems he was right: while trying to make my FFT code work, I have found myself dealing with all sorts of number theoretical stuff, ranging from simple ones, like factorization, to arcane, as the Chinese remainder theorem. And since I have this thing about complete expositions, I’m afraid that this interleaving of apparently unrelated number theory posts with the FFT ones is going to be the norm. As a matter of fact the odd post on the Euclidean algorithm was supposed to be a stepping stone to an improved radix 2 FFT, which I finally managed to do without.

But lets get going with today’s subject: computing the divisors of a positive integer. May not seem a very sexy topic at first glance, but there’s beauty in it, as the striking Divisor Plot website put together by Jeffrey Ventrella proves…

Getting into the algorithmic treatment of the problem, there is of course the dumber approach, on the line of “to compute all divisors of a number, do trial division by all numbers smaller than it.” I will not pursue that line of thinking though, and will focus on computing the divisors once the prime factorization of the number is known, which is the efficient way of dealing with this problem.

So I will assume that we know the expansion of the number into its prime factors as

\displaystyle n = \prod_{i=1}^{\omega(n)} p_i^{a_i},

where the function ω(n) represents the number of distinct prime factors, pi, each appearing with the corresponding exponent, ai.

A divisor will have a subset of the prime factors of n as its own factors, so we could represent any divisor as a list of ω(n) numbers, the i-th ranging from 0 up to and including ai. This scheme makes it very easy to calculate the total number of divisors as

\displaystyle d(n)=\prod_{i=1}^{\omega(n)} (a_i+1).

Unordered Generation of Divisors
Furthermore, it points out that the computation of all the divisors could be achieved using ω(n) nested loops. Of course the number of loops to nest will vary from input to input, so recursion is probably the best of choices, or so thinks Tim Peters, as the code below has only cosmetic modifications to what he posted on the Python-list@python.org mailing list:

def divisors(factors) :
    """
    Generates all divisors, unordered, from the prime factorization.
    """
    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

This approach, while extremely fast, is not much use if you want to produce the divisors sorted, as you will then have to first generate and store all divisors, then sort them.

Ordered Generation of Divisors
If what you are after is not the full set of divisors, or any arbitrary subset of them, but something like “the first divisor larger than x” or “the last divisor smaller than x,” it may be worth trying other algorithms which, while being less efficient, manage to produce the divisors properly sorted.

There’s little literature on the subject, so its hard to find a valid benchmark . Even PARI/GP, which is the golden standard on computational number theory seems, from the little sense I have managed to make out of their Subversion repository, to be generating all divisors (in a similar manner to above) and then sorting them.

So here’s a calculation scheme that can be effectively modified to yield ordered divisors… We start with a pivot value of 1, and an ordered list of all prime factors. Before every factor we split our flow in two:

  • one branch skips the factor, and proceeds to the first factor different from the current one, keeping the pivot unchanged
  • the other branch first updates the pivot, multiplying it by the factor, yields the result, and then proceeds to the next factor, with the modified pivot.

The flow is probably easier to understand with the example below, showing the branching when computing the divisors of 540 = 22·33·5. Blue arrows skip to the first different factors, red ones go to the next factors, a line going through a factor implies multiplication by it…

The branching of the algorithm to generate all the divisors of 540 in order.

To keep track of this structure, I’m going to get fancy, and will build a tree-like data structure, where each item will be a list, holding the factor in the first position, and pointers, i.e. nested lists, to the next elements in the second and third ones. Again, to make things clear, here’s the data structure for 540, using the same color code as before…

The data structure holding the prime factors of 540, arranged for ordered calculation of the divisors.

Rather than traversing this structure recursively, we will approach it in a more structured way, by keeping a separate list of factors that can be visited in a next step. In the example above, when we start out, we can visit the first factor (a 2), or skip to the third (the first 3), so the list will hold two elements. If we select going to the first 2, we then have two new options, i.e. going on to the second 2, or taking the shortcut to the first three. Note that the two instances point at the first three are different, since once has a pivot value of 1 (and will thus generate the divisor 3), while the other has it set to 2 (an will generate the divisor 6).

If we don’t care about the other of divisors, this method can even be implemented by hand, by keeping a list of pivots and remaining factors, and on each step erasing the first element of the list and generating a divisor, plus adding up to two new items at the bottom. The first steps, plus the resulting list of divisors, are worked out below, again for the case of n = 540.

Working out the divisors of 540

Now, while when doing it on paper there is not much more choice other than processing the items sequentially, we could choose to process always the element that produces the smallest divisor, and thus produce them in order. The simplest implementation of this idea would simply use a sort on the list of processable items before yielding the next divisor…

def divisors_ordered(factors) :
    """
    Generates all divisors, ordered, from the prime factorization
    """
    factors =  [(j,factors.count(j)) for j in sorted(set(factors))]
    omega = len(factors)
    fcts =[factors[0][0],[],[]]
    aux = [fcts, None]
    for j in xrange(omega) :
        p,a = factors[j]
        if omega - j > 1 :
            aux[1] = [factors[j+1][0],[],[]]
        else :
            aux[1] = None
        for k in xrange(1, a) :
            b = [p,[],[]]
            aux[0][1] = b
            aux[0][2] = aux[1]
            aux[0] = b
        aux[0][1] = aux[1]
        aux[0][2] = aux[1]
        aux[0] = aux[1]
    div_tree = [[fcts[0],1, fcts]]
    yield 1
    while div_tree :
        div_tree.sort()
        divisor, pivot, pointer = div_tree.pop(0)
        yield divisor
        if pointer[1] :
            div_tree += [[divisor * pointer[1][0] , divisor, pointer[1]]]
        if pointer[2] :
            div_tree += [[pivot * pointer[2][0], pivot, pointer[2]]]

Of course doing a full sort before yielding every divisor kills the efficiency of this algorithm. The fancy (or should I simply say correct?) way of keeping an ordered list, popping the smallest item only, and pushing arbitrary elements, is with a priority queue, which I will implement as a binary heap. Since we will be building our heap from scratch, we only need to worry about push and pop functions.

def push_heap(item, heap = None, compare = cmp) :
    """
    Pushes an item into a min-heap, or creates one.
    """
    heap_size = None
    if not heap :
        heap = [item, 1]
    else :
        heap_size = heap.pop()
        heap += [item, heap_size + 1]
    i = heap[-1] - 1
    val = heap[i]
    heap_ordered = False
    while i > 0 and not heap_ordered :
        n_i = i - 1
        if n_i & 1 : # faster than n_i % 2
            n_i -= 1
        n_i >>= 1 # faster than n_i // 2
        n_val = heap[n_i]
        if compare(val, n_val) < 0 :
            heap[i], heap[n_i] = heap[n_i], heap[i]
            i, val = n_i, heap[n_i]
        else :
            heap_ordered = True
    return heap

def pop_heap(heap, compare = cmp) :
    """
    Pops the smallest item from a min-heap.
    """
    heap_size = heap.pop()
    if not heap_size:
        return None
    elif heap_size == 1 :
        return heap.pop()
    ret, heap[0] = heap[0], heap.pop()
    heap_size -= 1
    heap += [heap_size]
    i = 0
    heap_ordered = False
    while not heap_ordered :
        n_i = 2 * i + 1
        if n_i < heap_size :
            if n_i + 1 < heap_size and compare(heap[n_i],heap[n_i + 1]) > 0 :
                n_i += 1
            if compare(heap[i], heap[n_i]) > 0 :
                heap[i], heap[n_i] = heap[n_i], heap[i]
                i = n_i
            else :
                heap_ordered = True
        else :
            heap_ordered = True
    return ret

def divisors_ordered(factors) :
    """
    Generates an ordered list of all divisors
    """
    factors =  [(j,factors.count(j)) for j in sorted(set(factors))]
    omega = len(factors)
    fcts =[factors[0][0],[],[]]
    aux = [fcts, None]
    for j in xrange(omega) :
        p,a = factors[j]
        if omega - j > 1 :
            aux[1] = [factors[j+1][0],[],[]]
        else :
            aux[1] = None
        for k in xrange(1, a) :
            b = [p,[],[]]
            aux[0][1] = b
            aux[0][2] = aux[1]
            aux[0] = b
        aux[0][1] = aux[1]
        aux[0][2] = aux[1]
        aux[0] = aux[1]
    yield 1
    div_heap = push_heap([fcts[0],1, fcts])
    while div_heap :
        divisor, pivot, pointer = pop_heap(div_heap)
        yield divisor
        if pointer[1] :
            div_heap = push_heap([divisor * pointer[1][0] , divisor, pointer[1]], div_heap)
        if pointer[2] :
            div_heap = push_heap([pivot * pointer[2][0], pivot, pointer[2]], div_heap)

The following graph shows the time required by each of the three implementations presented to produce an ordered list of all divisors of pn#, i.e. the product of the first n prime numbers. Do notice how the lower computational complexity of the heap sorting eventually beats the initial speed advantage of the sorting version. Still, for full lists of divisors, generating all of them and then sorting beats any other approach. It is also worth pointing out how the time required for the heap version of the ordered algorithm eventually settles around being 14 times slower than the unordered generator with delayed sorting, i.e. the computational complexities of both are the same, which probably means that sorting dominates the time of the unordered algorithm.

Time to compute all ordered divisors of the first few primorials.

Just to show why we could want to use an algorithm that is more than ten times slower, the following graph shows the time needed to compute the first divisor larger than 1000 for the first few primorials. While using the heap version of the ordered algorithm time is close to constant, the unordered algorithm shows the same behavior as in the previous test.

Time to compute the first divisor larger than 1000 for the first few primorials.

It is nevertheless worth noting that the break even point is p11# = 200,560,490,130, so unless you are after the first few divisors of a huge number, the first algorithm presented should always be the choice.

10 Responses to Divisors

  1. Steve says:

    In the first graph ‘unordered then sort’ is faster than ‘ordered with heap’. Could that be because Python sort is implemented in C whereas your heap is implemented in Python, so it is telling us that C is faster than Python – which we knew already!

    Python supplies a heap / priority queue in the ‘heapq’ module. It is implemented in C and is, I assume, quite fast.

  2. Jaime says:

    I don’t think so, Steve… The algorithm to generate divisors in a sortable manner is slower than the first one, even if you avoid the sorting altogether. If you extract the potential divisors from the end and not the beginning of the list, it seems to have the same computational complexity as the first one, but runs somewhere between 1.5 and 2 times slower.

    Using heapq instead of my own, things improve, specially at the lower runtime end, but overall for large inputs the performance is still 3 to 5 times slower than sorting after generation, not during.

    I nevertheless feel that using data structures without building them with your bare hands is close to cheating… ;-)

    Despite the better performance, it’s funny that the break even point to generate the first divisor larger than 1000 is not altered much: it is now around p9# = 223,092,870, which is still a huge number for the practical purposes I have in mind.

  3. Beetle B. says:

    Python supplies a heap / priority queue in the ‘heapq’ module. It is implemented in C and is, I assume, quite fast.

    Apparently not on all systems. I once read a blog post where someone used the heapq module for sorting, and I couldn’t duplicate his results. I emailed the author, and he discovered that it worked well on one computer of his, but the results were quite different in the other (which was also running a different OS).

    Turned out one had the module compiled into C. The other was plain Python.

  4. Steve says:

    Beetle B,
    I presume it was because the versions of Python were different.

    Python 2.3 introduced the new heapq module, implemented in Python.
    In Python 2.4 the heapq module was converted to C resulting in a tenfold improvement in speed.

  5. Beetle B. says:

    I doubt that was the reason, as we were both using 2.5. I was on Linux, and he was, I believe on a Mac.

  6. [...] repeating) list of primes. If I’d been paying slightly more attention, he’s got a lengthy post about computing the number of divisors for an integer as [...]

  7. This is a very useful algorithm, and should come in handy in solving some projecteuler problems. Thanks for the great explanation!

    Could you tell us what program you use to draw those beautiful-looking diagrams (of the linked structures)? They look very natural. Or did you draw them by hand?

  8. behindthefall says:

    Is there a mixup in the coloring and labeling of the lines in the two plots, that is, does BLUE (ordered with sort) in the next-to-last become ORANGE (unordered then sort) in the last, while GREEN (ordered with heap) stays the same? I ask this because the blue line in the next-to-last plot and the orange line in the last both have the hockey-stick shape and both cross the green line at about p-sub-11-#. I see that the elapsed times are different in the two problems. Just surprised that orange should have acquired such a pronounced, late upwards break … Your text is consistent with the labeling in the last plot, so I guess if there is a mixup it occurred in the next-to-last plot.

  9. wjjw says:

    Nice and interesting post! Both recursive and priority queue based version are beutifull. Will definietly use it from time to time :)

  10. wjjw says:

    Hello again.

    I rewritten (unsorted version) it into C (and added few other optimalisation). I needed divisors for determining an order of element a in finite field Z/pZ. (order of element a is smallest positive number k, such that a^k mod p == 1, and I know it is true for k==p-1, via Fermat’s theorem, but it is not nacassarly smallest k. so need to check also divisors of p-1, because if a^k == a^(k1*k2) == (a^k1)^k2 == 1^k2 , then we see k1 (a divisor of k), or its divisor, can be an order of a )

    You can find it at http://pastebin.com/tMBTYSx7
    It is small, fast and interesting code I think.

    I will also implement ordered heap to check if it faster, but for my purpose it can be beneficial to also use identity a^(x*y) = (a^x)^y first, as it will be faster to perform powpod probably incrementally, than checking from smallest ones (but early termination in ordered heap will be really nice for me also).

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

%d bloggers like this: