#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include "Chebyshev.h"
#include "Exception.h"
#include "MathP.h"

Chebyshev::Chebyshev(double lower, double upper, int nCoeffs, 
		     double (*func)(double))
{
  _lower = lower;
  _upper = upper;
  _range = _upper - _lower;

  if (nCoeffs <= 0) {
    throw Exception("Chebyshev::Chebyshev() - nCoeffs must be positive");
  }

  _nCoeffs = nCoeffs;

  if ((_coeffs = (double *)malloc(_nCoeffs * sizeof(double))) == 0) {
    throw Exception("Chebyshev::Chebyshev() - malloc() of _coeffs failed");
  }

  double *vector;
  if ((vector = (double *)malloc(_nCoeffs * sizeof(double))) == 0) {
    throw Exception("Chebyshev::Chebyshev() - malloc() of vector failed");
  }

  double bma = 0.5 * (_upper - _lower);
  double bpa = 0.5 * (_upper + _lower);
  
  int j, k;

  for (k = 0; k < _nCoeffs; k++) {
    double y = cos(M_PI *(k + 0.5) / _nCoeffs);
    vector[k] = (*func)(y * bma + bpa);
  }

  double fac = 2.0 / _nCoeffs;
  for (j = 0; j < _nCoeffs; j++) {
    double sum = 0.0;
    for (k = 0; k < _nCoeffs; k++) {
      sum += vector[k] * cos(PI*j*(k+0.5)/_nCoeffs);
    } 
    _coeffs[j] = fac * sum;
  }

  free((void *)vector);
}


Chebyshev::~Chebyshev()
{
  free((void *)_coeffs);
}



double Chebyshev::evaluate(double x)
{
  double d = 0.0, dd = 0.0, sv, y, y2;
  int j;

  if ((x - _lower) * (x - _upper) > 0.0) {
    // Out of range
    throw Exception("Chebyshev::evaluate() - x out of range");
  }

  y = (2. * x + _range) / _range;
  y2 = 2.0 * y;
  for (j = _nCoeffs - 1; j >= 1; j--) {
    sv = d;
    d = y2 * d - dd + _coeffs[j];
    dd = sv;
  }

  // Last step is different
  return y * d - dd + 0.5 * _coeffs[0];
}