File indexing completed on 2025-01-18 10:10:14
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
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
0031
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 }
0053
0054
0055
0056
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
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
0103
0104
0105
0106
0107 size_t i;
0108 double d1 = 0.0;
0109 double d2 = 0.0;
0110
0111
0112
0113
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
0128
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 }
0143
0144 }
0145
0146
0147
0148 #endif