Index: trunk/npstat/nm/AbsClassicalOrthoPoly1D.hh
===================================================================
--- trunk/npstat/nm/AbsClassicalOrthoPoly1D.hh	(revision 907)
+++ trunk/npstat/nm/AbsClassicalOrthoPoly1D.hh	(revision 908)
@@ -1,630 +1,630 @@
 #ifndef NPSTAT_ABSCLASSICALORTHOPOLY1D_HH_
 #define NPSTAT_ABSCLASSICALORTHOPOLY1D_HH_
 
 /*!
 // \file AbsClassicalOrthoPoly1D.hh
 //
 // \brief Base class for classical continuous orthonormal polynomials
 //
 // Author: I. Volobouev
 //
 // May 2017
 */
 
 #include <utility>
 #include <algorithm>
 #include <climits>
 #include <vector>
 #include <stdexcept>
 
 #include "geners/CPP11_auto_ptr.hh"
 
 #include "npstat/nm/Matrix.hh"
 #include "npstat/nm/SimpleFunctors.hh"
 #include "npstat/nm/AbsIntervalQuadrature1D.hh"
 
 namespace npstat {
     class StorablePolySeries1D;
     class OrthoPoly1DDeg;
     class OrthoPoly1DWeight;
 
     class AbsClassicalOrthoPoly1D
     {
     public:
         typedef OrthoPoly1DDeg degree_functor;
         typedef OrthoPoly1DWeight weight_functor;
 
         inline virtual ~AbsClassicalOrthoPoly1D() {}
 
         /** Virtual copy constructor */
         virtual AbsClassicalOrthoPoly1D* clone() const = 0;
 
         /** The weight function should be normalized */
         virtual long double weight(long double x) const = 0;
 
         //@{
         /** Support of the polynomial */
         virtual double xmin() const = 0;
         virtual double xmax() const = 0;
         //@}
 
         /** Maximum polynomial degree supported */
         inline virtual unsigned maxDegree() const {return UINT_MAX - 1U;}
 
         /**
         // Generate principal minor of order n of the Jacobi matrix.
         // n must not exceed the output of "maxDegree()" method.
         */
         Matrix<long double> jacobiMatrix(unsigned n) const;
 
         /** Polynomial values */
         long double poly(unsigned deg, long double x) const;
 
         /** 
         // Values of two orthonormal polynomials.
         // Faster than calling "poly" two times.
         */
         std::pair<long double,long double> twopoly(
             unsigned deg1, unsigned deg2, long double x) const;
 
         /**
         // Values of all orthonormal polynomials up to some degree.
         // Faster than calling "poly" multiple times. The size of
         // the "values" array should be at least maxdeg + 1.
         */
         void allpoly(long double x, long double* values, unsigned maxdeg) const;
 
         /** Polynomial series */
         double series(const double* coeffs, unsigned maxdeg, double x) const;
 
         /** Integral of the series from xmin() to xmax() */
         double integrateSeries(const double* coeffs, unsigned maxdeg) 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, class Quadrature>
         void calculateCoeffs(const Functor& fcn, const Quadrature& quad,
                              double* coeffs, unsigned maxdeg) const;
 
         /**
         // Build the coefficients of the orthogonal polynomial series
         // for the given sample of points (empirical density 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 Numeric>
         void sampleCoeffs(const Numeric* coords, unsigned long lenCoords,
                           double* coeffs, unsigned maxdeg) const;
 
         /**
         // Estimate the variances of the coefficients returned by
         // the "sampleCoeffs" function. The "coeffs" array should be
         // at least maxdeg + 1 elements long and should be filled
         // by a previous call to "sampleCoeffs". The "variances" array
         // should have at least maxdeg + 1 elements. It will contain
         // the respective variances upon return.
         */
         template <class Numeric>
         void sampleCoeffVars(const Numeric* coords, unsigned long lenCoords,
                              const double* coeffs, unsigned maxdeg,
                              double* variances) const;
 
         /**
         // Estimate the covariances of the coefficients returned by
         // the "sampleCoeffs" function. The "coeffs" array should be
         // at least maxdeg + 1 elements long and should be filled
         // by a previous call to "sampleCoeffs". The returned matrix
         // will be dimensioned (maxdeg + 1) x (maxdeg + 1).
         */
         template <class Numeric>
         npstat::Matrix<double>
         sampleCoeffCovariance(const Numeric* coords, unsigned long lenCoords,
                               const double* coeffs, unsigned maxdeg) const;
 
         /**
         // Build the coefficients of the orthogonal polynomial series
         // for the given sample of points (empirical density 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).
         // Before calculating the coefficients, the coordinates are shifted
         // and scaled according to x_new = (x_original - location)/scale.
         // The resulting coefficients are also divided by scale.
         */
         template <class Numeric>
         void sampleCoeffs(const Numeric* coords, unsigned long lenCoords,
                           Numeric location, Numeric scale,
                           double* coeffs, unsigned maxdeg) const;
 
         /**
         // Estimate the variances of the coefficients returned by
         // the "sampleCoeffs" function. The "coeffs" array should be
         // at least maxdeg + 1 elements long and should be filled
         // by a previous call to "sampleCoeffs" with the same location
         // and scale. The "variances" array should have at least maxdeg + 1
         // elements. It will contain the respective variances upon return.
         */
         template <class Numeric>
         void sampleCoeffVars(const Numeric* coords, unsigned long lenCoords,
                              Numeric location, Numeric scale,
                              const double* coeffs, unsigned maxdeg,
                              double* variances) const;
 
         /**
         // Estimate the covariances of the coefficients returned by
         // the "sampleCoeffs" function. The "coeffs" array should be
         // at least maxdeg + 1 elements long and should be filled
         // by a previous call to "sampleCoeffs" with the same location
         // and scale. The returned matrix will be dimensioned
         // (maxdeg + 1) x (maxdeg + 1).
         */
         template <class Numeric>
         npstat::Matrix<double>
         sampleCoeffCovariance(const Numeric* coords, unsigned long lenCoords,
                               Numeric location, Numeric scale,
                               const double* coeffs, unsigned maxdeg) const;
 
         /**
         // Build the coefficients of the orthogonal polynomial series for
         // the given sample of weighted points (empirical density function).
         // The first element of the pair will be treated as the coordinate
         // and the second element will be treated as weight. Weights must
         // not be negative. 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 Pair>
         void weightedSampleCoeffs(const Pair* points, unsigned long numPoints,
                                   double* coeffs, unsigned maxdeg) const;
 
         /**
         // Estimate the variances of the coefficients returned by the
         // "weightedSampleCoeffs" function. The "coeffs" array should
         // be at least maxdeg + 1 elements long and should be filled
         // by a previous call to "weightedSampleCoeffs". The "variances"
         // array should have at least maxdeg + 1 elements. It will contain
         // the respective variances upon return. This code assumes that
         // weights and coordinates are statistically independent from
         // each other.
         */
         template <class Pair>
         void weightedSampleCoeffVars(const Pair* points, unsigned long nPoints,
                                      const double* coeffs, unsigned maxdeg,
                                      double* variances) const;
 
         /**
         // Estimate the covariances of the coefficients returned by the
         // "weightedSampleCoeffs" function. The "coeffs" array should
         // be at least maxdeg + 1 elements long and should be filled
         // by a previous call to "weightedSampleCoeffs". The returned
         // matrix will be dimensioned (maxdeg + 1) x (maxdeg + 1). This
         // code assumes that weights and coordinates are statistically
         // independent from each other.
         */
         template <class Pair>
         npstat::Matrix<double>
         weightedSampleCoeffCovariance(const Pair* points, unsigned long nPoints,
                                       const double* coeffs, unsigned maxdeg) const;
 
         /**
         // Build the coefficients of the orthogonal polynomial series for
         // the given sample of weighted points (empirical density function).
         // The first element of the pair will be treated as the coordinate
         // and the second element will be treated as weight. Weights must
         // not be negative. 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). Before calculating the coefficients, the
         // coordinates are shifted and scaled according to
         // x_new = (x_original - location)/scale. The resulting coefficients
         // are also divided by scale.
         */
         template <class Pair>
         void weightedSampleCoeffs(const Pair* points, unsigned long numPoints,
                                   typename Pair::first_type location,
                                   typename Pair::first_type scale,
                                   double* coeffs, unsigned maxdeg) const;
 
         /**
         // Estimate the variances of the coefficients returned by the
         // "weightedSampleCoeffs" function. The "coeffs" array should
         // be at least maxdeg + 1 elements long and should be filled
         // by a previous call to "weightedSampleCoeffs" with the same
         // location and scale. The "variances" array should have at least
         // maxdeg + 1 elements. It will contain the respective variances
         // upon return. This code assumes that weights and coordinates
         // are statistically independent from each other.
         */
         template <class Pair>
         void weightedSampleCoeffVars(const Pair* points, unsigned long nPoints,
                                      typename Pair::first_type location,
                                      typename Pair::first_type scale,
                                      const double* coeffs, unsigned maxdeg,
                                      double* variances) const;
 
         /**
         // Estimate the covariances of the coefficients returned by the
         // "weightedSampleCoeffs" function. The "coeffs" array should
         // be at least maxdeg + 1 elements long and should be filled
         // by a previous call to "weightedSampleCoeffs" with the same
         // location and scale. The returned matrix will be dimensioned
         // (maxdeg + 1) x (maxdeg + 1). This code assumes that weights
         // and coordinates are statistically independent from each other.
         */
         template <class Pair>
         npstat::Matrix<double>
         weightedSampleCoeffCovariance(const Pair* points, unsigned long nPoints,
                                       typename Pair::first_type location,
                                       typename Pair::first_type scale,
                                       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.
         */
         template <class Quadrature>
         double empiricalKroneckerDelta(const Quadrature& quad,
                                        unsigned deg1, unsigned deg2) const;
 
         /**
         // A measure-weighted average of a product of four orthonormal
         // poly values
         */
         template <class Quadrature>
         double jointAverage(const Quadrature& q, unsigned deg1, unsigned deg2,
                             unsigned deg3, unsigned deg4) const;
 
         /**
         // A measure-weighted average of a product of multiple orthonormal
         // poly values. "checkedForZeros" argument should be set to "true"
         // if we are sure that there are no zero elements in the "degrees"
         // array.
         */
         template <class Quadrature>
         double jointAverage(const Quadrature& q,
                             const unsigned* degrees, unsigned nDegrees,
                             bool checkedForZeros = false) const;
 
         /**
         // Direct unweighted integration of of a product of four orthonormal
         // poly values. The integration function should support the "integrate"
         // method without limits. Intended for use with GaussHermiteQuadrature
         // and such.
         */
         template <class Quadrature>
         double directQuadrature(const Quadrature& q, unsigned deg1, unsigned deg2,
                                 unsigned deg3, unsigned deg4) const;
 
         /**
         // Direct unweighted integration of of a product of multiple orthonormal
         // poly values. The integration function should support the "integrate"
         // method without limits. Intended for use with GaussHermiteQuadrature
         // and such. "checkedForZeros" argument should be set to "true" if we are
         // sure that there are no zero elements in the "degrees" array.
         */
         template <class Quadrature>
         double directQuadrature(const Quadrature& q,
                                 const unsigned* degrees, unsigned nDegrees,
                                 bool checkedForZeros = false) const;
 
         /**
         // An unweighted integral of the orthonormal polynomial with the
         // given degree to some power. For this method, it is assumed
         // that the polynomials are supported on a closed interval (without
         // such an assumption unweighted integrals do not make much sense)
         // and that Gauss-Legendre quadratures can be used.
         */
         double integratePoly(unsigned degree, unsigned power) const;
 
         /**
         // An unweighted integral of a product of multiple orthonormal
         // polynomials. For this method, it is assumed that the polynomials
         // are supported on a closed interval (without such an assumption
         // unweighted integrals do not make much sense) and that Gauss-Legendre
         // quadratures can be used.
         */
         double jointIntegral(const unsigned* degrees, unsigned nDegrees) const;
 
         /**
         // Unweighted averages of the polynomial values for the given sample.
         // The length of array "averages" should be at least maxdeg + 1.
         */
         template <class Numeric>
         void sampleAverages(const Numeric* coords, unsigned long lenCoords,
                             double* averages, unsigned maxdeg) const;
 
         /**
         // Averages of the polynomial values for the given sample of
         // weighted points (not weighting by the polynomial weight function).
         // The first element of the pair will be treated as the coordinate
         // and the second element will be treated as weight. Weights must not
         // be negative. The length of array "averages" should be at least
         // maxdeg + 1.
         */
         template <class Pair>
         void weightedPointsAverages(const Pair* points, unsigned long numPoints,
                                     double* averages, unsigned maxdeg) const;
 
         /**
         // Averages of the polynomial values weighted by an external
         // weight function. The length of array "averages" should be
         // at least maxdeg + 1.
         */
         template <class Functor>
         void extWeightAverages(const Functor& extWeight,
                                const AbsIntervalQuadrature1D& quad,
                                double* averages, unsigned maxdeg) const;
 
         /**
         // Unweighted averages of the pairwise products of the polynomial
         // values for the given sample. The returned matrix will be symmetric
         // and will have the dimensions (maxdeg + 1) x (maxdeg + 1).
         */
         template <class Numeric>
         Matrix<double> sampleProductAverages(const Numeric* coords,
                                              unsigned long lenCoords,
                                              unsigned maxdeg) const;
 
         /**
         // Pairwise products of the polynomial values weighted by
         // an external weight function. The returned matrix will be
         // symmetric and will have the dimensions (maxdeg + 1) x (maxdeg + 1).
         */
         template <class Functor>
         Matrix<double> extWeightProductAverages(const Functor& extWeight,
             const AbsIntervalQuadrature1D& quad, unsigned maxdeg) const;
 
         /**
         // Unweighted average of a product of polynomial values for the
         // given sample. "degrees" is the collection of polynomial degrees.
         // Polynomials of these degrees will be included in the product.
         */
         template <class Numeric>
         double jointSampleAverage(const Numeric* coords,
                                   unsigned long lenCoords,
                                   const unsigned* degrees,
                                   unsigned lenDegrees) const;
 
         // Similar to the previous method but calculates two averages
         // simultaneously
         template <class Numeric>
         std::pair<double, double> twoJointSampleAverages(
             const Numeric* coords, unsigned long lenCoords,
             const unsigned* degrees1, unsigned lenDegrees1,
             const unsigned* degrees2, unsigned lenDegrees2) const;
 
         CPP11_auto_ptr<StorablePolySeries1D> makeStorablePolySeries(
             unsigned maxPolyDeg,
             const double *coeffs=0, unsigned maxdeg=0, double shift=0.0) const;
 
         // Functions needed for compatibility with certain template codes
         inline long double location() const {return 0.0;}
         inline long double scale() const {return 1.0;}
 
     protected:
         inline static int kdelta(const unsigned i, const unsigned j)
             {return i == j ? 1 : 0;}
 
         // Recurrence relationship function should return alpha
         // as the first element of the pair and the square root
         // of beta as the second
         virtual std::pair<long double,long double>
         recurrenceCoeffs(unsigned deg) const = 0;
 
-    private:
         // Helper classes
         class ProdFcn : public Functor1<long double, long double>
         {
         public:
             inline ProdFcn(const AbsClassicalOrthoPoly1D& p,
                            const unsigned d1,
                            const unsigned d2)
                 : poly(p), deg1(d1), deg2(d2) {}
 
             inline virtual ~ProdFcn() {}
 
             inline long double operator()(const long double& x) const
             {
                 const std::pair<long double,long double>& p =
                     poly.twopoly(deg1, deg2, x);
                 return p.first*p.second*poly.weight(x);
             }
 
         private:
             const AbsClassicalOrthoPoly1D& poly;
             unsigned deg1;
             unsigned deg2;
         };
 
         class MultiProdFcn;
         friend class MultiProdFcn;
 
         class MultiProdFcn : public Functor1<long double, long double>
         {
         public:
             inline MultiProdFcn(const AbsClassicalOrthoPoly1D& p,
                                 const unsigned* degrees, unsigned lenDegrees,
                                 const bool includeWeight=true)
                 : poly(p),
                   degs(degrees, degrees+lenDegrees),
                   maxdeg(0),
                   weighted(includeWeight)
             {
                 if (lenDegrees)
                     maxdeg = *std::max_element(degrees, degrees+lenDegrees);
             }
 
             inline virtual ~MultiProdFcn() {}
 
             inline long double operator()(const long double& x) const
             {
                 long double w = 1.0L, p = 1.0L;
                 if (weighted)
                     w = poly.weight(x);
                 if (maxdeg)
                     p = poly.normpolyprod(&degs[0], degs.size(), maxdeg, x);
                 return w*p;
             }
 
         private:
             const AbsClassicalOrthoPoly1D& poly;
             std::vector<unsigned> degs;
             unsigned maxdeg;
             bool weighted;
         };
 
+    private:
         // For calling the function below, maxdeg should be calculated as
         // *std::max_element(degrees, degrees+nDegrees)
         long double normpolyprod(const unsigned* degrees, unsigned nDegrees,
                                  unsigned maxdeg, long double x) const;
 #ifdef SWIG
     public:
         inline double weight2(const double x) const
             {return weight(x);}
 
         inline double poly2(const unsigned deg, const double x) const
             {return poly(deg, x);}
 
         inline std::pair<double,double> twopoly2(const unsigned deg1,
                                                  const unsigned deg2,
                                                  const double x) const
         {
             const std::pair<long double,long double>& p = twopoly(deg1, deg2, x);
             return std::pair<double,double>(p.first, p.second);
         }
 
         inline std::vector<double> allpoly2(const unsigned maxdeg, const double x) const
         {
             std::vector<long double> p(maxdeg+1);
             allpoly(x, &p[0], p.size());
             return std::vector<double>(p.begin(), p.end());
         }
 
         inline StorablePolySeries1D* makeStorablePolySeries_2(
             unsigned maxPolyDeg, const std::vector<double>& coeffs, double shift=0.0) const
         {
             CPP11_auto_ptr<StorablePolySeries1D> ptr;
             if (coeffs.empty())
                 ptr = this->makeStorablePolySeries(maxPolyDeg, 0, 0, shift);
             else
                 ptr = this->makeStorablePolySeries(maxPolyDeg, &coeffs[0],
                                                    coeffs.size() - 1, shift);
             return ptr.release();
         }
 
         inline Matrix<double> jacobiMatrix_2(const unsigned n) const
             {return jacobiMatrix(n);}
 
         inline double location2() const {return location();}
 
         inline double scale2() const {return scale();}
 #endif // SWIG
     };
 
     /**
     // A functor for the weight function of the given ortho poly system.
     // The poly system is not copied, only a reference is used. It is
     // a responsibility of the user to make sure that the lifetime of
     // the poly system object exceeds the lifetime of the functor.
     */
     class OrthoPoly1DWeight : public Functor1<long double, long double>
     {
     public:
         inline OrthoPoly1DWeight(const AbsClassicalOrthoPoly1D& fcn,
                                  const long double normfactor)
             : fcn_(fcn), norm_(normfactor) {}
 
         // Separate constructor here because swig gets confused
         // by long doubles with defaults
         inline explicit OrthoPoly1DWeight(const AbsClassicalOrthoPoly1D& fcn)
             : fcn_(fcn), norm_(1.0L) {}
 
         inline virtual ~OrthoPoly1DWeight() {}
 
         inline virtual long double operator()(const long double& a) const
             {return norm_*fcn_.weight(a);}
 
     private:
         OrthoPoly1DWeight();
 
         const AbsClassicalOrthoPoly1D& fcn_;
         long double norm_;
     };
 
     /**
     // A functor for a particular degree polynomial of the given ortho
     // poly system. The poly system is not copied, only a reference is used.
     // It is a responsibility of the user to make sure that the lifetime of
     // the poly system object exceeds the lifetime of the functor.
     */
     class OrthoPoly1DDeg : public Functor1<long double, long double>
     {
     public:
         inline OrthoPoly1DDeg(const AbsClassicalOrthoPoly1D& fcn,
                               const unsigned degree,
                               const long double normfactor)
             : fcn_(fcn), norm_(normfactor), deg_(degree)
         {
             if (deg_ > fcn_.maxDegree()) throw std::invalid_argument(
                 "In npstat::OrthoPoly1DDeg::constructor: "
                 "degree argument is out of range");
         }
 
         // Separate constructor here because swig gets confused
         // by long doubles with defaults
         inline OrthoPoly1DDeg(const AbsClassicalOrthoPoly1D& fcn,
                               const unsigned degree)
             : fcn_(fcn), norm_(1.0L), deg_(degree)
         {
             if (deg_ > fcn_.maxDegree()) throw std::invalid_argument(
                 "In npstat::OrthoPoly1DDeg::constructor: "
                 "degree argument is out of range");
         }
 
         inline virtual ~OrthoPoly1DDeg() {}
 
         inline virtual long double operator()(const long double& a) const
             {return norm_*fcn_.poly(deg_, a);}
 
     private:
         OrthoPoly1DDeg();
 
         const AbsClassicalOrthoPoly1D& fcn_;
         long double norm_;
         unsigned deg_;
     };
 
     class UnweightedTwoPolyProd : public Functor1<long double, long double>
     {
     public:
         inline UnweightedTwoPolyProd(const AbsClassicalOrthoPoly1D& fcn,
                                      const unsigned deg1, const unsigned deg2)
             : fcn_(fcn), degree1_(deg1), degree2_(deg2) {}
 
         inline virtual ~UnweightedTwoPolyProd() {}
 
         inline virtual long double operator()(const long double& a) const
         {
             std::pair<long double,long double> p = fcn_.twopoly(
                 degree1_, degree2_, a);
             return p.first * p.second;
         }
 
     private:
         const AbsClassicalOrthoPoly1D& fcn_;
         unsigned degree1_;
         unsigned degree2_;
     };
 }
 
 #include "npstat/nm/AbsClassicalOrthoPoly1D.icc"
 
 #endif // NPSTAT_ABSCLASSICALORTHOPOLY1D_HH_
Index: trunk/npstat/nm/AbsClassicalOrthoPoly1D.icc
===================================================================
--- trunk/npstat/nm/AbsClassicalOrthoPoly1D.icc	(revision 907)
+++ trunk/npstat/nm/AbsClassicalOrthoPoly1D.icc	(revision 908)
@@ -1,691 +1,691 @@
 #include <cassert>
 #include <vector>
 #include <stdexcept>
 #include <algorithm>
 
 #include "npstat/nm/polyPrivateUtils.hh"
 
 namespace npstat {
     namespace Private {
         template <class SomePoly1D, class Functor>
         class ClassOrthoPolyHlp : public Functor1<long double, long double>
         {
         public:
             inline ClassOrthoPolyHlp(const SomePoly1D& p,
                                      const Functor& f, const unsigned d1)
                 : poly(p), fcn(f), deg(d1) {}
 
             inline long double operator()(const long double& x) const
                 {return fcn(x)*poly.poly(deg, x)*poly.weight(x);}
 
         private:
             const SomePoly1D& poly;
             const Functor& fcn;
             unsigned deg;
         };
     }
 
     template <class Functor, class Quadrature>
     void AbsClassicalOrthoPoly1D::calculateCoeffs(
         const Functor& fcn, const Quadrature& quad,
         double *coeffs, const unsigned maxdeg) const
     {
         assert(coeffs);
 
         const double a = xmin();
         const double b = xmax();
 
         for (unsigned i=0; i<=maxdeg; ++i)
         {
             Private::ClassOrthoPolyHlp<AbsClassicalOrthoPoly1D,Functor>
                 func(*this, fcn, i);
             coeffs[i] = quad.integrate(func, a, b);
         }
     }
 
     template <class Quadrature>
     double AbsClassicalOrthoPoly1D::empiricalKroneckerDelta(
         const Quadrature& quad, const unsigned deg1, const unsigned deg2) const
     {
-        ProdFcn fcn(*this, deg1, deg2);
+        const ProdFcn fcn(*this, deg1, deg2);
         return quad.integrate(fcn, xmin(), xmax());
     }
 
     template <class Quadrature>
     double AbsClassicalOrthoPoly1D::directQuadrature(
         const Quadrature& quad, const unsigned deg1, const unsigned deg2,
         const unsigned deg3, const unsigned deg4) const
     {
         unsigned degs[4];
         unsigned nDegs = 0;
         if (deg1)
             degs[nDegs++] = deg1;
         if (deg2)
             degs[nDegs++] = deg2;
         if (deg3)
             degs[nDegs++] = deg3;
         if (deg4)
             degs[nDegs++] = deg4;
         return directQuadrature(quad, degs, nDegs, true);
     }
 
     template <class Quadrature>
     double AbsClassicalOrthoPoly1D::directQuadrature(
         const Quadrature& quad, const unsigned* degrees,
         const unsigned nDegrees, const bool checkedForZeros) const
     {
         if (nDegrees)
         {
             assert(degrees);
             if (!checkedForZeros)
             {
                 bool allNonZero = true;
                 for (unsigned ideg=0; ideg<nDegrees && allNonZero; ++ideg)
                     allNonZero = degrees[ideg];
                 if (!allNonZero)
                 {
                     unsigned degBuf[256];
                     if (nDegrees - 1U <= sizeof(degBuf)/sizeof(degBuf[0]))
                     {
                         unsigned nonZeroDegs = 0;
                         for (unsigned ideg=0; ideg<nDegrees; ++ideg)
                             if (degrees[ideg])
                                 degBuf[nonZeroDegs++] = degrees[ideg];
                         return directQuadrature(quad, degBuf, nonZeroDegs, true);
                     }
                 }
             }
 
-            MultiProdFcn fcn(*this, degrees, nDegrees, false);
+            const MultiProdFcn fcn(*this, degrees, nDegrees, false);
             return quad.integrate(fcn);
         }
         else
         {
-            ConstValue1<double,double> one(1.0);
+            const ConstValue1<double,double> one(1.0);
             return quad.integrate(one);
         }
     }
 
     template <class Quadrature>
     double AbsClassicalOrthoPoly1D::jointAverage(
         const Quadrature& quad, const unsigned deg1, const unsigned deg2,
         const unsigned deg3, const unsigned deg4) const
     {
         unsigned degs[4];
         unsigned nDegs = 0;
         if (deg1)
             degs[nDegs++] = deg1;
         if (deg2)
             degs[nDegs++] = deg2;
         if (deg3)
             degs[nDegs++] = deg3;
         if (deg4)
             degs[nDegs++] = deg4;
         return jointAverage(quad, degs, nDegs, true);
     }
 
     template <class Quadrature>
     double AbsClassicalOrthoPoly1D::jointAverage(
         const Quadrature& quad, const unsigned* degrees,
         const unsigned nDegrees, const bool checkedForZeros) const
     {
         if (nDegrees)
         {
             assert(degrees);
 
             if (nDegrees == 1U)
                 return kdelta(degrees[0], 0U);
             if (nDegrees == 2U)
                 return kdelta(degrees[0], degrees[1]);
 
             if (!checkedForZeros)
             {
                 bool allNonZero = true;
                 for (unsigned ideg=0; ideg<nDegrees && allNonZero; ++ideg)
                     allNonZero = degrees[ideg];
                 if (!allNonZero)
                 {
                     unsigned degBuf[256];
                     if (nDegrees - 1U <= sizeof(degBuf)/sizeof(degBuf[0]))
                     {
                         unsigned nonZeroDegs = 0;
                         for (unsigned ideg=0; ideg<nDegrees; ++ideg)
                             if (degrees[ideg])
                                 degBuf[nonZeroDegs++] = degrees[ideg];
                         return jointAverage(quad, degBuf, nonZeroDegs, true);
                     }
                 }
             }
 
-            MultiProdFcn fcn(*this, degrees, nDegrees);
+            const MultiProdFcn fcn(*this, degrees, nDegrees);
             return quad.integrate(fcn, xmin(), xmax());
         }
         else
             return 1.0;
     }
 
     template <class Numeric>
     Matrix<double> AbsClassicalOrthoPoly1D::sampleProductAverages(
         const Numeric* coords, const unsigned long lenCoords,
         const unsigned maxdeg) const
     {
         if (!lenCoords) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::sampleProductAverages:"
             " empty sample");
         assert(coords);
         const unsigned long nsums = (maxdeg + 1)*(unsigned long)(maxdeg + 2)/2;
         std::vector<long double> membuf(nsums + (maxdeg + 1U), 0.0L);
         long double* poly = &membuf[0];
         long double* sums = poly + (maxdeg + 1U);
         for (unsigned long i=0; i<lenCoords; ++i)
         {
             allpoly(coords[i], poly, maxdeg);
             unsigned long isum = 0;
             for (unsigned k=0; k<=maxdeg; ++k)
                 for (unsigned j=0; j<=k; ++j, ++isum)
                     sums[isum] += poly[k]*poly[j];
         }
         Matrix<double> mat(maxdeg + 1U, maxdeg + 1U);
         unsigned long isum = 0;
         for (unsigned k=0; k<=maxdeg; ++k)
             for (unsigned j=0; j<=k; ++j, ++isum)
             {
                 const double v = sums[isum]/lenCoords;
                 mat[k][j] = v;
                 if (j != k)
                     mat[j][k] = v;
             }
         return mat;
     }
 
     template <class Functor>
     void AbsClassicalOrthoPoly1D::extWeightAverages(
         const Functor& fcn,
         const AbsIntervalQuadrature1D& quad,
         double* averages, const unsigned maxdeg) const
     {
         Private::extWeightAverages(
             *this, fcn, quad, maxdeg, xmin(), xmax(), averages);
     }
 
     template <class Functor>
     Matrix<double> AbsClassicalOrthoPoly1D::extWeightProductAverages(
         const Functor& fcn, const AbsIntervalQuadrature1D& quad,
         const unsigned maxdeg) const
     {
         return Private::extWeightProductAverages(
             *this, fcn, quad, maxdeg, xmin(), xmax());        
     }
 
     template <class Numeric>
     void AbsClassicalOrthoPoly1D::sampleAverages(
         const Numeric* coords, const unsigned long lenCoords,
         double *averages, const unsigned maxdeg) const
     {
         if (!lenCoords) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::sampleAverages: empty sample");
         assert(coords);
         assert(averages);
 
         std::vector<long double> membuf(2U*(maxdeg + 1U), 0.0L);
         long double* poly = &membuf[0];
         long double* sums = poly + (maxdeg + 1U);
 
         for (unsigned long i=0; i<lenCoords; ++i)
         {
             allpoly(coords[i], poly, maxdeg);
             for (unsigned k=0; k<=maxdeg; ++k)
                 sums[k] += poly[k];
         }
         const long double ldsz = lenCoords;
         for (unsigned k=0; k<=maxdeg; ++k)
             averages[k] = sums[k]/ldsz;
     }
 
     template <class Numeric>
     void AbsClassicalOrthoPoly1D::sampleCoeffs(
         const Numeric* coords, const unsigned long lenCoords,
         const Numeric location, const Numeric scale,
         double *coeffs, const unsigned maxdeg) const
     {
         if (!lenCoords) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::sampleCoeffs: empty sample");
         assert(coords);
         assert(coeffs);
         const Numeric zero = Numeric();
         if (scale == zero) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::sampleCoeffs: zero scale");
         std::vector<long double> membuf(2U*(maxdeg + 1U), 0.0L);
         long double* poly = &membuf[0];
         long double* sums = poly + (maxdeg + 1U);
         for (unsigned long i=0; i<lenCoords; ++i)
         {
             const long double x = (coords[i] - location)/scale;
             const long double w = weight(x);
             allpoly(x, poly, maxdeg);
             for (unsigned k=0; k<=maxdeg; ++k)
                 sums[k] += poly[k]*w;
         }
         for (unsigned k=0; k<=maxdeg; ++k)
             coeffs[k] = sums[k]/scale/lenCoords;
     }
 
     template <class Numeric>
     void AbsClassicalOrthoPoly1D::sampleCoeffVars(
         const Numeric* coords, const unsigned long lenCoords,
         const Numeric location, const Numeric scale,
         const double *coeffs, const unsigned maxdeg,
         double *variances) const
     {
         if (lenCoords < 2UL) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::sampleCoeffVars: "
             "insufficient sample size");
         assert(coords);
         assert(coeffs);
         assert(variances);
         const Numeric zero = Numeric();
         if (scale == zero) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::sampleCoeffVars: zero scale");
         std::vector<long double> membuf(2U*(maxdeg + 1U), 0.0L);
         long double* poly = &membuf[0];
         long double* sums = poly + (maxdeg + 1U);
         for (unsigned long i=0; i<lenCoords; ++i)
         {
             const long double x = (coords[i] - location)/scale;
             const long double w = weight(x)/scale;
             allpoly(x, poly, maxdeg);
             for (unsigned k=0; k<=maxdeg; ++k)
             {
                 const long double delta = poly[k]*w - coeffs[k];
                 sums[k] += delta*delta;
             }
         }
         long double lm1 = lenCoords - 1UL;
         lm1 *= lenCoords;
         for (unsigned k=0; k<=maxdeg; ++k)
             variances[k] = sums[k]/lm1;
     }
 
     template <class Numeric>
     npstat::Matrix<double>
     AbsClassicalOrthoPoly1D::sampleCoeffCovariance(
         const Numeric* coords, const unsigned long lenCoords,
         const Numeric location, const Numeric scale,
         const double *coeffs, const unsigned maxdeg) const
     {
         if (lenCoords < 2UL) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::sampleCoeffCovariance: "
             "insufficient sample size");
         assert(coords);
         assert(coeffs);
         const Numeric zero = Numeric();
         if (scale == zero) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::sampleCoeffCovariance: "
             "zero scale");
         std::vector<long double> membuf(maxdeg + 1U, 0.0L);
         long double* poly = &membuf[0];
         npstat::Matrix<long double> sums(maxdeg + 1U, maxdeg + 1U, 0);
         for (unsigned long i=0; i<lenCoords; ++i)
         {
             const long double x = (coords[i] - location)/scale;
             const long double w = weight(x)/scale;
             allpoly(x, poly, maxdeg);
             for (unsigned k=0; k<=maxdeg; ++k)
             {
                 const long double deltak = poly[k]*w - coeffs[k];
                 long double* tmp = sums[k];
                 for (unsigned m=0; m<=k; ++m)
                 {
                     const long double deltam = poly[m]*w - coeffs[m];
                     tmp[m] += deltak*deltam;
                 }
             }
         }
         for (unsigned k=0; k<=maxdeg; ++k)
             for (unsigned m=k+1; m<=maxdeg; ++m)
                 sums[k][m] = sums[m][k];
         long double lm1 = lenCoords - 1UL;
         lm1 *= lenCoords;
         sums /= lm1;
         return sums;
     }
 
     template <class Numeric>
     void AbsClassicalOrthoPoly1D::sampleCoeffs(
         const Numeric* coords, const unsigned long lenCoords,
         double *coeffs, const unsigned maxdeg) const
     {
         const Numeric zero = Numeric();
         const Numeric one = static_cast<Numeric>(1);
         sampleCoeffs(coords, lenCoords, zero, one, coeffs, maxdeg);
     }
 
     template <class Numeric>
     void AbsClassicalOrthoPoly1D::sampleCoeffVars(
         const Numeric* coords, const unsigned long lenCoords,
         const double *coeffs, const unsigned maxdeg,
         double *variances) const
     {
         const Numeric zero = Numeric();
         const Numeric one = static_cast<Numeric>(1);
         sampleCoeffVars(coords, lenCoords, zero, one,
                         coeffs, maxdeg, variances);
     }
 
     template <class Numeric>
     npstat::Matrix<double>
     AbsClassicalOrthoPoly1D::sampleCoeffCovariance(
         const Numeric* coords, unsigned long lenCoords,
         const double* coeffs, unsigned maxdeg) const
     {
         const Numeric zero = Numeric();
         const Numeric one = static_cast<Numeric>(1);
         return sampleCoeffCovariance(coords, lenCoords, zero,
                                      one, coeffs, maxdeg);
     }
 
     template <class Pair>
     void AbsClassicalOrthoPoly1D::weightedSampleCoeffs(
         const Pair* points, const unsigned long numPoints,
         const typename Pair::first_type location,
         const typename Pair::first_type scale,
         double *coeffs, const unsigned maxdeg) const
     {
         typedef typename Pair::first_type Numeric;
 
         if (!numPoints) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffs: "
             "empty sample");
         assert(points);
         assert(coeffs);
         const Numeric zero = Numeric();
         if (scale == zero) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffs: "
             "zero scale");
         std::vector<long double> membuf(2U*(maxdeg + 1U), 0.0L);
         long double* poly = &membuf[0];
         long double* sums = poly + (maxdeg + 1U);
         long double wsum = 0.0L;
         for (unsigned long i=0; i<numPoints; ++i)
         {
             const long double ptw = points[i].second;
             if (ptw < 0.0L) throw std::invalid_argument(
                 "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffs: "
                 "all weights must be non-negative");
             wsum += ptw;
             const long double x = (points[i].first - location)/scale;
             const long double w = weight(x);
             allpoly(x, poly, maxdeg);
             for (unsigned k=0; k<=maxdeg; ++k)
                 sums[k] += poly[k]*w*ptw;
         }
         if (wsum == 0.0L) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffs: "
             "at least one weight must be positive");
         for (unsigned k=0; k<=maxdeg; ++k)
             coeffs[k] = sums[k]/wsum/scale;
     }
 
     template <class Pair>
     void AbsClassicalOrthoPoly1D::weightedPointsAverages(
         const Pair* points, const unsigned long numPoints,
         double *averages, const unsigned maxdeg) const
     {
         if (!numPoints) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::weightedPointsAverages: "
             "empty sample");
         assert(points);
         assert(averages);
         std::vector<long double> membuf(2U*(maxdeg + 1U), 0.0L);
         long double* poly = &membuf[0];
         long double* sums = poly + (maxdeg + 1U);
         long double wsum = 0.0L;
         for (unsigned long i=0; i<numPoints; ++i)
         {
             const long double ptw = points[i].second;
             if (ptw < 0.0L) throw std::invalid_argument(
                 "In npstat::AbsClassicalOrthoPoly1D::weightedPointsAverages: "
                 "all weights must be non-negative");
             wsum += ptw;
             allpoly(points[i].first, poly, maxdeg);
             for (unsigned k=0; k<=maxdeg; ++k)
                 sums[k] += poly[k]*ptw;
         }
         if (wsum == 0.0L) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::weightedPointsAverages: "
             "at least one weight must be positive");
         for (unsigned k=0; k<=maxdeg; ++k)
             averages[k] = sums[k]/wsum;
     }
 
     template <class Pair>
     void AbsClassicalOrthoPoly1D::weightedSampleCoeffVars(
         const Pair* points, const unsigned long numPoints,
         typename Pair::first_type location,
         typename Pair::first_type scale,
         const  double *coeffs, const unsigned maxdeg,
         double *variances) const
     {
         typedef typename Pair::first_type Numeric;
 
         if (numPoints < 2UL) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffVars: "
             "insufficient sample size");
         assert(points);
         assert(coeffs);
         assert(variances);
         const Numeric zero = Numeric();
         if (scale == zero) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffVars: "
             "zero scale");
         std::vector<long double> membuf(2U*(maxdeg + 1U), 0.0L);
         long double* poly = &membuf[0];
         long double* sums = poly + (maxdeg + 1U);
         long double wsum = 0.0L, wsqsum = 0.0L;
         for (unsigned long i=0; i<numPoints; ++i)
         {
             const long double ptw = points[i].second;
             if (ptw < 0.0L) throw std::invalid_argument(
                 "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffVars: "
                 "all weights must be non-negative");
             wsum += ptw;
             wsqsum += ptw*ptw;
             const long double x = (points[i].first - location)/scale;
             const long double w = weight(x)/scale;
             allpoly(x, poly, maxdeg);
             for (unsigned k=0; k<=maxdeg; ++k)
             {
                 const long double delta = poly[k]*w - coeffs[k];
                 sums[k] += delta*delta*ptw;
             }
         }
         if (wsum == 0.0L) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffVars: "
             "at least one weight must be positive");
         long double lm1 = wsum*wsum/wsqsum - 1.0L;
         if (lm1 <= 0.0L) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffVars: "
             "insufficient effective sample size");
         lm1 *= wsum;
         for (unsigned k=0; k<=maxdeg; ++k)
             variances[k] = sums[k]/lm1;
     }
 
     template <class Pair>
     npstat::Matrix<double>
     AbsClassicalOrthoPoly1D::weightedSampleCoeffCovariance(
         const Pair* points, const unsigned long numPoints,
         typename Pair::first_type location,
         typename Pair::first_type scale,
         const double *coeffs, const unsigned maxdeg) const
     {
         typedef typename Pair::first_type Numeric;
 
         if (numPoints < 2UL) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffCovariance: "
             "insufficient sample size");
         assert(points);
         assert(coeffs);
         const Numeric zero = Numeric();
         if (scale == zero) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffCovariance: "
             "zero scale");
         std::vector<long double> membuf(maxdeg + 1U, 0.0L);
         long double* poly = &membuf[0];
         npstat::Matrix<long double> sums(maxdeg + 1U, maxdeg + 1U, 0);
         long double wsum = 0.0L, wsqsum = 0.0L;
         for (unsigned long i=0; i<numPoints; ++i)
         {
             const long double ptw = points[i].second;
             if (ptw < 0.0L) throw std::invalid_argument(
                 "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffCovariance: "
                 "all weights must be non-negative");
             wsum += ptw;
             wsqsum += ptw*ptw;
             const long double x = (points[i].first - location)/scale;
             const long double w = weight(x)/scale;
             allpoly(x, poly, maxdeg);
             for (unsigned k=0; k<=maxdeg; ++k)
             {
                 const long double deltak = poly[k]*w - coeffs[k];
                 long double* tmp = sums[k];
                 for (unsigned m=0; m<=k; ++m)
                 {
                     const long double deltam = poly[m]*w - coeffs[m];
                     tmp[m] += deltak*deltam*ptw;;
                 }
             }
         }
         if (wsum == 0.0L) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffCovariance: "
             "at least one weight must be positive");
         long double lm1 = wsum*wsum/wsqsum - 1.0L;
         if (lm1 <= 0.0L) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffCovariance: "
             "insufficient effective sample size");
         lm1 *= wsum;
         for (unsigned k=0; k<=maxdeg; ++k)
             for (unsigned m=k+1; m<=maxdeg; ++m)
                 sums[k][m] = sums[m][k];
         sums /= lm1;
         return sums;
     }
 
     template <class Pair>
     void AbsClassicalOrthoPoly1D::weightedSampleCoeffs(
         const Pair* points, const unsigned long numPoints,
         double *coeffs, const unsigned maxdeg) const
     {
         typedef typename Pair::first_type Numeric;
 
         const Numeric zero = Numeric();
         const Numeric one = static_cast<Numeric>(1);
         weightedSampleCoeffs(points, numPoints, zero, one, coeffs, maxdeg);
     }
 
     template <class Pair>
     void AbsClassicalOrthoPoly1D::weightedSampleCoeffVars(
         const Pair* points, const unsigned long nPoints,
         const  double *coeffs, const unsigned maxdeg,
         double *variances) const
     {
         typedef typename Pair::first_type Numeric;
 
         const Numeric zero = Numeric();
         const Numeric one = static_cast<Numeric>(1);
         weightedSampleCoeffVars(points, nPoints, zero, one,
                                 coeffs, maxdeg, variances);
     }
 
     template <class Pair>
     npstat::Matrix<double>
     AbsClassicalOrthoPoly1D::weightedSampleCoeffCovariance(
         const Pair* points, unsigned long nPoints,
         const double* coeffs, unsigned maxdeg) const
     {
         typedef typename Pair::first_type Numeric;
 
         const Numeric zero = Numeric();
         const Numeric one = static_cast<Numeric>(1);
         return weightedSampleCoeffCovariance(points, nPoints, zero,
                                              one, coeffs, maxdeg);
     }
 
     template <class Numeric>
     double AbsClassicalOrthoPoly1D::jointSampleAverage(
         const Numeric* coords, const unsigned long lenCoords,
         const unsigned* degrees, const unsigned lenDegrees) const
     {
         if (!lenCoords) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::jointSampleAverage: "
             "empty sample");
         assert(coords);
         if (lenDegrees)
         {
             assert(degrees);
             const unsigned maxdeg = *std::max_element(degrees, degrees+lenDegrees);
             std::vector<long double> membuf(maxdeg + 1U, 0.0L);
             long double* poly = &membuf[0];
             long double sum = 0.0L;
             for (unsigned long i=0; i<lenCoords; ++i)
             {
                 allpoly(coords[i], poly, maxdeg);
                 long double prod = 1.0L;
                 for (unsigned k=0; k<lenDegrees; ++k)
                     prod *= poly[degrees[k]];
                 sum += prod;
             }
             return sum/lenCoords;
         }
         else
             return 1.0;
     }
 
     template <class Numeric>
     std::pair<double, double> AbsClassicalOrthoPoly1D::twoJointSampleAverages(
         const Numeric* coords, const unsigned long lenCoords,
         const unsigned* degrees1, const unsigned lenDegrees1,
         const unsigned* degrees2, const unsigned lenDegrees2) const
     {
         if (!lenCoords) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::twoJointSampleAverages: empty sample");
         assert(coords);
 
         unsigned maxdeg1 = 0, maxdeg2 = 0;
         if (lenDegrees1)
         {
             assert(degrees1);
             maxdeg1 = *std::max_element(degrees1, degrees1+lenDegrees1);
         }
         if (lenDegrees2)
         {
             assert(degrees2);
             maxdeg2 = *std::max_element(degrees2, degrees2+lenDegrees2);
         }
         const unsigned maxdeg = std::max(maxdeg1, maxdeg2);
 
         if (lenDegrees1 || lenDegrees2)
         {
             std::vector<long double> membuf(maxdeg + 1U, 0.0L);
             long double* poly = &membuf[0];
             long double sum1 = 0.0L, sum2 = 0.0L;
             for (unsigned long i=0; i<lenCoords; ++i)
             {
                 allpoly(coords[i], poly, maxdeg);
 
                 long double prod1 = 1.0L;
                 for (unsigned k=0; k<lenDegrees1; ++k)
                     prod1 *= poly[degrees1[k]];
                 sum1 += prod1;
 
                 long double prod2 = 1.0L;
                 for (unsigned k=0; k<lenDegrees2; ++k)
                     prod2 *= poly[degrees2[k]];
                 sum2 += prod2;
             }
             return std::pair<double, double>(sum1/lenCoords, sum2/lenCoords);
         }
         else
             return std::pair<double, double>(1.0, 1.0);
     }
 }