Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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. *
6 * All rights reserved. *
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
26namespace 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
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
#define c(i)
Definition RSha256.hxx:101
#define N
ChebyshevPol(unsigned int n)
double operator()(const double *x, const double *coeff)
return c1
Definition legend1.C:41
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
return c2
Definition legend2.C:14
return c3
Definition legend3.C:15
Namespace for new Math classes and functions.
double Eval< 3 >(double x, const double *c)
double Eval(double x, const double *c)
double Eval< 2 >(double x, const double *c)
double T< 1 >(double x)
double T< 2 >(double x)
double T< 3 >(double x)
double Eval< 0 >(double, const double *c)
double T< 0 >(double)
double Eval< 1 >(double x, const double *c)
double T(double x)
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)
double ChebyshevN(unsigned int n, double x, const double *c)
double Chebyshev0(double, double c0)
double Chebyshev2(double x, double c0, double c1, double c2)
double Chebyshev3(double x, double c0, double c1, double c2, double c3)
double Chebyshev4(double x, double c0, double c1, double c2, double c3, double c4)
double Chebyshev1(double x, double c0, double c1)
double Chebyshev6(double x, double c0, double c1, double c2, double c3, double c4, double c5, double c6)
double Chebyshev7(double x, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7)
double Chebyshev5(double x, double c0, double c1, double c2, double c3, double c4, double c5)
double Chebyshev8(double x, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7, double c8)
double Chebyshev9(double x, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7, double c8, double c9)
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...