def qsnorm(p):
    """
    rational approximation for x where q(x)=d, q being the cumulative
    normal distribution function. taken from Abramowitz & Stegun p. 933
    |error(x)| < 4.5*10**-4
    """
    d = p
    if d < 0. or d > 1.:
        print('d not in (1,1) ')
        sys.exit()
    x = 0.
    if (d - 0.5) > 0:
        d = 1. - d
    if (d - 0.5) < 0:
        t2 = -2. * np.log(d)
        t = np.sqrt(t2)
        x = t - old_div((2.515517 + .802853 * t + .010328 * t2),
                        (1. + 1.432788 * t + .189269 * t2 + .001308 * t * t2))
        if p < 0.5:
            x = -x
    return x