Metoda interpolacji Langrange'a

#include<iostream>
#include<vector>
using namespace std;
struct P
{
    double x, y;
    P( double a, double b )
        : x( a )
        , y( b )
    { }
    P() { };
};
void druk( vector < double > wsp )
{
    unsigned r = wsp.size() - 1;
    for( unsigned i = 0; i <= r; i++ )
    {
        if( wsp[ i ] == 0 )
            ;
        else if( i != r )
        {
            if( wsp[ i ] == 1.0 )
            {
                if( r - i == 1 )
                     cout << "x";
                else
                     cout << "x^" << r - i;
               
            }
            else
            {
                if( r - i == 1 )
                     cout << wsp[ i ] << "x";
                else
                     cout << wsp[ i ] << "x^" << r - i;
               
            }
            //                        cout<<endl;
        }
        else
             cout << wsp[ i ];
       
    }
}
vector < double > funk( vector < double > wsp, P nawias )
{
    //        cout<<showpos<<"Mnozenie przez nawias (x"<<nawias.x<<")"<<endl<<endl;
    int stopien = wsp.size() - 1;
    vector < double > wynik, temp;
    wynik = wsp;
    temp.push_back( 0 );
    for( int i = 0; i <= stopien; i++ )
         temp.push_back( nawias.x * wynik[ i ] );
   
    wynik.insert( wynik.end(), 0 );
    stopien = wynik.size();
    for( int i = 1; i <= stopien; i++ )
         wynik[ i ] = wynik[ i ] + temp[ i ];
   
    return wynik;
}
vector < double > iloczyn( vector < double > wsp )
{
    vector < double > wynik;
    P nawias( wsp[ 3 ], 0 );
    wynik.push_back( wsp[ 0 ] ); wynik.push_back( wsp[ 1 ] );
    wynik = funk( wynik, nawias );
    //        druk(wynik);
    if( wsp.size() == 3 )
         return wynik;
   
    for( unsigned i = 5; i < wsp.size(); i += 2 )
    {
        nawias.x = wsp[ i ];
        wynik = funk( wynik, nawias );
        //                druk(wynik);
    }
    return wynik;
}
vector < double > Lagrange( vector < P > wezly )
{
    vector < double > zapis, funk, temp;
    int stopien = wezly.size() - 1;
    zapis.resize( stopien + 1, 0 );
    double mianownik, licznik;
    /*        for(int i=0;i<=stopien;i++)
            {
    */
    for( int j = 0; j <= stopien; j++ )
    {
        licznik = wezly[ j ].y;
        cout << licznik << "*";
        mianownik = 1;
        cout.setf( ios::showpos );
        for( int k = 0; k <= stopien; k++ )
        {
            if( k == j )
                 continue;
           
            mianownik = mianownik *( wezly[ j ].x - wezly[ k ].x );
            temp.push_back( 1 );
            temp.push_back( - wezly[ k ].x );
            cout << "(x" <<- wezly[ k ].x << ")";
        }
        cout.unsetf( ios::showpos );
        cout << "/" << mianownik << endl;
        //                        cout.setf(ios::showpos);
        //                        cout.unsetf(ios::showpos);
        temp = iloczyn( temp );
        for( unsigned i = 0; i < temp.size(); i++ )
        {
            temp[ i ] = temp[ i ] * licznik / mianownik;
            zapis[ i ] = zapis[ i ] + temp[ i ];
        }
        //                        druk(temp);
        temp.clear();
    }
    //        }
    return zapis;
}
int main()
{
    vector < P > wezly;
    vector < double > wsp, pom;
    wezly.push_back( P( - 2, 46 ) );
    wezly.push_back( P( - 1, 4 ) );
    wezly.push_back( P( 3, 156 ) );
    wezly.push_back( P( 4, 484 ) );
    pom = Lagrange( wezly );
    cout << showpos;
    cout << endl;
    druk( pom );
    cout << endl << endl;
    system( "PAUSE" );
}