I’ve never been very good at doing manual computations, and whenever I need to do a tedious computation for an assignment, I like to automate it by writing a computer program. Usually I implemented an ad-hoc solution using Haskell, either using a simple library or rolling my own implementation if the library didn’t have it. But I found this solution to be unsatisfactory: my Haskell programs worked with integers and floating numbers and I couldn’t easily generalize it to work with symbolic expressions. So I looked to learn a CAS (computer algebra system), so in the future I won’t have to hack together buggy code for common math operations.

I have no experience with symbolic computing, so it wasn’t clear to me where to begin. To start off, there are many different competing computer algebra systems, all incompatible with each other, and it’s far from clear which one is best for my needs. I began to experiment with several systems, but after a few days I still couldn’t decide which one was the winner.

I narrowed it down to 3 platforms. Here’s my setup (all running on Windows 7):

**Mathematica**8.0**Maxima**5.32 with wxMaxima 13.04**Maple**18.00

So I came up with a trial — I had a short (but nontrivial) problem representative of the type of problem I’d be looking at, and I would try to solve it in all 3 languages, to determine which one was easiest to work with.

### The Problem

This problem came up as a part of a recent linear algebra assignment.

Let the field be (so all operations are taken modulo 5). Find all 2×2 matrices such that

We can break this problem into several steps:

- Enumerate all lists of length 4 of values between 0 to 4, that is, [[0,0,0,0],[0,0,0,1],…,[4,4,4,4]]. We will probably do this with a cartesian product or list comprehension.
- Figure out how to convert a list into a 2×2 matrix form that the system can perform matrix operations on. For example, [1,2,3,4] might become matrix([1,2],[3,4])
- Figure out how to do control flow, either by looping over a list (procedural) or with a map and filter (functional)
- Finally, multiply the matrices modulo 5 and check if it equals the identity matrix, and output.

This problem encompasses a lot of the challenges I have with CAS software, that is, utilize mathematical functions (in this case, we only use matrix multiplication and transpose), yet at the same time express a nontrivial control flow. There are 5^4=625 matrices to check, so performance is not a concern; I am focusing on ease of use.

For reference, here is the answer to this problem:

These are the 8 matrices that satisfy the desired property.

I have no prior experience in programming in any of the 3 languages, and I will try to solve this problem with the most straightforward way possible with each of the languages. I realize that my solutions will probably be redundant and inefficient because of my inexperience, but it will balance out in the end because I’m equally inexperienced in all of the languages.

### Mathematica

I started with Mathematica, a proprietary system by Wolfram Research and the engine behind Wolfram Alpha. Mathematica is probably the most powerful out of the three, with capabilities with working with data well beyond what I’d expect from a CAS.

What I found most jarring about Mathematica is its syntax. I’ve worked with multiple procedural and functional languages before, and there are certain things that Mathematica simply does differently from everybody else. Here are a few I ran across:

- To use a pure function (equivalent of a lambda expression), you refer to the argument as #, and the function must end with the & character
- The preferred shorthand for Map is /@ (although you can write the longhand Map)
- To create a cartesian product of a list with itself n times, the function is called Tuples, which I found pretty counterintuitive

Initially I wanted to convert my flat list into a nested list by pattern matching Haskell style, ie f [a,b,c,d] = [[a,b],[c,d]], but I wasn’t sure how to do that, or if the language supports pattern matching on lists. However I ran across Partition[xs,2] which does the job, so I went with that.

Despite the language oddities, the functions are very well documented, so I was able to complete the task fairly quickly. The UI is fairly streamlined and intuitive, so I’m happy with that. I still can’t wrap my head around the syntax — I would like it more if it behaved more like traditional languages — but I suppose I’ll get the hang of it after a while.

Here’s the program I came up with:

SearchSpaceLists := Tuples[Range[0, 4], 4] SearchSpaceMatrices := Map[Function[xs, Partition[xs, 2]], SearchSpaceLists] Middle := {{2, 0}, {0, 3}} FilteredMatrices := Select[SearchSpaceMatrices, Mod[Transpose[#].Middle.#, 5] == IdentityMatrix[2] &] MatrixForm[#] & /@ FilteredMatrices

### Maxima

Maxima is a lightweight, open source alternative to Mathematica; I’ve had friends recommend it as being small and easy to use.

The syntax for Maxima is more natural, with things like lists and loops and lambda functions working more or less the way I expect. However, whenever I tried to do something with a function that isn’t the most common use case, I found the documentation lacking and often ended up combing through old forum posts.

Initially I tried to generate a list with a cartesian product like my Mathematica version, but I couldn’t figure out how to do that, eventually I gave up and used 4 nested for loops because that was better documented.

Another thing I had difficulty with was transforming a nested list into a matrix using the matrix command. Normally you would create a matrix with matrix([1,2],[3,4]), so by passing in two parameters. The function doesn’t handle passing in matrix([[1,2],[3,4]]), so to get around that you need to invoke a macro: funmake(‘matrix,[[1,2],[3,4]]).

Overall I found that the lack of documentation made the system frustrating to work with. I would however use it for simpler computations that fall under the common use cases — these are usually intuitive in Maxima.

Here’s the program I came up with:

Middle:matrix([2,0],[0,3]); Ident:identfor(Middle); for a:0 thru 4 do for b:0 thru 4 do for c:0 thru 4 do for d:0 thru 4 do (P:funmake('matrix,[[a,b],[c,d]]), P2:transpose(P).Middle.P, if matrixmap(lambda([x],mod(x,5)),P2) = Ident then print(P));

Shortly after writing this I realized I didn’t actually need the funmake macro, since there’s no need to generate a nested list in the first place, I could simply do matrix([a,b],[c,d]). Oh well, the point still stands.

### Maple

Maple is a proprietary system developed by Maplesoft, a company based in Waterloo. Being a Waterloo student, I’ve had some contact with Maple: professors used it for demonstrations, some classes used it for grading. Hence I felt compelled to give Maple a shot.

At first I was pleasantly surprised that matrix multiplication in a finite field was easy — the code to calculate A*B in is simply A.B mod 5. But everything went downhill after that.

The UI for Maple feels very clunky. Some problems I encountered:

- It’s not clear how to halt a computation that’s in a an infinite loop. It doesn’t seem to be possible within the UI, and the documentation suggests it’s not possible in all cases (it recommends manually terminating the process). Of course, this loses all unsaved work, so I quickly learned to save before every computation.
- I can’t figure out how to delete a cell without googling it. It turns out you have to select your cell and a portion of the previous cell, then hit Del.
- Copy and pasting doesn’t work as expected. When I tried to copy code written inside Maple to a text file, all the internal formatting and syntax highlighting information came with it.
- Not an UI issue, but error reporting is poor. For example, the = operator works for integers, but when applied to matrices, it silently returns false. You have to use Equals(a,b) to compare matrices (this is kind of like java).

In the end, I managed to complete the task but the poor UI made the whole process fairly unpleasant. I don’t really see myself using Maple in the future; if I had to, I would try the command line.

Here’s the program I came up with:

with(LinearAlgebra): with(combinat, cartprod): L := [seq(0..4)]: T := cartprod([L, L, L, L]): Middle := <2,0;0,3>: while not T[finished] do pre_matrix := T[nextvalue](); matr := Matrix(2,2,pre_matrix); if Equal(Transpose(matr).Middle.matr mod 5, IdentityMatrix(2)) then print(matr); end if end do:

### Conclusion

After the brief trial, there is still no clear winner, but I have enough data to form some personal opinions:

- Mathematica is powerful and complete, but has a quirky syntax. It has the most potential — definitely the one I would go with if I were to invest more time into learning a CAS.
- Maxima is lightweight and fairly straightfoward, but because of lack of documentation, it might not be the best tool to do complicated things with. I would keep it for simpler calculations though.
- Maple may or may not be powerful compared to the other two, I don’t know enough to compare it. But its UI is clearly worse and it would take a lot to compensate for that.

I’ve posted this discussion on some discussion boards for Mathematica and Maple, you might find them informative.

It’s generally bad form in Mathematica to use “:=” (SetDelayed) instead of “=” (Set) unless defining a function. There are exceptions to this but if you are assigning a value to a variable “=” is almost always the correct choice. Also you don’t need to make a pure function in the last line of the code, you can just use the function name. (Also end lines with “;” so suppress the output)

Here is a slightly cleaner version of your code:

searchSpaceLists = Tuples[Range[0, 4], 4];

searchSpaceMatrices =

Map[Function[xs, Partition[xs, 2]], searchSpaceLists];

middle = {{2, 0}, {0, 3}};

filteredMatrices =

Select[searchSpaceMatrices,

Mod[Transpose[#].middle.#, 5] == IdentityMatrix[2] &];

MatrixForm /@ filteredMatrices

LikeLike

Thanks for the quick response! I suspected I’d do things wrong the first time around.

LikeLike

Where does Mathcad fit in? I find it powerful and easy to use. Also, with full-screen editing, it can produce nice-looking documents.

TKB

LikeLike

OK. I have a question. Why don’t you include Mathcad in your comparisons? It has some distinctive advantages [entry much like one would write on a board; regular math notation; full-screen editing], but seems generally to be ignored in blogs on CAS. I understand that not much is free, and that it has a strong engineering flavor.

TKB

________________________________

LikeLike

I’ve never heard of Mathcad before; of course there are numerous CAS systems other than Mathematica / Maxima / Maple and I can’t cover all of them.

LikeLike

p={{a,b},{c,d}}; pt=Transpose[p]; middle={{2,0},{0,3}};

Map[MatrixForm[p/.#]&, {ToRules[Reduce[

pt.middle.p==IdentityMatrix[2], {a,b,c,d},Modulus->5]]}]

LikeLike

Simpler

p={{a,b},{c,d}}; pt=Transpose[p]; middle={{2,0},{0,3}};

Map[MatrixForm[p/.#]&, Solve[pt.middle.p==IdentityMatrix[2],

{a,b,c,d}, Modulus->5]]

LikeLike

You probably can shorten the Maple code to half its size:

restart;

L := [seq(0..4)]:

T := combinat:-cartprod([L, L, L, L]):

while not T[finished] do

z:= Matrix(2,2, T[nextvalue]() );

if LinearAlgebra:-Equal( z^%T..z mod 5,) then print(z) fi;

od:

LikeLike

Maple can solve the problem in a “functional” style too:

L := []:

T := combinat:-cartprod([seq([1, 2, 3, 4], i = 1 .. 4)]):

while `not`(T[finished]) do L := [op(L), Matrix(2, 2, T[nextvalue]())] end do:

select(proc (x) LinearAlgebra:-Equal(`mod`(Transpose(x).Matrix([[2, 0], [0, 3]]).x, 5), IdentityMatrix(2)) end proc, L);

LikeLike

Came across your blog when writing up this one from a somewhat different perspective http://thingwy.blogspot.de/

LikeLike

The Maxima solution can be simplified quite a bit..

There is no need for (arcane) functions like funmake or matrixmap.

Also there is a standard function for identity matrices (ident(2) for a 2×2 one). And I prefer curly brackets for control blocks.

Then we get:

for a:0 thru 4 do

for b:0 thru 4 do

for c:0 thru 4 do

for d:0 thru 4 do

{

M: matrix([2,0],[0,3]),

P: matrix([a,b],[c,d]),

if mod(transpose(P).M.P, 5) = ident(2)

then print(P)

};

LikeLike

you can have Maxima do all the arithmetic mod 5 by setting

modulus:5.

It chooses the representation 0 1 2 -2 -1 instead of 0,1,2,3,4.

vals:map(rat,[0,1,2,3,4];

for a in vals do for b in vals do for c in vals do for d in vals

do (p:matrix([a,b],[c,d]),

if equal(transpose(p).matrix([2,0],[0,3]).p,

ident(2) )

then print (p));

LikeLike

I have just come across this post while searching for something else. I like it. Certainly an expert in any of those systems will construct a program which is efficient and simple. But the point here is how a non-expert would manage, starting from scratch, and with no prior background. I’m interested though in why those particular three were chosen? For matrix computations many people would use Matlab or something similar, and of course there’s also Sage. But I thought that this post was highly informative.

LikeLike

I haven’t looked into Sage, but I have a little bit of experience with Matlab. Matlab seems to be best suited for simulations and stuff involving heavy floating-point number crunching. In this post I’m only interested in comparing symbolic computer algebra systems.

LikeLike

Try using sagemath for a bit. You can do some nice things with the object oriented design and the whole python functionality. Your problem could simply be solved like this:

R = ZZ.quo(5) # Saving the Ring for convenience

MS = MatrixSpace(R, 2) # Space of 2×2-Mat. over Z/5Z

M = matrix(R, 2, [2,0,0,3]) # Def. Matrix

[x for x in MS if (x.transpose()*M*x).is_one()] # List of results

LikeLike

By the way, due to object orientation this also works for different Rings – let’s say, quotient polynomial rings over finite fields:

P. = ZZ.quo(2)[] # Polynomial Ring over Z/2Z in X

p = X**3 + X**2 + 1 # some Polynom in P

R = P.quo(p) # quotient polynomial ring

MS = MatrixSpace(R, 2) # matrix space over R

M = MS.random_element() # some random matrix

[x for x in MS if (x.transpose()*M*x).is_one()] # List of results

LikeLike

Mathematica has come a long way since this article was written. The solution in Mathematica can be found just using the built in functions without coding any of your own (which to me is a major use of a CAS). Here is the solution for Mathematica 11 in 6 lines.

P = {{a, b}, {c, d}};

Mid = {{2, 0}, {0, 3}};

Id = IdentityMatrix[2];

max = 5^4;

sols = FindInstance[{Transpose[P].Mid.P == Id}, {a, b, c, d}, max, Modulus -> 5];

Table[MatrixForm[{{a, b}, {c, d}} /. sols[[i]] ], {i, 1, Length[sols]}]

The final line prints the solutions in matrix form. Note: I have never solved integer equations in mathematica or used “FindInstance” but was able to code this up in just a few minutes. I really do think Mathematica is the best followed by either Maxima or Maple.

LikeLike

Wow, that’s like magic! Do you know which version of Mathematica started being able to solve matrix equations like this?

LikeLike

I’m not sure, I always remember being able to solve matrix equations since Mathematica 8 but I’m not sure when the restricted (modular) domains were introduced. If you are interested enough to look into it, the history of Mathematicas documentation is linked below.

https://reference.wolfram.com/history/

LikeLike