Geometric Algebra for Free!
The first bit of functionality I'll sacrfice is some of the genericity of my code so let's start with fixing a base field. Everything I say below will be over the reals:
> type R = Float
Anyway, today's theme is about how writing your code generically allows you to combine classes together in unanticipated ways to give you useful stuff - more or less for free.
If you've done some mathematical coding you probably have a complex number class lying around somewhere. If it's written in a modern language it's probably generic allowing you to choose your base structure. For example in C++ you might be able to use complex<float> or complex<int>. If you don't have such a class you can have a Haskell one for free:
> data Complex a = C a a deriving (Eq,Show)
> instance Num a => Num (Complex a) where
> fromInteger a = C (fromInteger a) 0
> (C a b)+(C a' b') = C (a+a') (b+b')
> (C a b)-(C a' b') = C (a-a') (b-b')
> (C a b)*(C a' b') = C (a*a'-b*b') (a*b'+b*a')
The usual complex numbers are given by Complex R. But here's a question: does Complex (Complex R) mean anything. Complex a requires a to be a Num, but Complex R is a Num. So Complex (Complex R) compiles fine. But does it mean anything mathematically? It does, it corresponds to C⊗C. (Remember, I'm working over the reals, so by ⊗ I mean ⊗R.) If you're shaky on tensor products then I can give a one line summary for the case of working over the reals: A⊗B means opening up the object of type A and replacing every occurence of a real number in it with a copy of a B. For example, a complex number can be thought of as an object that contains two real numbers. So C⊗C is an object that contains two complex numbers and which multiplies in a certain way. So the type constructor Complex actually corresponds quite well with the functor _⊗C. It's even better, if ⊗ is thought of as the tensor product in the category of algebras, not just vector spaces, then the Complex type constructor computes a tensor product of algebras. This is kinda obvious but I've never seen anyone say this explicitly. (Clearly the same is true of a C++ template like complex<>.)
So let's look at some examples. Maybe you also have a matrix class lying around. If you don't, here's a 2 by 2 matrix class for you:
> data Matrix a = M a a a a deriving (Eq,Show)
> instance Num a => Num (Matrix a) where
> fromInteger a = M (fromInteger a) 0 0 (fromInteger a)
> (M a b c d)+(M a' b' c' d') = M (a+a') (b+b')
> (c+c') (d+d')
> (M a b c d)-(M a' b' c' d') = M (a-a') (b-b')
> (c-c') (d-d')
> (M a b c d)*(M a' b' c' d') = M (a*a'+b*c') (a*b'+b*d')
> (c*a'+d*c') (c*b'+d*d')
We can start putting our words into practice. Use the notation A(n) to mean the n by n matrices with entries in A. Then A(n)=R(n)⊗A. So the type constructor Matrix is the functor _(2). Another theorem is that R(m)⊗R(n)=R(mn). This means that Matrix (Matrix R) gives 4 by 4 matrices. Neat eh? We're already getting some classes for free. But we can get bigger stuff for free.
I'm also a fan of geometric algebra. But I'm not a zealot. Geometric algebra can do a lot of neat stuff but I'm not one of those people campaigning to get everyone to throw out their old vectors and tensors to replace them with geometric algebra.
Geometric algebra is all about Clifford algebras. The Clifford algebra, Cliff(n), is the real numbers extended by throwing in n elements ei, for i from 1 to n, and all of the things you can get my multiplying and adding these numbers. In addition we impose the rules that ei²=-1 and eiej=-ejei. You already know some Clifford algebras. Cliff(0) is just r. Cliff(1) is C because i serves perfectly well as e1. Cliff(2) is probably familiar to you as well, it's just the quaternions. Ah...but quaternions have three values whose square is -1, so surely it's Cliff(3). Well, no actually. e1=i, e2=j and k is reallly just ij, ie. e1e2, and so isn't something new. In general, Cliff(n) is 2n dimensional with a basis consisting of all of the possible products of the ei. Anyway, if you don't already have a Quaternion class, here's one for you:
> data Quaternion a = Q a a a a deriving (Eq,Show)
> instance Num a => Num (Quaternion a) where
> fromInteger a = Q (fromInteger a) 0 0 0
> (Q a b c d)+(Q a' b' c' d') = Q (a+a') (b+b')
> (c+c') (d+d')
> (Q a b c d)-(Q a' b' c' d') = Q (a-a') (b-b')
> (c-c') (d-d')
> (Q a b c d)*(Q a' b' c' d') = Q (a*a'-b*b'-c*c'-d*d')
There are also another bunch of Clifford algebras called Cliff'(n). Cliff'(n), is the real numbers extended by throwing in n elements ei, for i from 1 to n. With the properties ei²=1 and eiej=-ejei. Note the 1 instead of -1. Cliff'(1) is known as the split complex numbers. You might not have an implementation of these, so here you go, gratis:
> data Split a = S a a deriving (Eq,Show)
> instance Num a => Num (Split a) where
> fromInteger a = S (fromInteger a) (fromInteger a)
> (S a b)+(S a' b') = S (a+a') (b+b')
> (S a b)-(S a' b') = S (a-a') (b-b')
> (S a b)*(S a' b') = S (a*a') (b*b')
You already have an implementation of Cliff'(2). Believe it or not, it's isomorphic to the 2 by 2 matrices. e1=[[-1,0],[0,1]] and e2=[[0,1],[1,0]], using the obvious Haskell lists to represent matrices. (Try squaring those matrices and seeing what you get.)
That's all well and good, but what about Cliff(n) for larger n. The obvious thing is to represent the elements as arrays of length 2n. But to multiply these things we'd have to loop over each element of the first multiplicand and then over each element of the second meaning a loop that takes time 22n. But who wants to write code when you can get stuff for free? Check outh this paper on the classification of Clifford algebras. On page 5 there's an important theorem
Now we know how to implement tensor products this becomes trivial. In Haskell notation,
> type Cliff1 = Complex R
> type Cliff1' = Split R
> type Cliff2 = Quaternion R
> type Cliff2' = Matrix R
> type Cliff3 = Quaternion Cliff1'
> type Cliff3' = Matrix Cliff1
> type Cliff4 = Quaternion Cliff2'
> type Cliff4' = Matrix Cliff2
and so on.
That's it. We've built Cliff(n) and Cliff'(n) out of a bunch of simple classes that you might already have had lying around. Geometric algebra for free! It almost seems sad to see a beautiful classification theorem demystified in such a simple way. Still, you don't see pairs of mutually recursive towers of types like this every day.
But actually this is all a con. It's no good having something isomorphic to what we want without the actual isomorphism. For example where is e3 in Cliff4', say? So we now have to do a tiny bit of real work. Start with a type class:
> class Clifford b where
> i :: R -> b
> e :: Int -> R -> b
i is the function that embeds our base field, the reals. And e i a is intended to be aei.
And now we simply assert:
> instance Clifford R where
> i a = a
> e _ _ = error ""
> instance Clifford (Complex R) where
> i a = C a 0
> e 1 a = C 0 a
> instance Clifford (Split R) where
> i a = S a a
> e 1 a = S a (-a)
Now we can define i and e for all of the algebras. For the details look at the proof of Theorem 3.4. They translate almost directly as:
> instance (Num b,Clifford b) => Clifford (Quaternion b) where
> i a = Q (i a) 0 0 0
> e 1 a = Q 0 (i a) 0 0
> e 2 a = Q 0 0 (i a) 0
> e n a = Q 0 0 0 (e (n-2) a)
> instance (Num b,Clifford b) => Clifford (Matrix b) where
> i a = let b = i a in M b 0 0 b
> e 1 a = let b = i a in M (-b) 0 0 b
> e 2 a = let b = i a in M 0 b b 0
> e n a = let b = e (n-2) a in M 0 b (-b) 0
And now we really do have usable Clifford algebras. Here are some more defintions:
> type Cliff5 = Quaternion Cliff3'
> type Cliff5' = Matrix Cliff3
> type Cliff6 = Quaternion Cliff4'
> type Cliff6' = Matrix Cliff4
> type Cliff7 = Quaternion Cliff5'
> type Cliff7' = Matrix Cliff5
> type Cliff8 = Quaternion Cliff6'
> type Cliff8' = Matrix Cliff6
> type Cliff9 = Quaternion Cliff7'
> type Cliff9' = Matrix Cliff7
(In C++ you can define 'metafunctions' that map integers to classes. Doing something similar in Haskell is pretty ugly.)
And I've learnt my lesson from my monad tutorial. Here's a little test to make sure I'm actually telling the truth:
> t i = e i 1 :: Cliff9
> t' i = e i 1 :: Cliff9'
> test a t n = map (\(i,j) -> (t i)*(t j)+(t j)*(t i)==0) [(i,j) | i <- [1..n], j <- [1..(i-1)]] ++
> map (\i -> (t i)*(t i)==a) [1..n]
> fulltest = and $ (test (-1) t 9) ++ (test 1 t' 9)
If fulltest evaluates to True then we know that for Cliff9 and Cliff9' we have created Haskell objects that satisfy the defining properties for ei.
One last thing, let's consider the complexity of this code. The naive algorithm multiplies two elements of Cliff(n) in time O(22n). It's an easy exercise to see that the multiplication time is O(27n/4). Not bad for something that was free. We can acualy do quite a bit better with some small changes - for example by using Strassen-like tricks for all but the Split multiplications.
But if you really want to do this stuff properly check out this paper which I haven't yet read properly.
I'd also like to add that Cliff(n+8)=R(16)⊗Cliff(n), so we have a kind of repeat with period 8. This is known as Bott periodicity and it lies at the heart of many of the most beautiful things in mathematics including some fascinating properties of the number 8.