function return_y=NORMDISTHULL(random_var)

%This is a polynomial approximation that gives six-decimal-place accuracy
%Ref. Hull, Options, Futures ans other derivatives, fourth edition, p.252
    gamma=0.23163419;
    a=[0.31938153 -0.356563782 1.781477937 -1.821255978 1.330274429];
    k=@(x) 1.0./(1+gamma*x);
    N_prime=@(x) exp(-.5*x.*x)/sqrt(2*pi);
    N=@(x, k_x) 1-N_prime(x).*(k_x.*(a(1)+k_x.*(a(2)+k_x.*(a(3)+k_x.*(a(4)+k_x*a(5))))));
    if random_var >= 0
        return_y=N(random_var, k(random_var));
    else
        return_y=1-N(-random_var, k(-random_var));
    end