Skip to Content »

Michi’s blog » Computational Group Theory in Haskell (1 in a series)

 Computational Group Theory in Haskell (1 in a series)

  • October 18th, 2006
  • 3:37 pm

This term, I’m listening to a lecture course on Computational Group Theory. As a good exercise, I plan to implement everything we talk about as Haskell functions as well.

The first lecture was mainly an introduction to the area, ending with a very naïve algorithm to generate a permutation group from a set of generators. Next week will bring less naïve algorithms with not quite as horrible complexity.

Before the algorithm can be brought, we’d want some undergrowth: we’d want to be able to work with permutations at all. So, we’ll start with the basic group theory and permutation implementations. A lot of this is stolen or rewritten from this permutation group code.

Our code will make use of two libraries, so if you collate code snippets while reading this, you’ll want to use

import Array
import List
 

If you don’t want to bother with that, the code is available here.

One of the things I love the most with Haskell is the way that you just write down what you think it means, and end up with good code. So, what is a group? It is this type class:

class Group a where
  identity :: Int -> a
  isIdentity :: a -> Bool
  (*>) :: a -> a -> a     – right multiplication
  inverse :: a -> a
 

with some conditions on the functions given. This is geared towards implementing permutation groups in the first line — so the idea is that you give your group as a subgroup of Sn, and you can retrieve the identity in Sn by using the identity function.

We tend to want to be able to conjugate elements and multiply a whole list of elements, so let’s add a few functions:

conjugate g h = inverse h *> g *> h
prod gs = foldl1 (*>) gs                                                        
 

And now we can start looking at permutations. We will store a permutation as the list of images, in order. So for instance, the permutation given by (1 2)(1 3) will be stored as the list [2,3,1]. (note that in the good tradition of computational group theory, we act from the right – thus motivating this result and not [2,1,3])

The interesting bit about permutations, though, is that they form a group; so we instantiate our type class as well.

newtype Permutation = PL [Int] deriving (Eq, Ord, Read, Show)                  
x .^ (PL g) = g !! (x-1)                    
instance Group Permutation where
  identity n = PL [1..n]                                                          
  isIdentity (PL xs) = all (\(i,j) -> i==j) (zip [1..] xs)                        
  (PL g) *> h = PL (map (.^ h) g)                                              
  inverse (PL g) = PL (map snd (sort (zip g [1..])))                                          
 

We are used to being able to write permutations in cycle notation, so we’d want to be able to convert to and from that, too.

fromCycles cs = PL (map snd (sort pointActions))                                
  where                                                                        
    pointActions = concat (map fromCycle cs)                                                                                                                    fromCycle is = zip is (rotateL is)                                                                                                                              fromCycles’ n cs = PL (elems (array (1,n) [(i,i) | i <- [1..n]] // (concat (map fromCycle cs))))                                                                
findOrbit (PL permutation) elt =                                                
  let buildOrbit os ps e =                                                          
    if (elem e os)                                                              
      then os                                                                  
      else let ei = elemIndex e ps in                                                
        case ei of                                                                
          Nothing -> os                                                        
          Just e’ -> buildOrbit ([e] ++ os) ps (e’+1)                  
  in buildOrbit [] permutation elt                                                                                                                              toCycles (PL permutation) =                                                    
  let
    getOrbits out [] per = out                                                  
    getOrbits out els per =                                                          
      let
        orbit = findOrbit per (head els)
        els’ = els \\ orbit                                                        
      in getOrbits (out ++ [orbit]) els’ per
  in getOrbits [] permutation (PL permutation)                                  
rotateL (x:xs) = xs ++ [x]                                                      
saturateCycles :: [[Int]] -> Int -> [[Int]]                                     saturateCycles cs n = cs ++ ( group ( [1..n] \\ foldl1 (++) cs ) )              
desaturateCycles cs = filter (\x -> length x > 1) cs
 

We can use saturateCycles and desaturateCycles in order to convert to and from the normal way of writing permutations. The toCycles and fromCycles functions both expect fixpoints to be listed as well. We need, however, to tell the saturateCycles which symmetric group we want to live in.

Now, the most naïve algorithm to list all elements in a group starts with the one element we know will live in the group: the identity. And for each iteration, all generators are multiplied onto the first known candidate in a list of group element candidates, and all results that are not yet known are added to this list.

In Haskell code:

generateGroup generators =
  let g0 = head generators
      i = g0 *> inverse g0
      addElts group [] gens = group
      addElts group cands gens =
        let cand = take 1 cands
            group’ = group ++ cand
            newcands = map (*> (head cand)) gens
            cands’ = (tail cands) ++ (newcands \\ (group’ ++ cands))
        in  
            addElts group’ cands’ gens
  in addElts [] [i] generators
 

We keep on adding tested candidates to the first list, and anything not yet seen that we get from our multiplication with generators gets added to the middle list. Once the middle list is empty, we’re done.

4 People had this to say...

Gravatar

Are you sure you want identity :: n -> a rather than identity :: a

Gravatar
  • Michi
  • October 19th, 2006
  • 14:56

Actually, I do, since I view every group as a subgroup of some symmetric group. The identity is a function taking as argument which symmetric group I operate in, and returning the appropriate identity element.

Sure, translating group axioms straight off would make identity :: a be much saner; but on the other hand, I write the group type class with a specific use in mind.

Gravatar

Hang on, I’m sure there was meant to be more of that comment, to the effect that what you really want is for each symmetric group to be a type in its own right. But then you’d need types parametrized by integers, which you can do in C++ or Ada, but not (IIRC) in Haskell. I remember running into this problem a while ago while implementing the algorithms from my Computational Methods in Finite Fields course.

Gravatar
  • Michi
  • October 19th, 2006
  • 15:38

Parametrizing the type “Subgroup of a symmetric group” would make a lot of things much easier; handling the cycle representation would look different, as would the multiplication code. But I’m not very much inclined to start writing things like
newtype Permutation3 = ..., and Haskell does not, at the current stage, seem to allow parametrized types.

There may possibly be some extension floating around in an alpha version capable of it though…

Want your say?

* Required fields. Your e-mail address will not be published on this site

You can use the following XHTML tags:
<a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>