Symbolic Computation with Python

In the beginning, computers were used for numerical computations, and it was good. But it was soon discovered that computers could perform symbolic computations as well. The first computer algebra system was Macsyma, which was developed by MIT in 1968. Today, the leading computer algebra systems are Mathematica and Maple. These systems are very powerful; but they are closed-source commercial products, and they can be quite expensive.

Fortunately, there are some excellent free alternatives. SymPy is an open-source Python library for symbolic computation. In this note, I will demonstrate the solution of two algebraic problems using SymPy.

A problem from Crux

The following problem is from Crux Mathematicorum, a journal for secondary and undergraduate students published by the Canadian Mathematical Society.

M379. The integers $27 + C$, $555 + C$, and $1371 + C$ are all perfect squares, the square roots of which form an arithmetic sequence. Determine all possible values of C. (Crux v35 n1, Feb 2009)

We need to solve the following nonlinear system of equations.

(x-y)^2 &= 27 + C \\
x^2 &= 555 + C \\
(x+y)^2 &= 1371 + C

The system can be solved manually by combining the equations in clever ways. (Hint: $(x-y)^2 + (x+y)^2 = 2x^2 + 2y^2$.) We present a solution using Python and SymPy.

from sympy.solvers import solve
from sympy import symbols

x, y, C = symbols('x y C')
system = [(x - y)**2 - (27 + C),
          x**2 - (555 + C),
          (x + y)**2 - (1371 + C)]
print solve(system)

The code is almost self-explanatory. Line 4 creates three symbols named x, y, and C. In lines 5-7 we define an array of equations (the right side of each equation is assumed to be zero). We print the solution in line 8. The program produces the following output:

[{C: 229, x: -28, y: -12}, {C: 229, x: 28, y: 12}]

Euler’s four-square Identity

Our second example is Euler’s four-square identity. This identity shows that if two numbers $A$ and $B$ can be expressed as the sum of four squares, then the product $AB$ can also be so expressed. This is a key step in the proof of Lagrange’s four square theorem, which states that every positive integer can be expressed as the sum of four squares.

Behold Euler’s four-square identity in all its glory.
(a_1^2+a_2^2+a_3^2+a_4^2)(b_1^2+b_2^2+b_3^2+b_4^2)= \\
(a_1 b_1 + a_2 b_2 + a_3 b_3 + a_4 b_4)^2 +\\
(a_1 b_2 – a_2 b_1 + a_3 b_4 – a_4 b_3)^2 +\\
(a_1 b_3 – a_2 b_4 – a_3 b_1 + a_4 b_2)^2 +\\
(a_1 b_4 + a_2 b_3 – a_3 b_2 – a_4 b_1)^2.

This identity can be verified by hand. The calculation is straightforward, but tedious. Why not let the computer do the work?

from sympy import symbols, expand

a1, a2, a3, a4 = symbols("a1 a2 a3 a4")
b1, b2, b3, b4 = symbols("b1 b2 b3 b4")

c1 = a1*b1 + a2*b2 + a3*b3 + a4*b4
c2 = a1*b2 - a2*b1 - a3*b4 + a4*b3
c3 = a1*b3 - a3*b1 + a2*b4 - a4*b2
c4 = a1*b4 - a4*b1 - a2*b3 + a3*b2

aa = a1**2 + a2**2 + a3**2 + a4**2
bb = b1**2 + b2**2 + b3**2 + b4**2
cc = c1**2 + c2**2 + c3**2 + c4**2

print expand(aa*bb - cc)

The result is 0, which confirms the identity.

Getting started with SymPy

I recommend the Anaconda distribution of Python, because it includes SymPy and most other tools that you will need for scientific computation. There is also a “live” version of SymPy that runs in a web browser and does not require installing additional software.