FejerQuadrature.hh

Go to the documentation of this file.

#ifndef NPSTAT_FEJERQUADRATURE_HH_
#define NPSTAT_FEJERQUADRATURE_HH_
 
/*!
// \file FejerQuadrature.hh
//
// \brief Fejer quadratures in long double precision
//
// Author: I. Volobouev
//
// July 2018
*/
 
#include <map>
#include <vector>
#include <utility>
#include <memory>
 
#include "npstat/nm/SimpleFunctors.hh"
#include "npstat/nm/AbsIntervalQuadrature1D.hh"
 
namespace npstat {
    /** 
    // Fejer quadrature. Internally, operations are performed
    // in long double precision.
    */
    class FejerQuadrature : public AbsIntervalQuadrature1D
    {
    public:
        enum {interval_quadrature = 1};
 
        explicit FejerQuadrature(unsigned npoints);
 
        inline virtual ~FejerQuadrature() {}
 
        /** Return the number of quadrature points */
        inline unsigned npoints() const {return npoints_;}
 
        /**
        // The abscissae of the integration rule.
        // The buffer length should be at least npoints.
        */
        void getAllAbscissae(long double* abscissae, unsigned len) const;
 
        /**
        // The weights of the integration rule.
        // The buffer length should be at least npoints.
        */
        void getAllWeights(long double* weights, unsigned len) const;
 
        /** Perform the quadrature on [a, b] interval */
        template <typename FcnResult, typename FcnArg>
        long double integrate(const Functor1<FcnResult,FcnArg>& fcn,
                              long double a, long double b) const;
 
        /**
        // Minimum number of points which integrates a polynomial
        // with the given degree exactly
        */
        static unsigned minimalExactRule(unsigned polyDegree);
 
    private:
        FejerQuadrature();
 
        std::vector<long double> a_;
        std::vector<long double> w_;
        mutable std::vector<std::pair<long double, long double> > buf_;
        unsigned npoints_;
 
#ifdef SWIG
    public:
        inline std::vector<double> abscissae2() const
        {
            return std::vector<double>(a_.begin(), a_.end());
        }
 
        inline std::vector<double> weights2() const
        {
            return std::vector<double>(w_.begin(), w_.end());
        }
 
        template <typename FcnResult, typename FcnArg>
        inline double integrate2(const Functor1<FcnResult,FcnArg>& fcn,
                                 const double a, const double b) const
        {
            return integrate(fcn, a, b);
        }
#endif // SWIG
    };
 
    /**
    // Static collection of Fejer quadratures. Call FejerQuadratureSet::quad
    // function in order to get the quadrature. The quadrature object will be
    // either recalled or created if it does not yet exist.
    */
    class FejerQuadratureSet
    {
    public:
        static std::shared_ptr<const FejerQuadrature> quad(unsigned npoints);
 
    private:
        typedef std::map<unsigned,std::shared_ptr<const FejerQuadrature> > FejerMap;
        static FejerMap& getFejerMap();
    };
}
 
#include "npstat/nm/FejerQuadrature.icc"
 
#endif // NPSTAT_FEJERQUADRATURE_HH_