Index: trunk/npstat/nm/SemiInfGaussianQuadrature.cc
===================================================================
--- trunk/npstat/nm/SemiInfGaussianQuadrature.cc	(revision 585)
+++ trunk/npstat/nm/SemiInfGaussianQuadrature.cc	(revision 586)
@@ -1,96 +1,115 @@
 #include <cassert>
 #include <stdexcept>
 
 #include "npstat/nm/SemiInfGaussianQuadrature.hh"
 
-static const unsigned allowed[] = {16U};
+static const unsigned allowed[] = {8U, 16U};
+
+static const long double x8[8] = {
+    0.07492311676455893922276557L, 0.3781584044761421077832691L, 
+    0.8715838973404047790594523L, 1.505071568300313929075073L, 
+    2.246981509078310537962837L, 3.088530913861651684254462L, 
+    4.04908276922991142368765L, 5.212801320529213982672174L
+};
+
+static const long double w8[8] = {
+    0.189659033149590363152723L, 0.3794769921770111393122674L, 
+    0.3902570380181665052619002L, 0.2226654966518453595846021L, 
+    0.0633767220309546943177233L, 0.007591407548118500187833945L, 
+    0.0002857611530656537853784407L, 1.686586748035605454239761e-6L
+};
 
 static const long double x16[16] = {
     0.02793589170072729925339991L, 0.1453843294578677844847097L, 
     0.3498731394391579314726867L, 0.6317238610158890058605988L, 
     0.980154254975464873819141L, 1.385086660726487105694736L, 
     1.838179686377704491754979L, 2.333246242708812302354253L, 
     2.866339577125563022190148L, 3.435761834337734010983741L, 
     4.042199289542951928362704L, 4.689205878459732580115888L, 
     5.384444355585162035569975L, 6.142787003373753071986793L, 
     6.995233726578337465407051L, 8.025687329032628646366107L
 };
 
 static const long double w16[16] = {
     0.07145261983854258606277078L, 0.1606667619613437734824571L, 
     0.2304055011905809966800292L, 0.2595970028857905512858743L, 
     0.2339655652545129730947059L, 0.164858397140373883445356L, 
     0.08768081108643750163479226L, 0.03382757775076844777847568L, 
     0.009064987912142623805170969L, 0.001606115711099596042573022L, 
     0.0001771814733827445172841926L, 0.00001124369887427931952844601L, 
     3.662823902026437858170459e-7L, 5.107501774601126771692312e-9L, 
     2.174604963637269706105677e-11L, 1.226717757951498517118791e-14L
 };
 
 namespace npstat {
     SemiInfGaussianQuadrature::SemiInfGaussianQuadrature(const unsigned npoints)
         : a_(0),
           w_(0),
           npoints_(npoints)
     {
         switch (npoints)
         {
+        case 8U:
+            a_ = x8;
+            w_ = w8;
+            break;
+
         case 16U:
             a_ = x16;
             w_ = w16;
             break;
 
         default:
             npoints_ = 0;
             throw std::invalid_argument(
                 "In npstat::SemiInfGaussianQuadrature constructor: "
                 "unsupported number of abscissae requested");
         }
     }
 
     std::vector<unsigned> SemiInfGaussianQuadrature::allowedNPonts()
     {
         std::vector<unsigned> v(allowed, allowed+sizeof(allowed)/sizeof(allowed[0]));
         return v;
     }
 
     unsigned SemiInfGaussianQuadrature::minimalExactRule(const unsigned polyDegree)
     {
         for (unsigned i=0; i<sizeof(allowed)/sizeof(allowed[0]); ++i)
             if (polyDegree < 2U*allowed[i])
                 return allowed[i];
         return 0U;
     }
 
     bool SemiInfGaussianQuadrature::isAllowed(const unsigned npoints)
     {
         for (unsigned i=0; i<sizeof(allowed)/sizeof(allowed[0]); ++i)
             if (npoints == allowed[i])
                 return true;
         return false;
     }
 
     void SemiInfGaussianQuadrature::getAbscissae(
         long double* abscissae, const unsigned len) const
     {
         if (len < npoints_) throw std::invalid_argument(
             "In npstat::SemiInfGaussianQuadrature::getAbscissae: "
             "unsifficient length of the output buffer");
         assert(abscissae);
         assert(a_);
         for (unsigned i=0; i<npoints_; ++i)
             abscissae[i] = a_[i];
     }
 
     void SemiInfGaussianQuadrature::getWeights(
         long double* weights, const unsigned len) const
     {
         if (len < npoints_) throw std::invalid_argument(
             "In npstat::SemiInfGaussianQuadrature::getWeights: "
             "unsifficient length of the output buffer");
         assert(weights);
         assert(w_);
         for (unsigned i=0; i<npoints_; ++i)
             weights[i] = w_[i];
     }
 }
Index: trunk/npstat/nm/SemiInfGaussianQuadrature.hh
===================================================================
--- trunk/npstat/nm/SemiInfGaussianQuadrature.hh	(revision 585)
+++ trunk/npstat/nm/SemiInfGaussianQuadrature.hh	(revision 586)
@@ -1,92 +1,92 @@
 #ifndef NPSTAT_SEMIINFGAUSSIANQUADRATURE_HH_
 #define NPSTAT_SEMIINFGAUSSIANQUADRATURE_HH_
 
 /*!
 // \file SemiInfGaussianQuadrature.hh
 //
 // \brief Gaussian quadrature with weight exp(-x^2/2) on [0, Infinity]
 //
 // Author: I. Volobouev
 //
 // June 2019
 */
 
 #include <vector>
 #include <utility>
 
 #include "npstat/nm/SimpleFunctors.hh"
 
 namespace npstat {
     /** 
     // Quadrature for functions of the type f(x) exp(-x^2/2) on [0, Infinity]
     */
     class SemiInfGaussianQuadrature
     {
     public:
         /**
         // At the moment, the following numbers of points are supported:
-        // 16.
+        // 8, 16.
         //
         // If an unsupported number of points is given in the
         // constructor, std::invalid_argument exception will be thrown.
         */
         explicit SemiInfGaussianQuadrature(unsigned npoints);
 
         /** Return the number of quadrature points */
         inline unsigned npoints() const {return npoints_;}
 
         /**
         // Get the abscissae for the integration points.
         // The buffer length should be at least npoints.
         */
         void getAbscissae(long double* abscissae, unsigned len) const;
 
         /**
         // Get the weights for the integration points.
         // The buffer length should be at least npoints.
         */
         void getWeights(long double* weights, unsigned len) const;
 
         /** Perform the quadrature */
         template <typename FcnResult, typename FcnArg>
         long double integrate(const Functor1<FcnResult,FcnArg>& fcn) const;
 
         /** Check if the rule with the given number of points is supported */
         static bool isAllowed(unsigned npoints);
 
         /** The complete set of allowed rules, in the increasing order */
         static std::vector<unsigned> allowedNPonts();
 
         /**
         // Minimum number of points, among the supported rules, which
         // integrates a polynomial with the given degree exactly (with
         // the appropriate weight). Returns 0 if the degree is out of range.
         */
         static unsigned minimalExactRule(unsigned polyDegree);
 
     private:
         SemiInfGaussianQuadrature();
 
         const long double* a_;
         const 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_, a_+npoints_);
         }
 
         inline std::vector<double> weights2() const
         {
             return std::vector<double>(w_, w_+npoints_);
         }
 #endif
     };
 }
 
 #include "npstat/nm/SemiInfGaussianQuadrature.icc"
 
 #endif // NPSTAT_SEMIINFGAUSSIANQUADRATURE_HH_