ROOT   6.10/09 Reference Guide
ChebyshevPol.h
Go to the documentation of this file.
1 // @(#)root/mathcore:$Id$
2 // Author: L. Moneta, 11/2012
3
4 /*************************************************************************
5  * Copyright (C) 1995-2012, Rene Brun and Fons Rademakers. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. * 9 * For the list of contributors see$ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11
12 //////////////////////////////////////////////////////////////////////////
13 // //
14 // Header file declaring functions for the evaluation of the Chebyshev //
15 // polynomials and the ChebyshevPol class which can be used for //
16 // creating a TF1. //
17 // //
18 //////////////////////////////////////////////////////////////////////////
19
20 #ifndef ROOT_Math_ChebyshevPol
21 #define ROOT_Math_ChebyshevPol
22
23 #include <sys/types.h>
24 #include <cstring>
25
26 namespace ROOT {
27
28  namespace Math {
29
30  /// template recursive functions for defining evaluation of Chebyshev polynomials
31  /// T_n(x) and the series S(x) = Sum_i c_i* T_i(x)
32  namespace Chebyshev {
33
34  template<int N> double T(double x) {
35  return (2.0 * x * T<N-1>(x)) - T<N-2>(x);
36  }
37
38  template<> double T<0> (double );
39  template<> double T<1> (double x);
40  template<> double T<2> (double x);
41  template<> double T<3> (double x);
42
43  template<int N> double Eval(double x, const double * c) {
44  return c[N]*T<N>(x) + Eval<N-1>(x,c);
45  }
46
47  template<> double Eval<0> (double , const double *c);
48  template<> double Eval<1> (double x, const double *c);
49  template<> double Eval<2> (double x, const double *c);
50  template<> double Eval<3> (double x, const double *c);
51
52  } // end namespace Chebyshev
53
54
55  // implementation of Chebyshev polynomials using all coefficients
56  // needed for creating TF1 functions
57  inline double Chebyshev0(double , double c0) {
58  return c0;
59  }
60  inline double Chebyshev1(double x, double c0, double c1) {
61  return c0 + c1*x;
62  }
63  inline double Chebyshev2(double x, double c0, double c1, double c2) {
64  return c0 + c1*x + c2*(2.0*x*x - 1.0);
65  }
66  inline double Chebyshev3(double x, double c0, double c1, double c2, double c3) {
67  return c3*Chebyshev::T<3>(x) + Chebyshev2(x,c0,c1,c2);
68  }
69  inline double Chebyshev4(double x, double c0, double c1, double c2, double c3, double c4) {
70  return c4*Chebyshev::T<4>(x) + Chebyshev3(x,c0,c1,c2,c3);
71  }
72  inline double Chebyshev5(double x, double c0, double c1, double c2, double c3, double c4, double c5) {
73  return c5*Chebyshev::T<5>(x) + Chebyshev4(x,c0,c1,c2,c3,c4);
74  }
75  inline double Chebyshev6(double x, double c0, double c1, double c2, double c3, double c4, double c5, double c6) {
76  return c6*Chebyshev::T<6>(x) + Chebyshev5(x,c0,c1,c2,c3,c4,c5);
77  }
78  inline double Chebyshev7(double x, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7) {
79  return c7*Chebyshev::T<7>(x) + Chebyshev6(x,c0,c1,c2,c3,c4,c5,c6);
80  }
81  inline double Chebyshev8(double x, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7, double c8) {
82  return c8*Chebyshev::T<8>(x) + Chebyshev7(x,c0,c1,c2,c3,c4,c5,c6,c7);
83  }
84  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) {
85  return c9*Chebyshev::T<9>(x) + Chebyshev8(x,c0,c1,c2,c3,c4,c5,c6,c7,c8);
86  }
87  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) {
88  return c10*Chebyshev::T<10>(x) + Chebyshev9(x,c0,c1,c2,c3,c4,c5,c6,c7,c8,c9);
89  }
90
91
92  // implementation of Chebyshev polynomial with run time parameter
93  inline double ChebyshevN(unsigned int n, double x, const double * c) {
94
95  if (n == 0) return Chebyshev0(x,c[0]);
96  if (n == 1) return Chebyshev1(x,c[0],c[1]);
97  if (n == 2) return Chebyshev2(x,c[0],c[1],c[2]);
98  if (n == 3) return Chebyshev3(x,c[0],c[1],c[2],c[3]);
99  if (n == 4) return Chebyshev4(x,c[0],c[1],c[2],c[3],c[4]);
100  if (n == 5) return Chebyshev5(x,c[0],c[1],c[2],c[3],c[4],c[5]);
101
102  /* do not use recursive formula
103  (2.0 * x * Tn(n - 1, x)) - Tn(n - 2, x) ;
104  which is too slow for large n
105  */
106
107  size_t i;
108  double d1 = 0.0;
109  double d2 = 0.0;
110
111  // if not in range [-1,1]
112  //double y = (2.0 * x - a - b) / (b - a);
113  //double y = x;
114  double y2 = 2.0 * x;
115
116  for (i = n; i >= 1; i--)
117  {
118  double temp = d1;
119  d1 = y2 * d1 - d2 + c[i];
120  d2 = temp;
121  }
122
123  return x * d1 - d2 + c[0];
124  }
125
126
127  // implementation of Chebyshev Polynomial class
128  // which can be used for building TF1 classes
129  class ChebyshevPol {
130  public:
131  ChebyshevPol(unsigned int n) : fOrder(n) {}
132
133  double operator() (const double *x, const double * coeff) {
134  return ChebyshevN(fOrder, x[0], coeff);
135  }
136  private:
137  unsigned int fOrder;
138  };
139
140
141
142  } // end namespace Math
143
144 } // end namespace ROOT
145
146
147
148 #endif // ROOT_Math_Chebyshev
double Chebyshev3(double x, double c0, double c1, double c2, double c3)
Definition: ChebyshevPol.h:66
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
double Eval< 1 >(double x, const double *c)
return c1
Definition: legend1.C:41
double T(double x)
Definition: ChebyshevPol.h:34
double Chebyshev2(double x, double c0, double c1, double c2)
Definition: ChebyshevPol.h:63
double T< 2 >(double x)
#define N
double Chebyshev8(double x, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7, double c8)
Definition: ChebyshevPol.h:81
double Chebyshev5(double x, double c0, double c1, double c2, double c3, double c4, double c5)
Definition: ChebyshevPol.h:72
TRObject operator()(const T1 &t1) const
double Eval< 2 >(double x, const double *c)
double Chebyshev6(double x, double c0, double c1, double c2, double c3, double c4, double c5, double c6)
Definition: ChebyshevPol.h:75
Double_t x[n]
Definition: legend1.C:17
double Chebyshev4(double x, double c0, double c1, double c2, double c3, double c4)
Definition: ChebyshevPol.h:69
double Chebyshev9(double x, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7, double c8, double c9)
Definition: ChebyshevPol.h:84
double Eval< 0 >(double, const double *c)
double ChebyshevN(unsigned int n, double x, const double *c)
Definition: ChebyshevPol.h:93
double Chebyshev7(double x, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7)
Definition: ChebyshevPol.h:78
return c2
Definition: legend2.C:14
ChebyshevPol(unsigned int n)
Definition: ChebyshevPol.h:131
double T< 0 >(double)
Namespace for new Math classes and functions.
double Eval(double x, const double *c)
Definition: ChebyshevPol.h:43
double Eval< 3 >(double x, const double *c)
double T< 3 >(double x)
double Chebyshev0(double, double c0)
Definition: ChebyshevPol.h:57
double Chebyshev1(double x, double c0, double c1)
Definition: ChebyshevPol.h:60
double T< 1 >(double x)
return c3
Definition: legend3.C:15
const Int_t n
Definition: legend1.C:16
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)
Definition: ChebyshevPol.h:87