#include "chebyshev.h"

#include "sbasis.h"

#include "sbasis-poly.h"

#include <vector>

using std::vector;

#include <gsl/gsl_math.h>

#include <gsl/gsl_chebyshev.h>

namespace Geom{

SBasis cheb(unsigned n) {

static std::vector<SBasis> basis;

if(basis.empty()) {

basis.push_back(Linear(1,1));

basis.push_back(Linear(0,1));

}

for(unsigned i = basis.size(); i <= n; i++) {

basis.push_back(Linear(0,2)*basis[i-1] - basis[i-2]);

}

return basis[n];

}

SBasis cheb_series(unsigned n, double* cheb_coeff) {

SBasis r;

for(unsigned i = 0; i < n; i++) {

double cof = cheb_coeff[i];

//if(i == 0)

//cof /= 2;

r += cheb(i)*cof;

}

return r;

}

SBasis clenshaw_series(unsigned m, double* cheb_coeff) {

/** b_n = a_n

double a = -1, b = 1;

SBasis d, dd;

SBasis y = (Linear(0, 2) - (a+b)) / (b-a);

SBasis y2 = 2*y;

for(int j = m - 1; j >= 1; j--) {

SBasis sv = d;

d = y2*d - dd + cheb_coeff[j];

dd = sv;

}

return y*d - dd + 0.5*cheb_coeff[0];

}

SBasis chebyshev_approximant (double (*f)(double,void*), int order, Interval in, void* p) {

gsl_cheb_series *cs = gsl_cheb_alloc (order+2);

gsl_function F;

F.function = f;

F.params = p;

gsl_cheb_init (cs, &F, in[0], in[1]);

SBasis r = compose(clenshaw_series(order, cs->c), Linear(-1,1));

gsl_cheb_free (cs);

return r;

}

struct wrap {

double (*f)(double,void*);

void* pp;

double fa, fb;

Interval in;

};

double f_interp(double x, void* p) {

struct wrap *wr = (struct wrap *)p;

double z = (x - wr->in[0]) / (wr->in[1] - wr->in[0]);

return (wr->f)(x, wr->pp) - ((1 - z)*wr->fa + z*wr->fb);

}

SBasis chebyshev_approximant_interpolating (double (*f)(double,void*),

int order, Interval in, void* p) {

double fa = f(in[0], p);

double fb = f(in[1], p);

struct wrap wr;

wr.fa = fa;

wr.fb = fb;

wr.in = in;

printf("%f %f\n", fa, fb);

wr.f = f;

wr.pp = p;

return compose(Linear(in[0], in[1]), Linear(fa, fb)) + chebyshev_approximant(f_interp, order, in, &wr) + Linear(fa, fb);

}

SBasis chebyshev(unsigned n) {

static std::vector<SBasis> basis;

if(basis.empty()) {

basis.push_back(Linear(1,1));

basis.push_back(Linear(0,1));

}

for(unsigned i = basis.size(); i <= n; i++) {

basis.push_back(Linear(0,2)*basis[i-1] - basis[i-2]);

}

return basis[n];

}

};