Logo ROOT  
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 "complex_quartic.h"
39
40#define DEBUG
41#ifdef DEBUG
42#include <iostream>
43#endif
44
45namespace ROOT {
46namespace Math {
47
48
50 // number of par is order + 1
51 ParFunc( n+1 ),
52 fOrder(n),
53 fDerived_params(std::vector<double>(n) )
54{
55}
56
57 //order 1
58Polynomial::Polynomial(double a, double b) :
59 ParFunc( 2 ),
60 fOrder(1),
61 fDerived_params(std::vector<double>(1) )
62{
63 fParams[0] = b;
64 fParams[1] = a;
65}
66
67// order 2
68Polynomial::Polynomial(double a, double b, double c) :
69 ParFunc( 3 ),
70 fOrder(2),
71 fDerived_params(std::vector<double>(2) )
72{
73 fParams[0] = c;
74 fParams[1] = b;
75 fParams[2] = a;
76}
77
78// order 3 (cubic)
79Polynomial::Polynomial(double a, double b, double c, double d) :
80 // number of par is order + 1
81 ParFunc( 4 ),
82 fOrder(3),
83 fDerived_params(std::vector<double>(3) )
84{
85 fParams[0] = d;
86 fParams[1] = c;
87 fParams[2] = b;
88 fParams[3] = a;
89}
90
91// order 3 (quartic)
92Polynomial::Polynomial(double a, double b, double c, double d, double e) :
93 // number of par is order + 1
94 ParFunc( 5 ),
95 fOrder(4),
96 fDerived_params(std::vector<double>(4) )
97{
98 fParams[0] = e;
99 fParams[1] = d;
100 fParams[2] = c;
101 fParams[3] = b;
102 fParams[4] = a;
103}
104
105
106
107// Polynomial::Polynomial(const Polynomial &)
108// {
109// }
110
111// Polynomial & Polynomial::operator = (const Polynomial &rhs)
112// {
113// if (this == &rhs) return *this; // time saving self-test
114
115// return *this;
116// }
117
118
119double Polynomial::DoEvalPar (double x, const double * p) const {
120
121 return gsl_poly_eval( p, fOrder + 1, x);
122
123}
124
125
126
127double Polynomial::DoDerivative(double x) const{
128
129 for ( unsigned int i = 0; i < fOrder; ++i )
130 fDerived_params[i] = (i + 1) * Parameters()[i+1];
131
132 return gsl_poly_eval( &fDerived_params.front(), fOrder, x);
133
134}
135
136double Polynomial::DoParameterDerivative (double x, const double *, unsigned int ipar) const {
137
138 return gsl_pow_int(x, ipar);
139}
140
141
142
144 Polynomial * f = new Polynomial(fOrder);
145 f->fDerived_params = fDerived_params;
146 f->SetParameters( Parameters() );
147 return f;
148}
149
150
151const std::vector< std::complex <double> > & Polynomial::FindRoots(){
152
153
154 // check if order is correct
155 unsigned int n = fOrder;
156 while ( Parameters()[n] == 0 ) {
157 n--;
158 }
159
160 fRoots.clear();
161 fRoots.reserve(n);
162
163
164 if (n == 0) {
165 return fRoots;
166 }
167 else if (n == 1 ) {
168 if ( Parameters()[1] == 0) return fRoots;
169 double r = - Parameters()[0]/ Parameters()[1];
170 fRoots.push_back( std::complex<double> ( r, 0.0) );
171 }
172 // quadratic equations
173 else if (n == 2 ) {
174 gsl_complex z1, z2;
175 int nr = gsl_poly_complex_solve_quadratic(Parameters()[2], Parameters()[1], Parameters()[0], &z1, &z2);
176 if ( nr != 2) {
177#ifdef DEBUG
178 std::cout << "Polynomial quadratic ::- FAILED to find roots" << std::endl;
179#endif
180 return fRoots;
181 }
182 fRoots.push_back(std::complex<double>(z1.dat[0],z1.dat[1]) );
183 fRoots.push_back(std::complex<double>(z2.dat[0],z2.dat[1]) );
184 }
185 // cubic equations
186 else if (n == 3 ) {
187 gsl_complex z1, z2, z3;
188 // renormmalize params in this case
189 double w = Parameters()[3];
190 double a = Parameters()[2]/w;
191 double b = Parameters()[1]/w;
192 double c = Parameters()[0]/w;
193 int nr = gsl_poly_complex_solve_cubic(a, b, c, &z1, &z2, &z3);
194 if (nr != 3) {
195#ifdef DEBUG
196 std::cout << "Polynomial cubic::- FAILED to find roots" << std::endl;
197#endif
198 return fRoots;
199
200 }
201 fRoots.push_back(std::complex<double> (z1.dat[0],z1.dat[1]) );
202 fRoots.push_back(std::complex<double> (z2.dat[0],z2.dat[1]) );
203 fRoots.push_back(std::complex<double> (z3.dat[0],z3.dat[1]) );
204 }
205 // quartic equations
206 else if (n == 4 ) {
207 // quartic eq.
208 gsl_complex z1, z2, z3, z4;
209 // renormalize params in this case
210 double w = Parameters()[4];
211 double a = Parameters()[3]/w;
212 double b = Parameters()[2]/w;
213 double c = Parameters()[1]/w;
214 double d = Parameters()[0]/w;
215 int nr = gsl_poly_complex_solve_quartic(a, b, c, d, &z1, &z2, &z3, & z4);
216 if (nr != 4) {
217#ifdef DEBUG
218 std::cout << "Polynomial quartic ::- FAILED to find roots" << std::endl;
219#endif
220 return fRoots;
221 }
222 fRoots.push_back(std::complex<double> (z1.dat[0],z1.dat[1]) );
223 fRoots.push_back(std::complex<double> (z2.dat[0],z2.dat[1]) );
224 fRoots.push_back(std::complex<double> (z3.dat[0],z3.dat[1]) );
225 fRoots.push_back(std::complex<double> (z4.dat[0],z4.dat[1]) );
226 }
227 else {
228 // for higher order polynomial use numerical fRoots
229 FindNumRoots();
230 }
231
232 return fRoots;
233
234 }
235
236
237std::vector< double > Polynomial::FindRealRoots(){
238 FindRoots();
239 std::vector<double> roots;
240 roots.reserve(fOrder);
241 for (unsigned int i = 0; i < fOrder; ++i) {
242 if (fRoots[i].imag() == 0)
243 roots.push_back( fRoots[i].real() );
244 }
245 return roots;
246}
247const std::vector< std::complex <double> > & Polynomial::FindNumRoots(){
248
249
250 // check if order is correct
251 unsigned int n = fOrder;
252 while ( Parameters()[n] == 0 ) {
253 n--;
254 }
255 fRoots.clear();
256 fRoots.reserve(n);
257
258
259 if (n == 0) {
260 return fRoots;
261 }
262
263 gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc( n + 1);
264 std::vector<double> z(2*n);
265 int status = gsl_poly_complex_solve ( Parameters(), n+1, w, &z.front() );
266 gsl_poly_complex_workspace_free(w);
267 if (status != GSL_SUCCESS) return fRoots;
268 for (unsigned int i = 0; i < n; ++i)
269 fRoots.push_back(std::complex<double> (z[2*i],z[2*i+1] ) );
270
271 return fRoots;
272}
273
274
275} // namespace Math
276} // namespace ROOT
double
Definition: Converters.cxx:939
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:66
IGenFunction * Clone() const
Clone a function.
Definition: Polynomial.cxx:143
double DoDerivative(double x) const
function to evaluate the derivative with respect each coordinate.
Definition: Polynomial.cxx:127
const std::vector< std::complex< double > > & FindRoots()
Find the polynomial roots.
Definition: Polynomial.cxx:151
std::vector< std::complex< double > > fRoots
Definition: Polynomial.h:167
Polynomial(unsigned int n=0)
Construct a Polynomial function of order n.
Definition: Polynomial.cxx:49
unsigned int fOrder
Definition: Polynomial.h:160
double DoParameterDerivative(double x, const double *p, unsigned int ipar) const
Evaluate the gradient, to be implemented by the derived classes.
Definition: Polynomial.cxx:136
double DoEvalPar(double x, const double *p) const
Implementation of the evaluation function using the x value and the parameters.
Definition: Polynomial.cxx:119
std::vector< double > FindRealRoots()
Find the only the real polynomial roots.
Definition: Polynomial.cxx:237
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:247
std::vector< double > fDerived_params
Definition: Polynomial.h:163
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.
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Definition: RNumpyDS.hxx:30
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)