I didn't peruse the source code. I just read the linked article in its entirety and it says

> The computation is quite slow. In order to stay as flexible as possible, I'm using the Monte Carlo method. Which means the calculator is running about 250K AST-based computations for every calculation you put forth.

So therefore I conclude Monte Carlo is being used.

It's dead simple. Here is the simplified version that returns the quantiles for '100 / 2 ~ 4'.

  import numpy as np
  
  def monte_carlo(formula, iterations=100000):
    res = [formula() for _ in range(iterations)]
    return np.percentile(res, [0, 2.5, \*range(10, 100, 10), 
    97.5, 100])

  def uncertain_division():
    return 100 / np.random.uniform(2, 4)

  monte_carlo(uncertain_division, iterations=100000)

Line 19 to 21 should be the Monte-Carlo sampling algorithm. The implementation is maybe a bit unintuitive but apparently he creates a function from the expression in the calculator, calling that function gives a random value from that function.