// @(#)root/mathcore:$Id$
// Author: L. Moneta,    11/2012

/*************************************************************************
 * Copyright (C) 1995-2012, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

//////////////////////////////////////////////////////////////////////////
//                                                                      //
// Header file declaring functions for the evaluation of the Chebyshev  //
// polynomials and the ChebyshevPol class which can be used for         //
// creating a TF1.                                                      //
//                                                                      //
//////////////////////////////////////////////////////////////////////////

#ifndef ROOT_Math_ChebyshevPol
#define ROOT_Math_ChebyshevPol

#include <sys/types.h>
#include <cstring>

namespace ROOT {

   namespace Math {

      /// template recursive functions for defining evaluation of Chebyshev polynomials
      ///  T_n(x) and the series S(x) = Sum_i c_i* T_i(x)
      namespace Chebyshev {

         template<int N> double T(double x) {
            return  (2.0 * x * T<N-1>(x)) - T<N-2>(x);
         }

         template<> double T<0> (double );
         template<> double T<1> (double x);
         template<> double T<2> (double x);
         template<> double T<3> (double x);

         template<int N> double Eval(double x, const double * c) {
            return c[N]*T<N>(x) + Eval<N-1>(x,c);
         }

         template<> double Eval<0> (double , const double *c);
         template<> double Eval<1> (double x, const double *c);
         template<> double Eval<2> (double x, const double *c);
         template<> double Eval<3> (double x, const double *c);

      } // end namespace Chebyshev


      // implementation of Chebyshev polynomials using all coefficients
      // needed for creating TF1 functions
      inline double Chebyshev0(double , double c0) {
         return c0;
      }
      inline double Chebyshev1(double x, double c0, double c1) {
         return c0 + c1*x;
      }
      inline double Chebyshev2(double x, double c0, double c1, double c2) {
         return c0 + c1*x + c2*(2.0*x*x - 1.0);
      }
      inline double Chebyshev3(double x, double c0, double c1, double c2, double c3) {
         return c3*Chebyshev::T<3>(x) + Chebyshev2(x,c0,c1,c2);
      }
      inline double Chebyshev4(double x, double c0, double c1, double c2, double c3, double c4) {
         return c4*Chebyshev::T<4>(x) + Chebyshev3(x,c0,c1,c2,c3);
      }
      inline double Chebyshev5(double x, double c0, double c1, double c2, double c3, double c4, double c5) {
         return c5*Chebyshev::T<5>(x) + Chebyshev4(x,c0,c1,c2,c3,c4);
      }


      // implementation of Chebyshev polynomial with run time parameter
      inline double ChebyshevN(unsigned int n, double x, const double * c) {

         if (n == 0) return Chebyshev0(x,c[0]);
         if (n == 1) return Chebyshev1(x,c[0],c[1]);
         if (n == 2) return Chebyshev2(x,c[0],c[1],c[2]);
         if (n == 3) return Chebyshev3(x,c[0],c[1],c[2],c[3]);
         if (n == 4) return Chebyshev4(x,c[0],c[1],c[2],c[3],c[4]);
         if (n == 5) return Chebyshev5(x,c[0],c[1],c[2],c[3],c[4],c[5]);

         /* do not use recursive formula
            (2.0 * x * Tn(n - 1, x)) - Tn(n - 2, x) ;
            which is too slow for large n
         */

         size_t i;
         double d1 = 0.0;
         double d2 = 0.0;

         // if not in range [-1,1]
         //double y = (2.0 * x - a - b) / (b - a);
         //double y = x;
         double y2 = 2.0 * x;

         for (i = n; i >= 1; i--)
         {
            double temp = d1;
            d1 = y2 * d1 - d2 + c[i];
            d2 = temp;
         }

         return x * d1 - d2 + c[0];
      }


      // implementation of Chebyshev Polynomial class
      // which can be used for building TF1 classes
      class ChebyshevPol {
      public:
         ChebyshevPol(unsigned int n) : fOrder(n) {}

         double operator() (const double *x, const double * coeff) {
            return ChebyshevN(fOrder, x[0], coeff);
         }
      private:
         unsigned int fOrder;
      };



   } // end namespace Math

} // end namespace ROOT



#endif  // ROOT_Math_Chebyshev
 ChebyshevPol.h:1
 ChebyshevPol.h:2
 ChebyshevPol.h:3
 ChebyshevPol.h:4
 ChebyshevPol.h:5
 ChebyshevPol.h:6
 ChebyshevPol.h:7
 ChebyshevPol.h:8
 ChebyshevPol.h:9
 ChebyshevPol.h:10
 ChebyshevPol.h:11
 ChebyshevPol.h:12
 ChebyshevPol.h:13
 ChebyshevPol.h:14
 ChebyshevPol.h:15
 ChebyshevPol.h:16
 ChebyshevPol.h:17
 ChebyshevPol.h:18
 ChebyshevPol.h:19
 ChebyshevPol.h:20
 ChebyshevPol.h:21
 ChebyshevPol.h:22
 ChebyshevPol.h:23
 ChebyshevPol.h:24
 ChebyshevPol.h:25
 ChebyshevPol.h:26
 ChebyshevPol.h:27
 ChebyshevPol.h:28
 ChebyshevPol.h:29
 ChebyshevPol.h:30
 ChebyshevPol.h:31
 ChebyshevPol.h:32
 ChebyshevPol.h:33
 ChebyshevPol.h:34
 ChebyshevPol.h:35
 ChebyshevPol.h:36
 ChebyshevPol.h:37
 ChebyshevPol.h:38
 ChebyshevPol.h:39
 ChebyshevPol.h:40
 ChebyshevPol.h:41
 ChebyshevPol.h:42
 ChebyshevPol.h:43
 ChebyshevPol.h:44
 ChebyshevPol.h:45
 ChebyshevPol.h:46
 ChebyshevPol.h:47
 ChebyshevPol.h:48
 ChebyshevPol.h:49
 ChebyshevPol.h:50
 ChebyshevPol.h:51
 ChebyshevPol.h:52
 ChebyshevPol.h:53
 ChebyshevPol.h:54
 ChebyshevPol.h:55
 ChebyshevPol.h:56
 ChebyshevPol.h:57
 ChebyshevPol.h:58
 ChebyshevPol.h:59
 ChebyshevPol.h:60
 ChebyshevPol.h:61
 ChebyshevPol.h:62
 ChebyshevPol.h:63
 ChebyshevPol.h:64
 ChebyshevPol.h:65
 ChebyshevPol.h:66
 ChebyshevPol.h:67
 ChebyshevPol.h:68
 ChebyshevPol.h:69
 ChebyshevPol.h:70
 ChebyshevPol.h:71
 ChebyshevPol.h:72
 ChebyshevPol.h:73
 ChebyshevPol.h:74
 ChebyshevPol.h:75
 ChebyshevPol.h:76
 ChebyshevPol.h:77
 ChebyshevPol.h:78
 ChebyshevPol.h:79
 ChebyshevPol.h:80
 ChebyshevPol.h:81
 ChebyshevPol.h:82
 ChebyshevPol.h:83
 ChebyshevPol.h:84
 ChebyshevPol.h:85
 ChebyshevPol.h:86
 ChebyshevPol.h:87
 ChebyshevPol.h:88
 ChebyshevPol.h:89
 ChebyshevPol.h:90
 ChebyshevPol.h:91
 ChebyshevPol.h:92
 ChebyshevPol.h:93
 ChebyshevPol.h:94
 ChebyshevPol.h:95
 ChebyshevPol.h:96
 ChebyshevPol.h:97
 ChebyshevPol.h:98
 ChebyshevPol.h:99
 ChebyshevPol.h:100
 ChebyshevPol.h:101
 ChebyshevPol.h:102
 ChebyshevPol.h:103
 ChebyshevPol.h:104
 ChebyshevPol.h:105
 ChebyshevPol.h:106
 ChebyshevPol.h:107
 ChebyshevPol.h:108
 ChebyshevPol.h:109
 ChebyshevPol.h:110
 ChebyshevPol.h:111
 ChebyshevPol.h:112
 ChebyshevPol.h:113
 ChebyshevPol.h:114
 ChebyshevPol.h:115
 ChebyshevPol.h:116
 ChebyshevPol.h:117
 ChebyshevPol.h:118
 ChebyshevPol.h:119
 ChebyshevPol.h:120
 ChebyshevPol.h:121
 ChebyshevPol.h:122
 ChebyshevPol.h:123
 ChebyshevPol.h:124
 ChebyshevPol.h:125
 ChebyshevPol.h:126
 ChebyshevPol.h:127
 ChebyshevPol.h:128
 ChebyshevPol.h:129
 ChebyshevPol.h:130
 ChebyshevPol.h:131
 ChebyshevPol.h:132
 ChebyshevPol.h:133