Rotate, Sort, Last Column: The Burrows Wheeler Transform in SPARK

Introduction

It’s been a real joy using modern AI models with Ada and SPARK, in particular SPARK. The agents can easily avoid the non-supported language constructs and respect the other language constraints. Proofs of absence of runtime errors are then usually a formality. It is very easy now to write real-world code that’s entirely in the SPARK subset, maybe except some boundary code, and that is proved to be free of all the issues that SPARK checks for.

Full functional correctness (platinum) proofs are somewhat harder and become more of a subproject in itself. But in many cases, they are quite a bit more achievable than I initially thought. In this blog post, I present how I tasked an agent with implementing the Burrows-Wheeler Transform and a more advanced variant, prove it at platinum level, and improve its runtime performance until it is within a factor of 2 or 3 of the reference implementations. If you want to jump right in, you can find the code in this github repository.

Speaking about AI, while AI was used heavily, or even exclusively, for the implementation and SPARK proof of the algorithms, AI was not used for the writing of this article, except for the animations and minor corrections such as typos.

Now let’s move on to some background about the transform itself.

The Burrows-Wheeler Transform

The Burrows-Wheeler Transform (BWT) is a deceptively simple string manipulation technique. It is used in compression and string searching. It starts by taking the input string and generating all rotations of it: removing the last character and adding it to the front produces one rotation; do this n times where n is the length of the string, and you have obtained all rotations and are back to the original string. Write all these rotations in a matrix (they are all of the same length, so the matrix is n times n). Now, sort the rotations in alphabetical order, and finally take the last letter of each string to obtain the result.

Why would I do that?

If we had taken the first letter of each rotation after sorting, we would simply have the letters of the input string, sorted alphabetically. By taking the last column instead, what we get is the letters sorted by the character that follows. Given that natural languages like English have a lot of patterns, we can expect this last column to have a lot of structure. For example, the letter t is often followed by h in English. The h ‘s are all sorted together in the first column, so we can expect to find a lot of t ‘s together in the last column. This structure can then be exploited by other processing techniques such as run-length encoding or move-to-front. By itself, the BWT does not compress anything as the output string is of exactly the same length.

Can I get my input string back?

Yes, and no. The BWT gives us one column of the matrix (the last column). But it’s easy to get another column: just sort input letters and we have the first column back. The sorting gives us a permutation (mapping) between neighboring columns, and we can use this permutation to get the second column from the first, and rebuild the entire matrix in this way.

The question is then, which row was the input string? Without any extra information, we don’t know. If we look at the construction of the BWT above, then any rotation of the input string would yield the same output. This is because we take all rotations of the input string, then sort them, obtaining the exact same matrix. So at that point we don’t know anymore which original rotation we started from. In other words, the BWT maps all rotations of a string to the output string, rather than mapping one string to another.

There are two simple ways to “fix” the procedure to identify the input string. One is to append a unique character that wasn’t in the input string, at the end. Then we know which row was the original: the one with the unique character at the end. Another way is to make the output of the BWT a tuple (string, i), where i is the index of the row of the sorted matrix that corresponds to the input string. Implementations of this classical BWT typically choose the second solution, as adding a character to the alphabet is often more awkward than having more than one output.

Can’t we have something nicer?

The situation is a bit unsatisfactory. To really get our input string back, we either need to add something to the input string, or add something to the output. Also, given that many input strings of length n map to the same output of length n, it means many strings of length n can never be the result of the BWT. Wouldn’t it be nice if every input string mapped to a different output string of length n, while still keeping the nice properties of the BWT which were the reason we use it? Enter the bijective BWT.

Lyndon words

The crux of the problem is identifying the “correct” string among all its rotations. What mathematicians do in this case is simply define what they want and go from there. So they defined a Lyndon word to be a word that’s strictly smaller than all its rotations. If the input to the BWT was a Lyndon word, then, after reconstructing the matrix from the last and first column, we would know the input must be the first row.

Lyndon factorization

Now remains the problem that the input is most likely not a Lyndon word. However, a theorem in word theory tells us that each string can be uniquely cut into a sequence of Lyndon words that get smaller (or equal) at each boundary. For example, the word GNATPROVE can be decomposed into the two Lyndon words GN and ATPROVE (you can check for yourself that both words are indeed strictly smaller than all their rotations, and also that GN is greater than ATPROVE). A technique called Duval’s algorithm is used to obtain the Lyndon factorization of the input string.

And now?

The bijective BWT (BBWT) works as follows. Cut the input string into non-increasing Lyndon words, then produce all rotations of each Lyndon word. To put all these strings into the same matrix, we need to “complete” the shorter ones by simply repeating them until they fill their row. Then we sort the matrix as before, and take the last letter of each word (not the last column this time).

The parts I left out

There is a fair amount of hand-waving in this text, in particular I didn’t explain properly how to obtain the input string back, and why exactly it works. If you want the full mathematical version, I invite you to read this excellent paper, which also attributes each idea to the researchers that came up with it.

Can we put all this into SPARK?

It turns out that we can. Here is the full implementation of classical and bijective BWT in SPARK, with full correctness proofs and roundtrip proofs (that decoding recovers the input string).

One noteworthy part is the signature of the main functions, which corresponds precisely to my explanation before (contracts omitted):

  type Classical_Result (Length : Natural) is record
      Last    : String (1 .. Length);
      Primary : Natural;
   end record;

  function Classical_Encode (S : String) return Classical_Result;

  function Bijective_Encode (S : String) return String;

As stated previously, the traditional (classical) BWT needs to produce the index of the row in the matrix alongside the string output, or it can’t be fully recovered. The bijective encode simply returns the string.

The big highlight is the roundtrip theorem (the fact that encoding a string, then decoding it, yields the input string) for both functions (this is the actual code from the files, no omissions this time):

   procedure Classical_Round_Trip (S : String)
   with
     Global => null,
     Pre    => Supported (S),
     Post   =>
       Classical_Decode
         (Classical_Encode (S).Last, Classical_Encode (S).Primary)
       = S;

   procedure Bijective_Round_Trip (S : String)
   with
     Global => null,
     Pre    => Supported (S),
     Post   => Bijective_Decode (Bijective_Encode (S)) = S;

   Max_Length : constant := 2**24;

   function Supported (S : String) return Boolean
   is (S'First = 1 and then S'Length <= Max_Length);

The precondition Supported is just there to fix the ’First of strings to 1, and cap their length, to simplify handling of string bounds and overflow checks.

Note that the two Encode functions have elaborate postconditions (omitted here) which make sure the implementation of Encode can’t just be the identity function, and that it actually corresponds to the BWT.

How big of an effort is this?

Pretty big. The theorems above hide a huge amount of helper lemmas. Also, the bijective variant is a much bigger proof effort.

Excluding sorting, the SPARK implementation of the classical BWT encoder and decoder fits in about 100 lines of executable SPARK (without contracts), but if you add contracts, ghost code and lemmas, this grows to ~700 lines.

The sorting routine, which is shared between the two encoders, is highly specialized for this algorithm, requiring 900 lines of code, and 2000 lines of ghost code and lemmas.

The bijective variant of the BWT is only ~150 lines of SPARK code, but close to 4000 lines of lemmas and ghost code. This part is so big because it also includes Duval’s algorithm for the Lyndon word factorization. However, the agent went a bit overboard here, as it also specified things that aren’t strictly needed, such as that the bijective BWT is actually bijective (any string can appear as the output of the encoder). This took a lot of extra space.

On my machine with 16 cores, the full proof takes a little under 17 minutes to prove all 10,406 checks. Maybe surprisingly, all checks could be proved with a 5s timeout, though some are quite tight so the Makefile still sets a 20s timeout.

How much handholding did the AI need for this?

Very little besides telling the agent to keep going. The code was 100% generated by AI, including type definitions, contracts, lemmas, loop invariants, everything. I only intervened to give the general direction of the effort, such as “now attempt the proof of the bijective variant“, or tell the agent to continue when it came back with a new lemma, but not the full proof.

The performance of the first fully proved implementation was not good. So then I told the agent to improve the performance while maintaining the proof status. It chose the optimizations and algorithmic improvements itself, and managed to keep the proofs green. As a result, this BWT implementation is roughly within a factor of 2 or 3 of state-of the art implementations in tools such as gzip and libsais, which is very good, given how optimized these tools are, and that some optimizations are still not implemented in this experiment. Some remaining optimization ideas are: algorithmic improvements, improving cache behavior, and tuning the generated assembly for specific processors.

Previous Formalizations

Just last year, in 2025, researchers implemented the classical BWT in the Isabelle proof assistant with full proofs.

I couldn’t find any proofs of the bijective BWT, nor Duval’s algorithm, that were of the implementation rather than only the algorithm (paper proof, or formalization of the paper proof). It seems to be the first time this is done, though I would be happy to be corrected on this one.