Learning Python with Advent of Code Walkthroughs

Dazbo's Advent of Code solutions, written in Python

The Python Journey - Maths with SymPy


Useful Links

SymPySymPy Tutorial

Page Contents

Work in progress…

What is SymPy?

SymPy is so cool! I came across when I was trying to solve some equations as part of 2023 AoC puzzles. It is computer alegebra systems (CAS): a library that allows us to perfrom algebra and mathematical computations.

You can use it for things like:

Basic Examples

Simplifying and Expanding

We need to define variables before we can use them in SymPy. We do this by creating symbols.

One interesting observation here is that when you evaluate a SymPy expression a Jupyter notebook, it is rendered in pretty-printed (mathematical) format. But when we print using print(), it is printed in conventional Python style.

If you want to explicitly pretty-print from a Python notebook, you can use IPython.display.

Alternatively, you can optionally run init_printing() to enable the best printer for your environment. Or run init_session() to automatically import everything in SymPy, create some common Symbols, setup plotting, and run init_printing().

import sympy
from sympy import latex
from IPython.display import display, Markdown

sympy.init_session() # set up default behavour and printing - not essential

a, b, x, y = sympy.symbols("a b x y") # define symbols (variables)
a = (x + 1)**2 # set a to be an expresssion
b = x**2 + 2*x + 1 # set b to be an equivalent expression
display(Markdown(f"$a = {latex(a)}$")) # pretty-print a
print(f"{latex('a = ')}{latex(a)}") # If we just want the raw latex

display(Markdown(f"$b = {latex(b)}$")) # pretty-print b

# Expand a...
display(Markdown("Expanding $a$..."))
expanded = sympy.expand(a)
display(Markdown(f"$a = {latex(expanded)}$"))

# Show a-b without simplifying...
a_minus_b = a-b
display(Markdown(f"$a-b = {latex(a_minus_b)}$"))

# Simplify...
display(Markdown(f"$a-b$ $= {latex(sympy.simplify(a_minus_b))}$"))

Simplifying and expanding

We can also test if two expressions are the same, e.g.


Assigning Variables and Evaluating

You can’t simply assign values to the Python variable. If you did so, they would cease to be a symbol. Instead, you need to use subs() to associate a value with a SymPy symbol.

x, y = sympy.symbols("x y")
expr = x+y
result = expr.subs({x: 2, y: 5})
print(f"{expr=}, {result=}")

SymPy variables

Differentiation and Integration

import sympy

x = sympy.symbols("x")
expr = x**2

deriv = sympy.diff(expr)

integral = sympy.integrate(deriv)

SymPy basics

Solving Quadratics

Here, I use solve() to determine the values of $x$. I’ve specified dict=True such that $x$ can be retrieved by passing its symbol as a dictionary key to the returned solutions dictionary.

expr = x**2 + 3*x - 10
solutions = sympy.solve(expr, x, dict=True)

for solution in solutions:

SymPy quadratic

Note that the expressions still need to be written in valid Python. So we can write 3*x, but we can’t write 3x.

Including Only Real Solutions

What if you want to ignore complex solutions and only include real solutions? In this case, you can use solveset(), and specify the domain as domain=sympy.S.Reals. (Whereas domain=S.Complexes is the default.)

expressions = []
expressions.append(x**2 - 9) # there are only real solutions
expressions.append(x**2 + 9) # there are only complex solutions

for expr in expressions:
    solutions = sympy.solveset(expr, x)
    print(f"There are {len(solutions)} solutions for x...")

solutions = sympy.solveset(expressions[-1], x, domain=sympy.S.Reals)
print(f"There are {len(solutions)} real solutions for x.")


Evaluating to Obtain Numeric Values

We can evaluate to obtain the float value of a symbol, and display to an arbitrary level of precision:

expr = sympy.sqrt(8)
display(expr.evalf(4)) # to 4 digits of precisions

SymPy Numeric Evaluation

AoC Examples

Boat Racing Quadratic

In this example taken from 2023 Day 6, I’m solving a quadratic:

\[h^2 - th + d = 0\]
def solve_part2_sympy_quadratic(data):
    """ h^2 - th + d = 0 """
    race_duration = int("".join(x for x in data[0].split(":")[1].split()))
    distance = int("".join(x for x in data[1].split(":")[1].split()))
    logger.debug(f"{race_duration=}, {distance=}")  

    # solve using quadratic with SymPy       
    h = sympy.symbols("h", real=True)
    equation = sympy.Eq(h**2 - race_duration*h + distance, 0)
    solutions = sympy.solve(equation, h, dict=True)
    answers = [solution[h].evalf() for solution in solutions] # there should be two

Trajectory Intersection

In this example taken from 2023 Day 24, I have a series of equations that are necessary to find several unknown variables. Specifically:

\[\begin{align} t &= \frac{x_{r} - x_{h}}{v_{x_{h}} - v_{x_{r}}} = \frac{y_{r} - y_{h}}{v_{y_{h}} - v_{y_{r}}} = \frac{z_{r} - z_{h}}{v_{z_{h}} - v_{z_{r}}} \\ \notag \\ (x_{r} - x_{h})(v_{y_{h}} - v_{y_{r}}) &= (y_{r} - y_{h})(v_{x_{h}} - v_{x_{r}}) \\ (y_{r} - y_{h})(v_{z_{h}} - v_{z_{r}}) &= (z_{r} - z_{h})(v_{y_{h}} - v_{y_{r}}) \\ (z_{r} - z_{h})(v_{x_{h}} - v_{x_{r}}) &= (x_{r} - x_{h})(v_{z_{h}} - v_{z_{r}}) \\ \end{align}\]

And here’s the code:

def solve_part2(data: list[str]):
    Determine the sum of the rock's (x,y,z) coordinate at t=0, for a rock that will hit every hailstone
    in our input data. The rock has constant velocity and is not affected by collisions.
    stones = parse_stones(data)
    logger.debug(f"We have {len(stones)} stones.")
    # define SymPy rock symbols - these are our unknowns representing:
    # initial rock location (xr, yr, zr)
    # rock velocity (vxr, vyr, vzr)
    xr, yr, zr, vxr, vyr, vzr = sympy.symbols("xr yr zr vxr vyr vzr")
    equations = [] # we assemble a set of equations that must be true
    for stone in stones[:10]: # we don't need ALL the stones to find a solution. We need just enough.
        x, y, z = stone.posn
        vx, vy, vz = stone.velocity
        equations.append(sympy.Eq((xr-x)*(vy-vyr), (yr-y)*(vx-vxr)))
        equations.append(sympy.Eq((yr-y)*(vz-vzr), (zr-z)*(vy-vyr)))

    solutions = sympy.solve(equations, dict=True) # SymPy does the hard work
    if solutions:
        solution = solutions[0]
        return sum([solution[xr], solution[yr], solution[zr]])
    logger.info("No solutions found.")
    return None