# Learning Python with Advent of Code Walkthroughs

Dazbo's Advent of Code solutions, written in Python

# The Python Journey - Maths with SymPy

## 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:

• Simplifying, rearranging, factoring and expanding algebraic expressions and equations.
• Solving equations of various types, including linear, polynomial, and differential.
• Performing calculus, including differentiation and integration.
• Plotting graphs.
• Generating LaTeX output.

## 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("Simplifying..."))
display(Markdown(f"$a-b$ $= {latex(sympy.simplify(a_minus_b))}$"))

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

a.equals(b)

### 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=}")

### Differentiation and Integration

import sympy

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

deriv = sympy.diff(expr)
print(deriv)

integral = sympy.integrate(deriv)
print(integral)

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)
solutions

for solution in solutions:
print(f"x={solution[x]}")

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:
display(expr)
solutions = sympy.solveset(expr, x)
print(f"There are {len(solutions)} solutions for x...")
display(solutions)

display(expressions[-1])
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)
display(expr.evalf(4)) # to 4 digits of precisions

## AoC Examples

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

$h^2 - th + d = 0$
""" 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]
logger.info(solution)
return sum([solution[xr], solution[yr], solution[zr]])

logger.info("No solutions found.")
return None