ROOT   Reference Guide
Searching...
No Matches
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 }
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
#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 a(i)
Definition RSha256.hxx:99
#define e(i)
Definition RSha256.hxx:103
#define GSL_SUCCESS
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition IFunction.h:112
virtual const double * Parameters() const
Access the parameter values.
Parametric Function class describing polynomials of order n.
Definition Polynomial.h:66
const std::vector< std::complex< double > > & FindRoots()
Find the polynomial roots.
std::vector< std::complex< double > > fRoots
Definition Polynomial.h:167
Polynomial(unsigned int n=0)
Construct a Polynomial function of order n.
double DoEvalPar(double x, const double *p) const override
Implementation of the evaluation function using the x value and the parameters.
IGenFunction * Clone() const override
Clone a function.
double DoParameterDerivative(double x, const double *p, unsigned int ipar) const override
Evaluate the gradient, to be implemented by the derived classes.
std::vector< double > FindRealRoots()
Find the only the real polynomial roots.
const std::vector< std::complex< double > > & FindNumRoots()
Find the polynomial roots using always an iterative numerical methods The numerical method used is fr...
double DoDerivative(double x) const override
Function to evaluate the derivative with respect each coordinate. To be implemented by the derived cl...
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.
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
TCanvas * roots()
Definition roots.C:1
int gsl_poly_complex_solve_cubic(double a, double b, double c, gsl_complex *z0, gsl_complex *z1, gsl_complex *z2)