#define     zero  0.00
#define     MAXITER    20
#define     EPS        1.E-6
#define     MAXORDER    15


# include <iostream>  
using namespace std;

# include <math.h>	 


class Polinom 
{
public:
        int         degree ;    
        double      *coefs ;
		double      *nule ;

Polinom(int d=0)
{
    int  j ;     
    degree = d ;
    coefs = new double[degree + 1] ; 
     for( j=0 ; j <= degree ; j++ ) 
        coefs[j] = 0.0 ;
	 if ( degree==0 )
	    {	nule = new double[1] ; 
	        nule[0] = 0.0;
	    }
	 else
	    {  for( j=0 ; j < degree ; j++ ) 
		  nule[j] = 0.0 ;
        } 

}
Polinom(int d, double *niz)
{
    int  j ;     
    degree = d ;
    coefs = new double[degree + 1] ; 
	nule = new double[degree] ;
	
    for( j=0 ; j < degree ; j++ ) 
	{
        coefs[j] = niz[j] ;
		nule[j] = 0.0 ;
    }
	coefs[j] = niz[j] ;
      Bairstow();
} 

Polinom( const Polinom &orig )
    {
    int     j ;     
    degree = orig.degree ;
    coefs = new double[degree + 1] ; 
	nule = new double[degree] ; 
    
    for( j=0 ; j <degree ; j++ ) {
        coefs[j] = orig.coefs[j] ;
		nule[j] = orig.nule[j] ;
    }
	coefs[j] = orig.coefs[j] ;
}  
    
~Polinom( void )

{
    delete[] coefs ;
}   
        
       
Polinom operator+( const Polinom &right )
{
    int     j, max_degree ;

    if( (*this).degree > right.degree ) {
        max_degree = (*this).degree ;
    }
    else {
        max_degree = right.degree ;
    }
 
    Polinom  sum( max_degree ) ;

    for( j=0 ; j <= (*this).degree ; j++ ) 
        sum.coefs[j] += (*this).coefs[j] ;

    for( j=0 ; j <= right.degree ; j++ ) 
        sum.coefs[j] += right.coefs[j] ;
	sum.Bairstow();
    
    return( sum ) ;
}


Polinom operator-( const Polinom &right )
{
    int     j, max_degree ;

    if( (*this).degree > right.degree ) {
        max_degree = (*this).degree ;
    }
    else {
        max_degree = right.degree ;
    }

    Polinom  diff( max_degree ) ;

    for( j=0 ; j <= (*this).degree ; j++ ) 
        diff.coefs[j] += (*this).coefs[j] ;

    for( j=0 ; j <= right.degree ; j++ ) 
        diff.coefs[j] -= right.coefs[j] ;
	diff.Bairstow();
    
    return( diff ) ;
}


Polinom operator*( const Polinom &right )
{
    int     j, k ;

    Polinom  product( (*this).degree + right.degree ) ;
        
//****************************************************
//	Work out the result.
//
    for( j=0 ; j <= (*this).degree ; j++ ) {
        for( k=0 ; k <= right.degree ; k++ )
            product.coefs[j+k] += (*this).coefs[j] * right.coefs[k] ;
    }
	product.Bairstow();
    
    return( product ) ;
}

Polinom operator*( const double &db )
{
   Polinom  product( (*this).degree ) ;
   for(int j=0 ; j <= (*this).degree ; j++ ) 
            product.coefs[j] = (*this).coefs[j] * db ;
            
		 product.Bairstow();
    
    return( product ) ;
}
                
        
Polinom operator=( const Polinom &right ) 
{
    int     j ;
    
    delete[]  coefs ;
	delete[]  nule ;
    coefs = new double[right.degree + 1] ;
    nule = new double[right.degree] ;
	for( j=0 ; j < right.degree ; j++ )
	  {
        coefs[j] = right.coefs[j] ;	
	    nule[j] = right.nule[j] ;	
	  }
        coefs[j] = right.coefs[j] ;	
        degree = right.degree ;
    
    return( *this ) ;
}
        

double operator[]( int n )
{
    if( (n >= 0) && (n <= degree) ) return( coefs[n] ) ;
    else                            return( zero ) ;
}
        
 
double operator()( double x )
{
    int         j ;
    double      result ;

    result = 0.0 ;
    for( j=degree ; j >= 0 ; j-- )
	{
        result *= x ;
        result += coefs[j] ;
    }
    
    return( result ) ;
}     
         

friend ostream &operator <<( ostream &to, const Polinom p )
{
    int     j ;
    int     n_printed ; 
    
    n_printed = 0 ;
    for( j = p.degree ; j > 0 ; j-- ) {
    //*****************************************
    //	Print only the nonzero coefficients
        if( p.coefs[j] != 0.0 ) {
            if( n_printed != 0 ) {
            //*********************************
            //	This is not the leading tern:
            // print a + or a - sign.
            //
                if( p.coefs[j] > 0 )    
                    to << " + " ;
                else    // p.coefs[j] < 0 
                    to << " - " ;
            }
            else if( p.coefs[j] < 0 ) {
                to << " - " ;
            }
            
        //*************************************
        // 	Print the coefficient itself,
        // unless it's one. Even then, make
        // an exception for the constant term.
        //
            if( (fabs( p.coefs[j]) != 1.0) || (j == 0) ) 
                to << fabs( p.coefs[j] ) ;
                
        //**************************************
        //	Print an appropriate power of x.
        //
            if( j > 1 ) 		to << " x^" << j ;
            else if( j == 1 )   to << " x" ;
            
            n_printed++ ;
        }
    }
    
//*************************************************
//	Treat the constant term specially
//
    if( n_printed > 0 ) {
        if( p.coefs[0] > 0 )    
            to << " + " << p.coefs[0] ;
        else    // coefs[0] < 0 
            to << " - " << fabs( p.coefs[0] ) ;
    }
    else {  // n_printed == 0 
        to << p.coefs[0] ;
    }
    
    return( to ) ;
}
//*****************************************************************************
//**************************************************************************
void Bairstow(double p_init=0.0, double q_init=0.0)
{


    int         numanswer;
    double      a[MAXORDER],b[MAXORDER],c[MAXORDER];
    double      ximag[MAXORDER];
    int         i, k, n;
    double      p, q, d, dp, dq;
	double      det;

	n=degree;

    for(k=0; k<MAXORDER; k++)
	{ 
		a[k] = b[k] = c[k] = 0.;
		nule[k] = ximag[k] = 0.;
	}


    for(k=0; k<=n; k++) a[k] = coefs[n-k];

	for(k=1; k<=n; k++)	a[k] /= a[0];
	a[0] = 1.;

    p = p_init;
    q = q_init;
    numanswer = 0;

    while(n>2)
	{

        for(i=0; i<MAXITER; i++)
		{
            b[0] = 1.;
            b[1] = a[1] - p*b[0];
            for(k=2; k<=n; k++) b[k] = a[k] - p*b[k-1] - q*b[k-2];

            c[0] = 1.;
            c[1] = b[1] - p*c[0];
            for(k=2; k<n; k++) c[k] = b[k] - p*c[k-1] - q*c[k-2];

            det=c[n-2]*c[n-2] - c[n-3]*(c[n-1]-b[n-1]);

            dp = (b[n-1]*c[n-2]-b[n]*c[n-3])/det;
            dq = (b[n]*c[n-2]-b[n-1]*(c[n-1]-b[n-1]))/det;

            if(fabs(dp) < EPS && fabs(dq) < EPS) break;

            p += dp;
            q += dq;
        }//..for.i
//*****************
		
          d=p*p-4.*q;
    if(d>=0.)
	{
        nule[numanswer  ] = (-p+sqrt(d))/2.;
		ximag[numanswer++] = 0.;
        nule[numanswer  ] = (-p-sqrt(d))/2.;
		ximag[numanswer++] = 0.;
    }
    else
	{
        nule[numanswer  ] = -p/2.;
        ximag[numanswer++] = sqrt(-d)/2.;
        nule[numanswer  ] = -p/2.;
        ximag[numanswer++] = -sqrt(-d)/2.;
    }
//****************
        n -= 2;
		for(k=0; k<=n; k++) a[k] = b[k];
    }/*..while(order>2)..*/

	if(n==2) 
	{
		p = a[1];
		q = a[2];
	//**************
		
    d=p*p-4.*q;
    if(d>=0.)
	{
        nule[numanswer  ] = (-p+sqrt(d))/2.;
		ximag[numanswer++] = 0.;
        nule[numanswer  ] = (-p-sqrt(d))/2.;
		ximag[numanswer++] = 0.;
    }
    else
	{
        nule[numanswer  ] = -p/2.;
        ximag[numanswer++] = sqrt(-d)/2.;
        nule[numanswer  ] = -p/2.;
        ximag[numanswer++] = -sqrt(-d)/2.;
    }
	
	}
	else     nule[numanswer++] = -a[1];
}/*..main..*/

} ;