# X ~ N(30, 4**2)
mu,sig = 30,4
INF = float('inf')
cases = [(-INF, 40), (21,INF), (30,35)]
def p(f: float): print(f'{f:.3f}')
`

`
### Method 1: bruteforce approx
import random
def sample(): return random.gauss(mu,sig)
# you can also sample randomly by taking 2 outputs of random.random() && using the Box-Muller formula.
# if `random` itself is banned, you can try making a simple LCG from scratch, but it may be too biased.
def brute(a, b, tries: int=99999):
  return sum(a<sample()<b for _ in range(tries)) / tries

for args in cases: p(brute(*args))
`

`
### Method 2: \Phi(z) = (1+erf(z/\sqrt{2}))/2
import math
def cdf_N(z: float): return (1 + math.erf(z/math.sqrt(2))) / 2
def cdf_X(z: float): return cdf_N((z-mu)/sig)

p(cdf_X(40))
p(1-cdf_X(21))
p(cdf_X(35)-cdf_X(30))
`

`
### Method 3: Riemann summation
# pdf(t) = \frac{e^{-t^2/2}}{\sqrt{2\pi}}
def pdf(t: float): return math.exp(-t*t/2) / math.sqrt(2*math.pi)
def phi(t): # approximate the CDF
    t_abs = abs(t)
    delta, eps = 1e-3, 1e-6
    s = 0
    while (d:=delta*pdf(t_abs)) >= eps:
      s += d
      t_abs += delta
    return s if t <= 0 else 1-s # \Phi(t) = 1-\Phi(t) if t > 0
def cdf(z): return phi((z-mu)/sig)
p(cdf(40))
p(1-cdf(21))
p(cdf(35)-cdf(30)