Precision and rounding errors

[Tarbal]

I wrote some code that will fit data from an ascii file to a polynomial. Ultimately, I will add regression (I've begun, and it is commented out in the below code). This is all to be incorporated into a much larger scientific program later. For now, my problem is this: my polynomial coefficients which are the current output do not match what I've calculated via matrices in wxMaxima. How can I increase accuracy and/or reduce rounding error? I'm still pretty new to c++.

Program output:

beta
231.165
-438.513
274.387
-75.3077
9.73383
-0.464325

wxMaxima solution:

226.8797363317376
−430.0542771046293
268.8091043714023
−73.70060682268274
9.523984400796692
−0.45415109784472

Input file data.dat:

1	1
2.5	8
3	27
4	64
5.1	125
6	216
7	343

Clearly, this would fit best to a 3rd order polynomial, but right now, I'm testing precision, so I want my result to match that of Maxima. Here is the code:

#include <iostream>
#include <fstream>
#include <istream>
#include <vector>
#include <stdio.h>
#include <cmath>
#include <string>
#include <sstream>
#include "../headers/tnt_array1d.h"
#include "../headers/tnt_array2d.h"
#include "../headers/jama_lu.h"


using namespace TNT;
using namespace JAMA;
using namespace std;


int poly=5;
int num_rows;
Array2D<long double> TX(poly+1,num_rows);
Array1D<long double> beta(poly+1);
Array1D<long double> Xbeta()
{
        LU<long double> lu(TX);
        return lu.solve(beta);
}

int main()
{
//      cout.precision(8);
//      cout.setf(ios::fixed);
        vector <long double> Tsigma;
        vector <long double> sigma;
        num_rows = 0;
        ifstream fin("data.dat");
        vector<vector <long double> > allData;
        string line;
        while(getline(fin,line))
        {
                vector<long double> lineData;
                long double val;
                istringstream lineStream(line);
                while (lineStream >> val)
                {
                        lineData.push_back(val);
                }
                allData.push_back(lineData);
                ++num_rows;
        }


        for (int i=0; i<num_rows; i++)
        {
                Tsigma.push_back(allData[i][0]);
                sigma.push_back(allData[i][1]);
//              cout << Tsigma[i] << " " << sigma[i] << endl;
        }



        Array1D<long double> Y(num_rows,1);
        for (int i = 0; i < num_rows; i++)
        {
                Y[i] = sigma[i];
        }

        Array2D <long double> X(num_rows,poly+1);
        for (int i = 0; i < num_rows; i++)
        {
                for (int j = 0; j <= poly; j++)
                {
                        X[i][j] = pow(Tsigma[i],j);
                }
        }


        LU<long double> lu(X);
        Array1D<long double> beta(poly+1,1);
        beta = lu.solve(Y);

/*      cout << "test2" << endl;
        cout << num_rows << endl;
        Array2D<long double> TX(poly+1,num_rows);
        for (int i = 0; i < num_rows; i++)
        {
                for (int j = 0; j < poly+1; j++)
                {
                        TX[j][i] = X[i][j];
                }
        }
        cout << "TX" << endl;
        cout << TX << endl;

        cout << "test3" << endl;
        Array1D<long double> epsilon(poly+1);
        Array1D<long double> Xbeta_var = Xbeta();
        cout << "Xbeta" << endl;
        cout << Xbeta_var << endl;
        cout << "epsilon" << endl;
        for (int i = 0; i < num_rows; i++)
        {
        epsilon[i] = Y[i] - Xbeta_var[i];
        }
        cout << epsilon << endl;
*/
        cout << "beta" << endl;
        for (int i = 0; i < poly+1; i++)
        {
                cout << beta[i] << endl;
        }
        return 0;
}

Why doesn't my beta value match that of the calculated value in Maxima? How can I increase precision and/or reduce rounding errors? Any help would be greatly appreciated.

[Tarbal]

Using a library called Eigen seemed to do the trick. Here is the code that gives high accuracy, in case you ever want to fit an x, y data file to a polynomial. However, if your data file has any blank x or y values or you have a blank line (including at the end of the file), then the program will crash.

#include <iostream>
#include <fstream>
#include <istream>
#include <vector>
#include <stdio.h>
#include <cmath>
#include <string>
#include <sstream>
#include "../headers/Eigen/Dense"
#include "../headers/Eigen/src/LU/Inverse.h"

using namespace std;
using namespace Eigen;

int poly=6;
int num_rows;

int main()
{
        cout.precision(8);
        cout.setf(ios::fixed);
        vector <long double> Tsigma;
        vector <long double> sigma;
        num_rows = 0;
        ifstream fin("data.dat");
        vector<vector <long double> > allData;
        string line;
        while(getline(fin,line))
        {
                vector<long double> lineData;
                long double val;
                istringstream lineStream(line);
                while (lineStream >> val)
                {
                        lineData.push_back(val);
                }
                allData.push_back(lineData);
                ++num_rows;
        }

        for (int i=0; i<num_rows; i++)
        {
                Tsigma.push_back(allData[i][0]);
                sigma.push_back(allData[i][1]);
        }



        MatrixXd Y(num_rows,1);
        for (int i = 0; i < num_rows; i++)
        {
                Y(i) = sigma[i];
        }

        MatrixXd X(num_rows,poly+1);
        for (int i = 0; i < num_rows; i++)
        {
                for (int j = 0; j <= poly; j++)
                {
                        X(i,j) = pow(Tsigma[i],j);
                }
        }

        MatrixXd TXX(poly+1,num_rows);
        TXX = X.transpose()*X;

        MatrixXd TXY(poly+1,1);
        TXY = X.transpose()*Y;

        MatrixXd beta(poly+1,1);
        beta = TXX.inverse()*TXY;

        cout << "beta\n" << scientific;
        cout << beta << endl;
        return 0;
}