Back to home page

EIC code displayed by LXR

 
 

    


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

0001 // @(#)root/mathcore:$Id$
0002 // Author: L. Moneta,    11/2012
0003 
0004 /*************************************************************************
0005  * Copyright (C) 1995-2012, Rene Brun and Fons Rademakers.               *
0006  * All rights reserved.                                                  *
0007  *                                                                       *
0008  * For the licensing terms see $ROOTSYS/LICENSE.                         *
0009  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
0010  *************************************************************************/
0011 
0012 //////////////////////////////////////////////////////////////////////////
0013 //                                                                      //
0014 // Header file declaring functions for the evaluation of the Chebyshev  //
0015 // polynomials and the ChebyshevPol class which can be used for         //
0016 // creating a TF1.                                                      //
0017 //                                                                      //
0018 //////////////////////////////////////////////////////////////////////////
0019 
0020 #ifndef ROOT_Math_ChebyshevPol
0021 #define ROOT_Math_ChebyshevPol
0022 
0023 #include <sys/types.h>
0024 #include <cstring>
0025 
0026 namespace ROOT {
0027 
0028    namespace Math {
0029 
0030       /// template recursive functions for defining evaluation of Chebyshev polynomials
0031       ///  T_n(x) and the series S(x) = Sum_i c_i* T_i(x)
0032       namespace Chebyshev {
0033 
0034          template<int N> double T(double x) {
0035             return  (2.0 * x * T<N-1>(x)) - T<N-2>(x);
0036          }
0037 
0038          template<> double T<0> (double );
0039          template<> double T<1> (double x);
0040          template<> double T<2> (double x);
0041          template<> double T<3> (double x);
0042 
0043          template<int N> double Eval(double x, const double * c) {
0044             return c[N]*T<N>(x) + Eval<N-1>(x,c);
0045          }
0046 
0047          template<> double Eval<0> (double , const double *c);
0048          template<> double Eval<1> (double x, const double *c);
0049          template<> double Eval<2> (double x, const double *c);
0050          template<> double Eval<3> (double x, const double *c);
0051 
0052       } // end namespace Chebyshev
0053 
0054 
0055       // implementation of Chebyshev polynomials using all coefficients
0056       // needed for creating TF1 functions
0057       inline double Chebyshev0(double , double c0) {
0058          return c0;
0059       }
0060       inline double Chebyshev1(double x, double c0, double c1) {
0061          return c0 + c1*x;
0062       }
0063       inline double Chebyshev2(double x, double c0, double c1, double c2) {
0064          return c0 + c1*x + c2*(2.0*x*x - 1.0);
0065       }
0066       inline double Chebyshev3(double x, double c0, double c1, double c2, double c3) {
0067          return c3*Chebyshev::T<3>(x) + Chebyshev2(x,c0,c1,c2);
0068       }
0069       inline double Chebyshev4(double x, double c0, double c1, double c2, double c3, double c4) {
0070          return c4*Chebyshev::T<4>(x) + Chebyshev3(x,c0,c1,c2,c3);
0071       }
0072       inline double Chebyshev5(double x, double c0, double c1, double c2, double c3, double c4, double c5) {
0073          return c5*Chebyshev::T<5>(x) + Chebyshev4(x,c0,c1,c2,c3,c4);
0074       }
0075       inline double Chebyshev6(double x, double c0, double c1, double c2, double c3, double c4, double c5, double c6) {
0076          return c6*Chebyshev::T<6>(x) + Chebyshev5(x,c0,c1,c2,c3,c4,c5);
0077       }
0078       inline double Chebyshev7(double x, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7) {
0079          return c7*Chebyshev::T<7>(x) + Chebyshev6(x,c0,c1,c2,c3,c4,c5,c6);
0080       }
0081       inline double Chebyshev8(double x, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7, double c8) {   
0082          return c8*Chebyshev::T<8>(x) + Chebyshev7(x,c0,c1,c2,c3,c4,c5,c6,c7);
0083       }
0084       inline double Chebyshev9(double x, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7, double c8, double c9) {
0085          return c9*Chebyshev::T<9>(x) + Chebyshev8(x,c0,c1,c2,c3,c4,c5,c6,c7,c8);
0086       }
0087       inline double Chebyshev10(double x, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7, double c8, double c9, double c10) {
0088          return c10*Chebyshev::T<10>(x) + Chebyshev9(x,c0,c1,c2,c3,c4,c5,c6,c7,c8,c9);
0089       }
0090 
0091 
0092       // implementation of Chebyshev polynomial with run time parameter
0093       inline double ChebyshevN(unsigned int n, double x, const double * c) {
0094 
0095          if (n == 0) return Chebyshev0(x,c[0]);
0096          if (n == 1) return Chebyshev1(x,c[0],c[1]);
0097          if (n == 2) return Chebyshev2(x,c[0],c[1],c[2]);
0098          if (n == 3) return Chebyshev3(x,c[0],c[1],c[2],c[3]);
0099          if (n == 4) return Chebyshev4(x,c[0],c[1],c[2],c[3],c[4]);
0100          if (n == 5) return Chebyshev5(x,c[0],c[1],c[2],c[3],c[4],c[5]);
0101 
0102          /* do not use recursive formula
0103             (2.0 * x * Tn(n - 1, x)) - Tn(n - 2, x) ;
0104             which is too slow for large n
0105          */
0106 
0107          size_t i;
0108          double d1 = 0.0;
0109          double d2 = 0.0;
0110 
0111          // if not in range [-1,1]
0112          //double y = (2.0 * x - a - b) / (b - a);
0113          //double y = x;
0114          double y2 = 2.0 * x;
0115 
0116          for (i = n; i >= 1; i--)
0117          {
0118             double temp = d1;
0119             d1 = y2 * d1 - d2 + c[i];
0120             d2 = temp;
0121          }
0122 
0123          return x * d1 - d2 + c[0];
0124       }
0125 
0126 
0127       // implementation of Chebyshev Polynomial class
0128       // which can be used for building TF1 classes
0129       class ChebyshevPol {
0130       public:
0131          ChebyshevPol(unsigned int n) : fOrder(n) {}
0132 
0133          double operator() (const double *x, const double * coeff) {
0134             return ChebyshevN(fOrder, x[0], coeff);
0135          }
0136       private:
0137          unsigned int fOrder;
0138       };
0139 
0140 
0141 
0142    } // end namespace Math
0143 
0144 } // end namespace ROOT
0145 
0146 
0147 
0148 #endif  // ROOT_Math_Chebyshev