// @(#)root/mathmore:$Name:  $:$Id: Polynomial.cxx,v 1.6 2006/11/23 17:28:57 moneta Exp $
// Authors: L. Moneta, A. Zsenei   08/2005 

 /**********************************************************************
  *                                                                    *
  * Copyright (c) 2004 ROOT Foundation,  CERN/PH-SFT                   *
  *                                                                    *
  * This library is free software; you can redistribute it and/or      *
  * modify it under the terms of the GNU General Public License        *
  * as published by the Free Software Foundation; either version 2     *
  * of the License, or (at your option) any later version.             *
  *                                                                    *
  * This library is distributed in the hope that it will be useful,    *
  * but WITHOUT ANY WARRANTY; without even the implied warranty of     *
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU   *
  * General Public License for more details.                           *
  *                                                                    *
  * You should have received a copy of the GNU General Public License  *
  * along with this library (see file COPYING); if not, write          *
  * to the Free Software Foundation, Inc., 59 Temple Place, Suite      *
  * 330, Boston, MA 02111-1307 USA, or contact the author.             *
  *                                                                    *
  **********************************************************************/

// Implementation file for class Polynomial
// 
// Created by: Lorenzo Moneta  at Wed Nov 10 17:46:19 2004
// 
// Last update: Wed Nov 10 17:46:19 2004
// 

#include "Math/Polynomial.h"


#include "gsl/gsl_math.h"
#include "gsl/gsl_errno.h"
#include "gsl/gsl_poly.h"
#include "gsl/gsl_poly.h"
#include "complex_quartic.h"

#define DEBUG
#ifdef DEBUG
#include <iostream>
#endif

namespace ROOT {
namespace Math {


Polynomial::Polynomial(unsigned int n) : 
  // number of par is order + 1
  ParamFunction( n+1, true, true), 
  fOrder(n), 
  fDerived_params(std::vector<double>(n) )
{
}

  //order 1
Polynomial::Polynomial(double a, double b) : 
  ParamFunction( 2, true, true), 
  fOrder(1), 
  fDerived_params(std::vector<double>(1) )
{
  fParams[0] = b; 
  fParams[1] = a; 
}

// order 2
Polynomial::Polynomial(double a, double b, double c) : 
  ParamFunction( 3, true, true), 
  fOrder(2), 
  fDerived_params(std::vector<double>(2) )
{
  fParams[0] = c; 
  fParams[1] = b; 
  fParams[2] = a; 
}

// order 3 (cubic)
Polynomial::Polynomial(double a, double b, double c, double d) : 
  // number of par is order + 1
  ParamFunction( 4, true, true), 
  fOrder(3), 
  fDerived_params(std::vector<double>(3) )
{
  fParams[0] = d; 
  fParams[1] = c; 
  fParams[2] = b; 
  fParams[3] = a; 
}

// order 3 (quartic)
Polynomial::Polynomial(double a, double b, double c, double d, double e) : 
  // number of par is order + 1
  ParamFunction( 5, true, true), 
  fOrder(4), 
  fDerived_params(std::vector<double>(4) )
{
  fParams[0] = e; 
  fParams[1] = d; 
  fParams[2] = c; 
  fParams[3] = b; 
  fParams[4] = a; 
}



Polynomial::~Polynomial() 
{
}

// Polynomial::Polynomial(const Polynomial &) 
// {
// }

// Polynomial & Polynomial::operator = (const Polynomial &rhs) 
// {
//    if (this == &rhs) return *this;  // time saving self-test

//    return *this;
// }


double  Polynomial::DoEval (double x) const { 
  
    return gsl_poly_eval( &fParams.front(), fOrder + 1, x); 

}


double  Polynomial::operator() (double x, const double *  p) { 
  
    return gsl_poly_eval( p, fOrder + 1, x); 

}


double  Polynomial::DoDerivative(double x) const{ 

   for ( unsigned int i = 0; i < fOrder; ++i ) 
    fDerived_params[i] =  (i + 1) * Parameters()[i+1]; 

   return gsl_poly_eval( &fDerived_params.front(), fOrder, x); 

}

void Polynomial::DoParameterGradient (double x, double * grad) const { 

   unsigned int npar = NPar(); 
   for (unsigned int i = 0; i < npar; ++i) 
      grad[i] = gsl_pow_int(x, i); 
      // return fParGradient; 

}



IGenFunction * Polynomial::Clone() const { 
    Polynomial * f =  new Polynomial(fOrder);
    f->fDerived_params = fDerived_params; 
    f->SetParameters( Parameters() ); 
    return f; 
}
	

const std::vector< std::complex <double> > &  Polynomial::FindRoots(){


    // check if order is correct
    unsigned int n = fOrder;
    while ( Parameters()[n] == 0 ) { 
      n--;
    } 

    fRoots.clear();  
    fRoots.reserve(n);  


    if (n == 0) {   
      return fRoots;
    } 
    else if (n == 1 ) { 
      if ( Parameters()[1] == 0) return fRoots;
      double r = - Parameters()[0]/ Parameters()[1]; 
      fRoots.push_back( std::complex<double> ( r, 0.0) ); 
    }
    // quadratic equations
    else if (n == 2 ) { 
      gsl_complex z1, z2;
      int nr = gsl_poly_complex_solve_quadratic(Parameters()[2], Parameters()[1], Parameters()[0], &z1, &z2); 
      if ( nr != 2) { 
#ifdef DEBUG
	std::cout << "Polynomial quadratic ::-  FAILED to find roots" << std::endl; 
#endif
	return fRoots;
      } 
      fRoots.push_back(std::complex<double>(z1.dat[0],z1.dat[1]) ); 
      fRoots.push_back(std::complex<double>(z2.dat[0],z2.dat[1]) );       
    }
    // cubic equations
    else if (n == 3 ) { 
      gsl_complex  z1, z2, z3;   
      // renormmalize params in this case
      double w = Parameters()[3]; 
      double a = Parameters()[2]/w;
      double b = Parameters()[1]/w;
      double c = Parameters()[0]/w;
      int nr = gsl_poly_complex_solve_cubic(a, b, c, &z1, &z2, &z3); 
      if (nr != 3) { 
#ifdef DEBUG
	std::cout << "Polynomial  cubic::-  FAILED to find roots" << std::endl; 
#endif
	return fRoots; 

      }
      fRoots.push_back(std::complex<double> (z1.dat[0],z1.dat[1]) ); 
      fRoots.push_back(std::complex<double> (z2.dat[0],z2.dat[1]) ); 
      fRoots.push_back(std::complex<double> (z3.dat[0],z3.dat[1]) );       
    }
    // quartic equations
    else if (n == 4 ) { 
      // quartic eq.  
      gsl_complex  z1, z2, z3, z4;   
      // renormalize params in this case
      double w = Parameters()[4]; 
      double a = Parameters()[3]/w;
      double b = Parameters()[2]/w;
      double c = Parameters()[1]/w;
      double d = Parameters()[0]/w;
      int nr = gsl_poly_complex_solve_quartic(a, b, c, d, &z1, &z2, &z3, & z4); 
      if (nr != 4) { 
#ifdef DEBUG
	std::cout << "Polynomial quartic ::-  FAILED to find roots" << std::endl; 
#endif
	return fRoots;
      } 
      fRoots.push_back(std::complex<double> (z1.dat[0],z1.dat[1]) ); 
      fRoots.push_back(std::complex<double> (z2.dat[0],z2.dat[1]) ); 
      fRoots.push_back(std::complex<double> (z3.dat[0],z3.dat[1]) );       
      fRoots.push_back(std::complex<double> (z4.dat[0],z4.dat[1]) );       
    }
    else { 
    // for higher order polynomial use numerical fRoots
      FindNumRoots();
    }      

    return fRoots; 

  }


std::vector< double >  Polynomial::FindRealRoots(){
  FindRoots(); 
  std::vector<double> roots; 
  roots.reserve(fOrder); 
  for (unsigned int i = 0; i < fOrder; ++i) { 
    if (fRoots[i].imag() == 0) 
      roots.push_back( fRoots[i].real() ); 
  }
  return roots; 
}
const std::vector< std::complex <double> > &  Polynomial::FindNumRoots(){


    // check if order is correct
    unsigned int n = fOrder;
    while ( Parameters()[n] == 0 ) { 
      n--;
    } 
    fRoots.clear();
    fRoots.reserve(n);  


    if (n == 0) {   
      return fRoots;
    } 

    gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc( n + 1); 
    std::vector<double> z(2*n);
    int status = gsl_poly_complex_solve ( Parameters(), n+1, w, &z.front() );  
    gsl_poly_complex_workspace_free(w);
    if (status != GSL_SUCCESS) return fRoots; 
    for (unsigned int i = 0; i < n; ++i) 
      fRoots.push_back(std::complex<double> (z[2*i],z[2*i+1] ) );      

    return fRoots; 
}


} // namespace Math
} // namespace ROOT