Logo ROOT   6.13/01
Reference Guide
Polynomial.cxx
Go to the documentation of this file.
1 // @(#)root/mathmore:$Id$
2 // Authors: L. Moneta, A. Zsenei 08/2005
3 
4  /**********************************************************************
5  * *
6  * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT *
7  * *
8  * This library is free software; you can redistribute it and/or *
9  * modify it under the terms of the GNU General Public License *
10  * as published by the Free Software Foundation; either version 2 *
11  * of the License, or (at your option) any later version. *
12  * *
13  * This library is distributed in the hope that it will be useful, *
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
16  * General Public License for more details. *
17  * *
18  * You should have received a copy of the GNU General Public License *
19  * along with this library (see file COPYING); if not, write *
20  * to the Free Software Foundation, Inc., 59 Temple Place, Suite *
21  * 330, Boston, MA 02111-1307 USA, or contact the author. *
22  * *
23  **********************************************************************/
24 
25 // Implementation file for class Polynomial
26 //
27 // Created by: Lorenzo Moneta at Wed Nov 10 17:46:19 2004
28 //
29 // Last update: Wed Nov 10 17:46:19 2004
30 //
31 
32 #include "Math/Polynomial.h"
33 
34 
35 #include "gsl/gsl_math.h"
36 #include "gsl/gsl_errno.h"
37 #include "gsl/gsl_poly.h"
38 #include "gsl/gsl_poly.h"
39 #include "complex_quartic.h"
40 
41 #define DEBUG
42 #ifdef DEBUG
43 #include <iostream>
44 #endif
45 
46 namespace ROOT {
47 namespace Math {
48 
49 
50 Polynomial::Polynomial(unsigned int n) :
51  // number of par is order + 1
52  ParFunc( n+1 ),
53  fOrder(n),
54  fDerived_params(std::vector<double>(n) )
55 {
56 }
57 
58  //order 1
59 Polynomial::Polynomial(double a, double b) :
60  ParFunc( 2 ),
61  fOrder(1),
62  fDerived_params(std::vector<double>(1) )
63 {
64  fParams[0] = b;
65  fParams[1] = a;
66 }
67 
68 // order 2
69 Polynomial::Polynomial(double a, double b, double c) :
70  ParFunc( 3 ),
71  fOrder(2),
72  fDerived_params(std::vector<double>(2) )
73 {
74  fParams[0] = c;
75  fParams[1] = b;
76  fParams[2] = a;
77 }
78 
79 // order 3 (cubic)
80 Polynomial::Polynomial(double a, double b, double c, double d) :
81  // number of par is order + 1
82  ParFunc( 4 ),
83  fOrder(3),
84  fDerived_params(std::vector<double>(3) )
85 {
86  fParams[0] = d;
87  fParams[1] = c;
88  fParams[2] = b;
89  fParams[3] = a;
90 }
91 
92 // order 3 (quartic)
93 Polynomial::Polynomial(double a, double b, double c, double d, double e) :
94  // number of par is order + 1
95  ParFunc( 5 ),
96  fOrder(4),
97  fDerived_params(std::vector<double>(4) )
98 {
99  fParams[0] = e;
100  fParams[1] = d;
101  fParams[2] = c;
102  fParams[3] = b;
103  fParams[4] = a;
104 }
105 
106 
107 
108 // Polynomial::Polynomial(const Polynomial &)
109 // {
110 // }
111 
112 // Polynomial & Polynomial::operator = (const Polynomial &rhs)
113 // {
114 // if (this == &rhs) return *this; // time saving self-test
115 
116 // return *this;
117 // }
118 
119 
120 double Polynomial::DoEvalPar (double x, const double * p) const {
121 
122  return gsl_poly_eval( p, fOrder + 1, x);
123 
124 }
125 
126 
127 
128 double Polynomial::DoDerivative(double x) const{
129 
130  for ( unsigned int i = 0; i < fOrder; ++i )
131  fDerived_params[i] = (i + 1) * Parameters()[i+1];
132 
133  return gsl_poly_eval( &fDerived_params.front(), fOrder, x);
134 
135 }
136 
137 double Polynomial::DoParameterDerivative (double x, const double *, unsigned int ipar) const {
138 
139  return gsl_pow_int(x, ipar);
140 }
141 
142 
143 
145  Polynomial * f = new Polynomial(fOrder);
147  f->SetParameters( Parameters() );
148  return f;
149 }
150 
151 
152 const std::vector< std::complex <double> > & Polynomial::FindRoots(){
153 
154 
155  // check if order is correct
156  unsigned int n = fOrder;
157  while ( Parameters()[n] == 0 ) {
158  n--;
159  }
160 
161  fRoots.clear();
162  fRoots.reserve(n);
163 
164 
165  if (n == 0) {
166  return fRoots;
167  }
168  else if (n == 1 ) {
169  if ( Parameters()[1] == 0) return fRoots;
170  double r = - Parameters()[0]/ Parameters()[1];
171  fRoots.push_back( std::complex<double> ( r, 0.0) );
172  }
173  // quadratic equations
174  else if (n == 2 ) {
175  gsl_complex z1, z2;
176  int nr = gsl_poly_complex_solve_quadratic(Parameters()[2], Parameters()[1], Parameters()[0], &z1, &z2);
177  if ( nr != 2) {
178 #ifdef DEBUG
179  std::cout << "Polynomial quadratic ::- FAILED to find roots" << std::endl;
180 #endif
181  return fRoots;
182  }
183  fRoots.push_back(std::complex<double>(z1.dat[0],z1.dat[1]) );
184  fRoots.push_back(std::complex<double>(z2.dat[0],z2.dat[1]) );
185  }
186  // cubic equations
187  else if (n == 3 ) {
188  gsl_complex z1, z2, z3;
189  // renormmalize params in this case
190  double w = Parameters()[3];
191  double a = Parameters()[2]/w;
192  double b = Parameters()[1]/w;
193  double c = Parameters()[0]/w;
194  int nr = gsl_poly_complex_solve_cubic(a, b, c, &z1, &z2, &z3);
195  if (nr != 3) {
196 #ifdef DEBUG
197  std::cout << "Polynomial cubic::- FAILED to find roots" << std::endl;
198 #endif
199  return fRoots;
200 
201  }
202  fRoots.push_back(std::complex<double> (z1.dat[0],z1.dat[1]) );
203  fRoots.push_back(std::complex<double> (z2.dat[0],z2.dat[1]) );
204  fRoots.push_back(std::complex<double> (z3.dat[0],z3.dat[1]) );
205  }
206  // quartic equations
207  else if (n == 4 ) {
208  // quartic eq.
209  gsl_complex z1, z2, z3, z4;
210  // renormalize params in this case
211  double w = Parameters()[4];
212  double a = Parameters()[3]/w;
213  double b = Parameters()[2]/w;
214  double c = Parameters()[1]/w;
215  double d = Parameters()[0]/w;
216  int nr = gsl_poly_complex_solve_quartic(a, b, c, d, &z1, &z2, &z3, & z4);
217  if (nr != 4) {
218 #ifdef DEBUG
219  std::cout << "Polynomial quartic ::- FAILED to find roots" << std::endl;
220 #endif
221  return fRoots;
222  }
223  fRoots.push_back(std::complex<double> (z1.dat[0],z1.dat[1]) );
224  fRoots.push_back(std::complex<double> (z2.dat[0],z2.dat[1]) );
225  fRoots.push_back(std::complex<double> (z3.dat[0],z3.dat[1]) );
226  fRoots.push_back(std::complex<double> (z4.dat[0],z4.dat[1]) );
227  }
228  else {
229  // for higher order polynomial use numerical fRoots
230  FindNumRoots();
231  }
232 
233  return fRoots;
234 
235  }
236 
237 
238 std::vector< double > Polynomial::FindRealRoots(){
239  FindRoots();
240  std::vector<double> roots;
241  roots.reserve(fOrder);
242  for (unsigned int i = 0; i < fOrder; ++i) {
243  if (fRoots[i].imag() == 0)
244  roots.push_back( fRoots[i].real() );
245  }
246  return roots;
247 }
248 const std::vector< std::complex <double> > & Polynomial::FindNumRoots(){
249 
250 
251  // check if order is correct
252  unsigned int n = fOrder;
253  while ( Parameters()[n] == 0 ) {
254  n--;
255  }
256  fRoots.clear();
257  fRoots.reserve(n);
258 
259 
260  if (n == 0) {
261  return fRoots;
262  }
263 
264  gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc( n + 1);
265  std::vector<double> z(2*n);
266  int status = gsl_poly_complex_solve ( Parameters(), n+1, w, &z.front() );
267  gsl_poly_complex_workspace_free(w);
268  if (status != GSL_SUCCESS) return fRoots;
269  for (unsigned int i = 0; i < n; ++i)
270  fRoots.push_back(std::complex<double> (z[2*i],z[2*i+1] ) );
271 
272  return fRoots;
273 }
274 
275 
276 } // namespace Math
277 } // namespace ROOT
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition: IFunction.h:135
Namespace for new ROOT classes and functions.
Definition: TFoamSampler.h:19
Polynomial(unsigned int n=0)
Construct a Polynomial function of order n.
Definition: Polynomial.cxx:50
std::vector< double > FindRealRoots()
Find the only the real polynomial roots.
Definition: Polynomial.cxx:238
std::vector< double > fDerived_params
Definition: Polynomial.h:162
const std::vector< std::complex< double > > & FindRoots()
Find the polynomial roots.
Definition: Polynomial.cxx:152
STL namespace.
virtual const double * Parameters() const
Access the parameter values.
Definition: ParamFunction.h:96
const std::vector< std::complex< double > > & FindNumRoots()
Find the polynomial roots using always an iterative numerical methods The numerical method used is fr...
Definition: Polynomial.cxx:248
double DoEvalPar(double x, const double *p) const
Implementation of the evaluation function using the x value and the parameters.
Definition: Polynomial.cxx:120
IGenFunction * Clone() const
Clone a function.
Definition: Polynomial.cxx:144
unsigned int fOrder
Definition: Polynomial.h:159
virtual void SetParameters(const double *p)
Set the parameter values.
std::vector< std::complex< double > > fRoots
Definition: Polynomial.h:166
Parametric Function class describing polynomials of order n.
Definition: Polynomial.h:63
* x
Deprecated and error prone model selection interface.
Definition: TRolke.cxx:630
int gsl_poly_complex_solve_cubic(double a, double b, double c, gsl_complex *z0, gsl_complex *z1, gsl_complex *z2)
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
Namespace for new Math classes and functions.
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
int gsl_poly_complex_solve_quartic(double a, double b, double c, double d, gsl_complex *z0, gsl_complex *z1, gsl_complex *z2, gsl_complex *z3)
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
double DoParameterDerivative(double x, const double *p, unsigned int ipar) const
Evaluate the gradient, to be implemented by the derived classes.
Definition: Polynomial.cxx:137
double DoDerivative(double x) const
function to evaluate the derivative with respect each coordinate.
Definition: Polynomial.cxx:128