AbsIntervalQuadrature1D.hh

Go to the documentation of this file.

#ifndef NPSTAT_ABSINTERVALQUADRATURE1D_HH_
#define NPSTAT_ABSINTERVALQUADRATURE1D_HH_
 
/*!
// \file AbsIntervalQuadrature1D.hh
//
// \brief Base class for quadratures on the [-1, 1] interval
//
// Author: I. Volobouev
//
// November 2020
*/
 
#include <vector>
#include <utility>
 
namespace npstat {
    class AbsIntervalQuadrature1D
    {
    public:
        inline virtual ~AbsIntervalQuadrature1D() {}
 
        virtual unsigned npoints() const = 0;
        virtual void getAllAbscissae(long double* absc, unsigned len) const = 0;
        virtual void getAllWeights(long double* weights, unsigned len) const = 0;
 
        //@{
        /**
        // Weighted integration points on the given interval, suitable
        // for constructing orthogonal polynomials w.r.t. the given weight
        // function (in particular, by ContOrthoPoly1D class). Naturally,
        // rule with "npoints" points must be able to calculate polynomial
        // normalization integrals exactly.
        */
        template <class Functor>
        inline std::vector<std::pair<double,double> > weightedIntegrationPoints(
            const Functor& weight, const long double left,
            const long double right) const
        {
            return weightedIntegrationPointsHlp<Functor,std::pair<double,double> >(
                weight, left, right);
        }
 
        template <class Functor>
        inline std::vector<std::pair<long double,long double> > weightedIntegrationPointsLD(
            const Functor& weight, const long double left,
            const long double right) const
        {
            return weightedIntegrationPointsHlp<Functor,std::pair<long double,long double> >(
                weight, left, right);
        }
        //@}
 
    private:
        template <class Functor, class Pair>
        inline std::vector<Pair> weightedIntegrationPointsHlp(
            const Functor& weight, const long double left,
            const long double right) const
        {
            const unsigned npt = this->npoints();
            std::vector<long double> buf(npt*2U);
            long double* absc = &buf[0];
            long double* w = absc + npt;
            this->getAllAbscissae(absc, npt);
            this->getAllWeights(w, npt);
            const long double midpoint = (left + right)/2.0L;
            const long double unit = (right - left)/2.0L;
 
            std::vector<Pair> result(npt);
            for (unsigned i=0; i<npt; ++i)
            {
                const long double x = midpoint + unit*absc[i];
                result[i] = Pair(x, unit*w[i]*weight(x));
            }
 
            return result;
        }
    };
}
 
#endif // NPSTAT_ABSINTERVALQUADRATURE1D_HH_