Balancing Chemical Equations Is Just Linear Algebra (Here's the Code)

If you've ever helped a friend with chemistry homework, you've seen it: they stare at C₃H₈ + O₂ → CO₂ + H₂O, add a coefficient, break something else, and erase the page down to the paper fibers.

👉 Try the balancer here

Most students are taught to guess and check. As a developer, you can see a cleaner way to look at it: balancing a chemical equation is a system of linear equations.

This post covers how that works, gives a working Python implementation, and links to a free web tool if you'd rather not write the parser yourself.

The problem

A balanced equation has the same number of each atom on both sides. You can change only the coefficients (the numbers in front of each compound), never the subscripts.

For propane combustion:

a C3H8 + b O2 → c CO2 + d H2O

Count atoms per element:

Element	Equation
C	3a = c
H	8a = 2d
O	2b = 2c + d

Move everything to one side and you have a homogeneous linear system:

3a + 0b - 1c + 0d = 0
8a + 0b + 0c - 2d = 0
0a + 2b - 2c - 1d = 0

In matrix form that's A · x = 0. The balanced coefficients are the null space of A, scaled to the smallest positive integers.

Step 1: Parse the formula

Formulas like Ca3(PO4)2 need nested parentheses handled. A small stack-based parser does it:

import re

def parse_formula(formula):
    tokens = re.findall(r"[A-Z][a-z]?|\(|\)|\d+", formula)
    stack, i = [{}], 0
    while i < len(tokens):
        t = tokens[i]
        if t == "(":
            stack.append({})
        elif t == ")":
            group = stack.pop()
            mult = 1
            if i + 1 < len(tokens) and tokens[i + 1].isdigit():
                mult = int(tokens[i + 1]); i += 1
            for el, n in group.items():
                stack[-1][el] = stack[-1].get(el, 0) + n * mult
        else:
            n = 1
            if i + 1 < len(tokens) and tokens[i + 1].isdigit():
                n = int(tokens[i + 1]); i += 1
            stack[-1][t] = stack[-1].get(t, 0) + n
        i += 1
    return stack[0]

print(parse_formula("Ca3(PO4)2"))
# {'Ca': 3, 'P': 2, 'O': 8}

Step 2: Build the matrix and find the null space

Reactant atoms count as positive, product atoms as negative. Then let SymPy do the exact rational arithmetic (floats would give you coefficients like 2.9999999):

python
from sympy import Matrix, ilcm, igcd

def balance(reactants, products):
    parsed = [parse_formula(s)