Symbolic computation with Clojure and Python

While putting together a course in the fundamentals of mathematics centered around the use of computer algebra systems (CAS), I've gone through several programming languages and feel like I've learned a lot.

As a Clojurian, I began by writing Clojure functions and publishing them as web pages with interactive code snippets. When it felt like I was beginning to hit a wall rolling my own solvers and such, I thought to learn other systems to see how they did it.

The first one I tried was Mathematica, which made it almost too easy. I did not like its proprietary nature, however I very much fell in love with the interactive notebook interface. This naturally led me to using Python in Jupyter Notebooks with the SymPy library, as it seemed to be the closest thing to an open source Wolfram notebook.

I was very pleased with the Python data science ecosystem as well as the convenient share-ability of the .ipynb format. But now I've found something even better, which lets me have all of the benefits I've mentioned while bringing Clojure back centerstage - nextjournal! The biggest advantage I've noticed so far is the ability to use multiple languages in the same notebook. As a demonstration, we will complete the square both in Clojure and Python.

This is a work in progress. Your feedback is much appreciated, and this new platform makes it super-easy to remix the notebook if anyone happens to feel so inspired! Enjoy using Clojure with nextjournal!

View or remix this notebook on nextjournal.

Rewrite the quadratic equation by completing the square.

Clojure

A quadratic equation consists of a quadratic coefficient :a, linear coefficient:b and a constant :c on the left, and 0 on the right. Let's write a function that will take these values and render it as a human-readable string:

(defn print-quadratic-equation [[a b c & r]]
  (str (if (not= 1 a) a) "x^2 "
       (if (not= 1 b)
         (if (pos? b)
           (str "+ " b "x ")
           (str "- " (- b) "x ")))
       (if (not= 0 c)
         (if (pos? c)
           (str "+ " c " ")
           (str "- " (- c) " ")))
       "= " (if (first r) (first r) 0)))

(-> [4 20 25]
  print-quadratic-equation)

The right side of the equation is assumed to be 0 if no r term is supplied. We've included conditionals for not rendering a constant of 0 or a coefficient of 1 because that would be redundant. If a number is negative, we print a minus - sign and its negation.

Now we can write functions that will return a new equation with the appropriate manipulation. Let's first move the constant :c over to the right by subtracting it from both sides:

(defn subtract-c [[a b c r]]
  [a b (- c c) (- 0 c)])

(-> [4 20 25]
  subtract-c
  print-quadratic-equation)

Now divide by the quadratic coefficient :a to get x^2 by itself:

(defn divide-by-a [[a b c r]]
  [1 (/ b a) 0 (/ r a)])

(-> [4 20 25]
  subtract-c
  divide-by-a
  print-quadratic-equation)

Now we want to complete the left side into a perfect square. To do that, we take half of :b, square it, and add it to both sides.

(defn complete-square [[a b c r]]
  (let [square (* (/ b 2) (/ b 2))]
    [a b square (+ square r)]))

(-> [4 20 25]
  subtract-c
  divide-by-a
  complete-square
  print-quadratic-equation)

We can now rewrite the left side of the equation as a squared term.

(defn rewrite-squared [[a b c r]]
  (str "(x " (if (pos? (/ b 2))
               (str "+ " (/ b 2))
               (str "- " (- (/ b 2))))
       ")^2 = " r))
       
(-> [4 20 25]
  subtract-c
  divide-by-a
  complete-square
  rewrite-squared)

Python

We'll be using SymPy, a Python library for symbolic computation.

Unlike other computer algebra systems like Mathematica, we need to first define the variable(s) we'll be using:

from sympy import *
x = Symbol('x')

In SymPy, an equation is a special data type created with the Eq function, where the arguments are the left and right sides of the equation separated by a comma:

Eq(4*x**2 + 20*x + 25, 0)

Also note that no implicit multiplication is allowed. That is,4x must be 4*x.

As before we will write a function to subtract the left constant from both sides:

def sub_const_left(a, b, c):
    return Eq(a*x**2 + b*x, -c)

sub_const_left(4, 20, 25)

Eq(4*x**2 + 20*x, -25)

Now we can divide by our leading coefficient:

def divide_by_a(a, b, c):
    return Eq(x**2 + (b*x / a), -Rational(c, a))

divide_by_a(4, 20, 25)

Eq(x**2 + 5*x, -25/4)

Now we want to complete the left side into a perfect square. To do that, we take half of :b, square it, and add it to both sides.

To be continued...