Happy (belated) π Day!

Happy belated π Day! π Day is the 14th of March, i.e. 03-14 both the U.S. or ISO 8601 date formats. Two times in the past, I have posted some kind of method (not necessarily efficient) to approximate π as a π Day celebration:

• 2021: Used a Monte Carlo method: randomly selecting points \((x, y) \in [0, 1] \times [0, 1]\), and seeing what proportion have \(x^2 + y^2 < 1\). Implemented in JavaScript.
• 2024: Used a power series approximation and an argument from Papa Rudin¹. Implemented in Julia

This year's entry, admittedly a little late, involves WebAssembly and differential equations.

A system of differential equations

This year, I decided to estimate π by approximating \(\sin t\) for values of \(t > 0\), using differential equations and Euler's method. Specifically, note that \(x = \sin t, y = \cos t\) is the only solution to the system of differential equations \[\begin{aligned} \frac{dx}{dt} &= y, \\ \frac{dy}{dt} &= -x \end{aligned}\] that satisfies the initial conditions \(x(0) = 0,\; y(0) = 1\).

Euler's method

Suppose we have a system of differential equations of the form \[\begin{aligned} \frac{dx}{dt} &= a x + b y,\\ \frac{dy}{dt} &= c x + d y. \end{aligned}\] We can use Euler's method to find numerical solutions (that is to say, approximations of the solutions) to this system. For any real number \(t_0\), and small real number \(\Delta t\), we can estimate the quantities \[\begin{aligned} \Delta x &= x(t_0 + \Delta t) - x(t_0),\\ \Delta y &= y(t_0 + \Delta t) - y(t_0), \end{aligned}\] with the approximations \[\begin{aligned} \Delta x &\approx \left(a x(t_0) + b y(t_0)\right) \Delta t,\\ \Delta y &\approx \left(c x(t_0) + d y(t_0)\right) \Delta t. \end{aligned}\]

Given initial values \(x(0) = x_0\) and \(y(0) = y_0\), we can come up with approximations \(x_k\) and \(y_k\) for \(x(k \Delta t)\) and \(y(k \Delta t)\), respectively, by recursively defining \[\begin{aligned} x_{k + 1} &= x_k + \left(a x_k + b y_k\right) \Delta t,\\ y_{k + 1} &= y_k + \left(c x_k + d y_k\right) \Delta t. \end{aligned}\]

Enter WASM

Lately, I have been learning WebAssembly (WASM), and thought it would be fun to implement this year's π Day demo in WASM. While traditional assembly language is used to specify processor instructions in order to program a computer directly, WebAssembly instructions are run on a stack-based virtual machine inside the web browser.² WASM in created by writing the instructions in a human-readable text format, WebAssembly text or WAT.

For example, here is some WAT code for my WASM implementation of Euler's method:

        ;; set delta_x = (a*x + b*y) * delta_t
        local.get $a
        local.get $x
        f64.mul
        local.get $b
        local.get $y
        f64.mul
        f64.add
        local.tee $dxdt
        local.get $delta_t
        f64.mul
        local.set $delta_x

Lines starting with ;; are comments. The virtual machine that executes WASM is a stack machine. In the above snippet³:

The overall effect of lines 44-54 is to:

• Compute $a × $x + $b × $y and store in $dxdt
• Compute $dxdt × $delta_t and store it in $delta_x.

This illustrates the downside to WASM: it takes 11 lines of code to accomplish what most programming languages can do in two.

Euler step demo

The total function to do one step of Euler's method takes about 40 lines of WAT code. This JavaScript on this page is able to access the Euler step function. To the browser, this function is called eulerStep. It takes 7 numbers as arguments (\(a, b, c, d, x(t), y(t),\) and \(\Delta T\)) and returns four numbers (estimates for \(x(t), y(t), x'(t),\) and \(y'(t)\). You can play with this function by opening your browser's JavaScript console and evaluating the function there, once the WASM is done loading:

eulerStep(1, 2, 3, 4, 1, 1, 0.0625)

[1.1875, 1.4375, 3, 7]

On the other hand, if mucking around in your browser's JavaScript console isn't your thing, you can also access eulerStep via the following demo:

Demo is not ready yet. Loading WASM...

Approximating π

One definition of π is that \(t = \pi\) smallest positive real solution to the equation \(\sin t = 0\). Earlier I mentioned that \(x(t) = \sin t, y(t) = \cos t\) is the unique solution to the system of differential equations \[\begin{aligned} \frac{dx}{dt} & = y, \\ \frac{dy}{dt} &= -x \end{aligned}\] that satisfies the initial conditions \(x(0) = 0, y(0) = 1\). My plan is to use Euler's method to approximate a solution to this system (with initial conditions), which means approximating \(\sin t\) and \(\cos t\), and looking for the first positive value of \(t\) for which \(\cos t \leq 0\), which will be our approximation of π. The smaller \(\Delta t\) used, the better approximation of π we get.

Once the WASM has loaded successfully, click on the Run button to come up with different estimates for π using smaller and smaller values of \(\Delta t\).

Our estimates are produced as IEEE 754 double-precision floating-point numbers. This 64-bit format represents numbers⁴ in the form \[(-1)^s \times \left(1 + f/2^{52}\right)\times 2^{e - 1023},\] using one bit for \(s\), 11 bits for \(e\), and 52 bits for \(f\). For example, the number 3.25 is represented as \[\begin{aligned} (-1)^0 \times \left(1 + 2814749767106560/2^{52}\right) \times 2^{1024 - 1023} &= \left(1 + \frac{10 \times 2^{48}}{2^{52}} \right) \times 2 \\ &= \left(1 + 10/16\right) \times 2 \\ &= 3.25. \end{aligned}\]

For each approximation produced, click on the Show me the IEEE! to see the approximation as it's represented internally.

Loading the WASM...

Run

• Standalone page with the interactive portions of this post

• Source code