cheby.h

Go to the documentation of this file.

#pragma once
/*
Chebyshev polynomial
*/
#include <stdexcept>
#include <iostream>
#include <limits>
#include <sstream>
#include <vector>

namespace tcc {

    class ChebyshevPolynomial {
    public:
        explicit ChebyshevPolynomial(
            std::vector<double> params, 
            double minX = -1,   
            double maxX = 1)    
        :
            params(params)
        {
            if (params.size() < 1) {
                throw std::runtime_error("ChebyshevPolynomial called with empty vector");
            }
            if (minX >= maxX) {
                std::ostringstream os;
                os << "ChebyshevPolynomial minX=" << minX << " >= " << maxX << "=maxX";
                throw std::runtime_error(os.str());
            }
            _initialize(minX, maxX);
        }
        
        explicit ChebyshevPolynomial()
        :
            params(1)
        {
            _initialize(-1.0, 1.0);
        }

        double getMinX() const { return _minX; };

        double getMaxX() const { return _maxX; };

        unsigned int getOrder() const { return this->params.size() - 1; };

        double operator() (double x) const {
            double xPrime = (x + _offset) * _scale;

            if ((xPrime < -1.0) || (xPrime > 1.0)) {
                return std::numeric_limits<double>::quiet_NaN();
            }
            
            // Clenshaw function for solving the Chebyshev polynomial
            // Non-recursive version from Kresimir Cosic
            int const order = getOrder();
            if (order == 0) {
                return this->params[0];
            } else if (order == 1) {
                return this->params[0] + (this->params[1] * xPrime);
            }
            double cshPrev = this->params[order];
            double csh = (2 * xPrime * this->params[order]) + this->params[order-1];
            for (int i = order - 2; i > 0; --i) {
                double cshNext = (2 * xPrime * csh) + this->params[i] - cshPrev;
                cshPrev = csh;
                csh = cshNext;
            }
            return (xPrime * csh) + this->params[0] - cshPrev;
        }

        std::string __repr__() const;
        
        std::vector<double> params;

    private:
        double _minX;    
        double _maxX;    
        double _scale;   
        double _offset;  

        void _initialize(double minX, double maxX) {
            _minX = minX;
            _maxX = maxX;
            _scale = 2.0 / (_maxX - _minX);
            _offset = -(_minX + _maxX) * 0.5;
        }
    };

    std::ostream &operator<<(std::ostream &os, ChebyshevPolynomial const &cheby);

}