1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:

from math import pi, exp, sqrt

def sgn(x):
    return -1 if x<0 else 0 if x==0 else 1

def bisection(f,x,a,b):
    for k in range(60):
        m = 0.5*(a+b);
        if f(m)-x<0: a = m
        else: b = m
    return m

def incomplete_Gamma_continued_fraction(a,x,n):
    y = 0
    for k in range(n,0,-1):
        y = k*(k-a)/(x-y-a+2*k+1)
    return exp(-x)*x**a/(x-y-a+1)

def incomplete_gamma_power_series(a,x,n):
    y = 0
    p = 1/a
    for k in range(1,n):
        y += p
        p = p*x/(a+k)
    return y*exp(-x)*x**a

def erf(x):
    if abs(x)>8: return sgn(x)
    if abs(x)>1.7:
        y = sqrt(pi)-incomplete_Gamma_continued_fraction(0.5,x*x,26)
    else:
        y = incomplete_gamma_power_series(0.5,x*x,26)
    return sgn(x)*y/sqrt(pi)

def cdf_normal(x,mu=0,sigma=1):
    return 0.5+0.5*erf((x-mu)/(sqrt(2)*sigma))

def cdf_normal_inv(x):
    return bisection(cdf_normal,x,-100,100)