Oriented forest TAoCP - Algorithm in python

Question

I try to implement Algorithm O (Oriented forests) from Donald E. Knuth: 'The Art of Computer Programming - Volume 4, Fascile 4, Generating All Trees' on page 24.

My Python solution is:

def generate_oriented_forest(n):
    """Algorithm O from Knuth TAoCP, Fascicle 4, p. 25. """
    p = range(-1, n)
    while True:
       yield p[1:]
       if p[n] > 0: p[n] = p[p[n]]
       else:
           k_largest =  0
           for k in range(1,n): 
               if p[k] != 0: k_largest = k
           k = k_largest
           if k == 0: return
           j =  p[k]
           d = k-j
           if p[k-d] == p[j]: p[k] = p[j]
           else: p[k] = p[k-d] + d
           while k != n:
               #print k, p
               k = k+1
               if p[k-d] == p[j]: p[k] = p[j]
               else: p[k] = p[k-d] + d

if __name__ == "__main__":
    for el in generate_oriented_forest(4):
        print el

    # According to page 23 and also Table 1 p.4 I would expect the sequence:
    # 0123, 0122, 0121, 0120, 0111, 0110, 0101, 0100, 0000

My Implementation gives me:

[0, 1, 2, 3], [0, 1, 2, 2], [0, 1, 2, 1], [0, 1, 2, 0], [0, 1, 1, 1], [0, 1, 1, 0],

[0, 1, 0, 3],

[0, 1, 0, 0], [0, 0, 0, 0].

I'm already looking too long for a bug. Hope someone can point me in the right direction of fix my code. Is my understanding of the algorithm correct? Any improvements on my python style are also appreciated. Thanks for your help.

Answer 1

I refactored your code a little, mostly to eliminate the duplication in step 5.
However the output is still the same as you are getting

def generate_oriented_forest(n):
    """Algorithm O from Knuth TAoCP, Fascicle 4, p. 25. """
    p = range(-1,n)
    while True:
        yield p[1:]
        if p[n] > 0:
            p[n] = p[p[n]]
            continue
        for k in range(n-1,0,-1):
            if p[k] != 0: break
        else:
            break
        j = p[k]
        d = k-j
        while True:
            if p[k-d] == p[j]:
                p[k] = p[j]
            else:
                p[k] = p[k-d] + d
            if k==n:
                break
            k+=1

if __name__ == "__main__":
    for el in generate_oriented_forest(4):
        print el

Answer 2

I do not understand the algorithm and have no idea if my suggested change to the algorithm (below) is correct or not. All I know is that it produces the expected output that you quote for n=4:

def generate_oriented_forest(n):
    """Algorithm O from Knuth TAoCP, Fascicle 4, p. 25. """
    p = range(-1,n)
    while True:
        yield p[1:]
        if p[n] > 0:
            p[n] = p[p[n]]
            continue
        for k in range(n-1,0,-1):
            if p[k] != 0: break
        else:
            break
        j = p[k]
        d = k-j
        while True:
            if p[k-d] == p[j]:
                p[k] = p[j]
            else:
                p[k] = p[k-d]
            if k==n:
                break
            k+=1

I used gnibbler's code as a starting point. I used traceit() and print statements to follow the code as it transitioned from [0, 1, 1, 0] --> [0, 1, 0, 3].

What I find is that this is state of the variables:

[0, 1, 1, 0]  # the last yield
k: 4
d: 2
j: 1
p: [-1, 0, 1, 0, 0]
[0, 1, 0, 3]  # the problem yield

and this is the only time that this code is executed:

__main__:19:             if p[k-d] == p[j]:
__main__:22:                 p[k] = p[k-d] + d

Since p[k-d]=p[2]=1, and you want p[k] to equal 1, I "guess" the correct formula should be p[k]=p[k-d].

I also changed

for k in range(n-1,-1,-1):

to

for k in range(n-1,0,-1):

to stop the code from yielding an extra [0, 0, 0, 0] at the end.

Answer 3

I digged up the original articel: SIAM Journal of Computing, T. Beyer and S.M. Hedetniemi, 'Constant Time generation of rooted trees". The paper explains the algorithm really well. Here's my implementation:

def generate_oriented_forest_2(n): 
    """SIAM J. COMPUT. Vol. 9, No. 4, November 1980
    T. Beyer and S.M. Hedetniemi: constant time generation of rooted trees.""" 

    L = range(-1, n) 
    while True: 
        yield L[1:] 
        if L[n] > 0: L[n] = L[L[n]] 
        else:
            for p in range(n-1, 0, -1): 
                if L[p] != 0: break
            else: break
            for q in range(p-1, 0, -1):
                if L[q] == L[p] - 1: break 
            i = p
            while i <= n:
                L[i] = L[i-(p-q)]
                i += 1

It gives me the expected result. But I'm still wondering why Knuth's algorithm doesn't work. Would be cool to find out.

Answer 4

Congratulations!

It sure looks like you just found an error in TAOCP. As you no doubt know, there's a reward of one hexadecimal dollar (drawn on the Bank of San Seriffe) for the first finder of such errors. I've got one, and I can tell you, it's a hell of thing to put on your wall.

It certainly appears to me that the "+d" in step O5 is in error; at least, I can't find a way to reconcile it with the description of the "cloning" step in the text description before the algorithm. I've checked the most recent Errata for V4f4, and this one is not in there, so it looks like you are the first to notice this.

To verify, I recommend that you calculate the values for n=5 both with and without the "+d", and see which match the expected results. If it goes the way I suspect, write it up and send it to Knuth by email (the address for TAOCP bugs is on his web site) along with your postal address, and you should receive a reply (by paper mail) within 6 months.