Реализация Андерсона Дарлинга на C++

Я уже много времени пытаюсь реализовать (ха-ха) Anderson-Darling проверка нормальности в C++. Вот мой код. Я знаю, что есть аналогичная тема здесь, но, к сожалению, это не решило мою проблему.

Дисперсия рассчитывается правильно, я бы предположил, что равномерное стандартизированное распределение тоже. Проблема в том, что b и c дают мне NAN для моих выборочных данных 1, 2....10.

Вы не знаете, где ошибка в формуле - Anderson_Darling(), см. код ниже?

Код вырезан из класса для большей ясности. Я не стал приводить здесь очевидные методы, такие как mean() и т.д.

double variance()
{
    double var_sum = 0.0;
    for(int i = 0; i < int(size); i++)            //(size is taken from a class)
        var_sum += pow(data.at(i)-mean(),2);
    return var_sum / (int(size)-1);
}

double phi(double x)
{
    double res =0.5 * erfc(-x * M_SQRT1_2);
    return res;
}

vector<double> tostdnormal()
{
    vector<double> Y (size);
    for(int i = 0; i < int(size); i++)
        Y.at(i) = (data.at(i) - mean())/(sqrt(variance()));
    return Y;
}

double Anderson_Darling()
{
    sort(data.begin(),data.end());
    int n = int(size);
    vector<double> Y = tostdnormal();

   double S = 0; double a = 0; double b = 0; double c = 0;
    for(int i = 0; i < n; i++)
    {
        a = 2.0 * (i+1) - 1;
        b = log(Y.at(i));
        c = log(1-Y.at(n-i-1));
        S += a * (b + c);
     }
return -n - S / n;


 }

Обновить. Я изменил b и c на это и получил ожидаемый результат.

b = log(phi(Y.at(i)));
c = log(1-phi(Y.at(n - i - 1)));