ClassicalOrthoPolys1D.hh

Go to the documentation of this file.

#ifndef NPSTAT_CLASSICALORTHOPOLYS1D_HH_
#define NPSTAT_CLASSICALORTHOPOLYS1D_HH_
 
/*!
// \file ClassicalOrthoPolys1D.hh
//
// \brief Orthonormal versions of some classical orthogonal polynomials
//
// Author: I. Volobouev
//
// July 2018
*/
 
#include <cmath>
#include <vector>
#include <algorithm>
 
#include "npstat/nm/AbsClassicalOrthoPoly1D.hh"
 
namespace npstat {
    /** Orthonormal Legendre polynomials on [-1, 1] */
    class LegendreOrthoPoly1D : public AbsClassicalOrthoPoly1D
    {
    public:
        inline virtual LegendreOrthoPoly1D* clone() const
            {return new LegendreOrthoPoly1D(*this);}
 
        inline virtual ~LegendreOrthoPoly1D() {}
 
        inline long double weight(const long double x) const
            {return (x >= -1.0L && x <= 1.0L) ? 0.5L : 0.0L;}
        inline double xmin() const {return -1.0L;}
        inline double xmax() const {return 1.0L;}
 
        /**
        // Derive the series coefficients for the integral of the
        // argument series from -1 to x. The array of coefficients
        // must be at least degree+1 long, and the buffer for the
        // integral coefficients must be at least degree+2 long.
        */
        void integralCoeffs(const double* coeffs, unsigned degree,
                            double* integCoeffs) const;
        /**
        // Derive the series coefficients for the derivative of the
        // argument series. The array of coefficients must be at least
        // degree+1 long, and the buffer for the derivative coefficients
        // must be at least degree long.
        */
        void derivativeCoeffs(const double* coeffs, unsigned degree,
                              double* derivCoeffs) const;
        /**
        // Derive the series coefficients for the monomial
        // expansion. The arrays of coefficients must be
        // at least degree+1 long.
        */
        void monomialCoeffs(const double* coeffs, unsigned degree,
                            long double* monoCoeffs) const;
        /**
        // Expand a squared polynomial into polynomial series.
        // The number of coefficients will be 2*degree + 1.
        */
        void squaredPolyCoeffs(unsigned degree,
                               long double* expansionCoeffs) const;
 
        inline virtual std::pair<long double,long double>
        monicRecurrenceCoeffs(const unsigned k) const
        {
            long double beta = 1.0L;
            if (k)
                beta = 1.0L/(4.0L - 1.0L/k/k);
            return std::pair<long double,long double>(0.0L, beta);
        }
 
        inline virtual std::pair<long double,long double>
        recurrenceCoeffs(const unsigned k) const
        {
            long double sqb = 1.0L;
            if (k)
                sqb = 1.0L/sqrtl(4.0L - 1.0L/k/k);
            return std::pair<long double,long double>(0.0L, sqb);
        }
 
#ifdef SWIG
    public:
        inline std::vector<double>
        integralCoeffs2(const double *c, const unsigned degreep1) const
        {
            std::vector<double> tmp(degreep1+1U);
            integralCoeffs(c, degreep1-1U, &tmp[0]);
            return tmp;
        }
 
        inline std::vector<double>
        derivativeCoeffs2(const double *c, const unsigned degreep1) const
        {
            std::vector<double> tmp(std::max(degreep1-1U, 1U));
            derivativeCoeffs(c, degreep1-1U, &tmp[0]);
            return tmp;
        }
 
        inline std::vector<double>
        monomialCoeffs2(const double *c, const unsigned degreep1) const
        {
            std::vector<long double> ldtmp(degreep1+1U);
            monomialCoeffs(c, degreep1-1U, &ldtmp[0]);
            return std::vector<double>(ldtmp.begin(), ldtmp.end());
        }
 
        inline std::vector<double>
        squaredPolyCoeffs2(const unsigned degree) const
        {
            std::vector<long double> tmp(2U*degree + 1U);
            squaredPolyCoeffs(degree, &tmp[0]);
            return std::vector<double>(tmp.begin(), tmp.end());
        }
#endif
    };
 
    /** Orthonormal Legendre polynomials on [0, 1] (that is, shifted) */
    class ShiftedLegendreOrthoPoly1D : public AbsClassicalOrthoPoly1D
    {
    public:
        inline virtual ShiftedLegendreOrthoPoly1D* clone() const
            {return new ShiftedLegendreOrthoPoly1D(*this);}
 
        inline virtual ~ShiftedLegendreOrthoPoly1D() {}
 
        inline long double weight(const long double x) const
            {return (x >= 0.0L && x <= 1.0L) ? 1.0L : 0.0L;}
        inline double xmin() const {return 0.0L;}
        inline double xmax() const {return 1.0L;}
 
        /**
        // Derive the series coefficients for the integral of the
        // argument series from 0 to x. The array of coefficients
        // must be at least degree+1 long, and the buffer for the
        // integral coefficients must be at least degree+2 long.
        */
        void integralCoeffs(const double* coeffs, unsigned degree,
                            double* integCoeffs) const;
        /**
        // Derive the series coefficients for the derivative of the
        // argument series. The array of coefficients must be at least
        // degree+1 long, and the buffer for the derivative coefficients
        // must be at least degree long.
        */
        void derivativeCoeffs(const double* coeffs, unsigned degree,
                              double* derivCoeffs) const;
        /**
        // Expand a squared polynomial into polynomial series.
        // The number of coefficients will be 2*degree + 1.
        */
        void squaredPolyCoeffs(unsigned degree,
                               long double* expansionCoeffs) const;
 
        inline virtual std::pair<long double,long double>
        monicRecurrenceCoeffs(const unsigned k) const
        {
            long double beta = 1.0L;
            if (k)
                beta = 0.25L/(4.0L - 1.0L/k/k);
            return std::pair<long double,long double>(0.5L, beta);
        }
 
        inline virtual std::pair<long double,long double>
        recurrenceCoeffs(const unsigned k) const
        {
            long double sqb = 1.0L;
            if (k)
                sqb = 0.5L/sqrtl(4.0L - 1.0L/k/k);
            return std::pair<long double,long double>(0.5L, sqb);
        }
 
#ifdef SWIG
    public:
        inline std::vector<double>
        integralCoeffs2(const double *c, const unsigned degreep1) const
        {
            std::vector<double> tmp(degreep1+1U);
            integralCoeffs(c, degreep1-1U, &tmp[0]);
            return tmp;
        }
 
        inline std::vector<double>
        derivativeCoeffs2(const double *c, const unsigned degreep1) const
        {
            std::vector<double> tmp(std::max(degreep1-1U, 1U));
            derivativeCoeffs(c, degreep1-1U, &tmp[0]);
            return tmp;
        }
 
        inline std::vector<double>
        squaredPolyCoeffs2(const unsigned degree) const
        {
            std::vector<long double> tmp(2U*degree + 1U);
            squaredPolyCoeffs(degree, &tmp[0]);
            return std::vector<double>(tmp.begin(), tmp.end());
        }
#endif
    };
 
    /**
    // The weight function for Jacobi polynomials is proportional
    // to (1 - x)^alpha (1 + x)^beta
    */
    class JacobiOrthoPoly1D : public AbsClassicalOrthoPoly1D
    {
    public:
        JacobiOrthoPoly1D(double alpha, double beta);
 
        inline virtual JacobiOrthoPoly1D* clone() const
            {return new JacobiOrthoPoly1D(*this);}
 
        inline virtual ~JacobiOrthoPoly1D() {}
 
        inline double alpha() const {return alpha_;}
        inline double beta() const {return beta_;}
 
        long double weight(long double x) const;
        inline double xmin() const {return -1.0L;}
        inline double xmax() const {return 1.0L;}
 
        /**
        // The following function generates coefficients for derivatives
        // of the Jacobi polynomial series. The derivative series should be
        // evaluated using a JacobiOrthoPoly1D(alpha_+1, beta_+1) object and
        // with maxdeg argument decreased by 1. The array "coeffs" should
        // have at least maxdeg+1 elements, while the "derivativeCoeffs"
        // buffer should have at least maxdeg elements.
        */
        void derivativeCoeffs(const double* coeffs,
                              unsigned maxdeg,
                              double* derivCoeffs) const;
 
        virtual std::pair<long double,long double>
        monicRecurrenceCoeffs(unsigned k) const;
 
        virtual std::pair<long double,long double>
        recurrenceCoeffs(unsigned k) const;
 
    private:
        double alpha_;
        double beta_;
        long double norm_;
 
#ifdef SWIG
    public:
        inline std::vector<double>
        derivativeCoeffs2(const double *c, const unsigned degreep1) const
        {
            std::vector<double> tmp(degreep1-1U);
            derivativeCoeffs(c, degreep1-1U,
                             tmp.empty() ? (double *)0 : &tmp[0]);
            return tmp;
        }
#endif
    };
 
    /** Orthonormal Chebyshev polynomials of the 1st kind */
    class ChebyshevOrthoPoly1st : public AbsClassicalOrthoPoly1D
    {
    public:
        inline virtual ChebyshevOrthoPoly1st* clone() const
            {return new ChebyshevOrthoPoly1st(*this);}
 
        inline virtual ~ChebyshevOrthoPoly1st() {}
 
        inline long double weight(const long double x) const
        {
            const long double Pi = 3.14159265358979323846264338L;
            return (x > -1.0L && x < 1.0L) ? 1.0/sqrtl(1.0L - x*x)/Pi : 0.0L;
        }
        inline double xmin() const {return -1.0L;}
        inline double xmax() const {return 1.0L;}
 
        inline virtual std::pair<long double,long double>
        monicRecurrenceCoeffs(const unsigned k) const
        {
            long double beta = 0.25L;
            if (k == 0U)
                beta = 3.14159265358979323846264338328L;
            else if (k == 1U)
                beta = 0.5L;
            return std::pair<long double,long double>(0.0L, beta);
        }
 
        inline virtual std::pair<long double,long double>
        recurrenceCoeffs(const unsigned k) const
        {
            long double sqb = 0.5L;
            if (k == 0U)
                sqb = 1.772453850905516027298167483L; // sqrt(Pi)
            else if (k == 1U)
                sqb = 0.7071067811865475244008443621L; // 1/sqrt(2)
            return std::pair<long double,long double>(0.0L, sqb);
        }
    };
 
    /** Orthonormal Chebyshev polynomials of the 2nd kind */
    class ChebyshevOrthoPoly2nd : public AbsClassicalOrthoPoly1D
    {
    public:
        inline virtual ChebyshevOrthoPoly2nd* clone() const
            {return new ChebyshevOrthoPoly2nd(*this);}
 
        inline virtual ~ChebyshevOrthoPoly2nd() {}
 
        inline long double weight(const long double x) const
        {
            const long double PiOver2 = 1.57079632679489661923132169L;
            return (x > -1.0L && x < 1.0L) ? sqrtl(1.0L - x*x)/PiOver2 : 0.0L;
        }
        inline double xmin() const {return -1.0L;}
        inline double xmax() const {return 1.0L;}
 
        inline virtual std::pair<long double,long double>
        monicRecurrenceCoeffs(const unsigned k) const
        {
            long double beta = 0.25L;
            if (!k)
                beta = 1.57079632679489661923132169164L; // Pi/2
            return std::pair<long double,long double>(0.0L, beta);
        }
 
        inline virtual std::pair<long double,long double>
        recurrenceCoeffs(const unsigned k) const
        {
            long double sqb = 0.5L;
            if (!k)
                sqb = 1.25331413731550025120788264L; // sqrt(Pi/2)
            return std::pair<long double,long double>(0.0L, sqb);
        }
    };
 
    /** Orthonormal(!) probabilist's Hermite polynomials */
    class HermiteProbOrthoPoly1D : public AbsClassicalOrthoPoly1D
    {
    public:
        inline virtual HermiteProbOrthoPoly1D* clone() const
            {return new HermiteProbOrthoPoly1D(*this);}
 
        inline virtual ~HermiteProbOrthoPoly1D() {}
 
        long double weight(long double x) const;
        inline double xmin() const {return -xmax();}
        inline double xmax() const {return 150.0;}
 
        virtual std::pair<long double,long double>
        monicRecurrenceCoeffs(unsigned k) const;
 
        virtual std::pair<long double,long double>
        recurrenceCoeffs(unsigned k) const;
    };
}
 
#endif // NPSTAT_CLASSICALORTHOPOLYS1D_HH_