Polynomial.cxx

Go to the documentation of this file.

// @(#)root/mathmore:$Id$
// 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 "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
  ParFunc( n+1 ),
  fOrder(n),
  fDerived_params(std::vector<double>(n) )
{
}
 
  //order 1
Polynomial::Polynomial(double a, double b) :
  ParFunc( 2 ),
  fOrder(1),
  fDerived_params(std::vector<double>(1) )
{
  fParams[0] = b;
  fParams[1] = a;
}
 
// order 2
Polynomial::Polynomial(double a, double b, double c) :
  ParFunc( 3 ),
  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
  ParFunc( 4 ),
  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
  ParFunc( 5 ),
  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(const Polynomial &)
// {
// }
 
// Polynomial & Polynomial::operator = (const Polynomial &rhs)
// {
//    if (this == &rhs) return *this;  // time saving self-test
 
//    return *this;
// }
 
 
double  Polynomial::DoEvalPar (double x, const double * p) const {
 
    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);
 
}
 
double Polynomial::DoParameterDerivative (double x, const double *, unsigned int ipar) const {
 
      return gsl_pow_int(x, ipar);
}
 
 
 
IGenFunction * Polynomial::Clone() const {
    Polynomial * f =  new Polynomial(fOrder);
    f->fDerived_params = fDerived_params;
    f->SetParameters( Parameters() );
    return f;
}
 
std::vector<std::complex<double>> Polynomial::FindRoots() const
{
 
   // check if order is correct
   unsigned int n = fOrder;
   while (Parameters()[n] == 0) {
      n--;
   }
 
   std::vector<std::complex<double>> roots;
 
   if (n == 0) {
      return roots;
   } else if (n == 1) {
      if (Parameters()[1] == 0)
         return roots;
      double r = - Parameters()[0]/ Parameters()[1];
      roots.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 roots;
      }
      roots.push_back(std::complex<double>(z1.dat[0], z1.dat[1]));
      roots.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 roots;
      }
      roots.push_back(std::complex<double>(z1.dat[0], z1.dat[1]));
      roots.push_back(std::complex<double>(z2.dat[0], z2.dat[1]));
      roots.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 roots;
      }
      roots.push_back(std::complex<double>(z1.dat[0], z1.dat[1]));
      roots.push_back(std::complex<double>(z2.dat[0], z2.dat[1]));
      roots.push_back(std::complex<double>(z3.dat[0], z3.dat[1]));
      roots.push_back(std::complex<double>(z4.dat[0], z4.dat[1]));
   } else {
      // for higher order polynomial use numerical roots
      return FindNumRoots();
   }
 
   return roots;
}
 
std::vector<double> Polynomial::FindRealRoots() const
{
   auto roots = FindRoots();
   std::vector<double> realRoots;
   for (const auto &r : roots)
      if (std::abs(r.imag()) < std::numeric_limits<double>::epsilon())
         realRoots.push_back(r.real());
 
   return realRoots;
}
std::vector<std::complex<double>> Polynomial::FindNumRoots() const
{
 
   // check if order is correct
   unsigned int n = fOrder;
   while (Parameters()[n] == 0) {
      n--;
   }
 
   std::vector<std::complex<double>> roots;
   if (n == 0) {
      return roots;
   }
 
   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 roots;
   for (unsigned int i = 0; i < n; ++i)
      roots.push_back(std::complex<double>(z[2 * i], z[2 * i + 1]));
 
   return roots;
}
 
} // namespace Math
} // namespace ROOT