Commutative Algebra in Haskell, part 1

Once again, it's been a little while since my last post, and once again, my excuse is partly that I've been too busy writing code.

I've just uploaded a new release, HaskellForMaths 0.3.4, which contains the following new modules:

Math.Core.Utils - this is a collection of utility functions used throughout the rest of the library. I've belatedly decided that it's better to put them all in one place rather than scattered here and there throughout other modules.

Math.Core.Field - this provides new, more efficient implementations of several finite fields. There already were implementations of these finite fields, in the Math.Algebra.Field.Base and ...Extension modules, as discussed here and here. However, that code was written to make the maths clear, rather than for speed. This new module is about speed. For the prime power fields in particular (eg F4, F8, F9), these implementations are significantly faster.

Math.Combinatorics.Matroid - Matroids are a kind of combinatorial abstraction of the concept of linear independence. They're something that I heard about years ago - both of my favourite combinatorics books have brief introductions - but I never bothered to follow up. Well anyway, so something finally piqued my curiosity, and I got Oxley's Matroid Theory. It turned out to be really interesting stuff, and this module is pretty much a translation of a large part of that book into Haskell code, written as I taught myself all about matroids.

Math.CommutativeAlgebra.Polynomial - Although I hadn't yet got around to discussing them in the blog, HaskellForMaths has always had modules for working with multivariate polynomials, namely Math.Algebra.Commutative.Monomial and ...MPoly. However, these were some of the earliest code I wrote, before my more recent free vector space and algebra code. So I saw an opportunity to simplify and improve this code, by building it on top of the free vector space code. Also, I'm trying to rationalise the module naming convention in HaskellForMaths, to more closely follow the categories used in arxiv.org or mathoverflow.net . In the long run, I expect this module to supercede the older modules.

Math.CommutativeAlgebra.GroebnerBasis - Again, there was already code for Groebner bases in Math.Algebra.Commutative.GBasis. This is pretty much the same code, ported to the new polynomial implementation, but I've also begun to build on this, with code to find the sum, product, intersection, and quotient of ideals.


So the matroid code was just new code that I wrote while teaching myself some new maths. But most of the other code comes from an ambition to organise and simplify the HaskellForMaths library. I've also been trying to improve the documentation.

My ultimate ambition is to get more people using the library. To do that, the structure of the library needs to be clearer, the documentation needs to be better, and I need to explain how to use it. So I thought I'd start by explaining how to use the new commutative algebra modules.

(So this is a bit of a digression from the series on quantum algebra that I've been doing the last few months. However, in terms of the cumulative nature of maths, it's probably better to do this first.)


Okay, so suppose we want to do some polynomial arithmetic. Well, first we need to create some variables to work with. How do we do that?

First, decide on a monomial ordering - that is, we need to decide in what order monomials are to be listed within a polynomial. For the moment, let's use "graded lexicographic" or Glex order. This says that you should put monomials of higher degree before those of lower degree (eg y^3 before x^2), and if two monomials have the same degree, you should use lexicographic (dictionary) order (eg xyz before y^3).

Next, decide on a field to work over. Most often, we'll want to work over Q, the rationals.

Then, our variables themselves can be of any Haskell type - but there are usually only two sensible choices:

The easiest way is to use String as the type for our variables.

Then we could make some variables like this:

> :l Math.CommutativeAlgebra.Polynomial
> let [x,y,z] = map glexvar ["x","y","z"]

And then we can do polynomial arithmetic:

> (x+y+z)^3
x^3+3x^2y+3x^2z+3xy^2+6xyz+3xz^2+y^3+3y^2z+3yz^2+z^3


If we want to use any other field besides Q, then we will have to use a type annotation to tell the compiler which field we're working over:

> let [x,y,z] = map var ["x","y","z"] :: [GlexPoly F3 String]
> (x+y+z)^3
x^3+y^3+z^3


The alternative to using String for our variables is to define our own type. For example

data Var = X | Y | Z | W deriving (Eq,Ord)



instance Show Var where

    show X = "x"

    show Y = "y"

    show Z = "z"

    show W = "w"



[x,y,z,w] = map glexvar [X,Y,Z,W]

So there you have it - now you can do polynomial arithmetic in Haskell.

So how does it work?

Well, fundamentally, k-polynomials are a free k-vector space on the basis of monomials. So we define a type to implement monomials:

data MonImpl v = M Int [(v,Int)] deriving (Eq)

-- The initial Int is the degree of the monomial. Storing it speeds up equality tests and comparisons



instance Show v => Show (MonImpl v) where

    show (M _ []) = "1"

    show (M _ xis) = concatMap (\(x,i) -> if i==1 then showVar x else showVar x ++ "^" ++ show i) xis

        where showVar x = filter ( /= '"' ) (show x) -- in case v == String

Notice that our monomial implementation is polymorphic in v, the type of the variables.

Next, monomials form a monoid, so we make them an instance of Mon (the HaskellForMaths class for monoids):

instance (Ord v) => Mon (MonImpl v) where

    munit = M 0 []

    mmult (M si xis) (M sj yjs) = M (si+sj) $ addmerge xis yjs

In principle, all we need to do now is define an Ord instance, and then an Algebra instance, using the monoid algebra construction.

However, for reasons that will become clear in future postings, we want to be able to work with various different orderings on monomials, such as Lex, Glex, or Grevlex. So we provide various newtype wrappers around this basic monomial implementation. Here's the code for the Glex ordering that we used above:

newtype Glex v = Glex (MonImpl v) deriving (Eq, Mon) -- GeneralizedNewtypeDeriving



instance Show v => Show (Glex v) where

    show (Glex m) = show m



instance Ord v => Ord (Glex v) where

    compare (Glex (M si xis)) (Glex (M sj yjs)) =

        compare (-si, [(x,-i) | (x,i) <- xis]) (-sj, [(y,-j) | (y,j) <- yjs])



type GlexPoly k v = Vect k (Glex v)



glexvar :: v -> GlexPoly Q v

glexvar v = return $ Glex $ M 1 [(v,1)]



instance (Num k, Ord v, Show v) => Algebra k (Glex v) where

    unit x = x *> return munit

    mult xy = nf $ fmap (\(a,b) -> a `mmult` b) xy

We also have similar newtypes for Lex and Grevlex orderings, which I'll discuss another time.

And that's pretty much it. Now that we have an instance of Algebra k (Glex v), we get a Num instance for free, so we get +, -, *, and fromInteger. That means we can enter expressions like the following:

> (2*x^2-y*z)^2
4x^4-4x^2yz+y^2z^2

Note that division is not supported: you can't write x/y, for example. However, as a convenience, I have defined a partial instance of Fractional, which does let you divide by scalars. That means that it's okay to write x/2, for example.

Next time, some more things you can do with commutative algebra.