So first off:*** This code is not being used in production software.It is a personal project of mine, trying to understand approximation theory and advanced curve fitting.  In other words, I'm trying to understand how it works, not trying to get a currently existing solution.So I have been trying to implement the Remez algorithm for polynomial approximation.I sort of/maybe/kind of have it working (not really).My current solution generate ok polynomials, but a)  the coefficients are not converging&b)  while monitoring the coefficients at each stage, I've noticed that the x-values seem to slip past each other.I'll give some examples to show what I mean.My base function I'm trying to model is the square root function with a 4 degree polynomial on the domain [0.25, 1]Round 1X[0] = 0.25X[1] = 0.4X[2] = 0.55X[3] = 0.7X[4] = 0.85X[5] = 1Round 2X[0] = 0.25X[1] = 0.595076928220583X[2] = 0.493988453622788X[3] = 0.714333640596557X[4] = 1.14135676154991X[5] = 1Round 3X[0] = 0.25X[1] = 0.638393337463021X[2] = 0.63752199068821X[3] = 0.538600997945798X[4] = 1.07101841739164X[5] = 1Round 4X[0] = 0.25X[1] = 0.559423143598625X[2] = 0.560673538304964X[3] = 0.580378820375143X[4] = 1.04454592077508X[5] = 1So here's a look at my actual code.template <typename func_t>type_t MiniMax(size_t Degree, const type_t& LowerLimit, const type_t& UpperLimit, unsigned char Iterations, func_t F0, func_t F1, func_t F2){    // (C) Jacob Wells 2015    // This code is licensed under the BSD 3-Clause License    // http://opensource.org/licenses/BSD-3-Clause    if((Degree < 1) || (Iterations < 1))    {        return (type_t)NAN;    }    const type_t ONE(1), NEG1(-1), ZERO(0);    matrix_t<type_t> M(Degree + 2, Degree + 3);    vector<type_t> XVal;    type_t Delta, Sign, Pow, Err;    type_t D1, D2;    size_t I, J;    Delta = (UpperLimit - LowerLimit) / (Degree + 1);    XVal.resize(Degree + 2);    Coef.resize(Degree + 1);    Sign = ONE;    for(I = 0; I < (Degree + 2); I++)               // Generate our initial x-values    {        XVal[I] = (I * Delta) + LowerLimit;    }    do    {        Sign = NEG1;        for(I = 0; I < (Degree + 2); I++)        {            Pow = ONE;            M[I][Degree + 1] = Sign;            // Enters the alternating error sign            M[I][Degree + 2] = F0(XVal[I]);     // Enters the f(x) value            Sign *= NEG1;                                   for(J = 0; J <= Degree; J++)        // Evaluates the polynomial for each power            {                M[I][J] = Pow;                Pow *= XVal[I];             }        }        rref(Degree + 2, M);                    // Use Row Reduction Echelon Form to find the polynomial coefficients        Err = M[Degree + 1][Degree + 2];        for(I = 0; I <= Degree; I++)            // Copy the coefficients into the polynomial class' array        {            Coef[I] = M[I][Degree + 2];        }        if(Iterations > 1)        {            for(I = 1; I <= Degree; I++)        // Use Newton's method to find our new x-values            {                D1 = nth_deriv(XVal[I], 1) - F1(XVal[I]);                D2 = nth_deriv(XVal[I], 2) - F2(XVal[I]);                if((D1 != ZERO) && (D2 != ZERO))                {                    XVal[I] -= (D1 / D2);                           }            }        }    }while(--Iterations != 0);    return Err;} A quick little guide to some of my code:  F0, F1, & F2 are, respectively, the square root functions, it's first derivative, and it's second derivative.nth_deriv calculate the Nth Derivative of the current polynomial.rref reduces the Matrix to Row Reduction Echelon Formmatrix_t is a bare bones matrix class I came up with.Now I have done a lot of testing on these functions, and I haven't found a single error with them, so I feel very confident that the problem is in the MiniMax Function.EDIT: My barebones matrix classtemplate <typename type_t>class matrix_t{    public:    ~matrix_t() {}    matrix_t(size_t ColL, size_t RowL)    {        Arr.resize(ColL * RowL);        RowLen = RowL;    }    type_t* operator [] (size_t I)    {        return &Arr[I * RowLen];    }    const type_t* operator [] (size_t I) const    {        return &Arr[I * RowLen];    }    type_t At(size_t I)    {        return Arr[I];    }    private:    size_t RowLen;    vector<type_t> Arr;    matrix_t();    matrix_t(const matrix_t& m);    void operator = (const matrix_t& m);};
Problems Implementing the Remez Algorithm
regression;approximation algorithms;polynomials;curve fitting