Stand-alone Haskell code for standard normal CDF

The standard normal cumulative distribution function (CDF), often written Φ(x), is closely related to the error function erf(x) but is more common than the latter in statistical applications. See these notes on how the two functions are related.

The following code computes Φ(x) in pure Haskell with no external dependencies.

phi :: Double -> Double
phi x = y 
    where
        a1 =  0.254829592
        a2 = -0.284496736
        a3 =  1.421413741
        a4 = -1.453152027
        a5 =  1.061405429
        p  =  0.3275911

        -- Abramowitz and Stegun formula 7.1.26
        sign = if x > 0
                   then  1
                   else -1
        t = 1.0/(1.0 + p * abs x / sqrt 2.0)
        e = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-x*x/2.0)
        y = 0.5*(sign*e + 1.0)

test_phi :: Bool
test_phi = maximum [ abs(phi x - y) | (x, y) <- zip xs ys ] < epsilon
    where
        epsilon = 1.5e-7 -- accuracy promised by A&S
        xs = [-3, -1, 0.0, 0.5, 2.1 ]
        ys = [0.00134989803163, 
              0.158655253931, 
              0.5, 
              0.691462461274, 
              0.982135579437]