The unsure calculator inspired by https://filiph.github.io/unsure/

unsure_calc.py

import numpy as np
import re

NUMBER = re.compile(r'[0-9]+(\.[0-9]+)?')

def tokenize(s):
    orig_str = s
    while len(s) > 0:
        m = NUMBER.match(s)
        if m:
            yield np.float64(float(m.group()))
            s = s[m.end():]
        elif s[0] in '+-*/~()':
            yield s[0]
            s = s[1:]
        elif s[0] in ' \t\n':
            s = s[1:]
        else:
            raise ValueError(f"Wrong syntax in expression {repr(orig_str)}: cannot parse {repr(s)}")

def shunting_yard(tokens):
    stack = []
    for tok in tokens:
        if isinstance(tok, np.float64):
            yield tok
        elif tok == '(':
            stack.append(tok)
        elif tok == ')':
            while len(stack) > 0 and stack[-1] != '(':
                yield stack.pop()
            if len(stack) == 0:
                raise ValueError("Found right paren without left paren")
            stack.pop()
        elif tok in '+-*/~':
            while len(stack) > 0 and stack[-1] != '(' and (
                    stack[-1] == '~' or
                    stack[-1] in '*/' and tok in '+-*/' or
                    stack[-1] in '+-' and tok in '+-'):
                yield stack.pop()
            stack.append(tok)
    while len(stack) > 0:
        yield stack.pop()


def eval_rpn(ops):
    stack = []
    rng = np.random.default_rng()
    for op in ops:
        if isinstance(op, np.float64):
            stack.append(op)
        elif op == '+':
            a, b = stack.pop(), stack.pop()
            stack.append(b + a)
        elif op == '-':
            a, b = stack.pop(), stack.pop()
            stack.append(b - a)
        elif op == '*':
            a, b = stack.pop(), stack.pop()
            stack.append(b * a)
        elif op == '/':
            a, b = stack.pop(), stack.pop()
            stack.append(b / a)
        elif op == '~':
            a, b = stack.pop(), stack.pop()
            assert isinstance(a, np.float64) and isinstance(b, np.float64)
            mean = (a + b) / 2.0
            s = np.abs(b - a) / 3.28970725
            stack.append(rng.normal(mean, s, size=1_000_000))
    if len(stack) == 0:
        return
    assert len(stack) == 1
    if isinstance(stack[0], np.ndarray):
        res = np.quantile(stack[0], [0.05, 0.95])
        print("%.6g~%.6g" % (res[0], res[1]))
    else:
        print("%.6g" % stack[0])

if __name__ == "__main__":
    eval_rpn(shunting_yard(tokenize("100~200")))
    eval_rpn(shunting_yard(tokenize("100~200*2~4")))
    eval_rpn(shunting_yard(tokenize("100~200*2~4/2~4")))