Критерий Шапиро-Уилка (Исправить программу)

// Shapiro–Wilk.cpp: определяет точку входа для консольного приложения.
//
 
#include "stdafx.h"
#include <stdio.h>
#include <math.h>
#include <stdlib.h> 
#include <iostream>
#include "Rmath.h"
 
 
using namespace std;
 
 
#ifndef min
# define min(a, b)      ((a) > (b) ? (b) : (a))
#endif
 
 
 
static double poly(const double *, int, double);
 
 
 
 
 
int sign(int x) 
{
    if (x == 0)
        return 0;
    return x > 0 ? 1 : -1;
}
 
 
 
static void swilk(double *x, int n, double *w, double *pw, int *ifault)
{
 
    int nn2 = n / 2;
    double *a = new double [nn2+1];
 
 
//Рассчитывает W тест Шапиро-Wilk и его уровень значимости
 
    double small = 1e-19;
 
    //Полиномиальные коэффициенты
    double g[2] = { -2.273,.459 };
    double c1[6] = { 0.,.221157,-.147981,-2.07119, 4.434685, -2.706056 };
    double c2[6] = { 0.,.042981,-.293762,-1.752461,5.682633, -3.582633 };
    double c3[4] = { .544,-.39978,.025054,-6.714e-4 };
    double c4[4] = { 1.3822,-.77857,.062767,-.0020322 };
    double c5[4] = { -1.5861,-.31082,-.083751,.0038915 };
    double c6[3] = { -.4803,-.082676,.0030302 };
 
    //Локальные переменные
    int i, j, i1;
 
    double ssassx, summ2, ssumm2, gamma, range;
    double a1, a2, an, m, s, sa, xi, sx, xx, y, w1;
    double fac, asa, an25, ssa, sax, rsn, ssx, xsx;
 
    *pw = 1.;
    if (n < 3) 
    { 
        *ifault = 1; 
        return;
    }
 
    an = (double) n;
 
    if (n == 3) 
    {
    a[1] = 0.70710678;/* = sqrt(1/2) */
    } 
    else 
    {
    an25 = an + .25;
    summ2 = 0.;
    for (i = 1; i <= nn2; i++) 
    {
        a[i] = qnorm((i - 0.375) / an25, 0., 1., 1, 0);
        double r__1 = a[i];
        summ2 += r__1 * r__1;
    }
    summ2 *= 2.;
    ssumm2 = sqrt(summ2);
    rsn = 1. / sqrt(an);
    a1 = poly(c1, 6, rsn) - a[1] / ssumm2;
 
 
    if (n > 5) 
    {
        i1 = 3;
        a2 = -a[2] / ssumm2 + poly(c2, 6, rsn);
        fac = sqrt((summ2 - 2. * (a[1] * a[1]) - 2. * (a[2] * a[2]))
               / (1. - 2. * (a1 * a1) - 2. * (a2 * a2)));
        a[2] = a2;
    } 
    else 
    {
        i1 = 2;
        fac = sqrt((summ2 - 2. * (a[1] * a[1])) /
               ( 1.  - 2. * (a1 * a1)));
    }
    a[1] = a1;
    for (i = i1; i <= nn2; i++) a[i] /= - fac;
    }
 
//Проверrка для нулевого диапазона
 
    range = x[n - 1] - x[0];
    if (range < small) 
    {
        *ifault = 6; 
        return;
    
    }
 
//Проверка правильность порядка сортировки на полигоне - в пересчете X
 
    // *ifault = 7; <-- так как он изменяется ниже, в любом случае!
    *ifault = 0;
    xx = x[0] / range;
    sx = xx;
    sa = -a[1];
    for (i = 1, j = n - 1; i < n; j--) 
    {
    xi = x[i] / range;
    if (xx - xi > small) 
    {
        *ifault = 7;
    }
    sx += xi;
    i++;
    if (i != j) sa += sign(i - j) * a[min(i, j)];
    xx = xi;
    }
    if (n > 5000) 
        *ifault = 2;
 
//Рассчитать W статистики как квадрат корреляции между данными и коэффициентами
 
    sa /= n;
    sx /= n;
    ssa = ssx = sax = 0.;
    for (i = 0, j = n - 1; i < n; i++, j--) 
    {
    if (i != j) 
        asa = sign(i - j) * a[1 + min(i, j)] - sa; 
    else asa = -sa;
    xsx = x[i] / range - sx;
    ssa += asa * asa;
    ssx += xsx * xsx;
    sax += asa * xsx;
    }
 
// W1 равно (1-W) рассчитывается, чтобы избежать чрезмерного ошибку округления для W самом ближайшем 1 (потенциальная проблема в очень больших выборках)
 
    ssassx = sqrt(ssa * ssx);
    w1 = (ssassx - sax) * (ssassx + sax) / (ssa * ssx);
    *w = 1. - w1;
 
//Расчет уровня значимости для W
 
    if (n == 3) 
    {
    //Точное значение P:
    double pi6 = 1.90985931710274, /* = 6/pi */
        stqr = 1.04719755119660; /* = asin(sqrt(3/4)) */
    *pw = pi6 * (asin(sqrt(*w)) - stqr);
    if(*pw < 0.) 
        *pw = 0.;
    return;
    }
    y = log(w1);
    xx = log(an);
    if (n <= 11) 
    {
    gamma = poly(g, 2, an);
    if (y >= gamma) 
    {
        *pw = 1e-99; //Очевидно значение 'маленькое', 1е-19f
        return;
    }
    y = -log(gamma - y);
    m = poly(c3, 4, an);
    s = exp(poly(c4, 4, an));
    }
    else 
    {// n >= 12
    m = poly(c5, 4, xx);
    s = exp(poly(c6, 3, xx));
    }
 
    *pw = pnorm(y, m, s, 0, 0);
 
    return;
} 
 
static double poly(const double *cc, int nord, double x)
{
 
/*
Рассчитывает алгебраический многочлен порядка nord-1 с 
Массивом коэффициентов  cc. Нулевой коэффициент cc(1) = cc[0]
*/
    double p, ret_val;
 
    ret_val = cc[0];
    if (nord > 1)
    {
    p = x * cc[nord-1];
    for (int j = nord - 2; j > 0; j--) p = (p + cc[j]) * x;
    ret_val += p;
    }
    return ret_val;
} 
 
//Shapiro-Wilk
 
int main ()
{
    swilk(x, n, w, pw, ifault);
}