code: def qsnorm(p):
      
      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
docs: 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