I've updated PADS, my Python algorithm implementation library, to include

• The condensation of a directed graph.


• A bit-parallel algorithm for testing reachability in a directed acyclic graph; after a linear number of bit-vector operations, one can test reachability between any pair of vertices in constant time.


• A 2-satisfiability solver that uses reachability in the condensation of an implication graph to determine which variables have values that are fixed to true or fixed to false in all solutions to a given 2SAT instance.


• An optional rule for my Sudoku solver that transforms the puzzle into a 2SAT instance with at least as many solutions, and uses the fixed variables of this instance to make inferences about the puzzle.

( Some details )

The new smarter Sudoku solver still needs to backtrack on some puzzles, but many fewer of them than before. Here's one that it still finds difficult:

 ----------------------------------- 
| .   .   . | 9   .   2 | .   .   . |
|           |           |           |
| 8   .   . | .   4   . | .   .   3 |
|           |           |           |
| .   4   . | 7   .   . | 1   .   . |
|-----------------------------------|
| .   .   6 | .   .   9 | .   .   1 |
|           |           |           |
| .   2   . | .   .   . | .   7   . |
|           |           |           |
| 5   .   . | 3   .   . | 4   .   . |
|-----------------------------------|
| .   .   7 | .   .   3 | .   6   . |
|           |           |           |
| 6   .   . | .   1   . | .   .   5 |
|           |           |           |
| .   .   . | 2   .   6 | .   .   . |
 ----------------------------------- 

Ars Technica: Firefox to get massive JavaScript performance boost.

Mozilla Blog: What can you do when your browser is 7 times faster?.

The technology behind these speedups is trace-based optimization, a method of performance prediction that can use the information about which parts of some code have been important already to concentrate compiler-optimization power on making the same code run faster the next time. The researchers behind this are UCI computer scientist Michael Franz and his post-doc Andreas Gal.

Personally, I'm not so thrilled about JavaScript, despite my browser using it all the time on various web sites. But I'm interested in this for another reason: I do a lot of my programming in Python, and sometimes when I do I feel like I need to apologize for Python's poor performance compared to C. But this technology should work just as well in Python as in JavaScript, and gives me hope that eventually my Python code will be as fast as if I had taken the extra effort to write in a compiled and optimized language.

I've updated my library of Python algorithm implementations to include another module: SortedSet. It maintains a set of items (like a Python set data type), but adds a method to report all the items in the set in sorted order.

The obvious ways of doing this are to sort everything at reporting time or to maintain the items in a balanced binary tree data structure. My solution combines the simplicity of the first solution with the asymptotic speed of the second: store, along with the set itself, the sorted output from the last report and smaller lists of changed items. A new report request can be handled by sorting the smaller lists and merging them with the previous output, in constant time for each previously listed item and logarithmic time for each changed item.

Some recent uploads to arxiv.org that have caught my attention:

Lens sequences, by J. Kocik, arXiv:0710.3226. ( Read more... )

Generalizations of Schöbi's tetrahedral dissection, by Sloane and Vaishampayan, arXiv:0710.3857. ( Read more... )

Simple, linear-time modular decomposition, by Tedder, Corneil, Habib, and Paul, arxiv:0710.3901. ( Read more... )

The lonely vertex problem, by D. Frettlöh and A. Glazyrin, arxiv:0710.4870. ( Read more... )

Triangular peg solitaire unlimited, by G. I. Bell, arxiv:0711.0486. ( Read more... )

Permutations defining convex permutominoes, by Bernini, Disanto, Pinzani, and Rinaldi, arxiv:0711.0582. ( Read more... )

Ed Pegg recently asked me whether various cubic Hamiltonian graphs were xyz graphs, and to test this out I added some code for LCF Notation to GraphExamples.py in PADS along with several specific graphs that can be described using this notation.

( Read more... )

I've put up a new version of PADS after making some bugfixes to the partial cube recognition code, which the two graphs below gave some trouble to. The left one was causing it to crash, and the right was incorrectly reported as being a partial cube; the correct behavior in both cases was to report that the graph is not a partial cube. Now fixed and added to the test cases within the code. Fortunately, both bugs involved minor implementation details rather than any issues with the algorithm itself.

Michael Mitzenmacher has another interesting post up, on programming in theory courses. I don't see a way to comment directly on his blog, so I'll put some of my own thoughts on the subject here:

( Read more... )

I've added an implementation of my new quadratic-time partial cube recognition algorithm to PADS, replacing a theoretically-slower algorithm for the same purpose that was already there.

The labeling phase of the algorithm, as described in my new paper, is recursive, possibly with Ω(n) levels of recursion, and Python doesn't necessarily handle deep recursion well, so I replaced that part of the algorithm with an iterative version using union-find to keep track of the edge equivalences. I also added some new modules: Medium.py handles the media-theoretic side of the algorithm (the part where we find shortest paths between all pairs of vertices), and BFS.py is a new breadth first search routine (I previously had LexBFS.py, which generates a breadth first traversal order of the vertices in a somewhat complicated way, but in this algorithm I needed the sequence of subgraphs connecting each BFS level to the next level). I also checked in a module I'd written some time ago: GrayGraph, a description of the Gray Graph, a cubic bipartite graph that looks a lot like a partial cube but isn't (so making an interesting test case for my new code).

The implementation ended up being rather more complex than I was hoping: including unit tests, comments, etc., Medium and PartialCube together now have 588 lines. The unit tests were very helpful in shaking the bugs out of the Medium side of the implementation, and (though surprisingly to me the PartialCube side worked correctly the first time) in making me confident of the correctness of that part as well.

Someone emailed me about a problem that involved reading graphs from files and running one of the algorithms in PADS on them, that led me to realize that I hadn't included any I/O functionality in PADS. So I dug up a module I'd written a few years ago for reading graphs in various formats, and added it to PADS as ReadUndirectedGraph.py. The formats it knows about are MALF, edge list, node list (possibly broken link? cached on Google), GraphML, graph6, sparse6, and LEDA.GRAPH. Fortunately, they all look different enough from each other that it can distinguish them without having to be told which format to use...

Mark-Jason Dominus writes about taking a too-sophisticated approach to a puzzler about generating triangular numbers with many divisors, ending up with a program that produced the answer in a minute and a half when someone else's much shorter code took longer to compute but was much easier to write. But it took me only five minutes or so (less than it took to write up this blog entry, anyway) to write an even more sophisticated version of Dominus' idea, that produces the correct answer instantly even running in Python. The trick: I had a generator for already-factored numbers prewritten in PADS and could just plug it in. And of course taking advantage of the analysis Mark wrote up for his solution kept my own thinking time low...here's the code, slightly cleaned up from my original version:

( Cut for Python code )

Størmer's theorem is really not so much a theorem but an algorithm, for finding pairs of consecutive smooth numbers. I implemented Lehmer's improved algorithm, in Python, as I needed an implementation to generate some examples while expanding the WP article (and maybe more importantly, to check my understanding of Lehmer's paper). It includes some careful bignum implementation of rounding quadratic irrationals, finding continued fraction expansions of them, and solving Pell's equation, which might be useful for other purposes.

Algorithm Education in Python. Pai Chou, a faculty member in a different department here at UCI, clarifies why Python is such a good match for the pseudocode commonly used in algorithms classes. An old paper from 2002, but still quite valid, I think. (Via.)

I just merged today two articles in Wikipedia on the same problem, regular numbers, that is, numbers of the form 2ⁱ3^j5^k. The ancient Babylonians liked these numbers because their reciprocals had terminating sexagesimal representations; more recently, "Hamming's problem" is the task of generating these numbers efficiently. Dijkstra described a very simple but inefficient functional solution: if H is the sequence of numbers to be generated, then H consists of the number 1 concatenated with the merger of 2H, 3H, and 5H. But there are two flaws with this. First, it generates the same numbers many times: e.g. 6 is generated both as 2*3 and as 3*2. This can be alleviated by letting A be the powers of 2, B be A merged with 3B, and H be B merged with 5H. But even if that improvement is made, and one does the merges in the most obvious and natural lazy-functional way, the sequence is generated inefficiently: the number of merges a generated number is involved in is O(j+k), and all these merges cause generating n numbers in this way to take Θ(n^4/3) arithmetic operations.

Below is an approach that instead generates n numbers in a linear number of operations. The basic principle is, when solving equations like H = merge(B, 5*H), instead of generating 5*H using a recursive instance of the same algorithm, to cache a single list of the whole sequence of H and reuse that stored sequence when computing 5*H.

( Cut for Python code )

However, the pure functional approach uses less space to generate the same sequence, O(n^2/3) instead of O(n) for the caching approach. I can find algorithms that generate the sequence in O(n^2/3) space and O(n log n) time (essentially, by turning it into a graph shortest path problem in an infinite grid graph, applying Dijkstra's algorithm for shortest paths, and only keeping track of nodes on the search frontier), but is it possible to solve optimally both in time and space?

ETA: Yes, it is, by a very simple change to the above code. Simply throw away the parts of the cached sequence that the back pointer has moved past. ( Read more... )

ETA2, March 14: Writing in a more imperative style actually gives tighter code, though perhaps it is not as readable: ( Read more... )

The sieve of Eratosthenes is normally used to generate prime numbers. But it's almost as easy and almost as efficient to make it generate the sequence of factorizations of all natural numbers: instead of eliminating all multiples of each prime p, add p to the list of factors of all its multiples. An implementation of this (with some more care about handling prime powers) is now in Eratosthenes.py in my PADS library.

As an application, I included in the same file a generator for practical numbers, based on a formula described in the link for testing whether a number is practical given its factorization. Despite being highly composite, practical numbers have a lot of similarities as a sequence to prime numbers; one of those similarities is that this formula is a lot like a sieve. But rather than implement it as a sieve, it seemed easier to use the simpler Eratosthenes sieve to generate the factorizations and then test each factorization independently.

Engel expansion! How did I overlook this in my previous investigations of how to generate Egyptian fractions? As I just finished including in that article, this technique was apparently known to Fibonacci, who also investigated the more well known greedy method for Egyptian fractions.

Anyway, my Egyptian fraction Python code now knows how to generate expansions of this type. If you run, e.g.,

egypt.py 7/17 engel

you will get the expansion

1/3 + 1/15 + 1/90 + 1/1530

where each denominator is a multiple of the previous one: 3, 3×5, 3×5×6, 3×5×6×17. By default it generates the standard expansion in which these multipliers 3, 5, 6, 17 are nondecreasing, but it can also generate infinitely many other expansions without this nondecreasing property, so the -a option should only be used in conjunction with -d or -l.

leonardo_m just posted a fine collection of links to algorithmic art. His link to nodebox attracted my attention first as it is based on Python on the Mac, but the art he shows in his post from chriscoyne is more visually appealing to me: of the fractal L-system type, but with very good control of color, line, and tonality. I also found the image below (from nodebox, created by Tom De Smedt) especially intriguing; it's a visualization of the pattern of motion of a flock of simulated birds, but I'm not sure whether it was intended more as a visual aid for understanding that motion or as a piece of abstract art.

( Cut for large image )

Dilworth's Theorem: in any partial ordering, the size of the largest possible antichain is equal to the smallest possible number of chains in a partition of the items into chains.

( Proof sketch )

Anyway, I've now implemented all this in PADS. It's in a new module PartialOrder.py, which also contains everything that used to be in TopologicalOrder.py. Along with the new MaximumAntichain and MinimumChainDecomposition functions, I included a TransitiveClosure function, as it was needed when the input order is given by a sparser DAG.

Golly is a cross-platform cellular automaton simulator that incorporates hashlife algorithms, allowing it to simulate most Game of Life patterns for astronomical numbers of generations quickly. Version 1.0 was just released, and now also incorporates Python scripting, allowing the construction of astronomically complex patterns via modular high level programming rather than bit-at-a-time copy and pasting.

For instance, here is bricklayer.py, a pattern first found by David Bell in 2002 that is formed by combining two p22 oscillators to make a gun, combining three p22 guns to make a p154 gun, and combining two p154 guns with a p7 glider reflector and a messy small Life pattern called "lumps of muck". By itself, lumps of muck forms two subpatterns that, individually, would themselves turn back into LoM again, but unfortunately interfere with each other creating a mess. The streams of gliders from the two guns and reflector repeatedly hit the back subpattern, creating a block and allowing the other subpattern to grow without interference, effectively pushing the LoM forward diagonally at each step. It ends up forming two long diagonal sequences of blocks, with the guns at one end of the sequence, the LoM on the other, and the gliders running down the middle.

( Python )

And here it is in the previously standard form for such patterns, a run-length-encoded text file:

( RLE )

If you want to see the pattern itself within Golly, either one works equally well, but I think you can see which one is more human readable and modifiable as source code.

This is making me think the push-pull replicator-based spaceships that should exist in HighLife and B368/S12578 may be within reach of construction now...

An antimatroid (also called a learning space) is a family of sets closed under union, such that for every nonempty set in the family there is an element that can removed to produce another set in the family. It turns out that, for any antimatroid A except powersets, there is a set S on the same elements that can be added to A to make a larger augmented antimatroid.

( Proof sketch )

This proof can be turned into a deterministic algorithm for choosing a specific augmentation out of all the possible ones. If one considers the augmented antimatroid to be the parent of the original antimatroid, this parent relation forms a tree with the powerset at its root:



We can then generate all the possible antimatroids using reverse search, which is just a fancy phrase for an algorithm that explores a tree like this one by reversing the parent relation. The time per antimatroid generated, over a k-element universe, is polynomial in 2^k, not bad since the number of antimatroids generated is double-exponential in k.

I implemented this in Python, from which I generated the 22 families in the tree of 3-element antimatroids above, matching some hand calculations I'd done previously. My program also found 485 antimatroids over 4 elements and 59386 over 5 elements, but I am a little distrustful of these numbers since the program is intricate and painful to debug. Unfortunately I couldn't find anything about this enumeration problem in the OEIS or with Google, so I'm at a bit of a loss for how to check these results; I suppose I could at least enumerate all 2¹⁶ 4-element set families and test which ones are antimatroids.

If the pattern of number of antimatroids being roughly 2^2k - 1 holds, the program as it stands should be able to list all 6-element antimatroids in about a month of laptop time, but I suspect reimplementation in a faster language and/or running it on a faster machine should speed that up significantly, to closer to a single day.

ETA 6/19: This much shorter program for brute-force testing all 2¹⁶ 4-element set families produces exactly the same output as my reverse search, so I am now much more confident in the reverse search results.

ETA2 6/20: Now up on the OEIS.

An st-orientation of an undirected graph is an assignment of a direction for each edge, turning it into a directed acyclic graph with one source and one sink. For any biconnected graph such an orientation may be found in linear time via depth first search: For each DFS tree edge (u,v) (where u is the parent), orient (u,v) oppositely to the DFS tree edge down from the low vertex of v into the subtree containing u and v. For each nontree edge (u,v) (where u is the top vertex), orient (u,v) the same as the DFS tree edge down from u into the subtree containing v. I've been meaning to add this to PADS for a while, after lecturing on it as part of the graph drawing section of my undergraduate graph algorithms class this quarter, but it's now there (in Biconnectivity.py) after I found motivation in a connection to xyz graph recognition:

Any biconnected 3-regular graph can be partitioned in at most 2^n/2-1 different ways into three matchings. The idea of the proof is to process the vertices in a topological ordering (also added to PADS, in TopologicalOrder.py) of an st-orientation, and when processing each vertex choose how to assign the adjacent edges to the three matchings, consistently with prior choices. We only have to make a binary choice at the n/2-1 vertices with one incoming and two outgoing edges; at all other vertices there is no choice to be made. This leads to an O(n 2^n/2) time algorithm for recognizing xyz graphs, which I implemented as part of PADS in xyzGraph.py.

( Read more... )