Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 10:10:20

0001 // @(#)root/mathmore:$Id$
0002 // Authors: L. Moneta, A. Zsenei   08/2005
0003 
0004  /**********************************************************************
0005   *                                                                    *
0006   * Copyright (c) 2004 ROOT Foundation,  CERN/PH-SFT                   *
0007   *                                                                    *
0008   * This library is free software; you can redistribute it and/or      *
0009   * modify it under the terms of the GNU General Public License        *
0010   * as published by the Free Software Foundation; either version 2     *
0011   * of the License, or (at your option) any later version.             *
0012   *                                                                    *
0013   * This library is distributed in the hope that it will be useful,    *
0014   * but WITHOUT ANY WARRANTY; without even the implied warranty of     *
0015   * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU   *
0016   * General Public License for more details.                           *
0017   *                                                                    *
0018   * You should have received a copy of the GNU General Public License  *
0019   * along with this library (see file COPYING); if not, write          *
0020   * to the Free Software Foundation, Inc., 59 Temple Place, Suite      *
0021   * 330, Boston, MA 02111-1307 USA, or contact the author.             *
0022   *                                                                    *
0023   **********************************************************************/
0024 
0025 // Header file for class Polynomial
0026 //
0027 // Created by: Lorenzo Moneta  at Wed Nov 10 17:46:19 2004
0028 //
0029 // Last update: Wed Nov 10 17:46:19 2004
0030 //
0031 #ifndef ROOT_Math_Polynomial
0032 #define ROOT_Math_Polynomial
0033 
0034 #include <complex>
0035 #include <vector>
0036 
0037 #include "Math/ParamFunction.h"
0038 
0039 // #ifdef _WIN32
0040 // #pragma warning(disable : 4250)
0041 // #endif
0042 
0043 namespace ROOT {
0044 namespace Math {
0045 
0046 //_____________________________________________________________________________________
0047   /**
0048      Parametric Function class describing polynomials of order n.
0049 
0050      <em>P(x) = p[0] + p[1]*x + p[2]*x**2 + ....... + p[n]*x**n</em>
0051 
0052      The class implements also the derivatives, \a dP(x)/dx and the \a dP(x)/dp[i].
0053 
0054      The class provides also the method to find the roots of the polynomial.
0055      It uses analytical methods up to quartic polynomials.
0056 
0057      Implements both the Parameteric function interface and the gradient interface
0058      since it provides the analytical gradient with respect to x
0059 
0060 
0061      @ingroup ParamFunc
0062   */
0063 
0064 class Polynomial : public ParamFunction<IParamGradFunction>,
0065                    public IGradientFunctionOneDim
0066 {
0067 
0068 
0069 public:
0070 
0071  typedef  ParamFunction<IParamGradFunction> ParFunc;
0072    /**
0073       Construct a Polynomial function of order n.
0074       The number of Parameters is n+1.
0075    */
0076 
0077    Polynomial(unsigned int n = 0);
0078 
0079    /**
0080       Construct a Polynomial of degree  1 : a*x + b
0081    */
0082    Polynomial(double a, double b);
0083 
0084    /**
0085       Construct a Polynomial of degree  2 : a*x**2 + b*x + c
0086    */
0087    Polynomial(double a, double b, double c);
0088 
0089    /**
0090       Construct a Polynomial of degree  3 : a*x**3 + b*x**2 + c*x + d
0091    */
0092    Polynomial(double a, double b, double c, double d);
0093 
0094    /**
0095       Construct a Polynomial of degree  4 : a*x**4 + b*x**3 + c*x**2 + dx  + e
0096    */
0097    Polynomial(double a, double b, double c, double d, double e);
0098 
0099 
0100    ~Polynomial() override {}
0101 
0102    // use default copy-ctor and assignment operators
0103 
0104 
0105 
0106 //   using ParamFunction::operator();
0107 
0108    /**
0109       Find the polynomial roots.
0110       For n <= 4, the roots are found analytically while for larger order an iterative numerical method is used
0111       The numerical method used is from GSL (see <A HREF="https://www.gnu.org/software/gsl/doc/html/poly.html">documentation</A> )
0112       For the case of n = 4 by default an analytical algorithm is used from an implementation by
0113        Andrew W. Steiner and Andy Buckley which is a translation from the original Cenrlib routine
0114        (< HREF="https://cds.cern.ch/record/2050876/files/c208.html">RRTEQ4</A> ).
0115       Note that depending on the coefficients the result could be not very accurate if the discriminant of the resolvent cubic
0116       equation is very small. In that case it might be more robust to use the numerical method, by calling directly FindNumRoots()
0117 
0118    */
0119    const std::vector<std::complex <double> > & FindRoots();
0120 
0121    /**
0122       Find the only the real polynomial roots.
0123       For n <= 4, the roots are found analytically while for larger order an iterative numerical method is used
0124       The numerical method used is from GSL (see <A HREF="https://www.gnu.org/software/gsl/doc/html/poly.html">documentation</A> )
0125    */
0126    std::vector<double > FindRealRoots();
0127 
0128    /**
0129       Find the polynomial roots using always an iterative numerical methods
0130       The numerical method used is from GSL (see <A HREF="https://www.gnu.org/software/gsl/doc/html/poly.html">documentation</A> )
0131    */
0132    const std::vector<std::complex <double> > & FindNumRoots();
0133 
0134    /**
0135       Order of Polynomial
0136    */
0137    unsigned int Order() const { return fOrder; }
0138 
0139 
0140    IGenFunction * Clone() const override;
0141 
0142    /**
0143        Optimized method to evaluate at the same time the function value and derivative at a point x.
0144        Implement the interface specified by ROOT::Math::IGradientOneDim.
0145        In the case of polynomial there is no advantage to compute both at the same time
0146    */
0147    void FdF (double x, double & f, double & df) const override {
0148       f = (*this)(x);
0149       df = Derivative(x);
0150    }
0151 
0152 
0153 private:
0154 
0155    double DoEvalPar ( double x, const double * p ) const override ;
0156 
0157    double DoDerivative (double x) const override ;
0158 
0159    double DoParameterDerivative(double x, const double * p, unsigned int ipar) const override;
0160 
0161 
0162    // cache order = number of params - 1)
0163    unsigned int fOrder;
0164 
0165    // cache Parameters for Gradient
0166    mutable std::vector<double> fDerived_params;
0167 
0168    // roots
0169 
0170    std::vector< std::complex < double > > fRoots;
0171 
0172 };
0173 
0174 } // namespace Math
0175 } // namespace ROOT
0176 
0177 
0178 #endif /* ROOT_Math_Polynomial */