Counting semi magic squares using generating functions and Sage.

+John Cook has written a couple of posts recently that caught my eye and they kind of start with this one: Rolling dice for normal sample: Python version.

Go check it out it's pretty cool but basically John Cook shows how to use the Sympy package to obtain coefficients from generating functions in python.

I thought I'd do something similar but with +Sage Mathematical Software System (which is built on python and would be using Sympy on this occasion) but also on a slightly different topic: Semi Magic Squares.

A semi magic square is defined as a square matrix whose row and columns sums are equal to some constant $r$.

So for example semi magic squares of size 3 will be of the form:

$\begin{pmatrix} a&b&r-a-b\\ c&d&r-c-d\\ r-a-c&r-b-d&a+b+c+d-r \end{pmatrix}$

(with $0\leq a, b, c, d, a+c, b+d, a+b, c+d \leq r \leq a+b+c+d$ and $a,b,c,d \in \mathbb{Z}$)

Using this it can be shown that $|\text{SMS}(3,r)|=\binom{r+4}{4}+\binom{r+3}{4}+\binom{r+2}{4}$ (where SMS(n,r) is the set of semi magic squares of size $n$ and line size $r$). A really nice and accessible proof of this is in this paper by Bona: A New Proof of the Formula for the Number of 3 by 3 Magic Squares.

That is however not the point of this blog post :) In fact I'm just going to show the Sage code needed to enumerate the set SMS(n,2).

There's a nice result obtained in A combinatorial distribution problem by Anand Dumir and Gupta in 1966 that expresses the number of elements in SMS(n,r) as a generating function:

$\frac{e^{\frac{x}{2}}}{\sqrt{1-x}}=\sum_{n=0}^{\infty}|\text{SMS}(n,2)|\frac{x^n}{n!^2}$

In the rest of this post I'll now show how to use +Sage Mathematical Software System to get the numbers of Semi Magic Squares of line sum 2 and size $n$.

First of all we need to define the above generating function. In Sage this is easy and as we're going to use the variable $x$ don't need to declare it as a symbolic variable (this is by default for $x$ in sage).

Here's the natural way of doing this:

sage: f(x) = e ^ (x / 2)/(sqrt(1 - x))

We can get a plot of $f(x)$ easily in Sage:

sage: plot(f(x),x,-10,1)

which gives:

To obtain the taylor series (up to $x^5$) expansion of this function we simply call:

sage: taylor(f(x),x,0,5)

Which returns:

69/160*x^5 + 47/96*x^4 + 7/12*x^3 + 3/4*x^2 + x + 1

Recalling the generating function formulae above we see that the size of $|\text{SMS}(5,2)|$ is: $(5!)^2\times \frac{69}{160}=6210$.

We can use usual python syntax to define a function that will return $|\text{SMS}(n,2)|$ for any $n$:

sage: def SMS(n,2):
....:    return taylor(f(x),x,0,n).coeff(x^n)*(factorial(n)^2)

This makes use of the coeff method on the polynomials in sage.

To get $|\text{SMS}(n,2)|$ for $0\leq n \leq 10$ we can use:

sage: print [SMS(n,2) for n in range(11)]

which gives:

[0, 1, 3, 21, 282, 6210, 202410, 9135630, 545007960, 41514583320]

If we throw that sequence in to the OEIS (my best friend during my PhD) we indeed get the correct sequence: A000681 (if you don't know about the OEIS it's awesome, be sure to click one of those two links).

In my previous post today I just used pure Python to carry out simulations to estimate $\pi$ but there are certain things that I find Sage to be better suited for (like generating functions :) ).

In my PhD I looked at a set of objects similar to Semi Magic Squares called: Alternating Sign Matrices.

You can find some of the arguments I briefly talk about here in my thesis (the thought of someone actually reading it is terrifying). I used Mathematica throughout my PhD but (without wanting to sound like a sales pitch) if I had known about Sage (not a sponsor) at the time I would have probably used Sage.

(All the math on this page is rendered using http://www.mathjax.org/ which I'm not entirely sure to have configured correctly for +Blogger; if some of the math hasn't rendered try reloading the page...)