Index: trunk/npstat/nm/FejerQuadrature.icc
===================================================================
--- trunk/npstat/nm/FejerQuadrature.icc	(revision 575)
+++ trunk/npstat/nm/FejerQuadrature.icc	(revision 576)
@@ -1,51 +1,51 @@
 #include <cmath>
 #include <algorithm>
 
 namespace npstat {
     template <typename FcnResult, typename FcnArg>
     inline long double FejerQuadrature::integrate(
         const Functor1<FcnResult,FcnArg>& function,
         const long double left, const long double right) const
     {
         if (buf_.size() != npoints_)
             buf_.resize(npoints_);
         std::pair<long double, long double>* results = &buf_[0];
         const long double midpoint = (left + right)/2.0L;
         const long double unit = (right - left)/2.0L;
         FcnArg a;
         long double v;
         for (unsigned i=0; i<npoints_; ++i)
         {
             a = static_cast<FcnArg>(midpoint + unit*a_[i]);
             v = w_[i]*static_cast<long double>(function(a));
             results[i].first = std::abs(v);
             results[i].second = v;
         }
         std::sort(buf_.begin(), buf_.end());
         v = 0.0L;
         for (unsigned i=0; i<npoints_; ++i)
             v += results[i].second;
         return v*unit;
     }
 
     template <class Functor>
     inline std::vector<std::pair<double,double> >
-    FejerQuadrature::weightedChebyshevPoints(
+    FejerQuadrature::weightedIntegrationPoints(
         const Functor& weight,
         const long double left, const long double right) const
     {
         typedef std::pair<double,double> DDPair;
 
         std::vector<DDPair> result;
         result.reserve(npoints_);
         const long double midpoint = (left + right)/2.0L;
         const long double unit = (right - left)/2.0L;
         for (unsigned i=0; i<npoints_; ++i)
         {
             const long double x = midpoint + unit*a_[i];
             const double w = w_[i]*weight(x);
             result.push_back(DDPair(x, w));
         }
         return result;
     }
 }
Index: trunk/npstat/nm/GaussLegendreQuadrature.icc
===================================================================
--- trunk/npstat/nm/GaussLegendreQuadrature.icc	(revision 575)
+++ trunk/npstat/nm/GaussLegendreQuadrature.icc	(revision 576)
@@ -1,56 +1,81 @@
 #include <cmath>
 #include <stdexcept>
 #include <algorithm>
 
 namespace npstat {
     template <typename FcnResult, typename FcnArg>
     inline long double GaussLegendreQuadrature::integrate(
         const Functor1<FcnResult,FcnArg>& function,
         const long double left, const long double right) const
     {
         if (buf_.size() != npoints_)
             buf_.resize(npoints_);
         std::pair<long double, long double>* results = &buf_[0];
         const unsigned halfpoints = npoints_/2;
         const long double midpoint = (left + right)/2.0L;
         const long double unit = (right - left)/2.0L;
         FcnArg a;
         long double v;
         for (unsigned i=0; i<halfpoints; ++i)
         {
             a = static_cast<FcnArg>(midpoint + unit*a_[i]);
             v = w_[i]*static_cast<long double>(function(a));
             results[2*i].first = std::abs(v);
             results[2*i].second = v;
             a = static_cast<FcnArg>(midpoint - unit*a_[i]);
             v = w_[i]*static_cast<long double>(function(a));
             results[2*i+1].first = std::abs(v);
             results[2*i+1].second = v;
         }
         std::sort(buf_.begin(), buf_.end());
         v = 0.0L;
         for (unsigned i=0; i<npoints_; ++i)
             v += results[i].second;
         return v*unit;
     }
 
     template <typename FcnResult, typename FcnArg>
     inline long double GaussLegendreQuadrature::integrate(
         const Functor1<FcnResult,FcnArg>& function,
         const long double left, const long double right,
         const unsigned nsplit) const
     {
         if (!nsplit) throw std::invalid_argument(
             "In npstat::GaussLegendreQuadrature::integrate: "
             "number of subintervals must be positive");
         const long double step = (right - left)/nsplit;
         long double acc = 0.0L;
         for (unsigned i=0; i<nsplit; ++i)
         {
             const long double a = left + step*i;
             const long double b = (i == nsplit - 1U ? right : a + step);
             acc += integrate(function, a, b);
         }
         return acc;
     }
+
+    template <class Functor>
+    inline std::vector<std::pair<double,double> >
+    GaussLegendreQuadrature::weightedIntegrationPoints(
+        const Functor& weight,
+        const long double left, const long double right) const
+    {
+        typedef std::pair<double,double> DDPair;
+
+        std::vector<DDPair> result;
+        result.reserve(npoints_);
+        const unsigned halfpoints = npoints_/2;
+        const long double midpoint = (left + right)/2.0L;
+        const long double unit = (right - left)/2.0L;
+        for (unsigned i=0; i<halfpoints; ++i)
+        {
+            long double x = midpoint + unit*a_[i];
+            double w = w_[i]*weight(x);
+            result.push_back(DDPair(x, w));
+            x = midpoint - unit*a_[i];
+            w = w_[i]*weight(x);
+            result.push_back(DDPair(x, w));
+        }
+        return result;
+    }
 }
Index: trunk/npstat/nm/FejerQuadrature.hh
===================================================================
--- trunk/npstat/nm/FejerQuadrature.hh	(revision 575)
+++ trunk/npstat/nm/FejerQuadrature.hh	(revision 576)
@@ -1,89 +1,92 @@
 #ifndef NPSTAT_FEJERQUADRATURE_HH_
 #define NPSTAT_FEJERQUADRATURE_HH_
 
 /*!
 // \file FejerQuadrature.hh
 //
 // \brief Fejer quadratures in long double precision
 //
 // Author: I. Volobouev
 //
 // July 2018
 */
 
 #include <vector>
 #include <utility>
 
 #include "npstat/nm/SimpleFunctors.hh"
 
 namespace npstat {
     /** 
     // Fejer quadrature. Internally, operations are performed
     // in long double precision.
     */
     class FejerQuadrature
     {
     public:
         explicit FejerQuadrature(unsigned npoints);
 
         /** 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 getAbscissae(long double* abscissae, unsigned len) const;
 
         /**
         // The weights of the integration rule.
         // The buffer length should be at least npoints.
         */
         void getWeights(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;
 
         /**
-        // Weighted integration points on the given interval, suitable
-        // for constructing orthogonal polynomials w.r.t. the given
-        // weight function (in particular, by ContOrthoPoly1D class).
+        // Weighted integration points (Chebyshev points in this case)
+        // 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>
-        std::vector<std::pair<double,double> > weightedChebyshevPoints(
+        std::vector<std::pair<double,double> > weightedIntegrationPoints(
             const Functor& weight, 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());
         }
 #endif
     };
 }
 
 #include "npstat/nm/FejerQuadrature.icc"
 
 #endif // NPSTAT_FEJERQUADRATURE_HH_
Index: trunk/npstat/nm/GaussLegendreQuadrature.hh
===================================================================
--- trunk/npstat/nm/GaussLegendreQuadrature.hh	(revision 575)
+++ trunk/npstat/nm/GaussLegendreQuadrature.hh	(revision 576)
@@ -1,104 +1,115 @@
 #ifndef NPSTAT_GAUSSLEGENDREQUADRATURE_HH_
 #define NPSTAT_GAUSSLEGENDREQUADRATURE_HH_
 
 /*!
 // \file GaussLegendreQuadrature.hh
 //
 // \brief Gauss-Legendre quadratures in long double precision
 //
 // Author: I. Volobouev
 //
 // April 2010
 */
 
 #include <vector>
 #include <utility>
 
 #include "npstat/nm/SimpleFunctors.hh"
 
 namespace npstat {
     /** 
     // Gauss-Legendre quadrature. Internally, operations are performed
     // in long double precision.
     */
     class GaussLegendreQuadrature
     {
     public:
         /**
         // At the moment, the following numbers of points are supported:
         // 2, 4, 6, 8, 10, 12, 16, 32, 64, 100, 128, 256, 512, 1024.
         //
         // If an unsupported number of points is given in the
         // constructor, std::invalid_argument exception will be thrown.
         */
         explicit GaussLegendreQuadrature(unsigned npoints);
 
         /** Return the number of quadrature points */
         inline unsigned npoints() const {return npoints_;}
 
         /**
         // The abscissae are returned for positive points only,
         // so the buffer length should be at least npoints/2.
         */
         void getAbscissae(long double* abscissae, unsigned len) const;
 
         /**
         // The weights are returned for positive points only,
         // so the buffer length should be at least npoints/2.
         */
         void getWeights(long double* weights, unsigned len) const;
 
-        /** Perform the quadrature */
+        /** 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;
 
         /**
         // This method splits the interval [a, b] into "nsplit"
         // subintervals of equal length, applies Gauss-Legendre
         // quadrature to each subinterval, and sums the results.
         */
         template <typename FcnResult, typename FcnArg>
         long double integrate(const Functor1<FcnResult,FcnArg>& fcn,
                               long double a, long double b,
                               unsigned nsplit) const;
 
+        /**
+        // 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>
+        std::vector<std::pair<double,double> > weightedIntegrationPoints(
+            const Functor& weight, long double a, long double b) 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.
         // Returns 0 if the degree is out of range.
         */
         static unsigned minimalExactRule(unsigned polyDegree);
 
     private:
         GaussLegendreQuadrature();
 
         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_/2U);
         }
 
         inline std::vector<double> weights2() const
         {
             return std::vector<double>(w_, w_+npoints_/2U);
         }
 #endif
     };
 }
 
 #include "npstat/nm/GaussLegendreQuadrature.icc"
 
 #endif // NPSTAT_GAUSSLEGENDREQUADRATURE_HH_
Index: trunk/tests/test_ContOrthoPoly1D.cc
===================================================================
--- trunk/tests/test_ContOrthoPoly1D.cc	(revision 575)
+++ trunk/tests/test_ContOrthoPoly1D.cc	(revision 576)
@@ -1,206 +1,230 @@
 #include "UnitTest++.h"
 #include "test_utils.hh"
 
 #include "npstat/nm/ContOrthoPoly1D.hh"
 #include "npstat/nm/ClassicalOrthoPolys1D.hh"
 #include "npstat/nm/FejerQuadrature.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);
 
             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);
                 CHECK_CLOSE(w, fvalue, 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)
+    TEST(ContOrthoPoly1D_jacobi_1)
     {
-        const double eps = 1.0e-12;
+        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.weightedChebyshevPoints(weight, -1.0, 1.0), OPOLY_LANCZOS);
+        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);
         }
     }
 }