polynomial.hpp

Go to the documentation of this file.

/*****************************************************************************
** Ifdefs
*****************************************************************************/

#ifndef ECL_GEOMETRY_POLYNOMIALS_HPP_
#define ECL_GEOMETRY_POLYNOMIALS_HPP_

/*****************************************************************************
** Includes
*****************************************************************************/

#include <cmath>
#include "cartesian_point.hpp"
#include "function_math.hpp"
#include "pascals_triangle.hpp"
#include <ecl/config/macros.hpp>
#include <ecl/concepts/macros.hpp>
#include <ecl/concepts/streams.hpp>
#include <ecl/containers/array.hpp>
#include <ecl/exceptions/standard_exception.hpp>
#include <ecl/errors/compile_time_assert.hpp>
#include <ecl/utilities/blueprints.hpp>
#include <ecl/formatters/floats.hpp>
#include <iostream>
#include "macros.hpp"

/*****************************************************************************
** Namespaces
*****************************************************************************/

namespace ecl {

/*****************************************************************************
** Forward declarations
*****************************************************************************/

template <unsigned int N> class Polynomial;

/*****************************************************************************
** BluePrintFactory
*****************************************************************************/

template<unsigned int N>
class BluePrintFactory< Polynomial<N> > {
public:
        BluePrintFactory() {};
        ~BluePrintFactory() {};
};

/*****************************************************************************
** Interface [Polynomial]
*****************************************************************************/
template <unsigned int N>
class ecl_geometry_PUBLIC Polynomial : public BluePrintFactory< Polynomial<N> >, public FunctionMath< Polynomial<N> > {
    public:
        /*********************
        ** Convenience
        **********************/
        typedef Array<double,N+1> Coefficients; 
        /*********************
        ** C&D
        **********************/
        Polynomial() : coeff(Coefficients::Constant(0.0)) {};

        template<typename Derived>
        Polynomial(const BluePrint< Derived > &blueprint) {
            blueprint.implementApply(*this);
        }
        virtual ~Polynomial() {}

        /*********************
        ** Transform
        **********************/
        void shift_horizontal(const double &shift);
//        void stretch_horizontal(double multiplier);

        /*********************
        ** Derivatives
        **********************/
        Polynomial<N-1> derivative() const {
            Polynomial<N-1> derivative_polynomial;
            typename Polynomial<N-1>::Coefficients &derivative_coefficients = derivative_polynomial.coefficients();
            for ( unsigned int i = 0; i < N; ++i ) {
                derivative_coefficients[i] = (i+1)*coeff[i+1];
            }
            return derivative_polynomial;
        }
        double derivative(const double &x) const;
        double dderivative(const double &x) const;

        /*********************
        ** Accessors
        **********************/
        Coefficients& coefficients() { return coeff; };
        const Coefficients& coefficients() const { return coeff; };

        double operator ()(const double &x) const;

        template <typename OutputStream, unsigned int Degree>
        friend OutputStream& operator << (OutputStream &ostream, const Polynomial<Degree> &polynomial);

    private:
        Coefficients coeff;
};

template <>
class ECL_PUBLIC Polynomial<0> {
    public:
        /*********************
        ** Convenience
        **********************/
        typedef Array<double,1> Coefficients; 
        /*********************
        ** C&D
        **********************/
        Polynomial() : coeff(Coefficients::Constant(0.0)) {};
        virtual ~Polynomial() {};

        /*********************
        ** Transform
        **********************/
        void shift_horizontal(const double& /* shift */) {}; // Does nothing.

        /*********************
        ** Derivatives
        **********************/
        Polynomial<0> derivative() const {
            return Polynomial<0>();
        }
        double derivative(const double& /* x */ ) const {
            return 0.0;
        };
        double dderivative(const double& /* x */) const {
            return 0.0;
        };

        /*********************
        ** Accessors
        **********************/
        Coefficients& coefficients() { return coeff; };
        const Coefficients& coefficients() const { return coeff; };

        double operator ()(const double& /* x */) const {
            return coeff[0];
        };

    private:
        Coefficients coeff;
};

/*****************************************************************************
** Typedefs
*****************************************************************************/

typedef Polynomial<1> LinearFunction;            
typedef Polynomial<2> QuadraticPolynomial;   
typedef Polynomial<3> CubicPolynomial;       
typedef Polynomial<5> QuinticPolynomial;     
/*****************************************************************************
 * Implementation [Polynomial - Transform]
 ****************************************************************************/
template <unsigned int N>
void Polynomial<N>::shift_horizontal(const double &shift)
{
    PascalsTriangle<N> pascals_triangle;
    typename PascalsTriangle<N>::const_iterator iter;
    double tmp;

    for ( unsigned int i = 0; i < N; ++i ) {
        tmp = -1*shift;
        int j = i+1;
        for ( iter = (pascals_triangle.begin(i)+1); iter != pascals_triangle.end(i); ++iter ) { // skip the first one
            coeff[i] += (*iter)*tmp*coeff[j];
            tmp *= (-1*shift);
            ++j;
        }
    }
}
/*****************************************************************************
 * Implementation [Polynomial - Access]
 ****************************************************************************/
template <unsigned int N>
double Polynomial<N>::operator()(const double &x) const
{
    double tmp = x;
    double value = coeff[0];
    for ( unsigned int i = 1; i <= N; ++i ) {
         value += coeff[i]*tmp;
         tmp *= x;
    }
    return value;
}
template <unsigned int N>
double Polynomial<N>::derivative(const double &x) const {

    if ( N > 0 ) {
        return derivative()(x);
    } else {
        return 0.0;
    }
}

template <unsigned int N>
double Polynomial<N>::dderivative(const double &x) const {

    if ( N > 1 ) {
        return derivative().derivative()(x);
    } else {
        return 0.0;
    }
}

/*****************************************************************************
** Streaming Operator [Polynomial]
*****************************************************************************/

template <typename OutputStream, unsigned int Degree>
OutputStream& operator << (OutputStream &ostream, const Polynomial<Degree> &polynomial)
{
        ecl_compile_time_concept_check(StreamConcept<OutputStream>);

        Format<double> format;
    format.precision(2);

    ostream << format(polynomial.coeff[0]);
    for (unsigned int i = 1; i <= Degree; ++i) {
        ostream << " + " << format(polynomial.coeff[i]) << "x^" << i;
    }
    ostream.flush();

    return ostream;
}

/*****************************************************************************
** BluePrints
*****************************************************************************/

namespace blueprints {

/*****************************************************************************
** Interface [LinearInterpolation]
*****************************************************************************/

class ecl_geometry_PUBLIC LinearInterpolation : public ecl::BluePrint<LinearInterpolation> {
    public:
        typedef ecl::LinearFunction base_type;
        LinearInterpolation(const double x_i, const double y_i, const double x_f, const double y_f) :
                    x_initial(x_i),
                    y_initial(y_i),
                    x_final(x_f),
                    y_final(y_f)
        {}

        virtual ~LinearInterpolation() {}

        ecl::LinearFunction instantiate();

        void apply(ecl::LinearFunction& function) const;

    private:
        double x_initial, y_initial;
        double x_final, y_final;
};

/*****************************************************************************
** Interface [LinearPointSlope]
*****************************************************************************/
class ecl_geometry_PUBLIC LinearPointSlopeForm : public ecl::BluePrint<LinearPointSlopeForm> {
    public:
        typedef ecl::LinearFunction base_type;
        LinearPointSlopeForm(const double x_f, const double y_f, const double s) :
                    slope(s),
                    x_final(x_f),
                    y_final(y_f)
        {}

        virtual ~LinearPointSlopeForm() {}

        ecl::LinearFunction instantiate();

        void apply(ecl::LinearFunction& function) const;

    private:
        double slope;
        double x_final, y_final;
};


/*****************************************************************************
** Interface [CubicDerivativeInterpolation]
*****************************************************************************/

class ecl_geometry_PUBLIC CubicDerivativeInterpolation : public ecl::BluePrint<CubicDerivativeInterpolation> {
    public:
        typedef ecl::CubicPolynomial base_type;
        CubicDerivativeInterpolation(const double x_i, const double y_i, const double ydot_i,
                const double x_f, const double y_f, const double ydot_f) :
                    x_initial(x_i),
                    y_initial(y_i),
                    ydot_initial(ydot_i),
                    x_final(x_f),
                    y_final(y_f),
                    ydot_final(ydot_f)
        {}

        virtual ~CubicDerivativeInterpolation() {}
        ecl::CubicPolynomial instantiate();

        void apply(ecl::CubicPolynomial& polynomial) const;

    private:
        double x_initial, y_initial, ydot_initial;
        double x_final, y_final, ydot_final;
};
/*****************************************************************************
** Interface [CubicSecondDerivativeInterpolation]
*****************************************************************************/

class ecl_geometry_PUBLIC CubicSecondDerivativeInterpolation : public ecl::BluePrint<CubicSecondDerivativeInterpolation>  {
    public:
        typedef ecl::CubicPolynomial base_type;
        CubicSecondDerivativeInterpolation(const double x_i, const double y_i, const double yddot_i,
                const double x_f, const double y_f, const double yddot_f) :
                    x_initial(x_i),
                    y_initial(y_i),
                    yddot_initial(yddot_i),
                    x_final(x_f),
                    y_final(y_f),
                    yddot_final(yddot_f)
        {}

        virtual ~CubicSecondDerivativeInterpolation() {}

        base_type instantiate();

        void apply(base_type& polynomial) const;

    private:
        double x_initial, y_initial, yddot_initial;
        double x_final, y_final, yddot_final;
};

/*****************************************************************************
** Interface [QuinticInterpolation]
*****************************************************************************/

class ecl_geometry_PUBLIC QuinticInterpolation : public ecl::BluePrint<QuinticInterpolation> {
    public:
        typedef ecl::QuinticPolynomial base_type;
        QuinticInterpolation(const double x_i, const double y_i, const double ydot_i, const double yddot_i,
                const double x_f, const double y_f, const double ydot_f, const double yddot_f) :
                    x_initial(x_i),
                    y_initial(y_i),
                    ydot_initial(ydot_i),
                    yddot_initial(yddot_i),
                    x_final(x_f),
                    y_final(y_f),
                    ydot_final(ydot_f),
                    yddot_final(yddot_f)
        {}

        virtual ~QuinticInterpolation() {}

        ecl::QuinticPolynomial instantiate();

        void apply(ecl::QuinticPolynomial& polynomial) const;

    private:
        double x_initial, y_initial, ydot_initial, yddot_initial;
        double x_final, y_final, ydot_final, yddot_final;
};

}; // namespace blueprints

/*****************************************************************************
** BluePrintFactory
*****************************************************************************/

using blueprints::LinearInterpolation;
using blueprints::LinearPointSlopeForm;
using blueprints::CubicDerivativeInterpolation;
using blueprints::CubicSecondDerivativeInterpolation;
using blueprints::QuinticInterpolation;

template<>
class BluePrintFactory< LinearFunction > {
public:
        static LinearInterpolation Interpolation(const double x_i, const double y_i, const double x_f, const double y_f);
        static LinearPointSlopeForm PointSlopeForm(const double x_f, const double y_f, const double slope);
        virtual ~BluePrintFactory() {}
};

template<>
class BluePrintFactory< CubicPolynomial > {
public:
        static CubicDerivativeInterpolation DerivativeInterpolation(const double x_i, const double y_i, const double ydot_i,
                        const double x_f, const double y_f, const double ydot_f);

        static CubicSecondDerivativeInterpolation SecondDerivativeInterpolation(const double x_i, const double y_i, const double yddot_i,
                        const double x_f, const double y_f, const double yddot_f);

        virtual ~BluePrintFactory() {}
};

template<>
class BluePrintFactory< QuinticPolynomial > {
    public:
        static QuinticInterpolation Interpolation(const double x_i, const double y_i, const double ydot_i, const double yddot_i,
                const double x_f, const double y_f, const double ydot_f, const double yddot_f);

        virtual ~BluePrintFactory() {}
};


/*****************************************************************************
** Interface [Maximum|Minimum][Polynomial]
*****************************************************************************/
template<>
class ecl_geometry_PUBLIC Maximum< LinearFunction > {
public:
        ECL_PUBLIC double operator()(const double& x_begin, const double& x_end, const LinearFunction &function);
        virtual ~Maximum() {}
};

template<>
class ecl_geometry_PUBLIC Maximum<CubicPolynomial> {
public:
        ECL_PUBLIC double operator()(const double& x_begin, const double& x_end, const CubicPolynomial& cubic);
        virtual ~Maximum() {}
};

template<>
class ecl_geometry_PUBLIC Minimum< LinearFunction > {
public:
        ECL_PUBLIC double operator()(const double& x_begin, const double& x_end, const LinearFunction &function);
        virtual ~Minimum() {}
};

template<>
class ecl_geometry_PUBLIC Minimum<CubicPolynomial> {
public:
        ECL_PUBLIC double operator()(const double& x_begin, const double& x_end, const CubicPolynomial& cubic);
        virtual ~Minimum() {}
};

/*****************************************************************************
** Interface [Intersection][LinearFunction]
*****************************************************************************/
template<>
class ecl_geometry_PUBLIC Intersection< LinearFunction > {
public:
        Intersection() : last_operation_failed(false) {}
        virtual ~Intersection() {}
        ECL_PUBLIC CartesianPoint2d operator()(const LinearFunction& f, const LinearFunction& g) ecl_throw_decl(StandardException);

        bool fail() const { return last_operation_failed; }

private:
        bool last_operation_failed;
};

/*****************************************************************************
** Interface [Synthetic Division]
*****************************************************************************/
template<>
class ecl_geometry_PUBLIC Division<QuadraticPolynomial> {
public:
        Division() {};
        virtual ~Division() {};
        ECL_PUBLIC LinearFunction operator()(const QuadraticPolynomial &p, const double &factor, double &remainder);
};

template<>
class ecl_geometry_PUBLIC Division<CubicPolynomial> {
public:
        Division() {};
        virtual ~Division() {};
        ECL_PUBLIC QuadraticPolynomial operator()(const CubicPolynomial &p, const double &factor, double &remainder);
};

/*****************************************************************************
** Interface [Roots]
*****************************************************************************/
template<>
class ecl_geometry_PUBLIC Roots<LinearFunction> {
public:
        Roots() {}
        virtual ~Roots() {}
        ECL_PUBLIC Array<double> operator()(const LinearFunction& p);
};

template<>
class ecl_geometry_PUBLIC Roots<QuadraticPolynomial> {
public:
        Roots() {}
        virtual ~Roots() {}
        ECL_PUBLIC Array<double> operator()(const QuadraticPolynomial& p);
};
template<>
class ecl_geometry_PUBLIC Roots<CubicPolynomial> {
public:
        Roots() {}
        virtual ~Roots() {}
        ECL_PUBLIC Array<double> operator()(const CubicPolynomial& p);
};

/*****************************************************************************
** Interfaces [FunctionMath]
*****************************************************************************/
template <>
class ECL_PUBLIC FunctionMath<LinearFunction> {
public:
        FunctionMath() {}; 
        virtual ~FunctionMath() {};
        static CartesianPoint2d Intersection(const LinearFunction &f, const LinearFunction &g) ecl_throw_decl(StandardException) {
                CartesianPoint2d point;
                ecl_try {
                        ecl::Intersection< LinearFunction > intersection;
                        point = intersection(f,g);
                } ecl_catch( const StandardException &e ) {
                        ecl_throw(StandardException(LOC,e));
                }
                return point;
        }
        static Array<double> Roots(const LinearFunction &function) {
                return ecl::Roots<LinearFunction>()(function);
        }
        static double Minimum(const double& x_begin, const double& x_end, const LinearFunction &function) {
                return ecl::Minimum<LinearFunction>()(x_begin, x_end, function);
        }
        static double Maximum(const double& x_begin, const double& x_end, const LinearFunction &function) {
                return ecl::Maximum<LinearFunction>()(x_begin, x_end, function);
        }
};
template <>
class ECL_PUBLIC FunctionMath<QuadraticPolynomial> {
public:
        FunctionMath() {}; 
        virtual ~FunctionMath() {};
        static Array<double> Roots(const QuadraticPolynomial &p) {
                return ecl::Roots<QuadraticPolynomial>()(p);
        }
        static LinearFunction Division(const QuadraticPolynomial &p, const double &factor, double &remainder) {
                LinearFunction f(ecl::Division<QuadraticPolynomial>()(p,factor,remainder));
                return f;
        }
};

template <>
class ECL_PUBLIC FunctionMath<CubicPolynomial> {
public:
        FunctionMath() {}; 
        virtual ~FunctionMath() {};
        static Array<double> Roots(const CubicPolynomial &p) {
                return ecl::Roots<CubicPolynomial>()(p);
        }
        static QuadraticPolynomial Division(const CubicPolynomial &p, const double &factor, double &remainder) {
                QuadraticPolynomial q(ecl::Division<CubicPolynomial>()(p,factor,remainder));
                return q;
        }
        static double Minimum(const double& x_begin, const double& x_end, const CubicPolynomial &function) {
                return ecl::Minimum<CubicPolynomial>()(x_begin, x_end, function);
        }
        static double Maximum(const double& x_begin, const double& x_end, const CubicPolynomial &function) {
                return ecl::Maximum<CubicPolynomial>()(x_begin, x_end, function);
        }
};

}; // namespace ecl

#endif /*ECL_GEOMETRY_POLYNOMIALS_HPP_*/