Hola, necesito que alguien me eche una ayuda con un algortimo matematico consistente en calcular los coeficientes de un polinomio mediante las identidades de newton. Por si sirve de ayuda, pondré aquí dichas identidades:
Página 1:
Página 2:
Página 3:
Aquí paso la función de código que tengo (lógicamente esta mal), donde se le pasa el vector dd que tiene las C0, c1, c2 .. , cn y el vector p que tiene los x0, x1, x2, ... xn de las imágenes que he puesto.
vector<long double> calculaCoeficientes(vector<long double>& p, vector<long double> & dd)
{
vector<vector<long double> > mat(p.size(),vector<long double> (p.size(), 0)) ;
long double suma = 0;
vector<long double> t;
vector<long double> b;
for (int i = 0; i< p.size(); i++)
{
t.push_back(1);
if (i != 0)
for (int k=1; k<=i; k++)
{
for (int j=0; j<i; j++)
{ suma+= pow(p[j], k); }
t.push_back(suma);
suma = 0;
}
//Aqui tenemos calculadas las Ti de las identidades de newton
b = calculaB (t);
for (int u=0; u< b.size(); u++)
{if (u%2 == 0) b=b;
else b*= (-1); }
//Aquí tenemos las Bi de las identidades de Newton con su correspondiete signo, es decir, lo que en las transparencias se llama an-i
for (int y=0; y < mat[0].size(); y++)
for (int x=0, h=b.size()-1; x<=y && h>=0; x++, h--)
mat[x][y] = b[h];
t.resize(0);
b.resize(0) ;
}
calculaSistema (dd, mat);
vector<long double> a;
//ahora ya tenemos cada fila multiplicada por las diferencia divididas, así que nos queda sumarlas para obtener cada coeficiente;
for (int i=0; i<mat.size(); i++)
{
long double suma=0;
for (int j=0; j<mat.size(); j++)
suma += mat[j];
a.push_back(suma);
}
// Esto calcula cada Ai del polinomio en potencias de X, resolviendo el sistema triangular de las transparencias, empezando por An
return a; // Vector con los coeficientes
}
vector<long double> calculaB (vector<long double> &te)
{
vector<long double> aux(te.size());
aux[0] = te[0];
aux[1] = te[1];
// Metemos b0=t0=1 y b1=t1
long double x=0;
long double y=0;
long double z=0;
if (te.size() != 1)
for (int i=2; i< te.size(); i++)
{
for (int j=1; j < i; j++)
y += pow(static_cast<long double> (-1), static_cast<long double> (j+1)) * (aux[j]*te[i-j]);
x += te-y;
z += (pow(static_cast<long double> (-1),static_cast<long double> (i+1)) / i) * x;
aux= z;
x=y=z=0;
}
return aux;
}
void calculaSistema (const vector<long double> &did, vector<vector<long double> > &matrix)
{
for (int i=0; i< matrix[0].size() ; i++)
for (int j=0; j< matrix.size() ; j++)
matrix[j] *= did;
for (int i=0; i<matrix.size() ; i++)
{
for (int j=0; j<matrix.size() ; i++)
cout << matrix[j]<< " ";
cout << endl;
}
}