template < int n > struct Tchebychev { double operator () ( double x ) { return 2.0 * x * Tchebychev < n-1>()(x) - Tchebychev < n-2>()(x); } }; template < > struct Tchebychev < 0 > { double operator () ( double x ) { return 1.0; } }; template < > struct Tchebychev < 1 > { double operator () ( double x ) { return x; } }; double TchebychevN(int n,double x) { switch(n) { case 0 : return Tchebychev<0>()(x); case 1 : return Tchebychev<1>()(x); case 2 : return Tchebychev<2>()(x); case 3 : return Tchebychev<3>()(x); case 4 : return Tchebychev<4>()(x); case 5 : return Tchebychev<5>()(x); case 6 : return Tchebychev<6>()(x); case 7 : return Tchebychev<7>()(x); case 8 : return Tchebychev<8>()(x); case 9 : return Tchebychev<9>()(x); NO_DEFAULT_CASE } return 0; } template < typename functor > void FindTchebychevCoefficients(functor f,double lo = 0.0, double hi = 1.0, int N = (1<<20)) { double PI_over_N = PI/N; #define NUM_T 6 double t[NUM_T] = { 0 }; for(int k=0;k < N;k++) { double pk = PI_over_N * (k + 0.5); double x_k = cos(pk); double farg = lo + (hi - lo) * 0.5 * (x_k+1.0); double fval = f( farg ); for(int j=0;j < NUM_T;j++) { t[j] += fval * cos( pk * j ); } } for(int j=0;j < NUM_T;j++) { t[j] *= 1.0 / N; if ( j != 0 ) t[j] *= 2.0; //lprintfvar(t[j]); lprintf("t %d: %16.9f , %16.8g\n",j,t[j],t[j]); } double errSqr[NUM_T] = { 0 }; double errMax[NUM_T] = { 0 }; for(int i=0;i < N;i++) { double xunit = (i+0.5)/N; double farg = lo + (hi - lo) * xunit; double xt = xunit*2.0 - 1.0; double y = f(farg); double p = 0.0; for(int j=0;j < NUM_T;j++) { p += t[j] * TchebychevN(j,xt); errSqr[j] += fsquare(y-p); errMax[j] = MAX( errMax[j] , fabs(y-p) ); } } for(int j=0;j < NUM_T;j++) { lprintf("%d : errSqr = %g errMax = %g\n",j,errSqr[j]/N,errMax[j]); } }

You can also very easily find inner products by using the Integrate<> template I posted earlier plus this simple product adaptor :

template < typename f1, typename f2 > struct FunctorProduct { double operator () ( double x ) { return f1()(x) * f2()(x); } };

(eg. you can trivially find Hermite or Legendre series this way by doing inner products).

Another handy helper is a range remapper :

template < typename functor > struct RemapFmTo { double m_fmLo; double m_fmHi; double m_toLo; double m_toHi; RemapFmTo( double fmLo, double fmHi, double toLo, double toHi ) : m_fmLo(fmLo), m_fmHi(fmHi), m_toLo(toLo), m_toHi(toHi) { } double operator () ( double x ) { double t = (x - m_fmLo) / (m_fmHi - m_fmLo); double y = t * (m_toHi - m_toLo) + m_toLo; return functor()(y); } }; template < typename functor > struct RemapUnitTo { double m_toLo; double m_toHi; RemapUnitTo( double toLo, double toHi ) : m_toLo(toLo), m_toHi(toHi) { } double operator () ( double x ) { double y = x * (m_toHi - m_toLo) + m_toLo; return functor()(y); } };

Now we can trivially do what we did before to find the optimal approximation in a known range :

struct SqrtOnePlusX { double operator () ( double x ) { return sqrt( 1 + x ); } }; RemapUnitTo< SqrtOnePlusX > tf(-0.075f,0.075f); FindTchebychevCoefficients( tf );and the output is :

t 0: 0.999647973 , 0.99964797 t 1: 0.037519818 , 0.037519818 t 2: -0.000352182 , -0.00035218223 t 3: 0.000006612 , 6.6121459e-006 t 4: -0.000000155 , -1.5518238e-007 t 5: 0.000000004 , 4.0790168e-009 0 : errSqr = 0.000469204 errMax = 0.0378787 1 : errSqr = 5.78832e-008 errMax = 0.000358952 2 : errSqr = 2.12381e-011 errMax = 6.77147e-006 3 : errSqr = 1.18518e-014 errMax = 1.59378e-007 4 : errSqr = 8.23726e-018 errMax = 4.19794e-009 5 : errSqr = 6.54835e-021 errMax = 1.19024e-010I guess Remez minimax polynomials are better, but this is so easy and it gives you a good starting point, then you can just numerically optimize after this anyway.

ADDENDUM :

obviously the TchebychevN dispatcher sucks and is silly, but you don't actually use it in practice; you know which polynomials you want and you use something like :

double approx = c0 * Tchebychev<0>(x) + c1 * Tchebychev<1>(x) + c2 * Tchebychev<2>(x) + c3 * Tchebychev<3>(x);which the compiler handles quite well.

I also did Legendre polynomials :

template < int n > struct Legendre { double operator () ( double x ) { double L_n1 = Legendre < n-1 >()(x); double L_n2 = Legendre < n-2 >()(x); const double B = (n-1)/(double)n; // n-1/n return (1+B)*x*L_n1 - B*L_n2; } }; template< > struct Legendre<0> { double operator () ( double x ) { return 1.0; } }; template< > struct Legendre<1> { double operator () ( double x ) { return x; } }; // but ryg's loop variants are mostly just better : __forceinline double legendre_ryg(int n, double x) { if (n == 0) return 1.0; if (n == 1) return x; double t = x, t1 = 1, t2; for (int i=1; i < n; i++) { t2 = t1; t1 = t; const double B = (i)/(double)(i+1); t = (1+B)*x*t1 - B*t2; } return t; } __forceinline double chebyshev_ryg(int n, double x) { if (n == 0) return 1.0; if (n == 1) return x; double t = x, t1 = 1, t2; for (int i=1; i < n; i++) { t2 = t1; t1 = t; t = 2.0 * x * t1 - t2; } return t; } // you can find Legendre coefficients like so : RemapUnitTo< Legendre<0> > ShiftedLegendre0(-1,1); double c0 = Integrate(0.0,1.0,steps,MakeFunctorProduct(functor(),ShiftedLegendre0)); // (don't forget normalization) // using new helper : template < typename f1, typename f2 > FunctorProduct< f1,f2 > MakeFunctorProduct( f1 _f1 , f2 _f2 ) { return FunctorProduct< f1,f2 >(_f1,_f2); }