Библиотека Eigen

VectorXd fit(MatrixXd data, VectorXd X)
{
    int
        points = data.cols(),
        nodes = points - 4,     // conection nodes
        neq = 4*(points-3);     // number equations (for piecewise function)
    MatrixXd 
        A(neq,neq);
    VectorXd
        x(points), y(points);
// get data
    x << data.col(0);
    y << data.col(1);
 
    VectorXd
        B(neq), coeffs(neq);
// initialization
    A.fill(0), B.fill(0);
 
    
 
// interpolation by function (line, parabola or hyperbole)
    if (points > 1 && points <= 4)
    {
        for(int i = 0; i < points; i++)
            for(int j = 0; j < points; j++)
                A(i,j) = pow(x(i), points-j-1);
        B = y;
    }
        
    
 
 
// = = = = C U B I C  S P L I N E = = = =
 
    if (points > 4)
    {
        // the first three equations
        for (int i = 0; i < 3; i++)
        {
            for (int j = 0; j < 4; j++)
                A(i,j) = pow(x(i), 3-j);
            B(i) = y(i);
        }
            
        // the last three equations
        for (int i = 1; i <= 3; i++)
        {
            for (int j = 0; j < 4; j++)
                A(neq-i,neq-4+j) = pow(x(points-i), 3-j);
            B(neq-i) = y(points-i);
        }
 
        // leftover equations
        int
            shift = 4,
            index = 2;
        for(int i = 5; i < neq-5; i+=4)
        {
            for(int j = 0; j < 4; j++)
            {
                A(i,shift+j) = pow(x(index), 3-j);
                A(i+1,shift+j) = pow(x(index+1), 3-j);
            }
            B(i) = y(index);
            B(i+1) = y(index+1);
            index ++;
            shift += 4;
        }
 
        // derivations 
        int 
            row = 3, 
            col = 0;
        for (int i = 2; i < nodes+2; i++)
        {
            A(row,col) = 3*pow(x(i),2);
            A(row,col+1) = 2*x(i);
            A(row,col+2) = 1;
 
            A(row+1,col) = 6*x(i);
            A(row+1,col+1) = 2;
            for (int j = col+4; j < col+7; j++)
            {
                A(row,j) = -A(row,j-4);
                A(row+1,j) = -A(row+1,j-4);
            }
            row += 4;
            col += 4;
        }
    }
    // returning the contants for the piecewise fucntion
    coeffs = A.lu().solve(B);
 
    return coeffs;
}