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); }
