Метод наименьших квадратов

#ifndef STRUCT_FLOATPOINT
#define STRUCT_FLOATPOINT
struct FLOATPOINT // the basic data type for the plot.
{  
   float    x;
   float    y;
};
#endif

void Regress(const FLOATPOINT* const pRange,DWORD dwSize,FLOATPOINT& Params)
{  
   double   SumX  = 0.0;
   double   SumY  = 0.0;
   double   SumX2 = 0.0;
   double   SumXY = 0.0;
 
   double   NN = (double)dwSize;
 
   for (DWORD ii = 0; ii < dwSize; ++ii)
   {
      double   U  = pRange[ii].x;
      double   V  = pRange[ii].y;
 
      SumX  += U;
      SumX2 += U * U;
 
      SumY  += V;
 
      SumXY += U * V;
   }
 
   double   C1 = (NN * SumXY - SumX * SumY) / (NN * SumX2 - SumX * SumX);
   double   C0 = (SumY - C1 * SumX) / NN;
 
   // LS approximation y = C1x + C0
   Params.x = (float)C1;
   Params.y = (float)C0;
}

template <class PointPtr> void createByOLSY(PointPtr p, const PointPtr end)
    {
        int n = end - p;
        ASSERT(n > 1);
        double sx = 0, sy = 0, sx2 = 0, sxy = 0;
        for (; p < end; p++)
        {
            sx += p->x;
            sy += p->y;
            sx2 += p->x*(double)p->x;
            sxy += p->x*(double)p->y;
        }
        double a0 = (sx*sx - n*sx2);
        if (a0 != 0)
        {
            double a1 = (sx*sy - n*sxy) / a0;
            a0 = (sx*sxy - sx2*sy) / a0;        //y = a1*x+a0, approximate OLS for y
            b = 1.0 / ((a0 < 0) ? sqrt(a1*a1 + 1) : -sqrt(a1*a1 + 1));
            a = a1*b;
            c = a0*b;
            b = -b;
        }
        else
        {
            a = (sx < 0) ? -1.0 : 1.0;
            b = 0.0;
            c = -a*sx / n;
        }
        ASSERT(fabs((a*sx + b*sy) / n + c) < o_z);
    }

    template <class PointPtr> void createByOLS(PointPtr p, const PointPtr end)
    {
        createByOLSY(p, end);
        int n = end - p;
        ASSERT(n > 1);
        double sx = 0, sy = 0, sx2 = 0, sy2 = 0, sxy = 0;
        for (; p < end; p++)
        {
            sx += p->x;
            sy += p->y;
            sx2 += p->x*(double)p->x;
            sy2 += p->y*(double)p->y;
            sxy += p->x*(double)p->y;
        }
        int n_2 = n*n;
        double sx_2 = sx*sx;
        double sy_2 = sy*sy;
        double sx_4 = sx_2*sx_2;
        double sy_4 = sy_2*sy_2;
        double sx2_2 = sx2*sx2;
        double sy2_2 = sy2*sy2;
        double sxy_2 = sxy*sxy;
 
        double sx_2sy_2 = sx_2*sy_2;
        double sx_2sx2 = sx_2*sx2;
        double m4sxsxysy = 4*sx*sxy*sy;
        double sx_2sy2 = sx_2*sy2;
        double sx2sy_2 = sx2*sy_2;
        double sy_2sy2 = sy_2*sy2;
        double sx2sy2 = sx2*sy2;
        
        //double zn = sx_4 + 2 * sx_2sy_2 + sy_4 + 2 * n * (sx_2sy2 - sx_2sx2 - m4sxsxysy + sx2sy_2 - sy_2sy2) + n_2*(sx2_2 - 2 * sx2sy2 + sy2_2 + 4 * sxy_2);
        double zn = sx_4 - 2*n*sx_2*sx2 + n_2*sx2_2 + 4*n_2*sxy_2 - 8*n*sx*sxy*sy + 2*sx_2*sy_2 + 2*n*sx2*sy_2 + sy_4 + 2*n*sx_2*sy2 - 2*n_2*sx2*sy2 - 2*n*sy_2*sy2 + n_2*sy2_2;
//      b = sqrt(((sx_4 + sy_4 + n_2*(sx2_2 + sy2_2)) / 2 + sx_2sy_2 + n*(sx2sy_2 - sx_2sx2 - m4sxsxysy + sx_2sy2 - sy_2sy2) + n_2*(2 * sxy_2 - sx2sy2)) / zn \
//          + (n*(sx2 - sy2) - sx_2 + sy_2) / (2 * sqrt(zn)));
        double sb0 = sx_4 / (2 * (zn)) - (n*sx_2*sx2) / (zn)+(n_2*sx2_2) / (2 * (zn)) + (2 * n_2*sxy_2) / (zn)-(4 * n*sx*sxy*sy) / (zn)+(sx_2*sy_2) / (zn)+(n*sx2*sy_2) / (zn)+sy_4 / (2 * (zn)) + (n*sx_2*sy2) / (zn)-(n_2*sx2*sy2) / (zn)-(n*sy_2*sy2) / (zn)+(n_2*sy2_2) / (2 * (zn));
        double sb1 = -sx_2 / (2 * sqrt(zn)) + (n*sx2) / (2 * sqrt(zn)) + sy_2 / (2 * sqrt(zn)) - (n*sy2) / (2 * sqrt(zn));
        double sb = 1.0, sa = 1.0;
        if (b != 0)
        {
            double a0 = (sx_2 - n*sx2);
            double a1 = (sx*sy - n*sxy) / a0;
            a0 = (sx*sxy - sx2*sy) / a0;        //y = a1*x+a0, approximate OLS for y
            if (a0 < 0)
            {
                sb = -1.0;
                if (a1 < 0)
                    sa = -1.0;
            }
            else if (a1 >= 0)
                sa = -1.0;
        }
        else if (sx < 0) 
            sa = -1.0;
        
        b = sb*sqrt(sb0+sb1);
        a = sa*sqrt(1 - (sb0 + sb1));
        c = -(a*sx + b*sy) / n;
 
    }