Prefix function and KMP

Every molecular biology laboratory has a tube labelled “λ/HindIII”. It holds the DNA of a bacterial virus, phage lambda, cut by the enzyme HindIII. On a gel this mixture separates into eight bands of known length: 23,130, 9416, 6557, 4361, 2322, 2027, 564 and 125 base pairs, and other molecules are measured against them like a ruler. Let us predict this ladder with a program. The enzyme cuts DNA wherever the motif AAGCTT occurs, so it is enough to find all occurrences of six letters among 48,502. Below, the Knuth–Morris–Pratt algorithm finds them, one letter at a time, from left to right. Choose another genome, another enzyme, or type your own motif.

The program found six sites, giving seven pieces, while the ladder has eight bands: instead of 6557 and 125 we have one piece of 6682. There is no error in the search. The NCBI database holds the genome of the wild phage, while the tubes sold contain the DNA of a laboratory strain, the mutant cI857ind1 Sam7. Its ind1 mutation is a change from C to T at position 37,589, and the sequence AAGCTC at that place becomes AAGCTT. Tick “λ from a tube”: the 6682 piece splits into 6557 and 125, and the ladder matches the marker down to the last band. One letter of 48,502 shows on the gel.

Now look at the counters. On lambda with HindIII the search made 60,083 letter comparisons, 1.24 per letter, and the naive try-every-position search 65,502, 1.35 per letter. The difference is small, and that is an honest result: in real DNA the very first letter of the pattern usually does not match, and the naive search has almost nothing to lose. The KMP algorithm is valued for something else. It never goes back in the text, so it suits a stream that cannot be rewound, and even on the most unfortunate text it makes fewer than 2n comparisons. The second scene will show what an unfortunate text looks like.

The idea

Suppose the pattern matched the text on its first k letters and differed on the next. We already know these k letters of the text without rereading them: they are the beginning of the pattern. So where to shift the pattern can be decided without looking at the text at all. Find the longest beginning of the pattern with which its first k letters end, and continue comparing from its end. These lengths are computed in advance, once for each k, and the pointer in the text never goes back again.

A beginning of a string that is also its end, but not the whole string, is called a border. The longest border of the word “abracadabra” is “abra”: the word begins and ends with it. “onion” has “on”, and GAATTC has no border. The prefix function of a string p is the array π, where π[i] is the length of the longest border of the beginning p[0..i]:

 a  b  r  a  c  a  d  a  b  r  a
 0  0  0  1  0  1  0  1  2  3  4

Under the last a stands 4: the whole word has the border abra. Under the second b, 2: “abracadab” has the border “ab”. The table depends only on the pattern, and the search, as we will see now, comes down to one number k and this table. The prefix function is built from scratch and step by step in the course chapter “A needle in a haystack”; here are the proofs, the behaviour on worst texts, code for a stream, and problems with no search in them at all.

Search without stepping back

Go through the text from left to right, keeping one number k: how many first letters of the pattern match the end of what has been read. Read a letter c:

1. if c equals p[k], the match is extended: k grows by one;
2. if not and k > 0, the pattern shifts: k becomes π[k − 1], and the same letter c is compared again, now with another letter of the pattern;
3. if not and k = 0, no letters have matched: move on to the next letter of the text.

When k reaches the length of the pattern m, an occurrence is found, and k becomes π[m − 1]: the next occurrence may start inside this one. Here are both searches on the same strings, one comparison per step. On the left of the chart is the number of the text letter the search is looking at, horizontally how many comparisons have been made so far.

On a telomere, as the ends of chromosomes are called, where in humans the unit TTAGGG repeats, a pattern of five units almost matches the text at every sixth position, and only the last letter C rejects it. The naive search each time goes through 29 matching letters, stumbles on the thirtieth, shifts by one position and rereads them: on the chart this is a saw. KMP after a mismatch shifts the pattern by exactly one unit: the longest border of the 29 matched letters is 23 letters, one unit shorter. And the comparison continues from the same letter of the text. Its line is a staircase that never goes down. On 300 letters of telomere the naive search makes 1650 comparisons, KMP 346.

Switch to “AAAA…”, a text of nothing but A and a pattern of nineteen A and one C. This is KMP’s hardest case: at each letter it compares twice, with C and, after a fallback, with A. That is 581 comparisons against 5620 for the naive search. And on a piece of real lambda DNA the two lines almost merge: on ordinary text the naive search is good, and the difference shows only in the worst cases.

How to compute the prefix function

Trying all possible border lengths for each beginning is on the order of m³ operations. But the prefix function is built from itself, with the same movement as the search. Suppose k = π[i − 1] is already known. If p[i] = p[k], the border extends: π[i] = k + 1. If not, a shorter border that does extend is needed. Which other borders a string has besides the longest is told by the lemma of the next section: the border of the border, then the border of the border of the border, and so on. Their lengths are already written in the table: π[k − 1], π[π[k − 1] − 1], …

Click any letter: under the string appears the chain of borders of the beginning that ends with this letter, and orange marks show its smallest period.

Look at “abaababaabaababaaba”. This is a Fibonacci string without its last two letters: f₁ = b, f₂ = a, and each next string is the previous one glued to the one before it, f_k = f_k−1f_k−2. It has an unusually long chain of borders for such a short string: 11, 6, 3, 1. And “aaaaaaaa” has the longest possible chain: each next border is exactly one letter shorter than the previous one.

There is also a shortcut, with no separate search function. Glue together the pattern, a separator that occurs in neither string, and the text: p + “#” + t. The prefix function of the glued string equals m exactly where an occurrence of the pattern ends in the text: a border cannot be longer than m, because it cannot step over “#”. The KMP search is this same computation, only without storing the glued string.

Why it is right

Correctness rests on one observation: when a comparison has failed, the only thing worth remembering about what has been read is its longest end that matches a beginning of the pattern. Everything else follows from the table. So the letters of the text need not be stored at all.

Code

The search is written as an object that is fed the text in pieces: between pieces it remembers only k and how many letters have been read. This is convenient for reading a big file line by line, a network packet by packet, a tape block by block. For comparison there is a naive search beside it, and both count letter comparisons, the same ones the scene above shows.

Look at the last lines of the output. On a telomere of 5000 units, the naive search checks the pattern “30 units and TTAGGC” almost a million times, KMP 34,970. And when the pattern is exactly 30 units, KMP finds 4971 occurrences and makes exactly one comparison per letter: after each occurrence found, the fallback by one unit leaves 174 letters matched, and the next letter extends them at once.

Now on a real file. The lambda genome is in FASTA format, 70 letters per line. The program does not glue the lines into one but feeds them in turn to nine searches at once, one per enzyme:

31 sites, split between two lines of the file: NdeI at 38357-38362

The NdeI site CATATG begins in one line of the file and ends in the next. No line contains it whole, and yet the search found it: the state k survived the line break. A search that falls back in the text on a mismatch would need a buffer. This is exactly the difficulty James Morris was solving when the algorithm was invented.

What it costs

This is an amortised bound: there can be many fallbacks on one letter of the text. The pattern of m − 1 letters A and a letter C on the text A…AB falls back m − 1 times in a row at the letter B, but these fallbacks are paid for by the m − 1 letters before them. On average there are at most two comparisons per letter, and an expensive letter cannot recur until k builds up again.

The benchmark: the text is a thousand telomere units, the pattern a few units with TTAGGC at the end, so that almost every sixth placement gets to the last letter. The pattern length grows, the text length stays the same.

The naive search gets dearer with the pattern: from 96 letters to 384 the pattern got four times longer, and the time grew about fourfold: this is the term n · m/6 from the placements at the start of each unit. On short patterns it is still overshadowed by the ordinary pass over the text, so the overall slope is less than one. The KMP line is horizontal: the pattern length hardly affects the search time. Python’s built-in search is written in C and is the fastest in this benchmark; for long patterns it has used, since version 3.10, the two-way algorithm of Crochemore and Perrin, which, like KMP, is linear in the worst case. Go’s built-in search for short patterns tries positions with machine instructions that compare many bytes at once, and for long ones, if the trying often fails, switches to the Rabin–Karp algorithm. In ordinary work write find, in and strings.Index; your own KMP is needed where the built-in search does not fit: in a stream, on sequences that are not characters, with fallbacks to count, as in the problems below.

The prefix function without search: the period

A number q is called a period of a string s if each letter repeats the one q places to the left: s[i] = s[i + q] for all i. “abracadabra” has period 7: “abracad|abra”. The length of the string need not be divisible by the period.

By the border lemma all the periods come out too: n minus the lengths of all the borders of the chain. If n is divisible by the smallest period q, the string is n/q copies of its beginning of length q. If not, it cannot be made of equal pieces shorter than the whole string; this follows from the Fine–Wilf lemma on two periods, which is also proved in the paper by Knuth, Morris and Pratt.

So the prefix function answers questions with no search in them at all: which unit a repeat in DNA is made of, whether a string compresses into “unit × number”, what is the shortest piece to repeat to get a given string. Try “TTAGGGTTAGGGTTA” in the scene: the period is 6, but 15 is not divisible by 6, the read broke off in the middle of a unit. This is the problem about the telomere below.

Where it works

KMP itself is rarely met in standard libraries. In glibc, the standard C library on Linux, long patterns are searched for by the two-way algorithm of Crochemore and Perrin: like KMP it is linear in the worst case, but instead of a table of m numbers it manages with a few variables. GNU grep searches for fixed strings with variants of the Boyer–Moore algorithm, which on ordinary text jumps over letters without reading them. But the idea of KMP, falling back by a table built in advance without touching the text, works wherever there are many patterns or the text cannot be rewound.

• Many patterns at once. The Aho–Corasick automaton (1975) is the prefix function carried over from one string to a trie of thousands of words: each node has a fallback link to the longest suffix that is also in the trie. The first version of fgrep in Unix was built on it, and so are the filters that look for thousands of signatures in traffic and files.
• An automaton. If the alphabet is small, as for DNA, one can compute in advance, for each k and each letter, where the search will end up after all the fallbacks. The result is a finite automaton that spends exactly one transition per letter, with no loop. Regular expression engines that build automata are made this way; more in the chapter “Automata and regular expressions”. You can assemble the KMP automaton for your own word in the course chapter.
• A stream. KMP reads each letter once and remembers one number, so it suits data that passes by and is not stored: network packets, sensor readings, a tape.
• The structure of a string. Periods, repeats, compression, the number of distinct substrings, comparing a string with its rotations: the prefix function and its relative, the Z-function, solve these in linear time, and in programming olympiads they are workhorses.

An editor, an automaton and a palindrome problem

In the summer of 1969 James Morris was writing a text editor for the CDC 6400 machine. Buffering the file was complicated, and the aim was to find a pattern without ever going back in the text. Morris devised a way based on the theory of finite automata, but in that form it was not evident that the time was linear. In the 1977 paper the co-authors tell how the procedure turned out too complicated for the system’s other implementers, and a few months later Morris found that uninvited “fixes” had reduced it to ruins.

Independently, at the beginning of 1970, Donald Knuth was working through Stephen Cook’s theorem on two-way pushdown automata. It implied that a problem Knuth could not solve faster than in about n² steps, recognising a string that begins with an even-length palindrome, can be solved in linear time. Knuth traced Cook’s construction by hand, extracted the mechanism from it and generalised it to a program that finds the longest beginning of one string occurring in another. In Knuth’s words, this was the first time the theory of automata had taught him to solve a practical problem better than he could before. Vaughan Pratt made the running time independent of the alphabet size, and when Pratt described the algorithm to Morris, Morris recognised it as the same one. Morris and Pratt described it in a Berkeley technical report in 1970, and the joint paper of the three authors, “Fast pattern matching in strings”, came out in SIAM Journal on Computing in 1977.

Earlier still, in 1969, a linear search for texts of zeros and ones was built by Yuri Matiyasevich, in the form of a Turing machine with two-dimensional memory; the course chapter tells about this.

Problems

In the first problem the genome is closed into a ring, and occurrences can overlap and pass through the place where the ring was cut when written to a file. In the second nothing has to be searched for: from a piece of a repeat that begins and breaks off in the middle of a unit you have to find the unit itself, and that is a problem about the period.