Connecting Elliptic Curves with Finite Fields

So here we are. We’ve studied the general properties of elliptic curves, written a program for elliptic curve arithmetic over the rational numbers, and taken a long detour to get some familiarity with finite fields (the mathematical background and a program that implements arbitrary finite field arithmetic).

And now we want to get back on track and hook our elliptic curve program up with our finite field program to make everything work. And indeed, for most cases it’s just that simple! For example, take the point P = (2,1) on the elliptic curve y = x^3 + x + 1 with coefficients in \mathbb{Z}/5. Using purely code produced in previous posts, we can do arithmetic:

>>> F5 = FiniteField(5, 1)
>>> C = EllipticCurve(a=F5(1), b=F5(1))
>>> P = Point(C, F5(2), F5(1))
>>> P
(2 (mod 5), 1 (mod 5))
>>> 2*P
(2 (mod 5), 4 (mod 5))
>>> 3*P
Ideal

Here’s an example of the same curve y^2 = x^3 + x + 1 with coefficients over the finite field of order 25 \mathbb{F}_{5^2}.

>>> F25 = FiniteField(5,2)
>>> F25.idealGenerator
3 + 0 t^1 + 1 t^2
>>> curve = EllipticCurve(a=F25([1]), b=F25([1]))
>>> x = F25([2,1])
>>> y = F25([0,2])
>>> y*y - x*x*x - x - 1
0 ∈ F_{5^2}
>>> curve.testPoint(x,y)
True
>>> P = Point(curve, x, y)
>>> -P
(2 + 1 t^1, 0 + 3 t^1)
>>> P+P
(3 + 1 t^1, 2)
>>> 4*P
(3 + 2 t^1, 4 + 4 t^1)
>>> 9*P
Ideal

There are some subtle issues, though, in that we shouldn’t use the code we have to work over any finite field. But we’ve come very far and covered a lot of technical details, so let’s briefly remember how we got here.

Taking a Step Back

At the beginning there was only \mathbb{Q}, the field of rational numbers. We had a really nice geometric picture of elliptic curves over this field, and using that picture we developed an algorithm for (geometrically) adding points.

add-points-exampleIf we assume the equation of the elliptic curve had this nice form (the so-called Weierstrass normal form, y^2 = x^3 + ax + b), then we were able to translate the geometric algorithm into an algebraic one. This made it possible to write a program to perform the additions, and this was our first programmatic milestone. Along the way, we learned about groups and projective geometry, which I explained was the proper mathematical setting for elliptic curves. In that setting, we saw that for most fields, every elliptic curve could be modified into one in Weierstrass normal form without changing the algebraic structure of the set of solutions. Moreover, we saw that you can replace the field \mathbb{Q} with the field of your choice. The set of solutions to an elliptic curve still forms a group and the same algebraic point-adding algorithm works. It’s just an interesting quirk of mathematics that one way to represent elements of finite fields are as polynomial remainders when dividing by a “prime” polynomial (analogous to modular arithmetic with integers). So we spent a while actually implementing finite fields in terms of this representation.

The reader has probably heard of this, but in practice one uses a (very large) finite field for the coefficients of their elliptic curve. Often this is \mathbb{Z}/p for some really large prime p, or the field of 2^m elements for some large integer m. But one would naturally complain: there are so many (infinitely many!) finite fields to choose from! Which one should we use, and how did they choose these?

As with most engineering problems the answer is a trade-off, in this case between efficiency and security. Arithmetic is faster in fields of characteristic 2 (and easy to implement at the hardware level!) but a lot is known about the finite field of 2^m elements. In fact, if you are sloppy in picking m you’ll get no security at all! One prominent example is the so-called Weil descent attack, which breaks security assumptions for elliptic curve cryptography when m is not prime. These attacks use some sophisticated machinery, but this is how it goes. An abstract mathematical breakthrough can immediately invalidate cryptography based on certain elliptic curves.

But before we get neck-deep in cryptography we have an even bigger problem: for some finite fields, not every elliptic curve has a Weierstrass normal form! So our program isn’t expressive enough to represent all elliptic curves we might want to. We could avoid these curves in our applications, but that would be unnecessarily limiting. With a bit more careful work, we can devise a more general algorithm (and a different normal form) that works for all fields. But let’s understand the problem first.

In general, you can have an elliptic curve of the form \sum_{i+j=3} a_{i,j}x^iy^j = 0. That is, it’s just a really general degree 3 polynomial in two variables. If we assume the discriminant of this polynomial is nonzero, we’ll get a smooth curve. And then to get to the Weierstrass normal form involves a bunch of changes of variables. The problem is that the algebraic manipulations you do require you to multiply and divide by 2 and 3. In a field of either characteristic, these operations are either destructive (multiplying by zero) or totally illegal (dividing by zero), and they ruin Weierstrass’s day.

So what can we do?

Well it turns out that there is a more general Weierstrass normal form, unsurprisingly called the generalized Weierstrass normal form. It looks like this

\displaystyle y^2 + a_1 xy + a_3y = x^3 + a_2x^2 + a_4x + a_6

The same geometric idea of drawing lines works for this curve as well. It’s just that now the formula is way more complicated. It involves computing a bunch of helper constants and computing far more arithmetic. My colleague Daniel Ngheim was kind enough to code up the algorithm, and here it is

    def __add__(self, Q):
        if isinstance(Q, Ideal):
            return Point(self.curve, self.x, self.y)

        a1,a2,a3,a4,a6 = (self.curve.a1, self.curve.a2, self.curve.a3, self.curve.a4, self.curve.a6)

        if self.x == Q.x:
            x = self.x
            if self.y + Q.y + a1*x + a3 == 0:
                return Ideal(self.curve)
            else:
                c = ((3*x*x + 2*a2*x + a4 - a1*self.y) / (2*self.y + a1*x + a3))
                d = (-(x*x*x) + a4*x + 2*a6 - a3*self.y) / (2*self.y + a1*x + a3)
                Sum_x = c*c + a1*c - a2 - 2*self.x
                Sum_y = -(c + a1) * Sum_x - d - a3
                return Point(self.curve, Sum_x, Sum_y)
        else:
            c =  (Q.y - self.y) / (Q.x - self.x)
            d =  (self.y*Q.x - Q.y*self.x) / (Q.x - self.x)
            Sum_x = c*c + a1*c - a2 - self.x - Q.x
            Sum_y = -(c + a1)*Sum_x - d - a3
            return Point(self.curve, Sum_x, Sum_y)

   def __neg__(self):
      return Point(self.curve, self.x, -self.y - self.curve.a1*self.x - self.curve.a3)

I trust that the devoted reader could derive this algorithm by hand, but for a more detailed derivation see the book of Silverman (it’s a graduate level text, but the point is that if you’re not really serious about implementing elliptic curve cryptography then you shouldn’t worry about this more general algorithm).

One might start to wonder: are there still other forms of elliptic curves that we could use to get around some of the difficulties of the Weierstrass normal form? The answer is yes, but we’ll defer their discussion to a future post. The brief explanation is that through a different choice of variable changes you can get to a different form of curve, and the algorithms you get from writing out the algebraic equations for adding points are slightly more efficient.

For the remainder of this series we’ll just work with one family of finite fields, those fields of the form \mathbb{Z}/p for some large p. There is one particularly famous elliptic curve over this field that is used in some of the most secure applications in existence, and this will roughly be our target. In either case, we have provided the combined elliptic curve and finite field code (and the generalized elliptic curve class) on this blog’s Github page.

So in the next post we’ll actually start talking about cryptography and how to use elliptic curves to do things like generate a shared secret key.

Until then!

About these ads

Want to make a great puzzle game? Get inspired by theoretical computer science.

Two years ago, Erik Demaine and three other researchers published a fun paper to the arXiv proving that most incarnations of classic nintendo games are NP-hard. This includes almost every Super Mario Brothers, Donkey Kong, and Pokemon title. Back then I wrote a blog post summarizing the technical aspects of their work, and even gave a talk on it to a room full of curious undergraduate math majors.

But while bad tech-writers tend to interpret NP-hard as “really really hard,” the truth is more complicated. It’s really a statement about computational complexity, which has a precise mathematical formulation. Sparing the reader any technical details, here’s what NP-hard implies for practical purposes:

You should abandon hope of designing an algorithm that can solve any instance of your NP-hard problem, but many NP-hard problems have efficient practical “good-enough” solutions.

The very definition of NP-hard means that NP-hard problems need only be hard in the worst case. For illustration, the fact that Pokemon is NP-hard boils down to whether you can navigate a vastly complicated maze of trainers, some of whom are guaranteed to defeat you. It has little to do with the difficulty of the game Pokemon itself, and everything to do with whether you can stretch some subset of the game’s rules to create a really bad worst-case scenario.

So NP-hardness has very little to do with human playability, and it turns out that in practice there are plenty of good algorithms for winning at Super Mario Brothers. They work really well at beating levels designed for humans to play, but we are highly confident that they would fail to win in the worst-case levels we can cook up. Why don’t we know it for a fact? Well that’s the P \ne NP conjecture.

Since Demaine’s paper (and for a while before it) a lot of popular games have been inspected under the computational complexity lens. Recently, Candy Crush Saga was proven to be NP-hard, but the list doesn’t stop with bad mobile apps. This paper of Viglietta shows that Pac-man, Tron, Doom, Starcraft, and many other famous games all contain NP-hard rule-sets. Games like Tetris are even known to have strong hardness-of-approximation bounds. Many board games have also been studied under this lens, when you generalize them to an n \times n sized board. Chess and checkers are both what’s called EXP-complete. A simplified version of Go fits into a category called PSPACE-complete, but with the general ruleset it’s believed to be EXP-complete [1]. Here’s a list of some more classic games and their complexity status.

So we have this weird contrast: lots of NP-hard (and worse!) games have efficient algorithms that play them very well (checkers is “solved,” for example), but in the worst case we believe there is no efficient algorithm that will play these games perfectly. We could ask, “We can still write algorithms to play these games well, so what’s the point of studying their computational complexity?”

I agree with the implication behind the question: it really is just pointless fun. The mathematics involved is the very kind of nuanced manipulations that hackers enjoy: using the rules of a game to craft bizarre gadgets which, if the player is to surpass them, they must implicitly solve some mathematical problem which is already known to be hard.

But we could also turn the question right back around. Since all of these great games have really hard computational hardness properties, could we use theoretical computer science, and to a broader extent mathematics, to design great games? I claim the answer is yes.

[1] EXP is the class of problems solvable in exponential time (where the exponent is the size of the problem instance, say n for a game played on an n \times n board), so we’re saying that a perfect Chess or Checkers solver could be used to solve any problem that can be solved in exponential time. PSPACE is strictly smaller (we think; this is open): it’s the class of all problems solvable if you are allowed as much time as you want, but only a polynomial amount of space to write down your computations. 

A Case Study: Greedy Spiders

Greedy spiders is a game designed by the game design company Blyts. In it, you’re tasked with protecting a set of helplessly trapped flies from the jaws of a hungry spider.

A screenshot from Greedy Spiders.

A screenshot from Greedy Spiders. Click to enlarge.

In the game the spider always moves in discrete amounts (between the intersections of the strands of spiderweb) toward the closest fly. The main tool you have at your disposal is the ability to destroy a strand of the web, thus prohibiting the spider from using it. The game proceeds in rounds: you cut one strand, the spider picks a move, you cut another, the spider moves, and so on until the flies are no longer reachable or the spider devours a victim.

Aside from being totally fun, this game is obviously mathematical. For the reader who is familiar with graph theory, there’s a nice formalization of this problem.

The Greedy Spiders Problem: You are given a graph G_0 = (V, E_0) and two sets S_0, F \subset V denoting the locations of the spiders and flies, respectively. There is a fixed algorithm A that the spiders use to move. An instance of the game proceeds in rounds, and at the beginning of each round we call the current graph G_i = (V, E_i) and the current location of the spiders S_i. Each round has two steps:

  1. You pick an edge e \in E_i to delete, forming the new graph G_{i+1} = (V, E_i).
  2. The spiders jointly compute their next move according to A, and each spider moves to an adjacent vertex. Thus S_i becomes S_{i+1}.

Your task is to decide whether there is a sequence of edge deletions which keeps S_t and F disjoint for all t \geq 0. In other words, we want to find a sequence of edge deletions that disconnects the part of the graph containing the spiders from the part of the graph containing the flies.

This is a slightly generalized version of Greedy Spiders proper, but there are some interesting things to note. Perhaps the most obvious question is about the algorithm A. Depending on your tastes you could make it adversarial, devising the smartest possible move at every step of the way. This is just as hard as asking if there is any algorithm that the spiders can use to win. To make it easier, A could be an algorithm represented by a small circuit to which the player has access, or, as it truly is in the Greedy Spiders game, it could be the greedy algorithm (and the player can exploit this).

Though I haven’t heard of the Greedy Spiders problem in the literature by any other name, it seems quite likely that it would arise naturally. One can imagine the spiders as enemies traversing a network (a city, or a virus in a computer network), and your job is to hinder their movement toward high-value targets. Perhaps people in the defense industry could use a reasonable approximation algorithm for this problem. I have little doubt that this game is NP-hard [2], but the purpose of this article is not to prove new complexity results. The point is that this natural theoretical problem is a really fun game to play! And the game designer’s job is to do what game designers love to do: add features and design levels that are fun to play.

Indeed the Greedy Spiders folks did just that: their game features some 70-odd levels, many with multiple spiders and additional tools for the player. Some examples of new tools are: the ability to delete a vertex of the graph and the ability to place a ‘decoy-fly’ which is (to the greedy-algorithm-following spiders) indistinguishable from a real fly. They player is usually given only one or two uses of these tools per level, but one can imagine that the puzzles become a lot richer.

[2]: In the adversarial case it smells like it’s PSPACE-complete, being very close to known PSPACE-hard problems like Cops and Robbers and Generalized Geography

Examples

I can point to a number of interesting problems that I can imagine turning into successful games, and I will in a moment, but before I want to make it clear that I don’t propose game developers study theoretical computer science just to turn our problems into games verbatim. No, I imagine that the wealth of problems in computer science can serve as inspiration, as a spring board into a world of interesting gameplay mechanics and puzzles. The bonus for game designers is that adding features usually makes problems harder and more interesting, and you don’t need to know anything about proofs or the details of the reductions to understand the problems themselves (you just need familiarity with the basic objects of consideration, sets, graphs, etc).

For a tangential motivation, I imagine that students would be much more willing to do math problems if they were based on ideas coming from really fun games. Indeed, people have even turned the stunningly boring chore of drawing an accurate graph of a function into a game that kids seem to enjoy. I could further imagine a game that teaches programming by first having a student play a game (based on a hard computational problem) and then write simple programs that seek to do well. Continuing with the spiders example they could play as the defender, and then switch to the role of the spider by writing the algorithm the spiders follow.

But enough rambling! Here is a short list of theoretical computer science problems for which I see game potential. None of them have, to my knowledge, been turned into games, but the common features among them all are the huge potential for creative extensions and interesting level design.

Graph Coloring

Graph coloring is one of the oldest NP-complete problems known. Given a graph G and a set of colors \{ 1, 2, \dots, k \}, one seeks to choose colors for the vertices of G so that no edge connects two vertices of the same color.

coloring

Now coloring a given graph would be a lame game, so let’s spice it up. Instead of one player trying to color a graph, have two players. They’re given a k-colorable graph (say, k is 3), and they take turns coloring the vertices. The first player’s goal is to arrive at a correct coloring, while the second player tries to force the first player to violate the coloring condition (that no adjacent vertices are the same color). No player is allowed to break the coloring if they have an option. Now change the colors to jewels or vegetables or something, and you have yourself an award-winning game! (Or maybe: Epic Cartographer Battles of History)

An additional modification: give the two players a graph that can’t be colored with k colors, and the first player to color a monochromatic edge is the loser. Add additional move types (contracting edges or deleting vertices, etc) to taste.

Art Gallery Problem

Given a layout of a museum, the art gallery problem is the problem of choosing the minimal number of cameras so as to cover the whole museum.

artgallery

This is a classic problem in computational geometry, and is well-known to be NP-hard. In some variants (like the one pictured above) the cameras are restricted to being placed at corners. Again, this is the kind of game that would be fun with multiple players. Rather than have perfect 360-degree cameras, you could have an angular slice of vision per camera. Then one player chooses where to place the cameras (getting exponentially more points for using fewer cameras), and the opponent must traverse from one part of the museum to the other avoiding the cameras. Make the thief a chubby pig stealing eggs from birds and you have yourself a franchise.

For more spice, allow the thief some special tactics like breaking through walls and the ability to disable a single camera.

This idea has of course been the basis of many single-player stealth games (where the guards/cameras are fixed by the level designer), but I haven’t seen it done as a multiplayer game. This also brings to mind variants like the recent Nothing to Hide, which counterintuitively pits you as both the camera placer and the hero: you have to place cameras in such a way that you’re always in vision as you move about to solve puzzles. Needless to say, this fruit still has plenty of juice for the squeezing.

Pancake Sorting

Pancake sorting is the problem of sorting a list of integers into ascending order by using only the operation of a “pancake flip.”

panackesortJust like it sounds, a pancake flip involves choosing an index in the list and flipping the prefix of the list (or suffix, depending on your orientation) like a spatula flips a stack of pancakes. Now I think sorting integers is boring (and it’s not NP-hard!), but when you forget about numbers and that one special configuration (ascending sorted order), things get more interesting. Instead, have the pancakes be letters and have the goal be to use pancake flips to arrive at a real English word. That is, you don’t know the goal word ahead of time, so it’s the anagram problem plus finding an efficient pancake flip to get there. Have a player’s score be based on the number of flips before a word is found, and make it timed to add extra pressure, and you have yourself a classic!

The level design then becomes finding good word scrambles with multiple reasonable paths one could follow to get valid words. My mother would probably play this game!

Bin Packing

Young Mikio is making sushi for his family! He’s got a table full of ingredients of various sizes, but there is a limit to how much he can fit into each roll. His family members have different tastes, and so his goal is to make everyone as happy as possible with his culinary skills and the options available to him.

Another name for this problem is bin packing. There are a collection of indivisible objects of various sizes and values, and a set of bins to pack them in. Your goal is to find the packing that doesn’t exceed the maximum in any bin and maximizes the total value of the packed goods.

binpacking

I thought of sushi because I recently played a ridiculously cute game about sushi (thanks to my awesome friend Yen over at Baking And Math), but I can imagine other themes that suggest natural modifications of the problem. The objects being packed could be two-dimensional, there could be bonuses for satisfying certain family members (or penalties for not doing so!), or there could be a super knife that is able to divide one object in half.

I could continue this list for quite a while, but perhaps I should keep my best ideas to myself in case any game companies want to hire me as a consultant. :)

Do you know of games that are based on any of these ideas? Do you have ideas for features or variations of the game ideas listed above? Do you have other ideas for how to turn computational problems into games? I’d love to hear about it in the comments.

Until next time!

Programming with Finite Fields

Back when I was first exposed to programming language design, I decided it would be really cool if there were a language that let you define your own number types and then do all your programming within those number types. And since I get excited about math, I think of really exotic number types (Boolean rings, Gaussian integers, Octonions, oh my!). I imagined it would be a language feature, so I could do something like this:

use gaussianintegers as numbers

x = 1 + i
y = 2 - 3i
print(x*y)

z = 2 + 3.5i     # error

I’m not sure why I thought this would be so cool. Perhaps I felt like I would be teaching a computer math. Or maybe the next level of abstraction in playing god by writing programs is to play god by designing languages (and I secretly satisfy a massive god complex by dictating the actions of my computer).

But despite not writing a language of my own, programming with weird number systems still has a special place in my heart. It just so happens that we’re in the middle of a long series on elliptic curves, and in the next post we’ll actually implement elliptic curve arithmetic over a special kind of number type (the finite field). In this post, we’ll lay the groundwork by implementing number types in Python that allow us to work over any finite field. This is actually a pretty convoluted journey, and to be totally rigorous we’d need to prove a bunch of lemmas, develop a bunch of ring theory, and prove the correctness of a few algorithms involving polynomials.

Instead of taking the long and winding road, we’ll just state the important facts with links to proofs, prove the easy stuff, and focus more heavily than usual on the particular Python implementation details. As usual, all of the code used in this post is available on this blog’s Github page.

Integers Modulo Primes

The simples kind of finite field is the set of integers modulo a prime. We’ve dealt with this number field extensively on this blog (in groups, rings, fields, with RSA, etc.), but let’s recall what it is. The modulo operator \mod (in programming it’s often denoted %) is a binary operation on integers such that x \mod y is the unique positive remainder of x when divided by y.

Definition: Let p be a prime number. The set \mathbb{Z}/p consists of the numbers \left \{ 0, 1, \dots, p-1 \right \}. If you endow it with the operations of addition (mod p) and multiplication (mod p), it forms a field.

To say it’s a field is just to say that arithmetic more or less behaves the way we expect it to, and in particular that every nonzero element has a (unique) multiplicative inverse. Making a number type for \mathbb{Z}/p in Python is quite simple.

def IntegersModP(p):
   class IntegerModP(FieldElement):
      def __init__(self, n):
         self.n = n % p
         self.field = IntegerModP

      def __add__(self, other): return IntegerModP(self.n + other.n)
      def __sub__(self, other): return IntegerModP(self.n - other.n)
      def __mul__(self, other): return IntegerModP(self.n * other.n)
      def __truediv__(self, other): return self * other.inverse()
      def __div__(self, other): return self * other.inverse()
      def __neg__(self): return IntegerModP(-self.n)
      def __eq__(self, other): return isinstance(other, IntegerModP) and self.n == other.n
      def __abs__(self): return abs(self.n)
      def __str__(self): return str(self.n)
      def __repr__(self): return '%d (mod %d)' % (self.n, self.p)

      def __divmod__(self, divisor):
         q,r = divmod(self.n, divisor.n)
         return (IntegerModP(q), IntegerModP(r))

      def inverse(self):
         ...?

   IntegerModP.p = p
   IntegerModP.__name__ = 'Z/%d' % (p)
   return IntegerModP

We’ve done a couple of things worth note here. First, all of the double-underscore methods are operator overloads, so they are called when one tries to, e.g., add two instances of this class together. We’ve also implemented a division algorithm via __divmod__ which computes a (quotient, remainder) pair for the appropriate division. The built in Python function divmod function does this for integers, and we can overload it for a custom type. We’ll write a more complicated division algorithm later in this post. Finally, we’re dynamically creating our class so that different primes will correspond to different types. We’ll come back to why this encapsulation is a good idea later, but it’s crucial to make our next few functions reusable and elegant.

Here’s an example of the class in use:

>>> mod7 = IntegersModP(7)
>>> mod7(3) + mod7(6)
2 (mod 7)

The last (undefined) function in the IntegersModP class, the inverse function, is our only mathematical hurdle. Luckily, we can compute inverses in a generic way, using an algorithm called the extended Euclidean algorithm. Here’s the mathematics.

Definition: An element d is called a greatest common divisor (gcd) of a,b if it divides both a and b, and for every other z dividing both a and b, z divides d. For \mathbb{Z}/p gcd’s and we denote it as \gcd(a,b). [1]

Note that we called it ‘a’ greatest common divisor. In general gcd’s need not be unique, though for integers one often picks the positive gcd. We’ll actually see this cause a tiny programmatic bug later in this post, but let’s push on for now.

Theorem: For any two integers a,b \in \mathbb{Z} there exist unique x,y \in \mathbb{Z} such that ax + by = \gcd(a,b).

We could beat around the bush and try to prove these things in various ways, but when it comes down to it there’s one algorithm of central importance that both computes the gcd and produces the needed linear combination x,y. The algorithm is called the Euclidean algorithm. Here is a simple version that just gives the gcd.

def gcd(a, b):
   if abs(a) < abs(b):
      return gcd(b, a)

   while abs(b) > 0:
      q,r = divmod(a,b)
      a,b = b,r

   return a

This works by the simple observation that \gcd(a, aq+r) = \gcd(a,r) (this is an easy exercise to prove directly). So the Euclidean algorithm just keeps applying this rule over and over again: take the remainder when dividing the bigger argument by the smaller argument until the remainder becomes zero. Then the \gcd(x,0) = x because everything divides zero.

Now the so-called ‘extended’ Euclidean algorithm just keeps track of some additional data as it goes (the partial quotients and remainders). Here’s the algorithm.

def extendedEuclideanAlgorithm(a, b):
   if abs(b) > abs(a):
      (x,y,d) = extendedEuclideanAlgorithm(b, a)
      return (y,x,d)

   if abs(b) == 0:
      return (1, 0, a)

   x1, x2, y1, y2 = 0, 1, 1, 0
   while abs(b) > 0:
      q, r = divmod(a,b)
      x = x2 - q*x1
      y = y2 - q*y1
      a, b, x2, x1, y2, y1 = b, r, x1, x, y1, y

   return (x2, y2, a)

Indeed, the reader who hasn’t seen this stuff before is encouraged to trace out a run for the numbers 4864, 3458. Their gcd is 38 and the two integers are 32 and -45, respectively.

How does this help us compute inverses? Well, if we want to find the inverse of a modulo p, we know that their gcd is 1. So compute the x,y such that ax + py = 1, and then reduce both sides mod p. You get ax + 0 = 1 \mod p, which means that x \mod p is the inverse of a. So once we have the extended Euclidean algorithm our inverse function is trivial to write!

def inverse(self):
   x,y,d = extendedEuclideanAlgorithm(self.n, self.p)
   return IntegerModP(x)

And indeed it works as expected:

>>> mod23 = IntegersModP(23)
>>> mod23(7).inverse()
10 (mod 23)
>>> mod23(7).inverse() * mod23(7)
1 (mod 23)

Now one very cool thing, and something that makes some basic ring theory worth understanding, is that we can compute the gcd of any number type using the exact same code for the Euclidean algorithm, provided we implement an abs function and a division algorithm. Via a chain of relatively easy-to-prove lemmas, if your number type has enough structure (in particular, if it has a division algorithm that satisfies some properties), then greatest common divisors are well-defined, and the Euclidean algorithm gives us that special linear combination. And using the same trick above in finite fields, we can use the Euclidean algorithm to compute inverses.

But in order to make things work programmatically we need to be able to deal with the literal ints 0 and 1 in the algorithm. That is, we need to be able to silently typecast integers up to whatever number type we’re working with. This makes sense because all rings have 0 and 1, but it requires a bit of scaffolding to implement. In particular, typecasting things sensibly is really difficult if you aren’t careful. And the problems are compounded in a language like Python that blatantly ignores types whenever possible. [2]

So let’s take a quick break to implement a tiny type system with implicit typecasting.

[1] The reader familiar with our series on category theory will recognize this as the product of two integers in a category whose arrows represent divisibility. So by abstract nonsense, this proves that gcd’s are unique up to multiplication by a unit in any ring.
[2] In the process of writing the code for this post, I was sorely missing the stronger type systems of Java and Haskell. Never thought I’d say that, but it’s true.

A Tiny Generic Type System

The main driving force behind our type system will be a decorator called @typecheck. We covered decorators toward the end of our primer on dynamic programming, but in short a decorator is a Python syntax shortcut that allows some pre- or post-processing to happen to a function in a reusable way. All you need to do to apply the pre/post-processing is prefix the function definition with the name of the decorator.

Our decorator will be called typecheck, and it will decorate binary operations on our number types. In its basic form, our type checker will work as follows: if the types of the two operands are the same, then the decorator will just pass them on through to the operator. Otherwise, it will try to do some typecasting, and if that fails it will raise exceptions with reckless abandon.

def typecheck(f):
   def newF(self, other):
      if type(self) is not type(other):
         try:
            other = self.__class__(other)
         except TypeError:
            message = 'Not able to typecast %s of type %s to type %s in function %s'
            raise TypeError(message % (other, type(other).__name__, type(self).__name__, f.__name__))
         except Exception as e:
            message = 'Type error on arguments %r, %r for functon %s. Reason:%s'
            raise TypeError(message % (self, other, f.__name__, type(other).__name__, type(self).__name__, e))

      return f(self, other)

   return newF

So this is great, but there are two issues. The first is that this will only silently typecast if the thing we’re casting is on the right-hand side of the expression. In other words, the following will raise an exception complaining that you can’t add ints to Mod7 integers.

>>> x = IntegersModP(7)(1)
>>> 1 + x

What we need are the right-hand versions of all the operator overloads. They are the same as the usual operator overloads, but Python gives preference to the left-hand operator overloads. Anticipating that we will need to rewrite these silly right-hand overloads for every number type, and they’ll all be the same, we make two common base classes.

class DomainElement(object):
   def __radd__(self, other): return self + other
   def __rsub__(self, other): return -self + other
   def __rmul__(self, other): return self * other

class FieldElement(DomainElement):
   def __truediv__(self, other): return self * other.inverse()
   def __rtruediv__(self, other): return self.inverse() * other
   def __div__(self, other): return self.__truediv__(other)
   def __rdiv__(self, other): return self.__rtruediv__(other)

And we can go ahead and make our IntegersModP a subclass of FieldElement. [3]

But now we’re wading into very deep waters. In particular, we know ahead of time that our next number type will be for Polynomials (over the integers, or fractions, or \mathbb{Z}/p, or whatever). And we’ll want to do silent typecasting from ints and IntegersModP to Polynomials! The astute reader will notice the discrepancy. What will happen if I try to do this?

>>> MyInteger() + MyPolynomial()

Let’s take this slowly: by our assumption both MyInteger and MyPolynomial have the __add__ and __radd__ functions defined on them, and each tries to typecast the other the appropriate type. But which is called? According to Python’s documentation if the left-hand side has an __add__ function that’s called first, and the right-hand sides’s __radd__ function is only sought if no __add__ function is found for the left operand.

Well that’s a problem, and we’ll deal with it in a half awkward and half elegant way. What we’ll do is endow our number types with an “operatorPrecedence” constant. And then inside our type checker function we’ll see if the right-hand operand is an object of higher precedence. If it is, we return the global constant NotImplemented, which Python takes to mean that no __add__ function was found, and it proceeds to look for __radd__. And so with this modification our typechecker is done. [4]

def typecheck(f):
   def newF(self, other):
      if (hasattr(other.__class__, 'operatorPrecedence') and
            other.__class__.operatorPrecedence > self.__class__.operatorPrecedence):
         return NotImplemented

      if type(self) is not type(other):
         try:
            other = self.__class__(other)
         except TypeError:
            message = 'Not able to typecast %s of type %s to type %s in function %s'
            raise TypeError(message % (other, type(other).__name__, type(self).__name__, f.__name__))
         except Exception as e:
            message = 'Type error on arguments %r, %r for functon %s. Reason:%s'
            raise TypeError(message % (self, other, f.__name__, type(other).__name__, type(self).__name__, e))

      return f(self, other)

   return newF

We add a default operatorPrecedence of 1 to the DomainElement base class. Now this function answers our earlier question of why we want to encapsulate the prime modulus into the IntegersModP class. If this typechecker is really going to be generic, we need to be able to typecast an int by passing the single int argument to the type constructor with no additional information! Indeed, this will be the same pattern for our polynomial class and the finite field class to follow.

Now there is still one subtle problem. If we try to generate two copies of the same number type from our number-type generator (in other words, the following code snippet), we’ll get a nasty exception.

>>> mod7 = IntegersModP(7)
>>> mod7Copy = IntegersModP(7)
>>> mod7(1) + mod7Copy(2)
... fat error ...

The reason for this is that in the type-checker we’re using the Python built-in ‘is’ which checks for identity, not semantic equality. To fix this, we simply need to memoize the IntegersModP function (and all the other functions we’ll use to generate number types) so that there is only ever one copy of a number type in existence at a time.

So enough Python hacking: let’s get on with implementing finite fields!

[3] This also compels us to make some slight modifications to the constructor for IntegersModP, but they’re not significant enough to display here. Check out the Github repo if you want to see.
[4] This is truly a hack, and we’ve considered submitting a feature request to the Python devs. It is conceivably useful for the operator-overloading aficionado. I’d be interested to hear your thoughts in the comments as to whether this is a reasonable feature to add to Python.

Polynomial Arithmetic

Recall from our finite field primer that every finite field can be constructed as a quotient of a polynomial ring with coefficients in \mathbb{Z}/p by some prime ideal. We spelled out exactly what this means in fine detail in the primer, so check that out before reading on.

Indeed, to construct a finite field we need to find some irreducible monic polynomial f with coefficients in \mathbb{Z}/p, and then the elements of our field will be remainders of arbitrary polynomials when divided by f. In order to determine if they’re irreducible we’ll need to compute a gcd. So let’s build a generic polynomial type with a polynomial division algorithm, and hook it into our gcd framework.

We start off in much the same way as with the IntegersModP:

# create a polynomial with coefficients in a field; coefficients are in
# increasing order of monomial degree so that, for example, [1,2,3]
# corresponds to 1 + 2x + 3x^2
@memoize
def polynomialsOver(field=fractions.Fraction):

   class Polynomial(DomainElement):
      operatorPrecedence = 2
      factory = lambda L: Polynomial([field(x) for x in L])

      def __init__(self, c):
         if type(c) is Polynomial:
            self.coefficients = c.coefficients
         elif isinstance(c, field):
            self.coefficients = [c]
         elif not hasattr(c, '__iter__') and not hasattr(c, 'iter'):
            self.coefficients = [field(c)]
         else:
            self.coefficients = c

         self.coefficients = strip(self.coefficients, field(0))

      def isZero(self): return self.coefficients == []

      def __repr__(self):
         if self.isZero():
            return '0'

         return ' + '.join(['%s x^%d' % (a,i) if i > 0 else '%s'%a
                              for i,a in enumerate(self.coefficients)])

      def __abs__(self): return len(self.coefficients)
      def __len__(self): return len(self.coefficients)
      def __sub__(self, other): return self + (-other)
      def __iter__(self): return iter(self.coefficients)
      def __neg__(self): return Polynomial([-a for a in self])

      def iter(self): return self.__iter__()
      def leadingCoefficient(self): return self.coefficients[-1]
      def degree(self): return abs(self) - 1

All of this code just lays down conventions. A polynomial is a list of coefficients (in increasing order of their monomial degree), the zero polynomial is the empty list of coefficients, and the abs() of a polynomial is one plus its degree. [5] Finally, instead of closing over a prime modulus, as with IntegersModP, we’re closing over the field of coefficients. In general you don’t have to have polynomials with coefficients in a field, but if they do come from a field then you’re guaranteed to get a sensible Euclidean algorithm. In the formal parlance, if k is a field then k[x] is a Euclidean domain. And for our goal of defining finite fields, we will always have coefficients from \mathbb{Z}/p, so there’s no problem.

Now we can define things like addition, multiplication, and equality using our typechecker to silently cast and watch for errors.

      @typecheck
      def __eq__(self, other):
         return self.degree() == other.degree() and all([x==y for (x,y) in zip(self, other)])

      @typecheck
      def __add__(self, other):
         newCoefficients = [sum(x) for x in itertools.zip_longest(self, other, fillvalue=self.field(0))]
         return Polynomial(newCoefficients)

      @typecheck
      def __mul__(self, other):
         if self.isZero() or other.isZero():
            return Zero()

         newCoeffs = [self.field(0) for _ in range(len(self) + len(other) - 1)]

         for i,a in enumerate(self):
            for j,b in enumerate(other):
               newCoeffs[i+j] += a*b

         return Polynomial(newCoeffs)

Notice that, if the underlying field of coefficients correctly implements the operator overloads, none of this depends on the coefficients. Reusability, baby!

And we can finish off with the division algorithm for polynomials.

      @typecheck
      def __divmod__(self, divisor):
         quotient, remainder = Zero(), self
         divisorDeg = divisor.degree()
         divisorLC = divisor.leadingCoefficient()

         while remainder.degree() >= divisorDeg:
            monomialExponent = remainder.degree() - divisorDeg
            monomialZeros = [self.field(0) for _ in range(monomialExponent)]
            monomialDivisor = Polynomial(monomialZeros + [remainder.leadingCoefficient() / divisorLC])

            quotient += monomialDivisor
            remainder -= monomialDivisor * divisor

         return quotient, remainder

Indeed, we’re doing nothing here but formalizing the grade-school algorithm for doing polynomial long division [6]. And we can finish off the function for generating this class by assigning the field member variable along with a class name. And we give it a higher operator precedence than the underlying field of coefficients so that an isolated coefficient is cast up to a constant polynomial.

@memoize
def polynomialsOver(field=fractions.Fraction):

   class Polynomial(DomainElement):
      operatorPrecedence = 2

      [... methods defined above ...]

   def Zero():
      return Polynomial([])

   Polynomial.field = field
   Polynomial.__name__ = '(%s)[x]' % field.__name__
   return Polynomial

We provide a modest test suite in the Github repository for this post, but here’s a sample test:

>>> Mod5 = IntegersModP(5)
>>> Mod11 = IntegersModP(11)
>>> polysOverQ = polynomialsOver(Fraction).factory
>>> polysMod5 = polynomialsOver(Mod5).factory
>>> polysMod11 = polynomialsOver(Mod11).factory
>>> polysOverQ([1,7,49]) / polysOverQ([7])
1/7 + 1 x^1 + 7 x^2
>>> polysMod5([1,7,49]) / polysMod5([7])
3 + 1 x^1 + 2 x^2
>>> polysMod11([1,7,49]) / polysMod11([7])
8 + 1 x^1 + 7 x^2

And indeed, the extended Euclidean algorithm works without modification, so we know our typecasting is doing what’s expected:

>>> p = polynomialsOver(Mod2).factory
>>> f = p([1,0,0,0,1,1,1,0,1,1,1]) # x^10 + x^9 + x^8 + x^6 + x^5 + x^4 + 1
>>> g = p([1,0,1,1,0,1,1,0,0,1])   # x^9 + x^6 + x^5 + x^3 + x^1 + 1
>>> theGcd = p([1,1,0,1]) # x^3 + x + 1
>>> x = p([0,0,0,0,1]) # x^4
>>> y = p([1,1,1,1,1,1]) # x^5 + x^4 + x^3 + x^2 + x + 1
>>> (x,y,theGcd) == extendedEuclideanAlgorithm(f, g)
True

[5] The mathematical name for the abs() function that we’re using is a valuation.
[6] One day we will talk a lot more about polynomial long division on this blog. You can do a lot of cool algebraic geometry with it, and the ideas there lead you to awesome applications like robot motion planning and automated geometry theorem proving.

Generating Irreducible Polynomials

Now that we’ve gotten Polynomials out of the way, we need to be able to generate irreducible polynomials over \mathbb{Z}/p of any degree we want. It might be surprising that irreducible polynomials of any degree exist [7], but in fact we know a lot more.

Theorem: The product of all irreducible monic polynomials of degree dividing m is equal to x^{p^m} - x.

This is an important theorem, but it takes a little bit more field theory than we have under our belts right now. You could summarize the proof by saying there is a one-to-one correspondence between elements of a field and monic irreducible polynomials, and then you say some things about splitting fields. You can see a more detailed proof outline here, but it assumes you’re familiar with the notion of a minimal polynomial. We will probably cover this in a future primer.

But just using the theorem we can get a really nice algorithm for determining if a polynomial f(x) of degree m is irreducible: we just look at its gcd with all the x^{p^k} - x for k smaller than m. If all the gcds are constants, then we know it’s irreducible, and if even one is a non-constant polynomial then it has to be irreducible. Why’s that? Because if you have some nontrivial gcd d(x) = \gcd(f(x), x^{p^k} - x) for k < m, then it’s a factor of f(x) by definition. And since we know all irreducible monic polynomials are factors of that this collection of polynomials, if the gcd is always 1 then there are no other possible factors to be divisors. (If there is any divisor then there will be a monic irreducible one!) So the candidate polynomial must be irreducible. In fact, with a moment of thought it’s clear that we can stop at k= m/2, as any factor of large degree will necessarily require corresponding factors of small degree. So the algorithm to check for irreducibility is just this simple loop:

def isIrreducible(polynomial, p):
   ZmodP = IntegersModP(p)
   poly = polynomialsOver(ZmodP).factory
   x = poly([0,1])
   powerTerm = x
   isUnit = lambda p: p.degree() == 0

   for _ in range(int(polynomial.degree() / 2)):
      powerTerm = powerTerm.powmod(p, polynomial)
      gcdOverZmodp = gcd(polynomial, powerTerm - x)
      if not isUnit(gcdOverZmodp):
         return False

   return True

We’re just computing the powers iteratively as x^p, (x^p)^p = x^{p^2}, \dots, x^{p^j} and in each step of the loop subtracting x and computing the relevant gcd. The powmod function is just there so that we can reduce the power mod our irreducible polynomial after each multiplication, keeping the degree of the polynomial small and efficient to work with.

Now generating an irreducible polynomial is a little bit harder than testing for one. With a lot of hard work, however, field theorists discovered that irreducible polynomials are quite common. In fact, if you just generate the coefficients of your degree n monic polynomial at random, the chance that you’ll get something irreducible is at least 1/n. So this suggests an obvious randomized algorithm: keep guessing until you find one.

def generateIrreduciblePolynomial(modulus, degree):
   Zp = IntegersModP(modulus)
   Polynomial = polynomialsOver(Zp)

   while True:
      coefficients = [Zp(random.randint(0, modulus-1)) for _ in range(degree)]
      randomMonicPolynomial = Polynomial(coefficients + [Zp(1)])

      if isIrreducible(randomMonicPolynomial, modulus):
         return randomMonicPolynomial

Since the probability of getting an irreducible polynomial is close to 1/n, we expect to require n trials before we find one. Moreover we could give a pretty tight bound on how likely it is to deviate from the expected number of trials. So now we can generate some irreducible polynomials!

>>> F23 = FiniteField(2,3)
>>> generateIrreduciblePolynomial(23, 3)
21 + 12 x^1 + 11 x^2 + 1 x^3

And so now we are well-equipped to generate any finite field we want! It’s just a matter of generating the polynomial and taking a modulus after every operation.

@memoize
def FiniteField(p, m, polynomialModulus=None):
   Zp = IntegersModP(p)
   if m == 1:
      return Zp

   Polynomial = polynomialsOver(Zp)
   if polynomialModulus is None:
      polynomialModulus = generateIrreduciblePolynomial(modulus=p, degree=m)

   class Fq(FieldElement):
      fieldSize = int(p ** m)
      primeSubfield = Zp
      idealGenerator = polynomialModulus
      operatorPrecedence = 3

      def __init__(self, poly):
         if type(poly) is Fq:
            self.poly = poly.poly
         elif type(poly) is int or type(poly) is Zp:
            self.poly = Polynomial([Zp(poly)])
         elif isinstance(poly, Polynomial):
            self.poly = poly % polynomialModulus
         else:
            self.poly = Polynomial([Zp(x) for x in poly]) % polynomialModulus

         self.field = Fq

      @typecheck
      def __add__(self, other): return Fq(self.poly + other.poly)
      @typecheck
      def __sub__(self, other): return Fq(self.poly - other.poly)
      @typecheck
      def __mul__(self, other): return Fq(self.poly * other.poly)
      @typecheck
      def __eq__(self, other): return isinstance(other, Fq) and self.poly == other.poly

      def __pow__(self, n): return Fq(pow(self.poly, n))
      def __neg__(self): return Fq(-self.poly)
      def __abs__(self): return abs(self.poly)
      def __repr__(self): return repr(self.poly) + ' \u2208 ' + self.__class__.__name__

      @typecheck
      def __divmod__(self, divisor):
         q,r = divmod(self.poly, divisor.poly)
         return (Fq(q), Fq(r))

      def inverse(self):
         if self == Fq(0):
            raise ZeroDivisionError

         x,y,d = extendedEuclideanAlgorithm(self.poly, self.idealGenerator)
         return Fq(x) * Fq(d.coefficients[0].inverse())

   Fq.__name__ = 'F_{%d^%d}' % (p,m)
   return Fq

And some examples of using it:

>>> F23 = FiniteField(2,3)
>>> x = F23([1,1])
>>> x
1 + 1 x^1 ∈ F_{2^3}
>>> x*x
1 + 0 x^1 + 1 x^2 ∈ F_{2^3}
>>> x**10
0 + 0 x^1 + 1 x^2 ∈ F_{2^3}
>>> 1 / x
0 + 1 x^1 + 1 x^2 ∈ F_{2^3}
>>> x * (1 / x)
1 ∈ F_{2^3}
>>> k = FiniteField(23, 4)
>>> k.fieldSize
279841
>>> k.idealGenerator
6 + 8 x^1 + 10 x^2 + 10 x^3 + 1 x^4
>>> y
9 + 21 x^1 + 14 x^2 + 12 x^3 ∈ F_{23^4}
>>> y*y
13 + 19 x^1 + 7 x^2 + 14 x^3 ∈ F_{23^4}
>>> y**5 - y
15 + 22 x^1 + 15 x^2 + 5 x^3 ∈ F_{23^4}

And that’s it! Now we can do arithmetic over any finite field we want.

[7] Especially considering that other wacky things happen like this: x^4 +1 is reducible over every finite field!

Some Words on Efficiency

There are a few things that go without stating about this program, but I’ll state them anyway.

The first is that we pay a big efficiency price for being so generic. All the typecasting we’re doing isn’t cheap, and in general cryptography needs to be efficient. For example, if I try to create a finite field of order 104729^{20}, it takes about ten to fifteen seconds to complete. This might not seem so bad for a one-time initialization, but it’s clear that our representation is somewhat bloated. We would display a graph of the expected time to perform various operations in various finite fields, but this post is already way too long.

In general, the larger and more complicated the polynomial you use to define your finite field, the longer operations will take (since dividing by a complicated polynomial is more expensive than dividing by a simple polynomial). For this reason and a few other reasons, a lot of research has gone into efficiently finding irreducible polynomials with, say, only three nonzero terms. Moreover, if we know in advance that we’ll only work over fields of characteristic two we can make a whole lot of additional optimizations. Essentially, all of the arithmetic reduces to really fast bitwise operations, and things like exponentiation can easily be implemented in hardware. But it also seems that the expense coming with large field characteristics corresponds to higher security, so some researchers have tried to meet in the middle an get efficient representations of other field characteristics.

In any case, the real purpose of our implementation in this post is not for efficiency. We care first and foremost about understanding the mathematics, and to have a concrete object to play with and use in the future for other projects. And we have accomplished just that.

Until next time!

Martingales and the Optional Stopping Theorem

This is a guest post by my colleague Adam Lelkes.

The goal of this primer is to introduce an important and beautiful tool from probability theory, a model of fair betting games called martingales. In this post I will assume that the reader is familiar with the basics of probability theory. For those that need to refresh their knowledge, Jeremy’s excellent primers (1, 2) are a good place to start.

The Geometric Distribution and the ABRACADABRA Problem

Before we start playing with martingales, let’s start with an easy exercise. Consider the following experiment: we throw an ordinary die repeatedly until the first time a six appears. How many throws will this take in expectation? The reader might recognize immediately that this exercise can be easily solved using the basic properties of the geometric distribution, which models this experiment exactly. We have independent trials, every trial succeeding with some fixed probability p. If X denotes the number of trials needed to get the first success, then clearly \Pr(X = k) = (1-p)^{k-1} p (since first we need k-1 failures which occur independently with probability 1-p, then we need one success which happens with probability p). Thus the expected value of X is

\displaystyle E(X) = \sum_{k=1}^\infty k P(X = k) = \sum_{k=1}^\infty k (1-p)^{k-1} p = \frac1p

by basic calculus. In particular, if success is defined as getting a six, then p=1/6 thus the expected time is 1/p=6.

Now let us move on to a somewhat similar, but more interesting and difficult problem, the ABRACADABRA problem. Here we need two things for our experiment, a monkey and a typewriter. The monkey is asked to start bashing random keys on a typewriter. For simplicity’s sake, we assume that the typewriter has exactly 26 keys corresponding to the 26 letters of the English alphabet and the monkey hits each key with equal probability. There is a famous theorem in probability, the infinite monkey theorem, that states that given infinite time, our monkey will almost surely type the complete works of William Shakespeare. Unfortunately, according to astronomists the sun will begin to die in a few billion years, and the expected time we need to wait until a monkey types the complete works of William Shakespeare is orders of magnitude longer, so it is not feasible to use monkeys to produce works of literature.

So let’s scale down our goals, and let’s just wait until our monkey types the word ABRACADABRA. What is the expected time we need to wait until this happens? The reader’s first idea might be to use the geometric distribution again. ABRACADABRA is eleven letters long, the probability of getting one letter right is \frac{1}{26}, thus the probability of a random eleven-letter word being ABRACADABRA is exactly \left(\frac{1}{26}\right)^{11}. So if typing 11 letters is one trial, the expected number of trials is

\displaystyle \frac1{\left(\frac{1}{26}\right)^{11}}=26^{11}

which means 11\cdot 26^{11} keystrokes, right?

Well, not exactly. The problem is that we broke up our random string into eleven-letter blocks and waited until one block was ABRACADABRA. However, this word can start in the middle of a block. In other words, we considered a string a success only if the starting position of the word ABRACADABRA was divisible by 11. For example, FRZUNWRQXKLABRACADABRA would be recognized as success by this model but the same would not be true for AABRACADABRA. However, it is at least clear from this observation that 11\cdot 26^{11} is a strict upper bound for the expected waiting time. To find the exact solution, we need one very clever idea, which is the following:

Let’s Open a Casino!

Do I mean that abandoning our monkey and typewriter and investing our time and money in a casino is a better idea, at least in financial terms? This might indeed be the case, but here we will use a casino to determine the expected wait time for the ABRACADABRA problem. Unfortunately we won’t make any money along the way (in expectation) since our casino will be a fair one.

Let’s do the following thought experiment: let’s open a casino next to our typewriter. Before each keystroke, a new gambler comes to our casino and bets $1 that the next letter will be A. If he loses, he goes home disappointed. If he wins, he bets all the money he won on the event that the next letter will be B. Again, if he loses, he goes home disappointed. (This won’t wreak havoc on his financial situation, though, as he only loses $1 of his own money.) If he wins again, he bets all the money on the event that the next letter will be R, and so on.

If a gambler wins, how much does he win? We said that the casino would be fair, i.e. the expected outcome should be zero. That means that it the gambler bets $1, he should receive $26 if he wins, since the probability of getting the next letter right is exactly \frac{1}{26} (thus the expected value of the change in the gambler’s fortune is \frac{25}{26}\cdot (-1) + \frac{1}{26}\cdot (+25) = 0.

Let’s keep playing this game until the word ABRACADABRA first appears and let’s denote the number of keystrokes up to this time as T. As soon as we see this word, we close our casino. How much was the revenue of our casino then? Remember that before each keystroke, a new gambler comes in and bets $1, and if he wins, he will only bet the money he has received so far, so our revenue will be exactly T dollars.

How much will we have to pay for the winners? Note that the only winners in the last round are the players who bet on A. How many of them are there? There is one that just came in before the last keystroke and this was his first bet. He wins $26. There was one who came three keystrokes earlier and he made four successful bets (ABRA). He wins \$26^4. Finally there is the luckiest gambler who went through the whole ABRACADABRA sequence, his prize will be \$26^{11}. Thus our casino will have to give out 26^{11}+26^4+26 dollars in total, which is just under the price of 200,000 WhatsApp acquisitions.

Now we will make one crucial observation: even at the time when we close the casino, the casino is fair! Thus in expectation our expenses will be equal to our income. Our income is T dollars, the expected value of our expenses is 26^{11}+26^4+26 dollars, thus E(T)=26^{11}+26^4+26. A beautiful solution, isn’t it? So if our monkey types at 150 characters per minute on average, we will have to wait around 47 million years until we see ABRACADABRA. Oh well.

Time to be More Formal

After giving an intuitive outline of the solution, it is time to formalize the concepts that we used, to translate our fairy tales into mathematics. The mathematical model of the fair casino is called a martingale, named after a class of betting strategies that enjoyed popularity in 18th century France. The gambler’s fortune (or the casino’s, depending on our viewpoint) can be modeled with a sequence of random variables. X_0 will denote the gambler’s fortune before the game starts, X_1 the fortune after one round and so on. Such a sequence of random variables is called a stochastic process. We will require the expected value of the gambler’s fortune to be always finite.

How can we formalize the fairness of the game? Fairness means that the gambler’s fortune does not change in expectation, i.e. the expected value of X_n, given X_1, X_2, \ldots, X_{n-1} is the same as X_{n-1}. This can be written as E(X_n | X_1, X_2, \ldots, X_{n-1}) = X_{n-1} or, equivalently, E(X_n - X_{n-1} | X_1, X_2, \ldots, X_{n-1}) = 0.

The reader might be less comfortable with the first formulation. What does it mean, after all, that the conditional expected value of a random variable is another random variable? Shouldn’t the expected value be a number? The answer is that in order to have solid theoretical foundations for the definition of a martingale, we need a more sophisticated notion of conditional expectations. Such sophistication involves measure theory, which is outside the scope of this post. We will instead naively accept the definition above, and the reader can look up all the formal details in any serious probability text (such as [1]).

Clearly the fair casino we constructed for the ABRACADABRA exercise is an example of a martingale. Another example is the simple symmetric random walk on the number line: we start at 0, toss a coin in each step, and move one step in the positive or negative direction based on the outcome of our coin toss.

The Optional Stopping Theorem

Remember that we closed our casino as soon as the word ABRACADABRA appeared and we claimed that our casino was also fair at that time. In mathematical language, the closed casino is called a stopped martingale. The stopped martingale is constructed as follows: we wait until our martingale X exhibits a certain behaviour (e.g. the word ABRACADABRA is typed by the monkey), and we define a new martingale X’ as follows: let X'_n = X_n if n < T and X'_n = X_T if n \ge T where T denotes the stopping time, i.e. the time at which the desired event occurs. Notice that T itself is a random variable.

We require our stopping time T to depend only on the past, i.e. that at any time we should be able to decide whether the event that we are waiting for has already happened or not (without looking into the future). This is a very reasonable requirement. If we could look into the future, we could obviously cheat by closing our casino just before some gambler would win a huge prize.

We said that the expected wealth of the casino at the stopping time is the same as the initial wealth. This is guaranteed by Doob’s optional stopping theorem, which states that under certain conditions, the expected value of a martingale at the stopping time is equal to its expected initial value.

Theorem: (Doob’s optional stopping theorem) Let X_n be a martingale stopped at step T, and suppose one of the following three conditions hold:

  1. The stopping time T is almost surely bounded by some constant;
  2. The stopping time T is almost surely finite and every step of the stopped martingale X_n is almost surely bounded by some constant; or
  3. The expected stopping time E(T) is finite and the absolute value of the martingale increments |X_n-X_{n-1}| are almost surely bounded by a constant.

Then E(X_T) = E(X_0).

We omit the proof because it requires measure theory, but the interested reader can see it in these notes.

For applications, (1) and (2) are the trivial cases. In the ABRACADABRA problem, the third condition holds: the expected stopping time is finite (in fact, we showed using the geometric distribution that it is less than 26^{12}) and the absolute value of a martingale increment is either 1 or a net payoff which is bounded by 26^{11}+26^4+26. This shows that our solution is indeed correct.

Gambler’s Ruin

Another famous application of martingales is the gambler’s ruin problem. This problem models the following game: there are two players, the first player has a dollars, the second player has b dollars. In each round they toss a coin and the loser gives one dollar to the winner. The game ends when one of the players runs out of money. There are two obvious questions: (1) what is the probability that the first player wins and (2) how long will the game take in expectation?

Let X_n denote the change in the second player’s fortune, and set X_0 = 0. Let T_k denote the first time s when X_s = k. Then our first question can be formalized as trying to determine \Pr(T_{-b} < T_a). Let t = \min \{ T_{-b}, T_a\}. Clearly t is a stopping time. By the optional stopping theorem we have that

\displaystyle 0=E(X_0)=E(X_t)=-b\Pr(T_{-b} < T_a)+a(1-\Pr(T_{-b} < T_a))

thus \Pr(T_{-b} < T_a)=\frac{a}{a+b}.

I would like to ask the reader to try to answer the second question. It is a little bit trickier than the first one, though, so here is a hint: X_n^2-n is also a martingale (prove it), and applying the optional stopping theorem to it leads to the answer.

A Randomized Algorithm for 2-SAT

The reader is probably familiar with 3-SAT, the first problem shown to be NP-complete. Recall that 3-SAT is the following problem: given a boolean formula in conjunctive normal form with at most three literals in each clause, decide whether there is a satisfying truth assignment. It is natural to ask if or why 3 is special, i.e. why don’t we work with k-SAT for some k \ne 3 instead? Clearly the hardness of the problem is monotone increasing in k since k-SAT is a special case of (k+1)-SAT. On the other hand, SAT (without any bound on the number of literals per clause) is clearly in NP, thus 3-SAT is just as hard as k-SAT for any k>3. So the only question is: what can we say about 2-SAT?

It turns out that 2-SAT is easier than satisfiability in general: 2-SAT is in P. There are many algorithms for solving 2-SAT. Here is one deterministic algorithm: associate a graph to the 2-SAT instance such that there is one vertex for each variable and each negated variable and the literals x and y are connected by a directed edge if there is a clause (\bar x \lor y). Recall that \bar x \lor y is equivalent to x \implies y, so the edges show the implications between the variables. Clearly the 2-SAT instance is not satisfiable if there is a variable x such that there are directed paths x \to \bar x and \bar x \to x (since x \Leftrightarrow \bar x is always false). It can be shown that this is not only a sufficient but also a necessary condition for unsatisfiability, hence the 2-SAT instance is satisfiable if and only if there is are no such path. If there are directed paths from one vertex of a graph to another and vice versa then they are said to belong to the same strongly connected component. There are several graph algorithms for finding strongly connected components of directed graphs, the most well-known algorithms are all based on depth-first search.

Now we give a very simple randomized algorithm for 2-SAT (due to Christos Papadimitriou in a ’91 paper): start with an arbitrary truth assignment and while there are unsatisfied clauses, pick one and flip the truth value of a random literal in it. Stop after O(n^2) rounds where n denotes the number of variables. Clearly if the formula is not satisfiable then nothing can go wrong, we will never find a satisfying truth assignment. If the formula is satisfiable, we want to argue that with high probability we will find a satisfying truth assignment in O(n^2) steps.

The idea of the proof is the following: fix an arbitrary satisfying truth assignment and consider the Hamming distance of our current assignment from it. The Hamming distance of two truth assignments (or in general, of two binary vectors) is the number of coordinates in which they differ. Since we flip one bit in every step, this Hamming distance changes by \pm 1 in every round. It also easy to see that in every step the distance is at least as likely to be decreased as to be increased (since we pick an unsatisfied clause, which means at least one of the two literals in the clause differs in value from the satisfying assignment).

Thus this is an unfair “gambler’s ruin” problem where the gambler’s fortune is the Hamming distance from the solution, and it decreases with probability at least \frac{1}{2}. Such a stochastic process is called a supermartingale — and this is arguably a better model for real-life casinos. (If we flip the inequality, the stochastic process we get is called a submartingale.) Also, in this case the gambler’s fortune (the Hamming distance) cannot increase beyond n. We can also think of this process as a random walk on the set of integers: we start at some number and in each round we make one step to the left or to the right with some probability. If we use random walk terminology, 0 is called an absorbing barrier since we stop the process when we reach 0. The number n, on the other hand, is called a reflecting barrier: we cannot reach n+1, and whenever we get close we always bounce back.

There is an equivalent version of the optimal stopping theorem for supermartingales and submartingales, where the conditions are the same but the consequence holds with an inequality instead of equality. It follows from the optional stopping theorem that the gambler will be ruined (i.e. a satisfying truth assignment will be found) in O(n^2) steps with high probability.

[1] For a reference on stochastic processes and martingales, see the text of Durrett .