А> Функция erf() возвращает интеграл вероятности ошибки для
А> аргумента x, который математически определяется как

А> erf(x) = 2/sqrt(pi)* интеграл от 0 до x от exp(-t*t) dt

А> Функция erfc() вычисляет дополнение к интегралу
А> вероятности ошибки для аргумента x, который представляет
А> собой 1.0 — erf(x). Введена дополнительно для уменьшения
А> вычислительных погрешностей при больших значениях x, когда
А> значение erf() мало отличается от 1.0.
А>

K>double erf(double x)
K>{
K>    double y,z,xnum,xden,del,result;
K>    int i;
K>    const double xbreak = 0.46875;
K>    const double pi=3.14159265358979;


K>    y = fabs(x);

K>//  evaluate  erf  for  |x| <= 0.46875
K>    if (y<xbreak)
K>    {
K>        double a[] = {3.16112374387056560e00, 1.13864154151050156e02, 3.77485237685302021e02,
K>                      3.20937758913846947e03, 1.85777706184603153e-1};
K>        double b[] = {2.36012909523441209e01, 2.44024637934444173e02,  1.28261652607737228e03,
K>                      2.84423683343917062e03};

K>        z = y * y;
K>        xnum = a[4] * z;
K>        xden = z;
K>        for (i=0; i<3; i++)
K>        {
K>           xnum = (xnum + a[i]) * z;
K>           xden = (xden + b[i]) * z;
K>        }
K>        result= x * (xnum + a[3]) / (xden + b[3]);
K>    }
K>//  evaluate  erfc  for 0.46875 <= |x| <= 4.0
K>    else if ((y > xbreak) && (y <= 4.))
K>    {
K>        double c[] = {5.64188496988670089e-1, 8.88314979438837594e00, 6.61191906371416295e01,
K>                      2.98635138197400131e02, 8.81952221241769090e02, 1.71204761263407058e03,
K>                      2.05107837782607147e03, 1.23033935479799725e03, 2.15311535474403846e-8};
K>        double d[] = {1.57449261107098347e01, 1.17693950891312499e02, 5.37181101862009858e02,
K>                      1.62138957456669019e03, 3.29079923573345963e03, 4.36261909014324716e03,
K>                      3.43936767414372164e03, 1.23033935480374942e03};

K>        xnum = c[8] * y;
K>        xden = y;
K>        for (i=0; i<7; i++)
K>        {
K>            xnum = (xnum + c[i]) * y;
K>            xden = (xden + d[i]) * y;
K>        }
K>        result = (xnum + c[7]) / (xden + d[7]);
K>        z = ((int)(y*16))/16.0;
K>        del = (y-z)*(y+z);
K>        result *= exp(-z*z-del);
K>    }
K>//   evaluate  erfc  for |x| > 4.0
K>    else
K>    {
K>        double p[] = {3.05326634961232344e-1, 3.60344899949804439e-1, 1.25781726111229246e-1,
K>                      1.60837851487422766e-2, 6.58749161529837803e-4, 1.63153871373020978e-2};
K>        double q[] = {2.56852019228982242e00, 1.87295284992346047e00, 5.27905102951428412e-1,
K>                      6.05183413124413191e-2, 2.33520497626869185e-3};
K>        z = 1 / (y * y);
K>        xnum = p[5] * z;
K>        xden = z;
K>        for (i=0; i<4; i++)
K>        {
K>           xnum = (xnum + p[i]) * z;
K>           xden = (xden + q[i]) * z;
K>        }
K>        result = z * (xnum + p[4]) / (xden + q[4]);
K>        result = (1/sqrt(pi) -  result) / y;
K>        z = ((int)(y*16))/16.0;
K>        del = (y-z) * (y+z);
K>        result *= exp(-z*z-del);
K>    }

K>//   fix up for negative argument, erf, etc.
K>    if (x > xbreak)
K>       result = 1 - result;
K>    if (x < -xbreak)
K>       result = result - 1;

K>    return result;
K>}

K>