Visualizing Galois Fields (Follow-up) May 21st, 2012
Patrick Stein

In my previous article, I should have finished by remapping the multiplication and addition tables so that the multiplication table was in cyclic order. In cyclic order, the zeroth column (or row) represents zero and the (i+1)-th column (or row) (for i \neq 0) represents a^i for some generator a. As such, the multiplication table is simply 0 in the zeroth row and column and 1 + ((i+j) \mod 255) for the spot (i+1,j+1).

Once resorted by the order of the powers of a generator a, the multiplication table look the same regardless of the a. The addition, however, looks different for different generators. Below are the addition tables for two of them: a = (1+x) on the right and a = x^3 + x^6 + x^7 on the left. They make as decent a stereo as any of the other pairs so far.

 

Here’s the code that I used to generate the remapping for a given generator.

(defun make-cyclic-remapping (fn generator &optional (limit 256))
  (let ((to (make-array (list limit) :initial-element 0))
        (from (make-array (list limit) :initial-element 0))
        (used (make-hash-table :test #'equal)))
    ;; fill up the lookup tables
    (setf (gethash 0 used) t)
    (nlet fill-tables ((exp 1) (acc 1))
      (when (< exp limit)
        (setf (aref to exp) acc
              (aref from acc) exp
              (gethash acc used) t)
        (fill-tables (1+ exp) (funcall fn acc generator))))
    ;; return a closure around the lookup tables
    (when (= (hash-table-count used) limit)
      (lambda (direction n)
        (ecase direction
          (:to (aref to n))
          (:from (aref from n)))))))

If you’ve read it, you can probably tell by my code that I’m still under the influence of Let Over Lambda. If you haven’t read it, it is quite worth the read.

Then, I used a modified version of the draw-heightmap function which also takes in the remapping function.

(defun draw-mapped-heightmap (op map filename &optional (limit 256))
  (let ((png (make-instance 'zpng:pixel-streamed-png
                            :color-type :grayscale
                            :width limit
                            :height limit)))
    (with-open-file (stream filename
                            :direction :output
                            :if-does-not-exist :create
                            :element-type '(unsigned-byte 8))
      (zpng:start-png png stream)
      (dotimes (xx limit)
        (dotimes (yy limit)
          (let* ((a (funcall map :to xx))
                 (b (funcall map :to yy))
                 (c (funcall map :from (funcall op a b))))
            (zpng:write-pixel (list c) png))))
      (zpng:finish-png png)))
  (values))

Visualizing Galois Fields May 17th, 2012
Patrick Stein

Galois fields are used in a number of different ways. For example, the AES encryption standard uses them.

Arithmetic in Galois Fields

The Galois fields of size 2^n for various n are convenient in computer applications because of how nicely they fit into bytes or words on the machine. The Galois field GF(2^n) has 2^n elements. These elements are represented as polynomials of degree less than n with all coefficients either 0 or 1. So, to encode an element of GF(2^n) as a number, you need an n-bit binary number.

For example, let us consider the Galois field GF(2^3). It has 2^3 = 8 elements. They are (as binary integers) 000, 001, 010, 011, 100, 101, 110, and 111. The element 110, for example, stands for the polynomial x^2 + x. The element 011 stands for the polynomial x + 1.

The coefficients add together just like the integers-modulo-2 add. In group theory terms, the coefficients are from \mathbb{Z}/2\mathbb{Z}. That means that 0 + 0 = 0, 0 + 1 = 1, 1 + 0 = 1, and 1 + 1 = 0. In computer terms, a + b is simply a XOR b. That means, to get the integer representation of the sum of two elements, we simply have to do the bitwise-XOR of their integer representations. Every processor that I’ve ever seen has built-in instructions to do this. Most computer languages support bitwise-XOR operations. In Lisp, the bitwise-XOR of a and b is (logxor a b).

Multiplication is somewhat trickier. Here, we have to multiply the polynomials together. For example, (x + 1) \cdot x = x^2 + x. But, if we did (x^2 + x) \cdot x, we’d end up with x^3 + x^2. We wanted everything to be polynomials of degree less than n and our n is 3. So, what do we do?

What we do, is we take some polynomial of degree n that cannot be written as the product of two polynomials of less than degree n. For our example, let’s use x^3 + x + 1 (which is 1011 in our binary scheme). Now, we need to take our products modulo this polynomial.

You may not have divided polynomials by other polynomials before. It’s a perfectly possible thing to do. When we divide a positive integer a by another positive integer b (with b bigger than 1), we get some answer strictly smaller than a with a remainder strictly smaller than b. When we divide a polynomial of degree m by a polynomial of degree n (with n greater than zero), we get an answer that is a polynomial of degree strictly less than m and a remainder that is a polynomial of degree strictly less than n.

Dividing proceeds much as long division of integers does. For example, if we take the polynomial (with integer coefficients) x^3 + 2x + 5 and divide it by the polynomial x + 3, we start by writing it as:

x + 3 \overline{\left) x^3 + 0x^2 + 2x + 5\right.}

We notice that to get (x+3) \cdot q(x) to start with x^3, we need q(x) to start with x^2. We then proceed to subtract (x+3) \cdot x^2 and then figure out that we need a -3x for the next term, and so on. We end up with x^2 - 3x + 11 with a remainder of -28 (a degree zero polynomial).

\begin{array}{ccrcrcrcr} & & & & x^2 & - & 3x & + & 11 \\ & \multicolumn{8}{l}{\hrulefill} \\(x+3) & ) & x^3 & + & 0x^2 & + & 2x & + & 5 \\ & & x^3 & + & 3x^2 & & & & \\ & & \multicolumn{3}{l}{\hrulefill} & & & & \\ & &     & - & 3x^2 & + & 2x & & \\ & &     & - & 3x^2 & - & 9x & & \\ & &     & \multicolumn{4}{l}{\hrulefill} & & \\ & &     &   &      &   & 11x & + & 5 \\ & &     &   &      &   & 11x & + & 33 \\ & &     &   &      &   & \multicolumn{3}{l}{\hrulefill} \\ & &     &   &      &   &     & - & 28\end{array}

For the example we cited earlier we had x^3 + x^2 which we needed to take modulo x^3 + x + 1. Well, dividing x^3 + x^2 by x^3 + x + 1, we see that it goes in one time with a remainder of x^2 + x + 1. [Note: addition and subtraction are the same in GF(2^n).]

\begin{array}{ccrcrcrcr} & & & & & & & & 1 \\ & \multicolumn{8}{l}{\hrulefill} \\(x^3+x+1) & ) & x^3 & + & x^2 & + & 0x & + & 0 \\ &            & x^3 & + & 0x^2 & + & x & + & 1 \\ & & \multicolumn{7}{l}{\hrulefill} \\ & &     &  & x^2 & + & x & + & 1\\\end{array}

For a reasonable way to accomplish this in the special case of our integer representations of polynomials in GF(2^n), see this article about Finite Field Arithmetic and Reed Solomon Codes. In (tail-call style) Lisp, that algorithm for GF(2^8) might look something like this to multiply a times b modulo m:

(flet ((next-a (a)
               (ash a -1))
       (next-b (b)
               (let ((overflow (plusp (logand b #x80))))
                 (if overflow
                     (mod (logxor (ash b 1) m) #x100)
                   (ash b 1)))))
  (labels ((mult (a b r)
             (cond
              ((zerop a) r)
              ((oddp a) (mult (next-a a) (next-b b) (logxor r b)))
              (t (mult (next-a a) (next-b b) r)))))
    (mult a b 0)))

How is the Galois field structured?

The additive structure is simple. Using our 8-bit representations of elements of GF(2^8), we can create an image where the pixel in the i-th row and j-th column is the sum (in the Galois field) of i and j (written as binary numbers). That looks like this:

Just before the above-mentioned article hit reddit, I got to wondering if the structure of the Galois field was affected at all by the choice of polynomial you used as the modulus. So, I put together some code to try out all of the polynomials of order 8.

Remember way back at the beginning of multiplication, I said that the modulus polynomial had to be one which couldn’t be written as the product of two polynomials of smaller degree? If you allow that, then you have two non-zero polynomials that when multiplied together will equal your modulus polynomial. Just like with integers, if you’ve got an exact multiple of the modulus, the remainder is zero. We don’t want to be able to multiply together two non-zero elements to get zero. Such elements would be called zero divisors.

Zero divisors would make these just be Galois rings instead of Galois fields. Another way of saying this is that in a field, the non-zero elements form a group under multiplication. If they don’t, but multiplication is still associative and distributes over addition, we call it a ring instead of a field.

Galois rings might be interesting in their own right, but they’re not good for AES-type encryption. In AES-type encryption, we’re trying to mix around the bits of the message. If our mixing takes us to zero, we can never undo that mixing—there is nothing we can multiply or divide by to get back what we mixed in.

So, we need a polynomial that cannot be factored into two smaller polynomials. Such a polynomial is said to be irreducible. We can just start building an image for the multiplication for a given modulus and bail out if it has two non-zero elements that multiply together to get zero. So, I did this for all elements which when written in our binary notation form odd numbers between (and including) 100000001 and 111111111 (shown as binary). These are the only numbers which could possibly represent irreducible polynomials of degree 8. The even numbers are easy to throw out because they can be written as x times a degree 7 polynomial.

The ones that worked were: 100011011, 100011101, 100101011, 100101101, 100111001, 100111111, 101001101, 101011111, 101100011, 101100101, 101101001, 101110001, 101110111, 101111011, 110000111, 110001011, 110001101, 110011111, 110100011, 110101001, 110110001, 110111101, 111000011, 111001111, 111010111, 111011101, 111100111, 111110011, 111110101, and 111111001. That first one that worked (100011011) is the one used in AES. Its multiplication table looks like:

Here’s it is again on the left with the multiplication table when 110011111 is the modulus on the right:

 

So, the addition image provides some insight into how addition works. The multiplication tables, at least for me, provide very little insight into anything. They don’t even make a good stereo pair.

To say two of the multiplication tables have the same structure, means there is some way to map back and forth between them so that the multiplication still works. If we have table X and table Y, then we need an invertible function f such that f(a \cdot_X b) = f(a) \cdot_Y f(b) for all a and b in the table X.

What’s next?

If there is an invertible map between two multiplication tables and there is some element a in the first table, you can take successive powers of it: a, a^2, a^3, \ldots. There are only 2^n elements in GF(2^n) no matter which polynomial we pick. So, somewhere in there, you have to start repeating. In fact, you have to get back to a. There is some smallest, positive integer k so that a^{k+1} \equiv a in GF(2^n). If we pick a = 0, then we simply have that k = 1. For non-zero a, we are even better off because a^k \equiv 1.

So, what if I took the powers of each element of GF(2^n) in turn? For each number, I would get a sequence of its powers. If I throw away the order of that sequence and just consider it a set, then I would end up with a subset of GF(2^n) for each a in GF(2^n). How many different subsets will I end up with? Will there be a different subset for each a?

I mentioned earlier that the non-zero elements of GF(2^n) form what’s called a group. The subset created by the powers of any fixed a forms what’s called a subgroup. A subgroup is a subset of a group such that given any two members of that subset, their product (in the whole group) is also a member of the subset. As it turns out, for groups with a finite number of elements, the number of items in a subgroup has to divide evenly into the number of elements in the whole group.

The element zero in GF(2^n) forms the subset containing only zero. The non-zero elements of GF(2^n) form a group of 2^n - 1 elements. The number 2^n - 1 is odd (for all n \ge 1). So, immediately, we know that all of the subsets we generate are going to have an odd number of items in them. For GF(2^8), there are 255 non-zero elements. The numbers that divide 255 are: 1, 3, 5, 15, 17, 51, 85, and 255.

It turns out that the non-zero elements of GF(2^n) form what’s called a cyclic group. That means that there is at least one element a whose subset is all 2^n - 1 of the non-zero elements. If take one of those a‘s in GF(2^8) whose subset is all 255 of the elements, we can quickly see that the powers of a^3 form a subset containing 85 elements, the powers of a^5 form a subset containing 51 elements, …, the powers of a^{85} form a subset containing 3 elements, and the powers of a^{255} form a subset containing 1 element. Further, if both a and b have all 255 elements in their subset, then a^k and b^k will have the same number of elements in their subsets for all k. We would still have to check to make sure that if a^i + a^j = a^k that b^i + b^j = b^k to verify the whole field structure is the same.

This means there are only 8 different subset of GF(2^8)‘s non-zero elements which form subgroups. Pictorially, if we placed the powers of a around a circle so that a^0 = 1 was at the top and the powers progressed around the circle and then drew a polygon connecting consecutive points, then consecutive points in the a^3 sequence and consecutive points in the a^5 sequence, etc…. we’d end up with this:

If we don’t reorder them and just leave them in the numerical order of their binary representation, the result isn’t as aesthetically pleasing as the previous picture. Here are the same two we used before 100011011 (on the left) and 110011111 (on the right). They are easier to look at. They do not lend much more insight nor make a better stereo pair.

 

*shrug* Here’s the source file that I used to generate all of these images with Vecto and ZPNG: group.lisp

If it quacks like a parabola… September 21st, 2011
Patrick Stein

I am working on a game idea that involves (special) relativistic mechanics instead of classical mechanics. Working out the details that I needed was easy enough if I assumed that:

  • ships had a maximum speed at which they could move relative to my base frame
  • ships could instantly go from stopped to maximum speed or vice-versa

I didn’t like those assumptions at all. So, I started playing with the equations for relativity. In classical mechanics, the rate-of-change of velocity equals the force you’re applying divided by your mass: \frac{dv}{dt} = \frac{F}{m}.

In special relativity, your mass increases with velocity. So, that equation becomes: \frac{d\left(\frac{v}{\sqrt{1-v^2}}\right)}{dt} = \frac{F}{m_0} (assuming units where the speed of light is 1 unit of distance per 1 unit of time and m_0 is your rest-mass).

For the purposes of this post, I’m going to assume the simplest initial conditions: you start motionless and at the origin. For ease of notation, let a = \frac{F}{m_0}. Solving the above differential equation to get a formula for velocity and solving the resulting differential equation to get the distance x you’ve travelled in my base frame by time t, the answer comes out to: x(t) = \frac{1}{a}\left(-1 + \sqrt{1+a^2t^2}\right).

I have solved this problem at least thirty times in the past two months. Sometimes I used the simple initial conditions as above. Sometimes I did it in all of its gory details (including the messy case where the applied force is not aligned with the current velocity).

I got the same answer (well, plus the extra mess when I did the full-on problem) every way that I tried it.

So, why did I do it over and over again?

If this were classical mechanics, the end equation would have been x(t) = \frac{1}{2}at^2. And, I know that for low velocities, the classical mechanics answer should be almost identical to the special relativity answer. And, there was no way that I thought \frac{1}{a}\left(-1 + \sqrt{1+a^2t^2}\right) \approx \frac{1}{2}at^2.

I knew what the graph x = \sqrt{a^2t^2} looked like when t \ge 0. It is a straight line. It doesn’t look much like the parabola x = \frac{1}{2}a t^2 at all.

My assumption was that since x = \sqrt{a^2t^2} was a straight line for t \ge 0, then x = \sqrt{1 + a^2t^2} would be a straight line shifted up one unit and bent (concave-down) a little bit like the graph of x = \sqrt{at} is bent.

Boy was I wrong. Here is a plot of the two together (created with fooplot). The red line is the classical mechanics version. The black line is the relativistic version. Here, the force is such that the body is accelerating at a rate of the speed of light per second so they’ve already gotten up to around 28,000 miles per second before you can see any separation in the graphs here.

distance (in light-seconds) vs. time (in seconds)

Definitely, I can see the resemblance. Now, to fix my intuition about square-roots.

Calculating the mean and variance with one pass February 15th, 2011
Patrick Stein

A friend showed me this about 15 years ago. I use it every time I need to calculate the variance of some data set. I always forget the exact details and have to derive it again. But, it’s easy enough to derive that it’s never a problem.

I had to derive it again on Friday and thought, I should make sure more people get this tool into their utility belts.

First, a quick refresher on what we’re talking about here. The mean \mu of a data set { x_1, x_2, \ldots, x_n } is defined to be \frac{1}{n} \sum_{i=1}^n x_i. The variance \sigma^2 is defined to be \frac{1}{n} \sum_{i=1}^n (x_i - \mu)^2.

A naïve approach to calculating the variance then goes something like this:

(defun mean-variance (data)
  (flet ((square (x) (* x x)))
    (let* ((n (length data))
           (sum (reduce #'+ data :initial-value 0))
           (mu (/ sum n))
           (vv (reduce #'(lambda (accum xi)
                           (+ accum (square (- xi mu))))
                       data :initial-value 0)))
      (values mu (/ vv n)))))

This code runs through the data list once to count the items, once to calculate the mean, and once to calculate the variance. It is easy to see how we could count the items at the same time we are summing them. It is not as obvious how we can calculate the sum of squared terms involving the mean until we’ve calculated the mean.

If we expand the squared term and pull the constant \mu outside of the summations it ends up in, we find that:

\frac{\sum (x_i - \mu)^2}{n} = \frac{\sum x_i^2}{n} - 2 \mu \frac{\sum x_i}{n} + \mu^2 \frac{\sum 1}{n}

When we recognize that \frac{\sum x_i}{n} = \mu and \sum_{i=1}^n 1 = n, we get:

\sigma^2 = \frac{\sum x_i^2}{n} - \mu^2 = \frac{\sum x_i^2}{n} - \left( \frac{\sum x_i}{n} \right)^2
.

This leads to the following code:

(defun mean-variance (data)
  (flet ((square (x) (* x x)))
    (destructuring-bind (n xs x2s)
        (reduce #'(lambda (accum xi)
                    (list (1+ (first accum))
                          (+ (second accum) xi)
                          (+ (third accum) (square xi))))
                data :initial-value '(0 0 0))
      (let ((mu (/ xs n)))
        (values mu (- (/ x2s n) (square mu)))))))

The code is not as simple, but you gain a great deal of flexibility. You can easily convert the above concept to continuously track the mean and variance as you iterate through an input stream. You do not have to keep data around to iterate through later. You can deal with things one sample at a time.

The same concept extends to higher-order moments, too.

Happy counting.

Edit: As many have pointed out, this isn’t the most numerically stable way to do this calculation. For my part, I was doing it with Lisp integers, so I’ve got all of the stability I could ever want. :) But, yes…. if you are intending to use these numbers for big-time decision making, you probably want to look up a really stable algorithm.

Making Fun Algebra Problems Funner October 29th, 2010
Patrick Stein

A month ago, a friend posted the following problem on Facebook. I just noticed it this week.
Drawing of a circle with radius r sitting on a line with two squares wedged against it.  One square has side-length 1, the other side-length 2.  There are six units between the squares.
The goal is to find the exact length of the radius r.

I love this kind of math problem.  It has a bunch of features that make it a great, toy math problem.

  • It looks like it’s going to be easy, but at the same time seems at a glance to not have enough information
  • It looks like a geometry problem but only requires that you know:
    • All radii of a circle have the same length
    • A radius to a point where a tangent line touches the circle is perpendicular to that tangent line
  • It requires only very basic algebra:
    • Pythagorean theorem
    • Solving quadratics
  • The numbers in the problem are small, non-zero integers

I spent the next 25 minutes and six pieces of paper working the problem. About 20% of the time that I spent was rechecking my work. Why did I bother rechecking my work on a toy problem?

Warning: Spoilers ahead. If you like this kind of problem, stop reading now and play with it first.


The problem with the problem

This problem failed to satisfy one of my other criterion for great, toy puzzle problems.

  • The answer is a small natural number (or the current year)

I am notorious for messing up signs when doing large algebra calculations. I had to check and recheck all of my work to make sure that I hadn’t done 5x - (3x - 2) and gotten 2x - 2 somewhere. If I had come up with an integer value for r, I’d have been done. Instead, the answer involved a half-integer plus a multiple of \sqrt{74}.

What? \sqrt{74}?!? Raise your hand if you’ve ever seen that written anywhere before this problem.

That’s what I thought.

So, I spent the next hour looking for a different set of numbers to put in for 1, 2, and 6 so that the radius r would come out to an integer. My approach was Brownian, at best. I threw Pythagorean triples at a wall until one stuck.

The sticky triple left me with [50,200,360] to put in place of [1,2,6]. So, I put the problem down for awhile.

The next time that I was in the car, I realized that the radius for the sticky triple was less than 200. That meant that the larger box was touching the upper half of the circle instead of the lower half. The circle had to overlap the box.

So, I was back to the drawing board. Of course, I’d have been back to the drawing board at some point anyway because that sticky triple violated both of my small numbers criteria.

Warning: If you want to play around with finding a set of numbers that work, do it before you read the following. There are even more spoilers ahead.

Infinitely many combinations to get an integer radius

When I next sat down to play (pronounced: /OBSESS/) with this problem, I quickly hit upon the following way to take any Pythagorean triple (a, b, c) and make a version of the above puzzle where the radius r is an integer.
Same sort of diagram as above except radius is c, boxes are height c-b and c-a, and some (a,b,c) triangles are drawn in

In fact, using a 3-4-5 triangle, it is obvious that had the distance between the blocks in the original problem been 7 instead of 6, the radius would have been a small natural number.

Now. Now, I had infinitely many ways to construct this problem so that the radius was an integer. But, did I have all of the ways?

All your configuration are belong to us

When looking at that last diagram, the first question that comes to mind (for me) is: Do I really need to use the same Pythagorean triple on the left and right triangles? The answer, of course, is No.

If I have two Pythagorean triples (a_1,b_1,c_1) and (a_2,b_2,c_2) and d is any (positive) divisor of both c_1 and c_2, then I can start with an a_1-b_1-c_1 triangle on the left and an a_2-b_2-c_2 triangle on the right. I can then scale up the left triangle by c_2/d and scale the right triangle up by c_1/d. Now, both triangles have the same hypotenuse.

This means that if the left block has height c_2 (c_1 - b_1)/d, the right block has height c_1(c_2 - b_2)/d, the distance between the two blocks is (c_2 a_1 + c_1 a_2)/d, and the resulting circle has radius c_1c_2/d.
A more general version of the previous diagrams

Does this cover all of the solutions? Could I possibly have non-integer a_1 and a_2 such that (c_2a_1 + c_1a_2)/d is still integer?

Well, for the problem to be formulated entirely with integers, we certainly need the hypotenuse of the right triangles to be an integer. Further, to make the blocks integer heights then both of the legs of the right triangles along the vertical radius must also be integer. So, if the triangle on the left had hypotenuse r, vertical leg b and horizontal leg a, then we know that a = \sqrt{r^2 - b^2} where both r and b are integers. Thus, a must be the square root of an integer.

It is easy to see that if a were a non-integer rational number, then its square would also be a non-integer rational number. So, either a is integer or it is the irrational square root of an integer.

So, can a be an irrational number? If a were an irrational number, then the horizontal leg of the other triangle would have to be some integer n minus the irrational a for the distance between the two blocks to be an integer. Let’s say the vertical leg of the triangle on the right is \beta. Just like b in the first triangle, \beta must be an integer. We also know that (n-a)^2 = r^2 - \beta^2. This, in turn, means that a = (\beta^2 + n^2 + a^2 - r^2) / (2n). The right hand side is clearly rational, so a could not have been irrational.

So, there you have it: all of the numbers that make this problem have all integers on the drawing and an integer answer.

Updates In Email

Email:

l