Extension fields

[New version HaskellForMaths 0.1.7 uploaded here.]

Last week, we looked at the finite fields of prime order, Fp. A field, remember, is a set in which you can do addition, subtraction, multiplication and division. For a given prime p, the field Fp consists of the set {0, 1, ..., p-1}, with arithmetic operations done modulo p. Why does p have to be a prime? We can do addition, subtraction, and multiplication modulo n for any n, prime or not. However, for division, we require p to be prime. (Why?)

You might think, therefore, that the fields Fp of prime order are the only finite fields. However, you would be wrong. We shall see that there are in fact fields Fq of order q (ie, with q elements) for every prime power q = p^n. However, before we can do that we need to understand about algebraic extensions of a field.

It can happen that a field is contained within another field. For example, we have Q

• At step 0, start with Q
• At step 1, add a
• At step n, add all elements that can be obtained by combining two elements from step n-1 using one of the arithmetic operators.

So for example, at step 2, we will have 1+a, 2*a, a*a, 1/a, etc. Anything which appears in the list at step n for some n is in Q(a).

The field Q(a), obtained by adjoining a single element, is called a simple extension of Q. (We could then go on to adjoin another element b, to get Q(a,b), and so on.) Among simple extensions, there are two ways it can go.

Suppose that a is a zero of some polynomial with coefficients in Q - for example x^2-2. In that case, we will construct the whole of Q(a) after only a finite number of steps. For example, in the case where a = sqrt 2, every element of Q(sqrt 2) can be expressed as c + d sqrt 2, for some c, d in Q. For example:
(sqrt 2)^2 = 2
1/(sqrt 2) = (sqrt 2) / 2
1/(1 + sqrt 2) = (1 - sqrt 2) / (1 - 2) = -1 + sqrt 2
This is called an algebraic extension.

Alternatively, if a is not the zero of a polynomial over Q (for example: e, pi), then the construction of Q(a) in steps will never finish. This is called a transcendental extension.

Algebraic extensions of Q are fascinating things. Initially, the main motivation for studying them came from number theory. Some other time, I'd like to take a closer look at them. However, for the moment, I just want to show how to construct them in Haskell, as what we're really aiming for is algebraic extensions of the finite fields Fp.

Let's suppose that we're trying to construct Q(sqrt2). The polynomial we're interested in is x^2-2. Then the basic idea is:

1. Form the polynomial ring Q[x], consisting of polynomials in x with coefficients in Q
2. Represent elements of Q(a) as polynomials in Q[x]
3. Do addition, subtraction, multiplication in Q(a) using the underlying operations in Q[x]
4. After any arithmetic operation, replace the result by its remainder on division by x^2-2.

Step 4 is the key step. Suppose that we end up with a polynomial f. We use division with remainder to write f = q(x^2-2)+r, and we replace f by r. In effect, this means that we're setting x^2-2 = 0, which is the same as setting x = sqrt 2. If we do this consistently, then the x ends up acting as sqrt2, and we end up with Q(sqrt2).

To construct Q(sqrt3), Q(i), etc, just replace x^2-2 by the appropriate polynomial.

Okay, rather belatedly, time for some code. First, we need to define a type for (univariate) polynomials:

newtype UPoly a = UP [a] deriving (Eq,Ord)

UP [c0,c1,...cn] is to be interpreted as the polynomial c0 + c1 x + ... + cn x^n.

Exercise: Define a Num instance for UPoly a.

If we then define

x = UP [0,1] :: UPoly Integer

together with a suitable Show instance, then we can do things like the following:

> (1+x)^3
1+3x+3x^2+x^3


Exercise: Write quotRemUP, to perform division with remainder in UPoly k, on the assumption that k is a field (that is, a Fractional instance).

So we now have a type, UPoly Q, representing Q[x]. Next, we want to wrap this in another type to represent extension fields Q(a). Rather than have to do this over again each time for Q(sqrt2), Q(sqrt3), Q(i), and so on, we use a little bit of phantom type trickery again. Last time, we used phantom types to represent integers; this time, we're going to use phantom types to represent the polynomials that define the fields (x^2-2, x^2-3, x^2+1, etc).

class PolynomialAsType k poly where
    pvalue :: (k,poly) -> UPoly k

data ExtensionField k poly = Ext (UPoly k) deriving (Eq,Ord)

Here, k represents the field we are extending - Q, to begin with - and if a is the element that we want to adjoin, then poly represents the polynomial over Q of which a is a zero.

From here, it's a short step to define some extension fields:

data Sqrt a = Sqrt a

-- n should be square-free
instance IntegerAsType n => PolynomialAsType Q (Sqrt n) where
    pvalue _ = convert $ x^2 - fromInteger (value (undefined :: n))

type QSqrt2 = ExtensionField Q (Sqrt T2)
sqrt2 = embed x :: QSqrt2

type QSqrt3 = ExtensionField Q (Sqrt T3)
sqrt3 = embed x :: QSqrt3

And now we can do arithmetic in Q(sqrt2), Q(sqrt3), etc. For example:

> :set +t
> (1+sqrt2)^2
3+2a :: QSqrt2
> (1+sqrt3)^2
4+2a :: QSqrt3

As you see, the show function of ExtensionField k poly shows the adjoined element as "a", regardless of which field we're working in. With a little more type hackery, we could have made that a parameter to the type too, so that it would show "sqrt2", "sqrt3", "i", depending which field we were in. However, that would have obscured the code even more. In practice, when we use this code we're only going to be working over one field at a time, so it's fine as it is.

Another limitation of this code is that it's going to be a bit unwieldy to construct Q(a,b) as the type ExtensionField (ExtensionField k poly1) poly2. Luckily, this doesn't matter, as any algebraic extension Q(a,b) is equal to Q(c) for some c.

Exercise: Show that Q(sqrt2, sqrt3) = Q(sqrt2 + sqrt3), and find the polynomial of which sqrt2 + sqrt3 is a zero.

I apologise that I've only really explained the bare bones of how this works. Hopefully you can fill in the gaps. (They're all in the HaskellForMaths source - see link at top of page.) This post has already taken me far too long to write, but there is just one other thing I ought to mention.

Field extensions can have automorphisms (symmetries). Recall that a symmetry is a change that leaves something looking the same. In the case of Q(sqrt2), there is a non-trivial symmetry that sends c + d sqrt 2 to c - d sqrt2 (conjugation). Field automorphisms have a very important place in the history of mathematics: Galois invented group theory in order to study these symmetries, during his investigations into the unsolvability (in general) of the quintic.

Finite fields, part 1

Okay, so you could think of this as the beginning of chapter two of the imaginary book that I'm writing in this blog. Chapter one was about graphs, groups and symmetry. In chapter two I want to talk about some other highly symmetric combinatorial objects: finite geometries and (probably) designs. But first, I need to talk about finite fields.

A field is roughly a set where you can add, subtract, multiply and divide. You're probably familiar with the following fields:

• Q, the rational numbers
• R, the real numbers
• C, the complex numbers

And also with the following non-fields:

• N, the natural numbers. Not a field, because it doesn't contain additive inverses (the negative numbers), so you can't always subtract.
  
• Z, the integers. Not a field, because it doesn't contain multiplicative inverses, so you can't always divide.

I can't resist mentioning a little curiosity. To a number theorist / discrete mathematician like myself, the most interesting of the above are Q and Z. I like to think that this is why they were given the most interesting letters. Indeed, the interestingness of the mathematical objects appears to be correlated with the Scrabble value of the letters: in Scrabble, Q and Z are worth 10, indicating highly interesting; R and N are worth 1, indicating a bit boring; C is worth 3 - moderately interesting. Going further, H, the quaternions, is worth 4 - slightly more interesting than C. However, O, the octonions, is only worth 1, which looks like an anomaly.

Anyway, back to fields. Before we go any further, here's the HaskellForMaths version of Q, the rationals:

import Data.Ratio

newtype Q = Q Rational deriving (Eq,Ord,Num,Fractional)

instance Show Q where
   show (Q x) | b == 1    = show a
              | otherwise = show a ++ "/" ++ show b
              where a = numerator x
                    b = denominator x

Silly, I know, but I just got a bit annoyed with those percentage signs all over the place when using Haskell's built-in rationals.

Anyway, in addition to the fields Q, R, and C, there are many other types of field. The ones I want to look at here are the finite fields of prime order. These are the fields you get if you do arithmetic modulo p, where p is a prime. For each different prime p, we have a field Fp, consisting of the set {0,1, ... p-1}, with arithmetic performed modulo p. For example, in the field F5, we have 3*3 = 4; while in the field F7, we have 3*3 = 2.

Okay, so how to represent these fields in Haskell?

Well, in an earlier version of the HaskellForMaths code, I tried this:

data Fp = F Integer Integer

The idea is that the first integer is p, the prime we are working over, and the second is the value. For example, F 5 3 would represent the value 3 in the field F5.

Next we either derive or define instances of Eq, Ord, Show, Num, Fractional. The Num instance looks something like this:

instance Num Fp where
    F p x + F q y | p == q  = F p $ (x+y) `mod` p
    F p x * F q y | p == q  = F p $ (x*y) `mod` p
    fromInteger _ = error "Fp.fromInteger: not well defined"

The drawbacks of this approach are probably pretty obvious:

• Whenever I have two Fp values, F p x and F q y, I first have to check that p = q. This is a run-time check of something that would be better handled at compile time by the type system
• When defining numerical types in Haskell, the fromInteger function gives us implicit type conversion from integers to our type, so that we can just write code like 3 * 3, and have the type system convert them into the correct type. However, in this case, we can't give a sensible definition of fromInteger, because we don't know which prime p the user wants to work over. This means that everything is going to have to be written longhand as F 5 3 * F 5 3, and so on.
• It feels like a waste of memory to be carrying that first Integer p around all the time, when if I'm using the code correctly, all the ps in a given calculation will be the same.
  

For these reasons, I decided to try another approach. The aim is to make F2, F3, F5, etc be different types, so that adding 3 :: F5 to 3 :: F7 is a type error, which can be spotted at compile time, rather than a runtime error.

It's fairly simple really.

First, define a whole bunch of phantom types representing the primes:

data T2
data T3
data T5
...

Now define a mechanism for retrieving the values from the types:

class IntegerAsType a where
   value :: a -> Integer

instance IntegerAsType T2 where value _ = 2
instance IntegerAsType T3 where value _ = 3
instance IntegerAsType T5 where value _ = 5
...

This enables us to do the following:

> value (undefined :: T2)
2
> value (undefined :: T3)
3

Now, we use the phantom types to parameterise a finite field data type:

newtype Fp n = Fp Integer

Here, Fp (the one on the left of the equals sign) is a type constructor which takes a type parameter n. This type parameter will be one of the phantom types we just defined. For example:

type F2 = Fp T2
type F3 = Fp T3
type F5 = Fp T5

As before, we need to derive or define Eq, Ord, Show, Num, and Fractional instances. The Num instance looks like this:

instance IntegerAsType n => Num (Fp n) where
    Fp x + Fp y = Fp $ (x+y) `mod` p where p = value (undefined :: n)
    Fp x * Fp y = Fp $ (x*y) `mod` p where p = value (undefined :: n)
    fromInteger m = Fp $ m `mod` p where p = value (undefined :: n)
    ...

This is all fairly standard phantom type trickery, but it now enables us to do things like the following:

> :load Math.Algebra.Field.Base
> 3*3 :: F5
4
> 3*3 :: F7
2
> (3 :: F5) * (3 :: F7)
#type error#

We still have to tell Haskell which prime p we're working over, but now we do it via a type annotation, rather than an argument to the constructor. This means that Haskell can catch errors where we mix values from different fields at compile time, rather than run time.

The module Math.Algebra.Field.Base provides types for the finite fields F2, F3, F5, ..., F97. The module also provides functions to list the elements of each field:

f2 = map fromInteger [0..1] :: [F2]
f3 = map fromInteger [0..2] :: [F3]
f5 = map fromInteger [0..4] :: [F5]
...

And that's about it.

However, we must recognise, in all honesty, that there are also drawbacks to this implementation of the fields Fp.

First, whenever we want to work over a new prime p, we have to define a new type. For example, if we want to work over F101, we are going to have to define a new type for it.

(For the avoidance of doubt: No, type level Peano arithmetic is not what I'm after.)

Second, I can't decide at run time what field to work over. For certain applications it would be very useful to choose my p at run time. I could do that with the first implementation, but I can't do it with the second.

What I really want then, is for Fp and Fq to be different types when p /= q, but to be able to decide at run-time which p to work over. For example, I might have calculated the p I want to work over as the result of a previous calculation. This is impossible in Haskell.

However, there is a technology called "dependent types", in which what I want is possible. Dependent types means roughly that types can depend upon values. So in a Haskell-like language that supported dependent types, I might be able to write code like the following (not valid Haskell):

> 3 * 3 :: F 5
4

Here, the value 5 (not the type T5) is being passed as a type parameter to the type constructor F.

For a while, I thought this was the holy grail. I even considered defecting from Haskell. I did look into languages that support dependent types, such as Coq. However, I concluded that, sadly, this is not the answer. The problem is that - I think - if you have types dependent on values, then you lose the distinction between compile time and run time, which means - I think - that you can't compile the language any more.

Or does someone else know better?


Where we've been, and where we're going

It feels like about time to take a step back and survey where we've been and where we're going in this blog.

What I'm trying to do in this blog is a kind of guided tour of my HaskellForMaths library. In turn, the library itself, although it contains efficient implementations of fundamental algorithms in computer algebra, is really an educational tool. I wrote it to help me understand maths I wanted to learn about, and especially, to help develop my intuitions about the objects of study. My hope is that by retracing my steps in this blog, I can help others to do the same.

Up to now, we've been looking at graphs and their symmetries, considered as permutation groups. Group theory is sometimes described as the study of symmetry, and I hope that by introducing groups through symmetries of graphs, I've made it clear why.

We've really only looked at quite a small amount of code. The main functions we've seen (forgetting a few stepping stones along the way) are:

• graph (vs,es) - for constructing a graph with given vertices and edges
• graphAuts g - returns a strong generating set for the group of symmetries / automorphisms of the graph
• orderSGS sgs - given a strong generating set, return the order of the group - for example, the total number of symmetries of the graph
• conjClassReps - given a generating set for a group, return representatives and sizes for the conjugacy classes - the classes of elements which are "the same, just viewed from a different angle"
• p [[1,2,3],[4,5]] - given a list of cycles, return the corresponding permutation
• g*h, 1, g^-1 - multiplication, identity and inverses for group elements
  
• v .^ g - the action of a permutation on a vertex
• e -^ g - the action of a permutation on an edge
• sgs gs - return a strong generating set, given a set of generators

This is the core of the HaskellForMaths code for graphs, graph automorphisms, and permutation groups. Incidentally, I think this code really shows off Haskell to advantage: code to do the same things would be much more complicated in C++, say. I'm especially pleased with the ability to define graphs and groups over arbitrary types, which enables very natural constructions of various classes of graphs and groups. (In this respect, though not in others, I think the library probably beats even dedicated packages such as GAP and MAGMA.)

However, the graph and permutation group code represents, I would say, only about one sixth of what's in the HaskellForMaths library. (Not to mention that I have quite a lot more code in various stages of development back in the lab.) So what else is there still to look at?

Well, what I want to do next is stick with permutation groups a little longer, but widen the set of combinatorial objects that we consider, to include finite geometries, designs, and perhaps other incidence structures.

One of my aims is to take a look at a few of the sporadic finite simple groups. Finite simple groups are the "atoms of symmetry". Most of them fall in one of several infinite families (the cyclic groups of prime order, the alternating groups, and the finite groups of Lie type, aka Chevalley groups). However, it turns out that there are also 26 "sporadic" finite simple groups, which don't belong to any of these families. The HaskellForMaths library contains code for investigating both the infinite families, and some of the smaller sporadic groups.

Then there are root systems and Coxeter groups, including string rewriting and the Knuth-Bendix algorithm. Root systems have some connection to Lie algebras, and are therefore important in physics, though I don't claim to fully understand all that stuff.

Another major part of the library, which I'll get round to eventually, is commutative algebra and Groebner bases. Groebner bases provide a way of calculating with ideals in polynomial rings, which is equivalent to doing algebraic geometry.

Finally, there's some code for working in non-commutative algebra. As an example, there's some code which uses Temperley-Lieb and Hecke algebras to calculate polynomial invariants of knots (the Jones polynomial and the HOMFLY polynomial). In due course, I'm hoping to expand this part of the library to do group algebras, representation theory, and other stuff.

In case I've managed to interest anyone, I should probably mention a few of the books that I used for inspiration:

• Ronan, Symmetry and the Monster
• Godsil and Royle, Algebraic Graph Theory
• Cameron, Combinatorics
• Seress, Permutation Group Algorithms
  
• Holt, Handbook of Computational Group Theory

Next time - back to the maths.

How to count the number of positions of Rubik's cube

Previously on this blog, we have been looking at the symmetries of graphs, which we defined as permutations of the vertices which leave the edges (collectively) in the same places. We think of permutations as active rearrangements of the vertices, rather than as static sequences. For example, we think of the static permutation [2,1,4,3] as the active rearrangement [[1,2],[3,4]], which swaps the 1 and 2 positions, and swaps the 3 and 4 positions. Given two permutations g and h, the active viewpoint enables us to ask what happens when you do g then h, which we think of as multiplication g*h; and how to undo g, which we think of as an inverse g^-1.

The symmetries of a graph form a group. This means, that if g and h are symmetries, then so are g*h, and g^-1. So a group is a set of permutations which is closed under multiplication and inverses.

We can use permutation groups to study all sorts of things besides graph symmetries. As an example, this week I'm going to look at Rubik's cube.

First, we need to work out a way to represent cube moves as permutations. That's easy - we just number the little squares as follows:

Then a clockwise rotation of the front face (the blue face, with an F in the centre), can be represented as the following permutation:

f = [[1,3,9,7],[2,6,8,4],[17,41,33,29],[18,44,32,26],[19,47,31,23]]

This says, that f sends the 1 square to the 3 position, 3 to 9, 9 to 7, 7 to 1 (the corners of the blue face); 2 to 6, 6 to 8, 8 to 4, 4 to 2 (the edges of the blue face); together with three other cycles (the corners and edges adjacent to the blue face).

It's then an easy matter to write down permutations corresponding to the clockwise rotations of the other faces:

b = p [[51,53,59,57],[52,56,58,54],[11,27,39,43],[12,24,38,46],[13,21,37,49]]
l = p [[21,23,29,27],[22,26,28,24],[ 1,31,59,11],[ 4,34,56,14],[ 7,37,53,17]]
r = p [[41,43,49,47],[42,46,48,44],[ 3,13,57,33],[ 6,16,54,36],[ 9,19,51,39]]
u = p [[11,13,19,17],[12,16,18,14],[ 1,21,51,41],[ 2,22,52,42],[ 3,23,53,43]]
d = p [[31,33,39,37],[32,36,38,34],[ 7,47,57,27],[ 8,48,58,28],[ 9,49,59,29]]

Note that I didn't bother to number the centre square in each face, or write down any permutations that move a middle slice. For example, we could have written down:

x = [[4,44,54,24],[5,45,55,25],[6,46,56,26]]

However, doing x is the same as doing u^-1 * d^-1 and then turning the whole cube round. What we're interested in is the possible positions of the 54 squares relative to one another, so we don't really want to count rotations of the whole cube. For this reason, we leave out rotations of the middle slices.

Okay, so we have written down six permutations, corresponding to clockwise rotations of the six faces. Hopefully it's clear that these six elements generate the Rubik's cube group. So we are now in a position to perform calculations within the Rubik's cube group.

Within a group, the order of an element g is defined as the least n>0 such that g^n = 1. In other words, the order of g is the number of times you have to do g to get back where you started. (Recall that the order of a group is the number of elements it contains. The order of an element is equal to the order of the group generated by the element, since this is {g, g^2, ... , g^n-1, g^n = 1}, having n elements.)

For permutations written in cycle notation, it is easy to see the order by inspection. For example:

• The order of g = [[1,2]] is 2 - you have to do g twice to get back where you started.
• The order of h = [[1,2,3]] is 3 - you have to do h three times to get back where you started.
• The order of k = [[1,2,3],[4,5]] is 6. At k^2, we have [[1,3,2]] - the 4 and 5 are in the correct positions, but the 1, 2, 3 are not. At k^3 we have [[4,5]] - now the 1, 2, 3 are in the correct positions, but the 4 and 5 are not. Only at k^6 does everything get back to the starting position.

Hopefully it's clear that the order of an element will always be the least common multiple of the lengths of the cycles:

orderElt g = foldl lcm 1 $ map length $ toCycles g


For example, in the Rubik's cube group we have:

> f*l
[[1,3,9,37,53,17,41,33,27,21,23,19,47,59,11],[2,6,8,34,56,14,4],[7,31,29],[18,44,32,28,24,22,26]]
> orderElt it
105

That is, if you do a front turn followed by a left turn, you will have to repeat 105 times before you get back to the starting position.

When we were looking at graph symmetries, we found that it was particularly useful to have our generators in the form of a strong generating set, since it was then straightforward, for example, to calculate the order of the group.

Rubik's cube isn't a graph. Luckily, the HaskellForMaths library has a function to calculate a strong generating set from any set of generators:

> mapM_ print $ sgs [f,b,l,r,u,d]
[[1,3,9,7],[2,6,8,4],[17,41,33,29],[18,44,32,26],[19,47,31,23]]
[[1,21,51,41],[2,22,52,42],[3,23,53,43],[11,13,19,17],[12,16,18,14]]
[[1,31,59,11],[4,34,56,14],[7,37,53,17],[21,23,29,27],[22,26,28,24]]
[[2,34,22,26,32,44,16,12,46,38,24],[3,31,47,13,49,37,21],[4,8,6,42,52,54,58,56,18,28,14],[7,9,43,39,27,11,41],[19,29,33,51,57,59,53]]
[[3,13,57,33],[6,16,54,36],[9,19,51,39],[41,43,49,47],[42,46,48,44]]
[[4,38,32,46,16,44,34,24,36,14,12],[6,28,56,48,22,52,26,58,8,54,42],[7,31,29],[9,59,57,51,53,33,37,49,13,21,47,27,39,43,11]]
[[6,52,56,54],[9,53,57,51],[11,39,43,47],[12,24,46,44],[13,33,21,49]]
[[7,47,57,27],[8,48,58,28],[9,49,59,29],[31,33,39,37],[32,36,38,34]]
[[8,56,16,52,54],[9,33,47],[12,46,32,24,42],[13,51,43]]
[[9,59],[12,38],[13,49],[24,46],[27,47],[33,37],[39,43],[51,57],[52,58],[54,56]]
[[11,27,39,43],[12,24,38,46],[13,21,37,49],[51,53,59,57],[52,56,58,54]]
[[11,37,39,43],[13,21,59,49],[16,56,54,58],[24,46,38,42],[27,57,51,53]]
[[12,38,46,42,24],[16,56,52,58,54],[27,37,59],[39,57,49]]
[[12,56,58,36],[16,54],[24,38,48,52],[27,37,59],[39,57,49],[42,46]]
[[13,39,43,57,51,49],[16,24,46,38,42,56,54,58],[27,37,59],[36,48]]
[[13,43,51],[27,37,59]]
[[14,28,38,48,54,16,56],[22,34,58,36,46,42,24],[27,59,37],[39,49,57]]
[[16,56,58],[24,38,42],[27,37,59],[39,57,49]]
[[24,36,56,48],[27,37,59],[38,54,58,46],[39,57,49]]
[[24,46,38],[27,37,59],[39,57,49],[54,58,56]]
[[27,39],[36,48],[37,49],[38,46,58,54],[57,59]]
[[27,59,37],[39,49,57]]
[[28,54,58],[34,46,38]]
[[36,38,54],[46,48,58]]
[[36,58,54],[38,46,48]]
[[38,58],[46,54]]

Okay, so strong generating sets aren't all that interesting in themselves - although the last few entries in this one are a little bit interesting, as they show that there are moves which just rotate a pair of corners, or a pair of edges.

Remember though, that SGS allow us to easily calculate the order of a group. In this case:

> orderSGS $ sgs [f,b,l,r,u,d]
43252003274489856000

That's about 4 * 10^19 different positions.

The sgs function, to calculate a strong generating set, uses the Schreier-Sims algorithm. I'll probably explain how it works some time, but not just yet.