double precision function cdfn(x)
        implicit real*8(a-h,o-z)
        y = x
        if (x.gt.20.0d0) y=20
        if (x.lt.-20.0d0) y = -20
        call cdfn2(y,a,b)
        if (a.lt.0.0000001d0) a = 0.0000001d0
        if (a.gt.0.9999999d0) a = 0.9999999d0
        cdfn = a
        return
        end


c***************************************************************
c
c   Subroutine to evaluate the cdf and pdf of the standard
c   normal distribution at the single point x using A&S, 26.2.17
c
c   x (input), cdf (output), and pdf (output) are double precision
c
c*****************************************************************
c
        subroutine cdfn2(x,cdf,pdf)
        implicit real*8(a-h,o-z)
        double precision x,one,pdf,t,t2,cdf,zpdf
        data one/1.0d0/
c
        pdf=zpdf(x)
        t=one/(one+.2316419d0*dabs(x))
        t2=t*t
        cdf=one-pdf*(.319381530*t-.356563782*t2+1.781477937*t*t2-
     1      1.821255978*t2*t2+1.330274429*t*t2*t2)
        if(x.lt.0.0d0) cdf=one-cdf
c
        return
        end

c******************************************************************
c
c   Function to evaluate the standard normal pdf fz at the
c   value z
c
c   z is double precision
c
c*******************************************************************
        double precision function zpdf(z)
        implicit real*8(a-h,o-z)
        double precision z
c
        zpdf=dexp(-z*z/2.0d0)/dsqrt(8.0d0*datan(1.0d0))
c
        return
        end