Logo ROOT   6.16/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
46namespace ROOT {
47namespace Math {
48
49
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
59Polynomial::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
69Polynomial::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)
80Polynomial::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)
93Polynomial::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
120double Polynomial::DoEvalPar (double x, const double * p) const {
121
122 return gsl_poly_eval( p, fOrder + 1, x);
123
124}
125
126
127
128double 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
137double 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);
146 f->fDerived_params = fDerived_params;
147 f->SetParameters( Parameters() );
148 return f;
149}
150
151
152const 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
238std::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}
248const 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
ROOT::R::TRInterface & r
Definition: Object.C:4
#define d(i)
Definition: RSha256.hxx:102
#define b(i)
Definition: RSha256.hxx:100
#define f(i)
Definition: RSha256.hxx:104
#define c(i)
Definition: RSha256.hxx:101
#define e(i)
Definition: RSha256.hxx:103
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition: IFunction.h:135
virtual const double * Parameters() const
Access the parameter values.
Definition: ParamFunction.h:96
Parametric Function class describing polynomials of order n.
Definition: Polynomial.h:65
IGenFunction * Clone() const
Clone a function.
Definition: Polynomial.cxx:144
double DoDerivative(double x) const
function to evaluate the derivative with respect each coordinate.
Definition: Polynomial.cxx:128
const std::vector< std::complex< double > > & FindRoots()
Find the polynomial roots.
Definition: Polynomial.cxx:152
std::vector< std::complex< double > > fRoots
Definition: Polynomial.h:166
Polynomial(unsigned int n=0)
Construct a Polynomial function of order n.
Definition: Polynomial.cxx:50
unsigned int fOrder
Definition: Polynomial.h:159
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 DoEvalPar(double x, const double *p) const
Implementation of the evaluation function using the x value and the parameters.
Definition: Polynomial.cxx:120
std::vector< double > FindRealRoots()
Find the only the real polynomial roots.
Definition: Polynomial.cxx:238
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
std::vector< double > fDerived_params
Definition: Polynomial.h:162
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)
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
Namespace for new Math classes and functions.
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
STL namespace.
TCanvas * roots()
Definition: roots.C:1
auto * a
Definition: textangle.C:12
int gsl_poly_complex_solve_cubic(double a, double b, double c, gsl_complex *z0, gsl_complex *z1, gsl_complex *z2)