Index: trunk/npstat/nm/SemiInfGaussianQuadrature.hh
===================================================================
--- trunk/npstat/nm/SemiInfGaussianQuadrature.hh	(revision 586)
+++ trunk/npstat/nm/SemiInfGaussianQuadrature.hh	(revision 587)
@@ -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:
-        // 8, 16.
+        // 4, 8, 16, 24, 32.
         //
         // 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_
Index: trunk/npstat/nm/SemiInfGaussianQuadrature.cc
===================================================================
--- trunk/npstat/nm/SemiInfGaussianQuadrature.cc	(revision 586)
+++ trunk/npstat/nm/SemiInfGaussianQuadrature.cc	(revision 587)
@@ -1,115 +1,210 @@
 #include <cassert>
 #include <stdexcept>
 
 #include "npstat/nm/SemiInfGaussianQuadrature.hh"
 
-static const unsigned allowed[] = {8U, 16U};
+static const unsigned allowed[] = {4U, 8U, 16U, 24U, 32U};
+
+static const long double x4[4] = {
+    0.1891884656679243290128979L, 0.8829284441871048265177099L,
+    1.898635201026033700523118L, 3.199890790487880003820947L
+};
+
+static const long double w4[4] = {
+    0.4600479141368865531562111L, 0.5955353746508146156968605L,
+    0.1887161938025806930202006L, 0.009014654725218389334610484L
+};
 
 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
 };
 
+static const long double x24[24] = {
+    0.01546572428284263967158309L, 0.08102467456890859053226466L, 
+    0.1971495597774370175831004L, 0.3610060868778612268925103L, 
+    0.5689668040943317804934842L, 0.8170238621296853354221345L, 
+    1.101140512570927519762387L, 1.417517883185255653127699L, 
+    1.762770389346924215078403L, 2.134023433270852872718858L, 
+    2.528955881522317545103751L, 2.945810315977091985315583L, 
+    3.383390800591344448153514L, 3.841064291962996178510113L, 
+    4.318779879592193935827666L, 4.817121166723755257882443L, 
+    5.337413129910965390503086L, 5.881919719503349848587692L, 
+    6.45420236959330413008426L, 7.059790264811718937141413L, 
+    7.707527148649140750285707L, 8.412628265821523231899599L, 
+    9.205170595246298124280356L, 10.16471787591427543551836L
+};
+
+static const long double w24[24] = {
+    0.03963177373583971623520682L, 0.09092634417200772682452021L, 
+    0.1378388705514305352209355L, 0.1747888048604865886698802L, 
+    0.1945140755122695262325578L, 0.191050388085933371932191L, 
+    0.1640804718514537606961942L, 0.1213237868670473137964304L, 
+    0.07585085352846821444096815L, 0.03933583490926127777881721L, 
+    0.01659145184837393220345831L, 0.005577558962540112735262391L, 
+    0.001462903752340742782512164L, 0.0002925377088122708645115568L, 
+    0.00004346182159496116522023304L, 4.654928147248950254973396e-6L, 
+    3.465516883188886818939015e-7L, 1.712916429499539738584348e-8L, 
+    5.291582829662127654312429e-10L, 9.396968158076974083736813e-12L, 
+    8.476746787737468334671944e-14L, 3.178625220390072186141808e-16L, 
+    3.40152384541126386684788e-19L, 4.065997163351630907255484e-23L
+};
+
+static const long double x32[32] = {
+    0.01012878189691214624350868L, 0.05319424244367037906119417L, 
+    0.1299787175951869406561715L, 0.2393691656090843915139899L, 
+    0.3798543436482732972094256L, 0.5496455979334752304054161L, 
+    0.7467901653542001362974571L, 0.9692758581539488690148097L, 
+    1.215119792756784593423518L, 1.482437552089186111423254L, 
+    1.769492731870456339905001L, 2.074729217030429624151045L, 
+    2.396789740726666517125893L, 2.734524585497434861613664L, 
+    3.086994052653475580640136L, 3.453467859430628176713422L, 
+    3.833424142404338574398058L, 4.226550390562291219746039L, 
+    4.632748495051340594863436L, 5.052146267879433987095272L, 
+    5.485118365418631013862636L, 5.932320765862010050938686L, 
+    6.394745207225159107082344L, 6.873804136784972385013229L, 
+    7.371464544977494715968631L, 7.890464600740635307583666L, 
+    8.434680157505573739773468L, 9.009785727262104274673841L, 
+    9.624559242192175035329744L, 10.29381932756854327573723L, 
+    11.04655434209949069750599L, 11.95902534091863383425986L
+
+};
+
+static const long double w32[32] = {
+    0.02597243110108903860819742L, 0.05998522294742457352814454L, 
+    0.09252557558551830268517582L, 0.12168023193887471387632L, 
+    0.1446335439357877876644259L, 0.1580326574713564676245316L, 
+    0.159011590601382923917026L, 0.1465931521001040721470653L, 
+    0.1227765781762241411489409L, 0.09246888730922383493713899L, 
+    0.06193718064832836683995573L, 0.0364772456420070742443799L, 
+    0.01867171317426992524623413L, 0.00821064244096799641031966L, 
+    0.003065342546002567706332197L, 0.000959945434125516832948747L, 
+    0.0002490149086547299311600449L, 0.00005279972103866911857442108L, 
+    9.01991491001127595483372e-6L, 1.221817490526120766201247e-6L, 
+    1.288842643783733817900534e-7L, 1.036826850826123748376688e-8L, 
+    6.20549587716442676346212e-10L, 2.681447130674892888385382e-11L, 
+    8.058145645779520023829929e-13L, 1.605094827988058779886767e-14L, 
+    1.987693079561206625936811e-16L, 1.399402923082974754596922e-18L, 
+    4.90394420250063462879642e-21L, 6.889367999769572485562911e-24L, 
+    2.581217575596987383656126e-27L, 9.281477150840712552103502e-32L
+};
+
+
 namespace npstat {
     SemiInfGaussianQuadrature::SemiInfGaussianQuadrature(const unsigned npoints)
         : a_(0),
           w_(0),
           npoints_(npoints)
     {
         switch (npoints)
         {
+        case 4U:
+            a_ = x4;
+            w_ = w4;
+            break;
+
         case 8U:
             a_ = x8;
             w_ = w8;
             break;
 
         case 16U:
             a_ = x16;
             w_ = w16;
             break;
 
+        case 24U:
+            a_ = x24;
+            w_ = w24;
+            break;
+
+        case 32U:
+            a_ = x32;
+            w_ = w32;
+            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/tests/test_GaussHermiteQuadrature.cc
===================================================================
--- trunk/tests/test_GaussHermiteQuadrature.cc	(revision 586)
+++ trunk/tests/test_GaussHermiteQuadrature.cc	(revision 587)
@@ -1,102 +1,121 @@
 #include <cmath>
 #include <stdexcept>
 
 #include "UnitTest++.h"
 
 #include "test_utils.hh"
 
 #include "npstat/nm/GaussHermiteQuadrature.hh"
 #include "npstat/nm/SemiInfGaussianQuadrature.hh"
 
 using namespace npstat;
 using namespace std;
 
 namespace {
     struct Xsquared : public Functor1<long double, long double>
     {
         long double operator()(const long double& x) const
             {return x*x;}
     };
 
     struct X11 : public Functor1<long double, long double>
     {
         long double operator()(const long double& x) const
             {return powl(x, 11);}
     };
 
+    struct X41 : public Functor1<long double, long double>
+    {
+        long double operator()(const long double& x) const
+            {return powl(x, 41);}
+    };
+
     struct Expl : public Functor1<long double, long double>
     {
         inline explicit Expl(long double k) : k_(k) {}
         long double operator()(const long double& x) const
             {return expl(k_*x);}
     private:
         Expl();
         long double k_;
     };
 
     TEST(GaussHermiteQuadrature)
     {
         const unsigned nsup[] = {16, 32, 64, 100, 128, 256};
 
         for (unsigned icycle=0; icycle<sizeof(nsup)/sizeof(nsup[0]); ++icycle)
         {
             GaussHermiteQuadrature quad(nsup[icycle]);
             CHECK_EQUAL(nsup[icycle], quad.npoints());
 
             long double v = quad.integrate(Xsquared());
             CHECK_CLOSE(sqrt(M_PI)/2.0, v, 1.e-15);
 
             v = quad.integrateProb(3.0, 2.0, Xsquared());
             CHECK_CLOSE(13.0, v, 1.e-15);
 
             v = quad.integrateProb(0.0, 1.0, Xsquared());
             CHECK_CLOSE(1.0, v, 1.e-15);
 
             v = quad.integrate(Expl(0.5));
             CHECK_CLOSE(expl(0.5*0.5/4.0)*sqrtl(3.141592653589793238462643L),
                         v, powl(10.0, -0.15*nsup[icycle]));
 
             v = quad.integrateProb(3.0, 2.0, Expl(0.5));
             CHECK_CLOSE(expl(2.0), v, powl(10.0, -0.1*nsup[icycle]));
         }
 
         bool thrown = false;
         try
         {
             GaussHermiteQuadrature quad(1234567U);
         }
         catch (std::invalid_argument& e)
         {
             thrown = true;
         }
         CHECK(thrown);
     }
 
     TEST(SemiInfGaussianQuadrature)
     {
-        const unsigned nsup[] = {16};
+        const unsigned nsup[] = {4U, 8U, 16U, 24U, 32U};
 
         for (unsigned icycle=0; icycle<sizeof(nsup)/sizeof(nsup[0]); ++icycle)
         {
             SemiInfGaussianQuadrature quad(nsup[icycle]);
             CHECK_EQUAL(nsup[icycle], quad.npoints());
 
             long double v = quad.integrate(Xsquared());
             CHECK_CLOSE(sqrt(M_PI/2.0), v, 1.e-15);
 
-            v = quad.integrate(X11());
-            CHECK_CLOSE(3840.0, v, 3840.0*1.e-15);
+            v = quad.integrate(ConstValue1<long double,long double>(1.0L));
+            CHECK_CLOSE(sqrt(M_PI/2.0), v, 1.e-15);
+
+            if (nsup[icycle] > 5)
+            {
+                v = quad.integrate(X11());
+                CHECK_CLOSE(3840.0, v, 3840.0*1.e-15);
+            }
+
+            if (nsup[icycle] > 20)
+            {
+                const double expected = 2.55108265612582846464e24;
+                v = quad.integrate(X41());
+                CHECK_CLOSE(expected, v, expected*1.e-15);
+            }
         }
 
         bool thrown = false;
         try
         {
             SemiInfGaussianQuadrature quad(1234567U);
         }
         catch (std::invalid_argument& e)
         {
             thrown = true;
         }
         CHECK(thrown);
     }
 }
Index: trunk/tests/test_ContOrthoPoly1D.cc
===================================================================
--- trunk/tests/test_ContOrthoPoly1D.cc	(revision 586)
+++ trunk/tests/test_ContOrthoPoly1D.cc	(revision 587)
@@ -1,244 +1,244 @@
 #include "UnitTest++.h"
 #include "test_utils.hh"
 
 #include "npstat/nm/ContOrthoPoly1D.hh"
 #include "npstat/nm/ClassicalOrthoPolys1D.hh"
 #include "npstat/nm/FejerQuadrature.hh"
 #include "npstat/nm/StorablePolySeries1D.hh"
 
 #include "npstat/rng/MersenneTwister.hh"
 
 using namespace npstat;
 using namespace std;
 
 inline static int kdelta(const unsigned i, const unsigned j)
 {
     return i == j ? 1 : 0;
 }
 
 namespace {
     TEST(ContOrthoPoly1D_orthonormalization)
     {
         const double eps = 1.0e-15;
 
         const OrthoPolyMethod method[] = {OPOLY_STIELTJES, OPOLY_LANCZOS};
         const unsigned nMethods = sizeof(method)/sizeof(method[0]);
             
         const unsigned npoints = 64;
         const unsigned maxdeg = 10;
         const unsigned ntries = 10;
 
         std::vector<double> points(2U*npoints);
         std::vector<ContOrthoPoly1D::MeasurePoint> measure;
         measure.reserve(npoints);
         std::vector<double> coeffs(maxdeg + 1);
 
         MersenneTwister rng;
 
         for (unsigned imeth=0; imeth<nMethods; ++imeth)
             for (unsigned itry=0; itry<ntries; ++itry)
             {
                 const double xmin = rng()*5.0 - 6.0;
                 const double xmax = rng()*2.0 + 1.0;
                 const double range = xmax - xmin;
                 rng.run(&points[0], 2U*npoints, 2U*npoints);
 
                 measure.clear();
                 for (unsigned i=0; i<npoints; ++i)
                 {
                     ContOrthoPoly1D::MeasurePoint p(points[2*i]*range + xmin, points[2*i+1]);
                     measure.push_back(p);
                 }
                 ContOrthoPoly1D poly(maxdeg, measure, method[imeth]);
 
                 for (unsigned i=0; i<=maxdeg; ++i)
                     for (unsigned j=0; j<=maxdeg; ++j)
                     {
                         const double d = poly.empiricalKroneckerDelta(i, j);
                         CHECK_CLOSE(static_cast<double>(kdelta(i, j)), d, eps);
                     }
 
                 measure.clear();
                 for (unsigned i=0; i<npoints; ++i)
                 {
                     ContOrthoPoly1D::MeasurePoint p(points[2*i]*range + xmin, 1.0);
                     measure.push_back(p);
                 }
                 ContOrthoPoly1D poly2(maxdeg, measure, method[imeth]);
 
                 poly2.weightCoeffs(&coeffs[0], maxdeg);
                 for (unsigned i=0; i<=maxdeg; ++i)
                     CHECK_CLOSE(static_cast<double>(kdelta(i, 0)), coeffs[i], eps);
             }
     }
 
     TEST(ContOrthoPoly1D_weightCoeffs)
     {
         const double eps = 1.0e-7;
 
         const unsigned maxdeg = 10;
         const unsigned ntries = 10;
         const unsigned npoints = maxdeg + 1U;
 
         std::vector<double> points(2U*npoints);
         std::vector<ContOrthoPoly1D::MeasurePoint> measure;
         measure.reserve(npoints);
         std::vector<double> coeffs(maxdeg + 1);
 
         MersenneTwister rng;
 
         for (unsigned itry=0; itry<ntries; ++itry)
         {
             const double xmin = rng()*5.0 - 6.0;
             const double xmax = rng()*2.0 + 1.0;
             const double range = xmax - xmin;
             rng.run(&points[0], 2U*npoints, 2U*npoints);
 
             measure.clear();
             for (unsigned i=0; i<npoints; ++i)
             {
                 ContOrthoPoly1D::MeasurePoint p(points[2*i]*range + xmin, points[2*i+1]);
                 measure.push_back(p);
             }
             ContOrthoPoly1D poly(maxdeg, measure, OPOLY_LANCZOS);
             poly.weightCoeffs(&coeffs[0], maxdeg);
             CPP11_auto_ptr<StorablePolySeries1D> poly2(
                 poly.makeStorablePolySeries(xmin, xmax));
 
             for (unsigned deg=0; deg<=maxdeg; ++deg)
                 for (unsigned ipow=0; ipow<4; ++ipow)
                 {
                     const double i1 = poly.integratePoly(deg, ipow, xmin, xmax);
                     const double i2 = poly2->integratePoly(deg, ipow);
                     if (fabs(i2) > 1.0e-10)
                         CHECK_CLOSE(1.0, i1/i2, 1.0e-10);
                 }
 
             for (unsigned i=0; i<npoints; ++i)
             {
                 const double x = measure[i].first;
                 const double w = measure[i].second;
                 const double fvalue = poly.series(&coeffs[0], maxdeg, x);
                 const double f2 = poly2->series(&coeffs[0], maxdeg, x);
-                CHECK_CLOSE(w, fvalue, eps);
+                CHECK_CLOSE(w, fvalue, 10*eps);
                 CHECK_CLOSE(fvalue, f2, eps);
             }
         }
     }
 
     TEST(ContOrthoPoly1D_cov8)
     {
         const double eps = 1.0e-12;
 
         MersenneTwister rng;
         const unsigned npoints = 64;
         const unsigned maxdeg = 6;
         const unsigned ntries = 5;
         const unsigned degtries = 100;
         std::vector<double> points(npoints);
 
         for (unsigned itry=0; itry<ntries; ++itry)
         {
             rng.run(&points[0], npoints, npoints);
             ContOrthoPoly1D poly(maxdeg, points);
 
             for (unsigned ideg=0; ideg<degtries; ++ideg)
             {
                 unsigned d[8];
                 bool has0 = true;
                 while (has0)
                 {
                     for (unsigned i=0; i<8; ++i)
                         d[i] = (maxdeg + 0.99)*rng();
                     has0 = (d[0] == 0 && d[1] == 0) ||
                            (d[2] == 0 && d[3] == 0) ||
                            (d[4] == 0 && d[5] == 0) ||
                            (d[6] == 0 && d[7] == 0);
                 }
                 const double cov8 = poly.cov8(d[0],d[1],d[2],d[3],d[4],d[5],d[6],d[7]);
                 const double check = poly.slowCov8(d[0],d[1],d[2],d[3],d[4],d[5],d[6],d[7]);
                 CHECK_CLOSE(check, cov8, eps);
             }
         }
     }
 
     TEST(ContOrthoPoly1D_epsCovariance)
     {
         const double eps = 1.0e-15;
         const unsigned npoints = 64;
         const unsigned maxdeg = 6;
         const unsigned ntries = 5;
 
         MersenneTwister rng;
         std::vector<double> points(npoints);
 
         for (unsigned itry=0; itry<ntries; ++itry)
         {
             rng.run(&points[0], npoints, npoints);
             ContOrthoPoly1D poly(maxdeg, points);
 
             for (unsigned i=0; i<=maxdeg; ++i)
                 for (unsigned j=0; j<=i; ++j)
                 {
                     CHECK(poly.epsExpectation(i, j, false) == 0.0);
 
                     for (unsigned k=0; k<=maxdeg; ++k)
                         for (unsigned m=0; m<=k; ++m)
                         {
                             const double eps1 = poly.empiricalKroneckerCovariance(i, j, k, m);
                             const double eps2 = poly.epsCovariance(i, j, k, m, false);
                             CHECK_CLOSE(eps1, eps2, eps);
                         }
                 }
         }
     }
 
     // Test that we can reproduce Jacobi polynomials using ContOrthoPoly1D
     TEST(ContOrthoPoly1D_jacobi_1)
     {
         const double eps = 1.0e-13;
 
         const unsigned alpha = 3;
         const unsigned beta = 4;
         const unsigned maxdeg = 5;
         const unsigned npoints = 64;
 
         JacobiOrthoPoly1D jac(alpha, beta);
         OrthoPoly1DWeight weight(jac);
 
         const unsigned npt = FejerQuadrature::minimalExactRule(2*maxdeg+alpha+beta);
         FejerQuadrature fej(npt);
         ContOrthoPoly1D poly(maxdeg, fej.weightedIntegrationPoints(weight, -1.0, 1.0), OPOLY_LANCZOS);
 
         for (unsigned i=0; i<npoints; ++i)
         {
             const double x = test_rng()*2.0 - 1.0;
             for (unsigned deg=0; deg<=maxdeg; ++deg)
                 CHECK_CLOSE(jac.poly(deg, x), poly.poly(deg, x), eps);
         }
     }
 
     TEST(ContOrthoPoly1D_jacobi_2)
     {
         const double eps = 1.0e-13;
 
         const unsigned alpha = 3;
         const unsigned beta = 4;
         const unsigned maxdeg = 5;
         const unsigned npoints = 64;
 
         JacobiOrthoPoly1D jac(alpha, beta);
         OrthoPoly1DWeight weight(jac);
 
         const unsigned npt = GaussLegendreQuadrature::minimalExactRule(2*maxdeg+alpha+beta);
         GaussLegendreQuadrature glq(npt);
         ContOrthoPoly1D poly(maxdeg, glq.weightedIntegrationPoints(weight, -1.0, 1.0), OPOLY_LANCZOS);
 
         for (unsigned i=0; i<npoints; ++i)
         {
             const double x = test_rng()*2.0 - 1.0;
             for (unsigned deg=0; deg<=maxdeg; ++deg)
                 CHECK_CLOSE(jac.poly(deg, x), poly.poly(deg, x), eps);
         }
     }
 }