Index: trunk/tests/jacobiEpsTest.cc
===================================================================
--- trunk/tests/jacobiEpsTest.cc	(revision 554)
+++ trunk/tests/jacobiEpsTest.cc	(revision 555)
@@ -1,148 +1,189 @@
 #include <string>
 #include <sstream>
 #include <iostream>
 
 #include "test_utils.hh"
 
 #include "npstat/nm/ClassicalOrthoPolys1D.hh"
+#include "npstat/nm/ContOrthoPoly1D.hh"
 #include "npstat/nm/GaussLegendreQuadrature.hh"
 
 #include "npstat/stat/Distributions1D.hh"
 #include "npstat/stat/SampleAccumulator.hh"
 #include "npstat/stat/orthoPoly1DVProducts.hh"
 
 #include "npstat/rng/MersenneTwister.hh"
 
 using namespace npstat;
 using namespace std;
 
 
+inline static double kdelta(const unsigned i,
+                            const unsigned j)
+{
+    return i == j ? 1.0 : 0.0;
+}
+
+
 static double calculateExpectation(
     const JacobiOrthoPoly1D& jac, const unsigned nPoints,
     const unsigned m, const unsigned n, const bool highOrder)
 {
     const unsigned sumdeg = (unsigned)(jac.alpha()) + (unsigned)(jac.beta()) +
                             4U*std::max(m, n);
     const unsigned nGood = GaussLegendreQuadrature::minimalExactRule(sumdeg);
     assert(nGood);
     GaussLegendreQuadrature glq(nGood);
     return oracleEpsExpectation(jac, glq, nPoints, m, n, highOrder);
 }
 
 
 static double calculateCovariance(
     const JacobiOrthoPoly1D& jac, const unsigned nPoints,
     const unsigned m1, const unsigned n1,
     const unsigned m2, const unsigned n2,
     const bool highOrder)
 {
     const unsigned sumdeg = (unsigned)(jac.alpha()) + (unsigned)(jac.beta()) +
         8U*std::max(std::max(m1, n1), std::max(m2, n2));
     const unsigned nGood = GaussLegendreQuadrature::minimalExactRule(sumdeg);
     assert(nGood);
     GaussLegendreQuadrature glq(nGood);
     return oracleEpsCovariance(jac, glq, nPoints, m1, n1, m2, n2, highOrder);
 }
 
 
-static string eps_name(const unsigned m, const unsigned n)
+static string eps_name(const unsigned m, const unsigned n, const bool isReal)
 {
     ostringstream os;
+    if (isReal)
+        os << "r_";
     os << "eps_" << m << '_' << n;
     return os.str();
 }
 
 
 template <class Acc>
 static void print_accumulator(const string& which,
                               const Acc& acc, const double expMean)
 {
     const double m = acc.mean();
     const double u = acc.meanUncertainty();
     const double nsig = (m - expMean)/u;
     
     cout << which << " average is " << m << " +- " << u
          << ", expected " << expMean << " (" << nsig << " sigma)" << endl;
 }
 
 
 int main(int argc, char const* argv[])
 {
     if (argc != 10)
     {
         cout << "Usage: " << parse_progname(argv[0])
              << " alpha beta pointsPerPseudo nPseudo m n p q highOrder" << endl;
         exit(1);
     }
 
     unsigned alpha, beta;
     if (parse_unsigned(argv[1], &alpha)) exit(1);
     if (parse_unsigned(argv[2], &beta)) exit(1);
 
     unsigned perpseudo;
     if (parse_unsigned(argv[3], &perpseudo)) exit(1);
     if (!perpseudo)
     {
         cerr << "pointsPerPseudo must be positive" << endl;
         exit(1);
     }
 
     unsigned npseudo;
     if (parse_unsigned(argv[4], &npseudo)) exit(1);
     if (npseudo < 2)
     {
         cerr << "nPseudo must be at least 2" << endl;
         exit(1);
     }
 
     unsigned m, n, p, q;
     if (parse_unsigned(argv[5], &m)) exit(1);
     if (parse_unsigned(argv[6], &n)) exit(1);
     if (parse_unsigned(argv[7], &p)) exit(1);
     if (parse_unsigned(argv[8], &q)) exit(1);
+    const unsigned maxdeg = std::max(std::max(m, n), std::max(p, q));
 
     bool highOrder;
     if (parse_bool(argv[9], &highOrder)) exit(1);
 
     vector<double> rn;
     rn.reserve(perpseudo);
 
     Beta1D distro(-1.0, 2.0, beta+1, alpha+1);
     MersenneTwister rng;
     JacobiOrthoPoly1D jac(alpha, beta);
+    OrthoPoly1DWeight jacWeight(jac);
 
-    SampleAccumulator<double> acc0, acc1;
+    const unsigned sumdeg = (unsigned)(jac.alpha()) + (unsigned)(jac.beta()) + 2U*maxdeg;
+    const unsigned integPoints = GaussLegendreQuadrature::minimalExactRule(sumdeg);
+
+    SampleAccumulator<double> acc0, acc1, acc2, acc3;
     acc0.reserve(npseudo);
     acc1.reserve(npseudo);
+    acc2.reserve(npseudo);
+    acc3.reserve(npseudo);
 
     for (unsigned ips=0; ips<npseudo; ++ips)
     {
         rn.clear();
         double tmp;
         for (unsigned i=0; i<perpseudo; ++i)
         {
             distro.random(rng, &tmp);
             rn.push_back(tmp);
         }
+        ContOrthoPoly1D samplePoly(maxdeg, rn);
 
         const double v0 = sampleEpsValue(jac, &rn[0], perpseudo, m, n, highOrder);
         acc0.accumulate(v0);
 
         const double v1 = sampleEpsValue(jac, &rn[0], perpseudo, p, q, highOrder);
         acc1.accumulate(v1);
+
+        const double winteg = samplePoly.integrateEWPolyProduct(
+            jacWeight, 0, 1, -1.0, 1.0, integPoints);
+        cout << "winteg is " << winteg << endl;
+
+        const double eps0 = samplePoly.integrateEWPolyProduct(
+            jacWeight, m, n, -1.0, 1.0, integPoints) - kdelta(m, n);
+        acc2.accumulate(eps0);
+
+        const double eps1 = samplePoly.integrateEWPolyProduct(
+            jacWeight, p, q, -1.0, 1.0, integPoints) - kdelta(m, n);
+        acc3.accumulate(eps1);
     }
 
-    const std::string& name0 = eps_name(m, n);
+    const std::string& name0 = eps_name(m, n, false);
+    const std::string& rname0 = eps_name(m, n, true);
     const double em0 = calculateExpectation(jac, perpseudo, m, n, highOrder);
     print_accumulator(name0, acc0, em0);
+    cout << endl;
+    print_accumulator(rname0, acc2, em0);
+    cout << endl;
 
-    const std::string& name1 = eps_name(p, q);
+    const std::string& name1 = eps_name(p, q, false);
+    const std::string& rname1 = eps_name(p, q, true);
     const double em1 = calculateExpectation(jac, perpseudo, p, q, highOrder);
     print_accumulator(name1, acc1, em1);
+    cout << endl;
+    print_accumulator(rname1, acc3, em1);
+    cout << endl;
 
     const double cov = calculateCovariance(jac, perpseudo, m, n, p, q, highOrder);
     cout << name0 << '/' << name1 << " covariance is " << acc0.cov(acc1)
          << ", expected " << cov << endl;
+    cout << endl;
+    cout << rname0 << '/' << rname1 << " covariance is " << acc2.cov(acc3)
+         << ", expected " << cov << endl;
 
     return 0;
 }
Index: trunk/npstat/nm/ContOrthoPoly1D.hh
===================================================================
--- trunk/npstat/nm/ContOrthoPoly1D.hh	(revision 554)
+++ trunk/npstat/nm/ContOrthoPoly1D.hh	(revision 555)
@@ -1,474 +1,474 @@
 #ifndef NPSTAT_CONTORTHOPOLY1D_HH_
 #define NPSTAT_CONTORTHOPOLY1D_HH_
 
 /*!
 // \file ContOrthoPoly1D.hh
 //
 // \brief Continuous orthogonal polynomial series in one dimension
 //        generated for a discrete measure
 //
 // Author: I. Volobouev
 //
 // May 2017
 */
 
 #include <vector>
 #include <utility>
 #include <cassert>
 #include <cmath>
 #include <map>
 
 #include "npstat/nm/OrthoPolyMethod.hh"
 #include "npstat/nm/Matrix.hh"
 #include "npstat/nm/SimpleFunctors.hh"
 
 namespace npstat {
     class ContOrthoPoly1D
     {
     public:
         // The first element of the pair is the coordinate and the
         // second measure is the weight (all weights must be non-negative)
         typedef std::pair<double,double> MeasurePoint;
 
         /**
         // Main constructor, with obvious arguments. Internally, the
         // measure will be sorted in the order of increasing weight
         // and the measure coordinates will be shifted so that their
         // weighted mean is at 0.
         */
         ContOrthoPoly1D(unsigned maxDegree,
                         const std::vector<MeasurePoint>& measure,
                         OrthoPolyMethod m = OPOLY_STIELTJES);
 
         /** Constructor which assumes that all weights are 1.0 */
         ContOrthoPoly1D(unsigned maxDegree,
                         const std::vector<double>& coords,
                         OrthoPolyMethod m = OPOLY_STIELTJES);
 
         //@{
         /** A simple inspector of object properties */
         inline unsigned maxDegree() const {return maxdeg_;}
         inline unsigned long measureLength() const {return measure_.size();}
         inline MeasurePoint measure(const unsigned index) const
             {return measure_.at(index);}
         inline double coordinate(const unsigned index) const
             {return measure_.at(index).first;}
         inline double weight(const unsigned index) const
             {return measure_.at(index).second;}
         inline double sumOfWeights() const {return wsum_;}
         inline double sumOfWeightSquares() const {return wsumsq_;}
         inline bool areAllWeightsEqual() const {return allWeightsEqual_;}
         inline double minCoordinate() const {return minX_;}
         inline double maxCoordinate() const {return maxX_;}
         inline double meanCoordinate() const {return meanX_;}
         //@}
 
         /** Kish's effective sample size for the measure */
         double effectiveSampleSize() const;
 
         /** Return the value of one of the normalized polynomials */
         double poly(unsigned deg, double x) const;
 
         /** Return the values of two of the normalized polynomials */
         std::pair<double,double> polyPair(
             unsigned deg1, unsigned deg2, double x) const;
 
         //@{
         /**
         // Return the values of all orthonormal polynomials up to some
         // maximum degree. Size of the "polyValues" array should be
         // at least maxdeg + 1.
         */
         void allPolys(double x, unsigned maxdeg, double *polyValues) const;
         void allPolys(double x, unsigned maxdeg, long double *polyValues) const;
         //@}
 
         /** Return the value of the orthonormal polynomial series */
         double series(const double *coeffs, unsigned maxdeg, double x) const;
 
         /**
         // Build the coefficients of the orthogonal polynomial series
         // for the given function. The length of the array "coeffs"
         // should be at least "maxdeg" + 1. Note that the coefficients
         // are returned in the order of increasing degree (same order
         // is used by the "series" function).
         */
         template <class Functor>
         void calculateCoeffs(const Functor& fcn,
                              double *coeffs, unsigned maxdeg) const;
 
         /**
         // Integrated squared error between the given function and the
         // polynomial series. Parameter "integrationRulePoints" specifies
         // the parameter "npoints" for the GaussLegendreQuadrature class.
         // If "integrationRulePoints" is 0, the rule will be picked
         // automatically in such a way that it integrates a polynomial
         // of degree maxdeg*2 exactly.
         */
         template <class Functor>
         double integratedSquaredError(const Functor& fcn,
                                       const double *coeffs, unsigned maxdeg,
                                       double xmin, double xmax,
                                       unsigned integrationRulePoints=0U) const;
 
         /**
         // Squared error between the given function and the polynomial
         // series, weighted according to the measure
         */
         template <class Functor>
         double weightedSquaredError(const Functor& fcn, const double *coeffs,
                                     unsigned maxdeg) const;
 
         /**
         // This method is useful for testing the numerical precision
         // of the orthonormalization procedure. It returns the scalar
         // products between various polynomials.
         */
         double empiricalKroneckerDelta(unsigned deg1, unsigned deg2) const;
 
         /**
         // If the Kronecker deltas were calculated from a sample, they
         // would be random and would have a covariance matrix. This
         // is an estimate of the terms in that covariance matrix. The
         // covariance is between delta_{deg1,deg2} and delta_{deg3,deg4}.
         */
         double empiricalKroneckerCovariance(unsigned deg1, unsigned deg2,
                                             unsigned deg3, unsigned deg4) const;
 
         /** 
         // Examine the recurrence coefficients. The first element of
         // the returned pair is alpha, and the second element is beta.
         */
         std::pair<double,double> recurrenceCoeffs(unsigned deg) const;
 
         /**
         // Generate principal minor of order n of the Jacobi matrix.
         // n must not exceed "maxDegree"
         */
         Matrix<double> jacobiMatrix(unsigned n) const;
 
         /** 
         // Roots of the polynomial with the given degree. Naturally,
         // the degree argument must not exceed "maxDegree".
         */
         void calculateRoots(double *roots, unsigned degree) const;
 
         /**
         // Integral of an orthonormal polynomials to some power without
         // the weight (so this is not an inner product with the poly of
         // degree 0).
         */
         double integratePoly(unsigned degree, unsigned power,
                              double xmin, double xmax) const;
 
         /**
         // Integral of the product of three orthonormal polynomials
         // without the weight (so this is not an inner product)
         */
         double integrateTripleProduct(unsigned deg1, unsigned deg2,
                                       unsigned deg3, double xmin,
                                       double xmax) const;
 
         /**
         // Integral of the product of two orthonormal polynomials
         // with some external weight function (e.g., oracle density).
         // "EW" in the method name stands for "externally weighted".
         //
         // "integrationRulePoints" argument will be passed to the
         // GaussLegendreQuadrature constructor and must be supported
         // by that class.
         */
         template <class Functor>
         double integrateEWPolyProduct(const Functor& weight,
                                       unsigned deg1, unsigned deg2,
                                       double xmin, double xmax,
                                       unsigned integrationRulePoints) const;
 
         /**
         // A measure-weighted average of orthonormal poly values
         // to some power
         */
         double powerAverage(unsigned deg, unsigned power) const;
 
         /**
         // A measure-weighted average of the product of two orthonormal
         // poly values to some powers
         */
         double jointPowerAverage(unsigned deg1, unsigned power1,
                                  unsigned deg2, unsigned power2) const;
 
         /** 
         // A measure-weighted average of the product of several
         // orthonormal poly values
         */
         double jointAverage(const unsigned* degrees, unsigned nDegrees,
                             bool degreesSortedIncreasingOrder = false) const;
 
         //@{
         /**
         // A measure-weighted average that will be remembered internally
         // so that another call to this function with the same arguments
         // will return the answer quickly
         */
         double cachedJointAverage(unsigned deg1, unsigned deg2,
                                   unsigned deg3, unsigned deg4) const;
         double cachedJointAverage(unsigned deg1, unsigned deg2,
                                   unsigned deg3, unsigned deg4,
                                   unsigned deg5, unsigned deg6) const;
         double cachedJointAverage(unsigned deg1, unsigned deg2,
                                   unsigned deg3, unsigned deg4,
                                   unsigned deg5, unsigned deg6,
                                   unsigned deg7, unsigned deg8) const;
         //@}
 
         /** Covariance between v_{ab} and v_{cd} */
         double cov4(unsigned a, unsigned b, unsigned c, unsigned d) const;
 
         /** Covariance between v_{ab}*v_{cd} and v_{ef} */
         double cov6(unsigned a, unsigned b, unsigned c, unsigned d,
                     unsigned e, unsigned f) const;
 
         /** Covariance between v_{ab}*v_{cd} and v_{ef}*v_{gh} */
         double cov8(unsigned a, unsigned b, unsigned c, unsigned d,
                     unsigned e, unsigned f, unsigned g, unsigned h) const;
 
         /**
         // Alternative implementation of "cov8". It is slower than "cov8"
         // but easier to program and verify. Useful mainly for testing "cov8".
         */
         double slowCov8(unsigned a, unsigned b, unsigned c, unsigned d,
                         unsigned e, unsigned f, unsigned g, unsigned h) const;
 
         double epsExpectation(unsigned m, unsigned n, bool highOrder) const;
 
         double epsCovariance(unsigned m1, unsigned n1,
                              unsigned m2, unsigned n2, bool highOrder) const;
 
     private:
         struct Recurrence
         {
             inline Recurrence(const long double a,
                               const long double b)
                 : alpha(a), beta(b)
             {
                 assert(beta > 0.0L);
                 sqrbeta = sqrtl(beta);
             }
 
             long double alpha;
             long double beta;
             long double sqrbeta;
         };
 
         class PolyFcn;
         friend class PolyFcn;
 
         class PolyFcn : public Functor1<long double, long double>
         {
         public:
             inline PolyFcn(const ContOrthoPoly1D& p,
                            const unsigned d1, const unsigned power)
                 : poly(p), deg1(d1), polypow(power) {}
 
             inline long double operator()(const long double& x) const
             {
                 long double p = 1.0L;
                 if (polypow)
                 {
                     p = poly.normpoly(deg1, x);
                     switch (polypow)
                     {
                     case 1U:
                         break;
                     case 2U:
                         p *= p;
                         break;
                     default:
                         p = powl(p, polypow);
                         break;
                     }
                 }
                 return p;
             }
 
         private:
             const ContOrthoPoly1D& poly;
             unsigned deg1;
             unsigned polypow;
         };
 
         class TripleProdFcn;
         friend class TripleProdFcn;
 
         class TripleProdFcn : public Functor1<long double, long double>
         {
         public:
             inline TripleProdFcn(const ContOrthoPoly1D& p, const unsigned d1,
                                  const unsigned d2, const unsigned d3)
                 : poly(p)
             {
                 degs[0] = d1;
                 degs[1] = d2;
                 degs[2] = d3;
                 degmax = std::max(d1, std::max(d2, d3));
             }
 
             inline long double operator()(const long double& x) const
                 {return poly.normpolyprod(degs, 3, degmax, x);}
 
         private:
             const ContOrthoPoly1D& poly;
             unsigned degs[3];
             unsigned degmax;
         };
 
         template<class Funct> class EWPolyProductFcn;
         template<class Funct> friend class EWPolyProductFcn;
 
         template<class Funct>
         class EWPolyProductFcn : public Functor1<long double, long double>
         {
         public:
             inline EWPolyProductFcn(const ContOrthoPoly1D& p, const Funct& w,
                                     const unsigned d1, const unsigned d2)
-                : poly(p), wfun(w)
+                : poly(p), wf(w)
             {
                 degs[0] = d1;
                 degs[1] = d2;
                 degmax = std::max(d1, d2);
             }
 
             inline long double operator()(const long double& x) const
-                {return poly.normpolyprod(degs, 2, degmax, x)*wfun(x);}
+                {return poly.normpolyprod(degs,2,degmax,x-poly.meanX_)*wf(x);}
 
         private:
             const ContOrthoPoly1D& poly;
-            const Funct& wfun;
+            const Funct& wf;
             unsigned degs[2];
             unsigned degmax;
         };
 
         class MemoKey
         {
         public:
             inline MemoKey(const unsigned d0, const unsigned d1,
                            const unsigned d2, const unsigned d3)
                 : nDegs_(4U)
             {
                 degs_[0] = d0;
                 degs_[1] = d1;
                 degs_[2] = d2;
                 degs_[3] = d3;                
                 std::sort(degs_, degs_+nDegs_);
             }
 
             inline MemoKey(const unsigned d0, const unsigned d1,
                            const unsigned d2, const unsigned d3,
                            const unsigned d4, const unsigned d5)
                 : nDegs_(6U)
             {
                 degs_[0] = d0;
                 degs_[1] = d1;
                 degs_[2] = d2;
                 degs_[3] = d3;                
                 degs_[4] = d4;                
                 degs_[5] = d5;                
                 std::sort(degs_, degs_+nDegs_);
             }
 
             inline MemoKey(const unsigned d0, const unsigned d1,
                            const unsigned d2, const unsigned d3,
                            const unsigned d4, const unsigned d5,
                            const unsigned d6, const unsigned d7)
                 : nDegs_(8U)
             {
                 degs_[0] = d0;
                 degs_[1] = d1;
                 degs_[2] = d2;
                 degs_[3] = d3;                
                 degs_[4] = d4;                
                 degs_[5] = d5;                
                 degs_[6] = d6;                
                 degs_[7] = d7;                
                 std::sort(degs_, degs_+nDegs_);
             }
 
             inline const unsigned* degrees() const {return degs_;}
             inline unsigned nDegrees() const {return nDegs_;}
 
             inline bool operator<(const MemoKey& r) const
             {
                 if (nDegs_ < r.nDegs_)
                     return true;
                 if (r.nDegs_ < nDegs_)
                     return false;
                 for (unsigned i=0; i<nDegs_; ++i)
                 {
                     if (degs_[i] < r.degs_[i])
                         return true;
                     if (r.degs_[i] < degs_[i])
                         return false;
                 }
                 return false;
             }
 
             inline bool operator==(const MemoKey& r) const
             {
                 if (nDegs_ != r.nDegs_)
                     return false;
                 for (unsigned i=0; i<nDegs_; ++i)
                     if (degs_[i] != r.degs_[i])
                         return false;
                 return true;
             }
 
             inline bool operator!=(const MemoKey& r) const
                 {return !(*this == r);}
 
         private:
             unsigned degs_[8];
             unsigned nDegs_;
         };
 
         void calcRecurrenceCoeffs(OrthoPolyMethod m);
 
         void calcRecurrenceStieltjes();
         void calcRecurrenceLanczos();
 
         long double monicpoly(unsigned degree, long double x) const;
 
         std::pair<long double,long double> monicInnerProducts(
             unsigned degree) const;
 
         long double normpoly(unsigned degree, long double x) const;
 
         std::pair<long double,long double> twonormpoly(
             unsigned deg1, unsigned deg2, long double x) const;
 
         long double normpolyprod(const unsigned* degrees, unsigned nDeg,
                                  unsigned degmax, long double x) const;
 
         long double normseries(const double *coeffs, unsigned maxdeg,
                                long double x) const;
 
         template <typename Numeric>
         void getAllPolys(double x, unsigned maxdeg, Numeric *polyValues) const;
 
         double getCachedAverage(const MemoKey& key) const;
 
         std::vector<MeasurePoint> measure_;
         std::vector<Recurrence> rCoeffs_;
         mutable std::map<MemoKey,double> cachedAverages_;
         long double wsum_;
         long double wsumsq_;
         double meanX_;
         double minX_;
         double maxX_;
         unsigned maxdeg_;
         bool allWeightsEqual_;
     };
 }
 
 #include "npstat/nm/ContOrthoPoly1D.icc"
 
 #endif // NPSTAT_CONTORTHOPOLY1D_HH_
Index: trunk/npstat/nm/ContOrthoPoly1D.icc
===================================================================
--- trunk/npstat/nm/ContOrthoPoly1D.icc	(revision 554)
+++ trunk/npstat/nm/ContOrthoPoly1D.icc	(revision 555)
@@ -1,121 +1,123 @@
 #include <stdexcept>
 
 #include "npstat/nm/GaussLegendreQuadrature.hh"
 
 namespace npstat {
     template <class Functor>
     void ContOrthoPoly1D::calculateCoeffs(const Functor& fcn,
                                           double *coeffs,
                                           const unsigned maxdeg) const
     {
         if (maxdeg > maxdeg_) throw std::invalid_argument(
             "In npstat::ContOrthoPoly1D::calculateCoeffs: "
             "maximum degree out of range");
         assert(coeffs);
 
         std::vector<long double> scalarProducts(maxdeg+1U, 0.0L);
 
         const unsigned long mSize = measure_.size();
         for (unsigned long i = 0; i < mSize; ++i)
         {
             const double x = measure_[i].first;
             const double w = measure_[i].second;
             const long double f = fcn(x + meanX_);
 
             long double polyk = 1.0L, polykm1 = 0.0L;
             for (unsigned k=0; k<maxdeg; ++k)
             {
                 scalarProducts[k] += polyk*f*w;
                 const long double poly = ((x - rCoeffs_[k].alpha)*polyk -
                                           rCoeffs_[k].sqrbeta*polykm1)/rCoeffs_[k+1].sqrbeta;
                 polykm1 = polyk;
                 polyk = poly;
             }
             scalarProducts[maxdeg] += polyk*f*w;
         }
 
         for (unsigned deg=0; deg<=maxdeg; ++deg)
             coeffs[deg] = scalarProducts[deg]/wsum_;
     }
 
     template <class Functor>
     double ContOrthoPoly1D::weightedSquaredError(const Functor& fcn,
                                                  const double *coeffs,
                                                  const unsigned maxdeg) const
     {
         if (maxdeg > maxdeg_) throw std::invalid_argument(
             "In npstat::ContOrthoPoly1D::weightedSquaredError: "
             "maximum degree out of range");
         assert(coeffs);
 
         long double sum = 0.0L;
         const unsigned long mSize = measure_.size();
         for (unsigned long i = 0; i < mSize; ++i)
         {
             const double x = measure_[i].first;
             const long double delta = normseries(coeffs, maxdeg, x) - fcn(x + meanX_);
             sum += delta*delta*measure_[i].second;
         }
         return sum/wsum_;
     }
 
     template <class Functor>
     double ContOrthoPoly1D::integratedSquaredError(
         const Functor& infcn,
         const double *coeffs, const unsigned maxdeg,
         const double xmin, const double xmax,
         unsigned integrationRulePoints) const
     {
         if (maxdeg > maxdeg_) throw std::invalid_argument(
             "In npstat::ContOrthoPoly1D::integratedSquaredError: "
             "maximum degree out of range");
         assert(coeffs);
 
         if (!integrationRulePoints)
         {
             integrationRulePoints = GaussLegendreQuadrature::minimalExactRule(2*maxdeg);
             if (!integrationRulePoints) throw std::invalid_argument(
                 "In npstat::ContOrthoPoly1D::integratedSquaredError: "
                 "order of requested integration rule is too high");
         }
         DeltaSquaredSerFcn<ContOrthoPoly1D,Functor> fcn(*this, infcn, coeffs, maxdeg);
         GaussLegendreQuadrature glq(integrationRulePoints);
         return glq.integrate(fcn, xmin, xmax);
     }
 
     template <class Functor>
     double ContOrthoPoly1D::integrateEWPolyProduct(
         const Functor& weight,
         const unsigned deg1, const unsigned deg2,
         const double xmin, const double xmax,
         const unsigned integrationRulePoints) const
     {
         if (deg1 > maxdeg_ || deg2 > maxdeg_) throw std::invalid_argument(
             "In npstat::ContOrthoPoly1D::integrateEWPolyProduct: "
             "degree argument out of range");
         GaussLegendreQuadrature glq(integrationRulePoints);
         EWPolyProductFcn<Functor> fcn(*this, weight, deg1, deg2);
+
+        // The shift of the poly argument will be performed inside "fcn"
         return glq.integrate(fcn, xmin, xmax);
     }
 
     template <typename Numeric>
     void ContOrthoPoly1D::getAllPolys(const double x, const unsigned maxdeg,
                                       Numeric *polyValues) const
     {
         if (maxdeg > maxdeg_) throw std::invalid_argument(
             "In npstat::ContOrthoPoly1D::getAllPolys: "
             "maximum degree out of range");
         assert(polyValues);
 
         long double polyk = 1.0L, polykm1 = 0.0L;
         polyValues[0] = polyk;
         for (unsigned k=0; k<maxdeg; ++k)
         {
             const long double p = ((x - rCoeffs_[k].alpha)*polyk -
                                    rCoeffs_[k].sqrbeta*polykm1)/rCoeffs_[k+1].sqrbeta;
             polyValues[k+1] = p;
             polykm1 = polyk;
             polyk = p;
         }
     }
 }