Index: trunk/tests/test_ContOrthoPoly1D.cc
===================================================================
--- trunk/tests/test_ContOrthoPoly1D.cc	(revision 568)
+++ trunk/tests/test_ContOrthoPoly1D.cc	(revision 569)
@@ -1,179 +1,179 @@
 #include "UnitTest++.h"
 #include "test_utils.hh"
 
 #include "npstat/nm/ContOrthoPoly1D.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-8;
+        const double eps = 2.0e-8;
 
         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);
             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);
                         }
                 }
         }
     }
 }
Index: trunk/npstat/stat/orthoPoly1DVProducts.hh
===================================================================
--- trunk/npstat/stat/orthoPoly1DVProducts.hh	(revision 568)
+++ trunk/npstat/stat/orthoPoly1DVProducts.hh	(revision 569)
@@ -1,163 +1,190 @@
 #ifndef NPSTAT_ORTHOPOLY1DVPRODUCTS_HH_
 #define NPSTAT_ORTHOPOLY1DVPRODUCTS_HH_
 
 /*!
 // \file orthoPoly1DVProducts.hh
 //
 // \brief Utility functions for calculating statistical properties of
 //        1-d orthogonal polynomials
 //
 // Author: I. Volobouev
 //
 // September 2017
 */
 
 #include "npstat/nm/ArrayND.hh"
 #include "npstat/nm/AbsClassicalOrthoPoly1D.hh"
 
 namespace npstat {
     typedef std::pair<unsigned,unsigned> UUPair;
 
     // See the "Empirical Chaos Polynomials" write-up for the
     // definition of v_{ab} and various related formulae for
     // expectations and covariances
     template <class Quadrature>
     double expectedVProduct(const AbsClassicalOrthoPoly1D& poly,
                             const Quadrature& quad, unsigned long nPoints,
                             UUPair ab, UUPair cd);
 
     template <class Quadrature>
     double expectedVProduct(const AbsClassicalOrthoPoly1D& poly,
                             const Quadrature& quad, unsigned long nPoints,
                             UUPair ab, UUPair cd, UUPair ef);
 
     template <class Quadrature>
     double expectedVProduct(const AbsClassicalOrthoPoly1D& poly,
                             const Quadrature& quad, unsigned long nPoints,
                             UUPair ab, UUPair cd, UUPair ef, UUPair gh);
 
     // Covariance between v_{ab} and v_{cd}
     template <class Quadrature>
     double expectedVCovariance(const AbsClassicalOrthoPoly1D& poly,
                                const Quadrature& quad, unsigned long nPoints,
                                UUPair ab, UUPair cd);
 
     // Covariance between v_{ab}*v_{cd} and v_{ef}
     template <class Quadrature>
     double expectedVCovariance(const AbsClassicalOrthoPoly1D& poly,
                                const Quadrature& quad, unsigned long nPoints,
                                UUPair ab, UUPair cd, UUPair ef);
 
     // Covariance between v_{ab}*v_{cd} and v_{ef}*v_{gh}
     template <class Quadrature>
     double expectedVCovariance(const AbsClassicalOrthoPoly1D& poly,
                                const Quadrature& quad, unsigned long nPoints,
                                UUPair ab, UUPair cd, UUPair ef, UUPair gh);
 
     // v_{ab} calculated for a given sample of unweighted points.
     // For samples generated using the polynomial weight function
     // as the density, the expectation value of v_{ab} is 0.
     template <typename Numeric>
     double sampleVProduct(const AbsClassicalOrthoPoly1D& poly,
                           const Numeric* coords, unsigned long nCoords,
                           UUPair ab);
 
     // v_{ab}*v_{cd} calculated for a given sample of unweighted points
     template <typename Numeric>
     double sampleVProduct(const AbsClassicalOrthoPoly1D& poly,
                           const Numeric* coords, unsigned long nCoords,
                           UUPair ab, UUPair cd);
 
     // v_{ab}*v_{cd}*v_{ef} calculated for a given sample
     template <typename Numeric>
     double sampleVProduct(const AbsClassicalOrthoPoly1D& poly,
                           const Numeric* coords, unsigned long nCoords,
                           UUPair ab, UUPair cd, UUPair ef);
 
     // v_{ab}*v_{cd}*v_{ef}*v_{gh} calculated for a given sample
     template <typename Numeric>
     double sampleVProduct(const AbsClassicalOrthoPoly1D& poly,
                           const Numeric* coords, unsigned long nCoords,
                           UUPair ab, UUPair cd, UUPair ef, UUPair gh);
 
+    // Expectation of v_{ab}*v_{cd} using sample averages of oracle polys
+    template <typename Numeric>
+    double sampleVProductExp(const AbsClassicalOrthoPoly1D& poly,
+                             const Numeric* coords, unsigned long nCoords,
+                             UUPair ab, UUPair cd);
+
+    // Expectation of v_{ab}*v_{cd}*v_{ef} using sample averages
+    // of oracle polys
+    template <typename Numeric>
+    double sampleVProductExp(const AbsClassicalOrthoPoly1D& poly,
+                             const Numeric* coords, unsigned long nCoords,
+                             UUPair ab, UUPair cd, UUPair ef);
+
+    // Expectation of v_{ab}*v_{cd}*v_{ef}*v_{gh} using sample averages
+    // of oracle polys
+    template <typename Numeric>
+    double sampleVProductExp(const AbsClassicalOrthoPoly1D& poly,
+                             const Numeric* coords, unsigned long nCoords,
+                             UUPair ab, UUPair cd, UUPair ef, UUPair gh);
+
     // Covariance between v_{ab} and v_{cd} estimated using
     // polynomial averages taken from the given sample
     template <typename Numeric>
     double sampleVCovariance(const AbsClassicalOrthoPoly1D& poly,
                              const Numeric* coords, unsigned long nCoords,
                              UUPair ab, UUPair cd);
 
     // Covariance between v_{ab}*v_{cd} and v_{ef} estimated using
     // polynomial averages taken from the given sample
     template <typename Numeric>
     double sampleVCovariance(const AbsClassicalOrthoPoly1D& poly,
                              const Numeric* coords, unsigned long nCoords,
                              UUPair ab, UUPair cd, UUPair ef);
 
     // Covariance between v_{ab}*v_{cd} and v_{ef}*v_{gh} estimated
     // using polynomial averages taken from the given sample
     template <typename Numeric>
     double sampleVCovariance(const AbsClassicalOrthoPoly1D& poly,
                              const Numeric* coords, unsigned long nCoords,
                              UUPair ab, UUPair cd, UUPair ef, UUPair gh);
 
     // Oracle counterpart of the "epsExpectation" method of the
     // ContOrthoPoly1D class. If "highOrder" parameter is "true",
     // the calculations will be performed to O(N^{-3/2}). If "highOrder"
     // is "false", the calculations will be to O(N^{-1}).
     template <class Quadrature>
     double oracleEpsExpectation(const AbsClassicalOrthoPoly1D& poly,
                                 const Quadrature& quad, unsigned long nPoints,
                                 unsigned m, unsigned n, bool highOrder);
 
     // Oracle counterpart of the "epsCovariance" method of the
     // ContOrthoPoly1D class
     template <class Quadrature>
     double oracleEpsCovariance(const AbsClassicalOrthoPoly1D& poly,
                                const Quadrature& quad, unsigned long nPoints,
                                unsigned m1, unsigned n1,
                                unsigned m2, unsigned n2, bool highOrder);
 
     // Fill a 4-d array with oracle covariances. The code knows that the
     // covariances are symmetric with respect to permuting some indices, and
     // does not recalculate them when permuting indices is sufficient. The
     // array span will be maxdeg+1 in each dimension.
     template <class Quadrature>
     ArrayND<double> oracleEpsCovarianceArray(
         const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad,
         unsigned long nPoints, unsigned maxdeg, bool highOrder);
 
     // Calculate the "eps_{mn}" approximation for a given sample
     // of unweighted points (up to O(N^{-1}) or O(N^{-3/2})). In
     // this approximation, all v_{ij} values are calculated using
     // sample averages of oracle polynomials. Note that the result
     // is not the exact "eps_{mn}" for the sample but just this
-    // particular special approximation. It is useful for checking
-    // the effect of replacing the empirical polys with the oracle
-    // ones in the "eps_{mn}" approximation.
+    // particular special approximation.
     template <typename Numeric>
     double sampleEpsValue(const AbsClassicalOrthoPoly1D& poly,
                           const Numeric* coords, unsigned long nCoords,
                           unsigned m, unsigned n, bool highOrder);
 
+    // Counterpart of the "epsExpectation" method of the ContOrthoPoly1D
+    // class in which chaos polys are replaced by oracle polys. It is useful
+    // for checking the effect of replacing the empirical polys with the
+    // oracle ones in the "eps_{mn}" approximation.
+    template <typename Numeric>
+    double sampleEpsExpectation(const AbsClassicalOrthoPoly1D& poly,
+                                const Numeric* coords, unsigned long nCoords,
+                                unsigned m, unsigned n, bool highOrder);
+
     // eps_{mn} covariance approximation using polynomial product
     // averages from the given sample of points
     template <typename Numeric>
     double sampleEpsCovariance(const AbsClassicalOrthoPoly1D& poly,
                                const Numeric* coords, unsigned long nCoords,
                                unsigned m1, unsigned n1,
                                unsigned m2, unsigned n2, bool highOrder);
 
     // Fill a 4-d array with sample covariances. The code knows that the
     // covariances are symmetric with respect to permuting some indices, and
     // does not recalculate them when permuting indices is sufficient. The
     // array span will be maxdeg+1 in each dimension.
     template <typename Numeric>
     ArrayND<double> sampleEpsCovarianceArray(
         const AbsClassicalOrthoPoly1D& poly, const Numeric* coords,
         unsigned long nCoords, unsigned maxdeg, bool highOrder);
 }
 
 #include "npstat/stat/orthoPoly1DVProducts.icc"
 
 #endif // NPSTAT_ORTHOPOLY1DVPRODUCTS_HH_
Index: trunk/npstat/stat/orthoPoly1DVProducts.icc
===================================================================
--- trunk/npstat/stat/orthoPoly1DVProducts.icc	(revision 568)
+++ trunk/npstat/stat/orthoPoly1DVProducts.icc	(revision 569)
@@ -1,522 +1,722 @@
 #include <string>
 #include <cassert>
 #include <stdexcept>
 #include <algorithm>
 #include <utility>
 
 namespace npstat {
     namespace Private {
         inline void validateNPoints(const char* where,
                                     const unsigned long nPoints)
         {
             if (!nPoints)
             {
                 std::string message("In npstat::");
                 message += where;
                 message += ": sample size must be positive";
                 throw std::invalid_argument(message);
             }
         }
 
         inline int kdelta(const UUPair ab)
         {
             return ab.first == ab.second ? 1 : 0;
         }
 
         inline bool isZeroPair(const UUPair pair)
         {
             return !(pair.first || pair.second);
         }
     }
 
     template <class Quadrature>
     double expectedVProduct(const AbsClassicalOrthoPoly1D& poly,
                             const Quadrature& quad, const unsigned long nPoints,
                             const UUPair ab, const UUPair cd)
     {
         Private::validateNPoints("expectedVProduct", nPoints);
         if (Private::isZeroPair(ab) || Private::isZeroPair(cd))
             return 0.0;
         else
             return (poly.jointAverage(quad, ab.first, ab.second, cd.first, cd.second) -
                     Private::kdelta(ab)*Private::kdelta(cd))/nPoints;
     }
 
+    template <typename Numeric>
+    double sampleVProductExp(const AbsClassicalOrthoPoly1D& poly,
+                             const Numeric* coords, const unsigned long nCoords,
+                             const UUPair ab, const UUPair cd)
+    {
+        Private::validateNPoints("sampleVProductExp", nCoords);
+        assert(coords);
+        if (Private::isZeroPair(ab) || Private::isZeroPair(cd))
+            return 0.0;
+        else
+        {
+            unsigned degs[4];
+            degs[0] = ab.first;
+            degs[1] = ab.second;
+            degs[2] = cd.first;
+            degs[3] = cd.second;
+            return (poly.jointSampleAverage(coords, nCoords, degs, 4) -
+                    Private::kdelta(ab)*Private::kdelta(cd))/nCoords;
+        }
+    }
+
     template <class Quadrature>
     inline double expectedVCovariance(const AbsClassicalOrthoPoly1D& poly,
                                       const Quadrature& quad,
                                       const unsigned long nPoints,
                                       const UUPair ab, const UUPair cd)
     {
         Private::validateNPoints("expectedVCovariance", nPoints);
         return expectedVProduct(poly, quad, nPoints, ab, cd);
     }
 
     template <class Quadrature>
     double expectedVProduct(const AbsClassicalOrthoPoly1D& poly,
                             const Quadrature& quad, const unsigned long nPoints,
                             const UUPair ab, const UUPair cd, const UUPair ef)
     {
         Private::validateNPoints("expectedVProduct", nPoints);
         if (Private::isZeroPair(ab) || Private::isZeroPair(cd) ||
             Private::isZeroPair(ef))
             return 0.0;
 
         const double pcdef = poly.jointAverage(quad, cd.first, cd.second, ef.first, ef.second);
         const double pabef = poly.jointAverage(quad, ab.first, ab.second, ef.first, ef.second);
         const double pabcd = poly.jointAverage(quad, ab.first, ab.second, cd.first, cd.second);
 
         unsigned degs[6];
         degs[0] = ab.first;
         degs[1] = ab.second;
         degs[2] = cd.first;
         degs[3] = cd.second;
         degs[4] = ef.first;
         degs[5] = ef.second;
         const double pabcdef = poly.jointAverage(quad, degs, 6);
 
         const int deltaprod = Private::kdelta(ab)*Private::kdelta(cd)*
                               Private::kdelta(ef);
         return (pabcdef - Private::kdelta(ab)*pcdef - Private::kdelta(cd)*pabef -
                 Private::kdelta(ef)*pabcd + 2.0*deltaprod)/nPoints/nPoints;
     }
 
+    template <typename Numeric>
+    double sampleVProductExp(const AbsClassicalOrthoPoly1D& poly,
+                             const Numeric* coords, const unsigned long nPoints,
+                             const UUPair ab, const UUPair cd, const UUPair ef)
+    {
+        Private::validateNPoints("sampleVProductExp", nPoints);
+        assert(coords);
+
+        if (Private::isZeroPair(ab) || Private::isZeroPair(cd) ||
+            Private::isZeroPair(ef))
+            return 0.0;
+
+        unsigned degs[6];
+        degs[0] = cd.first;
+        degs[1] = cd.second;
+        degs[2] = ef.first;
+        degs[3] = ef.second;
+        const double pcdef = poly.jointSampleAverage(coords, nPoints, degs, 4);
+
+        degs[0] = ab.first;
+        degs[1] = ab.second;
+        const double pabef = poly.jointSampleAverage(coords, nPoints, degs, 4);
+
+        degs[0] = ab.first;
+        degs[1] = ab.second;
+        degs[2] = cd.first;
+        degs[3] = cd.second;
+        const double pabcd = poly.jointSampleAverage(coords, nPoints, degs, 4);
+
+        degs[4] = ef.first;
+        degs[5] = ef.second;
+        const double pabcdef = poly.jointSampleAverage(coords, nPoints, degs, 6);
+
+        const int deltaprod = Private::kdelta(ab)*Private::kdelta(cd)*
+                              Private::kdelta(ef);
+        return (pabcdef - Private::kdelta(ab)*pcdef - Private::kdelta(cd)*pabef -
+                Private::kdelta(ef)*pabcd + 2.0*deltaprod)/nPoints/nPoints;
+    }
+
     template <class Quadrature>
     inline double expectedVCovariance(const AbsClassicalOrthoPoly1D& poly,
                                       const Quadrature& quad,
                                       const unsigned long nPoints,
                                       const UUPair ab, const UUPair cd,
                                       const UUPair ef)
     {
         Private::validateNPoints("expectedVCovariance", nPoints);
         return expectedVProduct(poly, quad, nPoints, ab, cd, ef);
     }
 
     template <class Quadrature>
     double expectedVProduct(const AbsClassicalOrthoPoly1D& poly,
                             const Quadrature& quad,
                             const unsigned long nPoints,
                             const UUPair ab, const UUPair cd,
                             const UUPair ef, const UUPair gh)
     {
         Private::validateNPoints("expectedVProduct", nPoints);
         if (Private::isZeroPair(ab) || Private::isZeroPair(cd) ||
             Private::isZeroPair(ef) || Private::isZeroPair(gh))
             return 0.0;
 
         const double pabcd = poly.jointAverage(quad, ab.first, ab.second, cd.first, cd.second);
         const double pefgh = poly.jointAverage(quad, ef.first, ef.second, gh.first, gh.second);
         const double pabef = poly.jointAverage(quad, ab.first, ab.second, ef.first, ef.second);
         const double pcdgh = poly.jointAverage(quad, cd.first, cd.second, gh.first, gh.second);
         const double pabgh = poly.jointAverage(quad, ab.first, ab.second, gh.first, gh.second);
         const double pcdef = poly.jointAverage(quad, cd.first, cd.second, ef.first, ef.second);
 
         unsigned degs[8];
         degs[0] = ab.first;
         degs[1] = ab.second;
         degs[2] = cd.first;
         degs[3] = cd.second;
         degs[4] = ef.first;
         degs[5] = ef.second;
         degs[6] = gh.first;
         degs[7] = gh.second;
         const double pabcdefgh = poly.jointAverage(quad, degs, 8);
         const double pabcdef = poly.jointAverage(quad, degs, 6);
 
         degs[4] = gh.first;
         degs[5] = gh.second;
         const double pabcdgh = poly.jointAverage(quad, degs, 6);
 
         degs[2] = ef.first;
         degs[3] = ef.second;
         const double pabefgh = poly.jointAverage(quad, degs, 6);
 
         degs[0] = cd.first;
         degs[1] = cd.second;
         const double pcdefgh = poly.jointAverage(quad, degs, 6);
 
         const int deltaprod = Private::kdelta(ab)*Private::kdelta(cd)*
                               Private::kdelta(ef)*Private::kdelta(gh);
 
         const double tmp = Private::kdelta(ef)*Private::kdelta(gh)*pabcd +
             Private::kdelta(cd)*Private::kdelta(gh)*pabef +
             Private::kdelta(cd)*Private::kdelta(ef)*pabgh +
             Private::kdelta(ab)*Private::kdelta(cd)*pefgh +
             Private::kdelta(ab)*Private::kdelta(ef)*pcdgh +
             Private::kdelta(ab)*Private::kdelta(gh)*pcdef;
 
         const double term1 = pabcd*pefgh + pabef*pcdgh +
             pabgh*pcdef + 3.0*deltaprod - tmp;
 
         const double term2 = pabcdefgh - 6.0*deltaprod -
             Private::kdelta(ab)*pcdefgh - Private::kdelta(cd)*pabefgh -
             Private::kdelta(ef)*pabcdgh - Private::kdelta(gh)*pabcdef -
             pabcd*pefgh - pabef*pcdgh - pabgh*pcdef + 2.0*tmp;
 
         return (term1 + term2/nPoints)/nPoints/nPoints;
     }
 
+    template <typename Numeric>
+    double sampleVProductExp(const AbsClassicalOrthoPoly1D& poly,
+                             const Numeric* coords, const unsigned long nPoints,
+                             const UUPair ab, const UUPair cd,
+                             const UUPair ef, const UUPair gh)
+    {
+        Private::validateNPoints("sampleVProductExp", nPoints);
+        assert(coords);
+
+        if (Private::isZeroPair(ab) || Private::isZeroPair(cd) ||
+            Private::isZeroPair(ef) || Private::isZeroPair(gh))
+            return 0.0;
+
+        std::pair<double, double> dp;
+        unsigned degs[8];
+        degs[0] = ab.first;
+        degs[1] = ab.second;
+        degs[2] = cd.first;
+        degs[3] = cd.second;
+        degs[4] = ef.first;
+        degs[5] = ef.second;
+        degs[6] = gh.first;
+        degs[7] = gh.second;
+        dp = poly.twoJointSampleAverages(coords, nPoints, degs, 4, degs+4, 4);
+        const double pabcd = dp.first;
+        const double pefgh = dp.second;
+        
+        degs[2] = ef.first;
+        degs[3] = ef.second;
+        degs[4] = cd.first;
+        degs[5] = cd.second;
+        dp = poly.twoJointSampleAverages(coords, nPoints, degs, 4, degs+4, 4);
+        const double pabef = dp.first;
+        const double pcdgh = dp.second;
+
+        degs[2] = gh.first;
+        degs[3] = gh.second;
+        degs[6] = ef.first;
+        degs[7] = ef.second;
+        dp = poly.twoJointSampleAverages(coords, nPoints, degs, 4, degs+4, 4);
+        const double pabgh = dp.first;
+        const double pcdef = dp.second;
+
+        degs[2] = cd.first;
+        degs[3] = cd.second;
+        degs[4] = ef.first;
+        degs[5] = ef.second;
+        degs[6] = gh.first;
+        degs[7] = gh.second;
+        dp = poly.twoJointSampleAverages(coords, nPoints, degs, 8, degs, 6);
+        const double pabcdefgh = dp.first;
+        const double pabcdef = dp.second;
+
+        degs[4] = gh.first;
+        degs[5] = gh.second;
+        const double pabcdgh = poly.jointSampleAverage(coords, nPoints, degs, 6);
+
+        degs[2] = ef.first;
+        degs[3] = ef.second;
+        const double pabefgh = poly.jointSampleAverage(coords, nPoints, degs, 6);
+
+        degs[0] = cd.first;
+        degs[1] = cd.second;
+        const double pcdefgh = poly.jointSampleAverage(coords, nPoints, degs, 6);
+
+        const int deltaprod = Private::kdelta(ab)*Private::kdelta(cd)*
+                              Private::kdelta(ef)*Private::kdelta(gh);
+
+        const double tmp = Private::kdelta(ef)*Private::kdelta(gh)*pabcd +
+            Private::kdelta(cd)*Private::kdelta(gh)*pabef +
+            Private::kdelta(cd)*Private::kdelta(ef)*pabgh +
+            Private::kdelta(ab)*Private::kdelta(cd)*pefgh +
+            Private::kdelta(ab)*Private::kdelta(ef)*pcdgh +
+            Private::kdelta(ab)*Private::kdelta(gh)*pcdef;
+
+        const double term1 = pabcd*pefgh + pabef*pcdgh +
+            pabgh*pcdef + 3.0*deltaprod - tmp;
+
+        const double term2 = pabcdefgh - 6.0*deltaprod -
+            Private::kdelta(ab)*pcdefgh - Private::kdelta(cd)*pabefgh -
+            Private::kdelta(ef)*pabcdgh - Private::kdelta(gh)*pabcdef -
+            pabcd*pefgh - pabef*pcdgh - pabgh*pcdef + 2.0*tmp;
+
+        return (term1 + term2/nPoints)/nPoints/nPoints;
+    }
+
     template <class Quadrature>
     inline double expectedVCovariance(const AbsClassicalOrthoPoly1D& poly,
                                       const Quadrature& quad,
                                       const unsigned long nPoints,
                                       const UUPair ab, const UUPair cd,
                                       const UUPair ef, const UUPair gh)
     {
         Private::validateNPoints("expectedVCovariance", nPoints);
         return expectedVProduct(poly, quad, nPoints, ab, cd, ef, gh) -
                expectedVProduct(poly, quad, nPoints, ab, cd)*
                expectedVProduct(poly, quad, nPoints, ef, gh);
     }
 
     template <typename Numeric>
     inline double sampleVProduct(const AbsClassicalOrthoPoly1D& poly,
                                  const Numeric* coords,
                                  const unsigned long nPoints,
                                  const UUPair ab)
     {
         Private::validateNPoints("sampleVProduct", nPoints);
         assert(coords);
         if (Private::isZeroPair(ab))
             return 0.0;
         else
         {
             unsigned degs[2];
             degs[0] = ab.first;
             degs[1] = ab.second;
 
             return poly.jointSampleAverage(coords, nPoints, degs, 2) -
                    Private::kdelta(ab);
         }
     }
 
     template <typename Numeric>
     inline double sampleVProduct(const AbsClassicalOrthoPoly1D& poly,
                                  const Numeric* coords,
                                  const unsigned long nPoints,
                                  const UUPair ab, const UUPair cd)
     {
         Private::validateNPoints("sampleVProduct", nPoints);
         assert(coords);
         if (Private::isZeroPair(ab) || Private::isZeroPair(cd))
             return 0.0;
         else
         {
             unsigned degs1[2];
             degs1[0] = ab.first;
             degs1[1] = ab.second;
 
             unsigned degs2[2];
             degs2[0] = cd.first;
             degs2[1] = cd.second;
 
             const std::pair<double, double>& dp = poly.twoJointSampleAverages(
                 coords, nPoints, degs1, 2, degs2, 2);
             return (dp.first - Private::kdelta(ab))*
                    (dp.second - Private::kdelta(cd));
         }
     }
 
     template <typename Numeric>
     inline double sampleVCovariance(const AbsClassicalOrthoPoly1D& poly,
                                     const Numeric* points,
                                     const unsigned long nPoints,
                                     const UUPair ab, const UUPair cd)
     {
         Private::validateNPoints("sampleVCovariance", nPoints);
-        return sampleVProduct(poly, points, nPoints, ab, cd);
+        return sampleVProductExp(poly, points, nPoints, ab, cd);
     }
 
     template <typename Numeric>
     inline double sampleVProduct(const AbsClassicalOrthoPoly1D& poly,
                                  const Numeric* coords,
                                  const unsigned long nPoints,
                                  const UUPair ab, const UUPair cd,
                                  const UUPair ef)
     {
         Private::validateNPoints("sampleVProduct", nPoints);
         assert(coords);
         if (Private::isZeroPair(ab) || Private::isZeroPair(cd) ||
             Private::isZeroPair(ef))
             return 0.0;
         else
             return sampleVProduct(poly, coords, nPoints, ab, cd)*
                    sampleVProduct(poly, coords, nPoints, ef);
     }
 
     template <typename Numeric>
     inline double sampleVCovariance(const AbsClassicalOrthoPoly1D& poly,
                                     const Numeric* points,
                                     const unsigned long nPoints,
                                     const UUPair ab, const UUPair cd,
                                     const UUPair ef)
     {
         Private::validateNPoints("sampleVCovariance", nPoints);
-        return sampleVProduct(poly, points, nPoints, ab, cd, ef);
+        return sampleVProductExp(poly, points, nPoints, ab, cd, ef);
     }
 
     template <typename Numeric>
     inline double sampleVProduct(const AbsClassicalOrthoPoly1D& poly,
                                  const Numeric* coords,
                                  const unsigned long nPoints,
                                  const UUPair ab, const UUPair cd,
                                  const UUPair ef, const UUPair gh)
     {
         Private::validateNPoints("sampleVProduct", nPoints);
         assert(coords);
         if (Private::isZeroPair(ab) || Private::isZeroPair(cd) ||
             Private::isZeroPair(ef) || Private::isZeroPair(gh))
             return 0.0;
         else
             return sampleVProduct(poly, coords, nPoints, ab, cd)*
                    sampleVProduct(poly, coords, nPoints, ef, gh);
     }
 
     template <typename Numeric>
     inline double sampleVCovariance(const AbsClassicalOrthoPoly1D& poly,
                                     const Numeric* points,
                                     const unsigned long nPoints,
                                     const UUPair ab, const UUPair cd,
                                     const UUPair ef, const UUPair gh)
     {
         Private::validateNPoints("sampleVCovariance", nPoints);
-        return sampleVProduct(poly, points, nPoints, ab, cd, ef, gh) -
-               sampleVProduct(poly, points, nPoints, ab, cd)*
-               sampleVProduct(poly, points, nPoints, ef, gh);
+        assert(points);
+        return sampleVProductExp(poly, points, nPoints, ab, cd, ef, gh) -
+               sampleVProductExp(poly, points, nPoints, ab, cd)*
+               sampleVProductExp(poly, points, nPoints, ef, gh);
     }
 
     template <class Quadrature>
     double oracleEpsExpectation(const AbsClassicalOrthoPoly1D& poly,
                                 const Quadrature& quad,
                                 const unsigned long nPoints,
                                 const unsigned m_in, const unsigned n_in,
                                 const bool highOrder)
     {
         if (highOrder)
         {
             Private::validateNPoints("oracleEpsExpectation", nPoints);
 
             long double sum = 0.0L;
             if (m_in || n_in)
             {
                 const unsigned m = std::min(m_in, n_in);
                 const unsigned n = std::max(m_in, n_in);
 
                 for (unsigned k=0; k<=n; ++k)
                 {
                     const double f = k == n ? 1.0 : (k > m ? 1.0 : 2.0);
                     sum += f*poly.jointAverage(quad, k, k, m, n);
                 }
                 if (m == n)
                     sum -= 1.0;
                 else
                     sum -= (poly.jointAverage(quad, m, m, m, n) +
                             poly.jointAverage(quad, m, n, n, n))/2.0;
             }
             return sum/nPoints;
         }
         else
             return 0.0;
     }
 
     template <typename Numeric>
+    double sampleEpsExpectation(const AbsClassicalOrthoPoly1D& poly,
+                                const Numeric* coords, const unsigned long nPoints,
+                                const unsigned m_in, const unsigned n_in,
+                                const bool highOrder)
+    {
+        if (highOrder)
+        {
+            Private::validateNPoints("sampleEpsExpectation", nPoints);
+            assert(coords);
+
+            long double sum = 0.0L;
+            if (m_in || n_in)
+            {
+                const unsigned m = std::min(m_in, n_in);
+                const unsigned n = std::max(m_in, n_in);
+
+                unsigned degs1[4];
+                degs1[2] = m;
+                degs1[3] = n;
+                
+                for (unsigned k=0; k<=n; ++k)
+                {
+                    const double f = k == n ? 1.0 : (k > m ? 1.0 : 2.0);
+                    degs1[0] = k;
+                    degs1[1] = k;
+                    sum += f*poly.jointSampleAverage(coords, nPoints, degs1, 4);
+                }
+
+                if (m == n)
+                    sum -= 1.0;
+                else
+                {
+                    degs1[0] = m;
+                    degs1[1] = m;
+
+                    unsigned degs2[4];
+                    degs2[0] = m;
+                    degs2[1] = n;
+                    degs2[2] = n;
+                    degs2[3] = n;
+
+                    const std::pair<double,double>& aves = poly.twoJointSampleAverages(
+                        coords, nPoints, degs1, 4, degs2, 4);
+                    sum -= (aves.first + aves.second)/2.0;
+                }
+            }
+            return sum/nPoints;
+        }
+        else
+            return 0.0;
+    }
+
+    template <typename Numeric>
     double sampleEpsValue(const AbsClassicalOrthoPoly1D& poly,
                           const Numeric* coords, unsigned long nPoints,
                           const unsigned m_in, const unsigned n_in,
                           const bool highOrder)
     {
         Private::validateNPoints("sampleEpsValue", nPoints);
         assert(coords);
 
         if (highOrder)
         {
             long double sum = 0.0L;
             if (m_in || n_in)
             {
                 const unsigned m = std::min(m_in, n_in);
                 const unsigned n = std::max(m_in, n_in);
 
                 const double vmn = sampleVProduct(poly, coords, nPoints,
                                                   UUPair(m, n));
                 if (m == n)
                     sum -= (vmn + vmn*vmn);
                 else
                 {
                     const double vnn = sampleVProduct(poly, coords, nPoints,
                                                       UUPair(n, n));
                     const double vmm = sampleVProduct(poly, coords, nPoints,
                                                       UUPair(m, m));
                     sum -= (vmn + 0.5*vmn*(vnn + vmm));
                 }
                 for (unsigned k=0; k<=n; ++k)
                 {
                     const double f = k > m ? 1.0 : 2.0;
                     const double vkmkn = sampleVProduct(poly, coords, nPoints,
                                                         UUPair(k, m), UUPair(k, n));
                     sum += f*vkmkn;
                 }
             }
             return sum;
         }
         else
             return -sampleVProduct(poly, coords, nPoints, UUPair(m_in, n_in));
     }
 
 #define npstat_EpsCovarianceCalculation_mainbody(calculator, argument) do {             \
         const bool has0 = (m1_in == 0 && n1_in == 0) ||                                 \
                           (m2_in == 0 && n2_in == 0);                                   \
         if (has0)                                                                       \
             return 0.0;                                                                 \
         if (highOrder)                                                                  \
         {                                                                               \
             const unsigned m1 = std::min(m1_in, n1_in);                                 \
             const unsigned n1 = std::max(m1_in, n1_in);                                 \
             const unsigned m2 = std::min(m2_in, n2_in);                                 \
             const unsigned n2 = std::max(m2_in, n2_in);                                 \
             long double sum = 0.0L;                                                     \
             /* Process the -v_{m1,n1} term */                                           \
             sum += calculator (poly, argument, nPoints,                                 \
                                UUPair(m2, n2), UUPair(m1, n1));                         \
             sum += calculator (poly, argument, nPoints,                                 \
                                UUPair(m2, n2), UUPair(n2, n2), UUPair(m1, n1))/2.0;     \
             sum += calculator (poly, argument, nPoints,                                 \
                                UUPair(m2, n2), UUPair(m2, m2), UUPair(m1, n1))/2.0;     \
             for (unsigned k=0; k<=n2; ++k)                                              \
             {                                                                           \
                 const double factor = k > m2 ? 1.0 : 2.0;                               \
                 sum -= factor*calculator (poly, argument, nPoints,                      \
                                           UUPair(k, m2), UUPair(k, n2), UUPair(m1, n1));\
             }                                                                           \
             /* Process the term -v_{m1,n1}/2 (v_{n1,n1} + v_{m1,m1}) */                 \
             unsigned idx[2];                                                            \
             idx[0] = n1;                                                                \
             idx[1] = m1;                                                                \
             for (unsigned ii=0; ii<2; ++ii)                                             \
             {                                                                           \
                 const unsigned mOrN = idx[ii];                                          \
                 sum += calculator (poly, argument, nPoints,                             \
                                    UUPair(m1, n1), UUPair(mOrN, mOrN),                  \
                                    UUPair(m2, n2))/2.0;                                 \
                 sum += calculator (poly, argument, nPoints,                             \
                                    UUPair(m1, n1), UUPair(mOrN, mOrN),                  \
                                    UUPair(m2, n2), UUPair(n2, n2))/4.0;                 \
                 sum += calculator (poly, argument, nPoints,                             \
                                    UUPair(m1, n1), UUPair(mOrN, mOrN),                  \
                                    UUPair(m2, n2), UUPair(m2, m2))/4.0;                 \
                 for (unsigned k=0; k<=n2; ++k)                                          \
                 {                                                                       \
                     const double factor = k > m2 ? 1.0 : 2.0;                           \
                     sum -= factor/2.0*calculator (poly, argument, nPoints,              \
                                                   UUPair(m1, n1), UUPair(mOrN, mOrN),   \
                                                   UUPair(k, m2), UUPair(k, n2));        \
                 }                                                                       \
             }                                                                           \
             /* Process the sum */                                                       \
             for (unsigned k1=0; k1<=n1; ++k1)                                           \
             {                                                                           \
                 const double f1 = k1 > m1 ? 1.0 : 2.0;                                  \
                 sum -= f1*calculator (poly, argument, nPoints,                          \
                                       UUPair(k1, m1), UUPair(k1, n1), UUPair(m2, n2));  \
                 sum -= f1*calculator (poly, argument, nPoints,                          \
                                       UUPair(k1, m1), UUPair(k1, n1),                   \
                                       UUPair(m2, n2), UUPair(n2, n2))/2.0;              \
                 sum -= f1*calculator (poly, argument, nPoints,                          \
                                       UUPair(k1, m1), UUPair(k1, n1),                   \
                                       UUPair(m2, n2), UUPair(m2, m2))/2.0;              \
                 for (unsigned k=0; k<=n2; ++k)                                          \
                 {                                                                       \
                     const double factor = k > m2 ? 1.0 : 2.0;                           \
                     sum += f1*factor*calculator (poly, argument, nPoints,               \
                                                  UUPair(k1, m1), UUPair(k1, n1),        \
                                                  UUPair(k, m2), UUPair(k, n2));         \
                 }                                                                       \
             }                                                                           \
             return sum;                                                                 \
         }                                                                               \
         else                                                                            \
             return calculator (poly, argument, nPoints,                                 \
                                UUPair(m2_in, n2_in), UUPair(m1_in, n1_in));             \
 } while(0);
 
     template <class Quadrature>
     double oracleEpsCovariance(const AbsClassicalOrthoPoly1D& poly,
                                const Quadrature& quad,
                                const unsigned long nPoints,
                                const unsigned m1_in, const unsigned n1_in,
                                const unsigned m2_in, const unsigned n2_in,
                                const bool highOrder)
     {
         Private::validateNPoints("oracleEpsCovariance", nPoints);
         npstat_EpsCovarianceCalculation_mainbody(expectedVCovariance, quad);
     }
 
     template <typename Numeric>
     double sampleEpsCovariance(const AbsClassicalOrthoPoly1D& poly,
                                const Numeric* coords,
                                const unsigned long nPoints,
                                const unsigned m1_in, const unsigned n1_in,
                                const unsigned m2_in, const unsigned n2_in,
                                const bool highOrder)
     {
         Private::validateNPoints("sampleEpsCovariance", nPoints);
         npstat_EpsCovarianceCalculation_mainbody(sampleVCovariance, coords);
     }
 
 #define npstat_EpsCovarianceArray_mainbody(calculator, argument) do {                    \
         const unsigned mdp1 = maxdeg + 1;                                                \
         ArrayND<double> arr(mdp1, mdp1, mdp1, mdp1);                                     \
         unsigned long nDefined = 0;                                                      \
         for (unsigned i=0; i<=maxdeg; ++i)                                               \
             for (unsigned j=0; j<=i; ++j)                                                \
                 for (unsigned k=0; k<=maxdeg; ++k)                                       \
                     for (unsigned m=0; m<=k; ++m)                                        \
                         if (std::make_pair(j,i) <= std::make_pair(m,k))                  \
                         {                                                                \
                             arr(j, i, m, k) = calculator (                               \
                                 poly, argument, nPoints, j, i, m, k, highOrder);         \
                             ++nDefined;                                                  \
                         }                                                                \
         for (unsigned i=0; i<=maxdeg; ++i)                                               \
             for (unsigned j=0; j<=maxdeg; ++j)                                           \
                 for (unsigned k=0; k<=maxdeg; ++k)                                       \
                     for (unsigned m=0; m<=maxdeg; ++m)                                   \
                        if (!(j<=i && m<=k && std::make_pair(j,i) <= std::make_pair(m,k)))\
                        {                                                                 \
                            unsigned i1 = i;                                              \
                            unsigned j1 = j;                                              \
                            unsigned k1 = k;                                              \
                            unsigned m1 = m;                                              \
                            if (!(j1 <= i1))                                              \
                                std::swap(i1, j1);                                        \
                            if (!(m1 <= k1))                                              \
                                std::swap(m1, k1);                                        \
                            if (!(std::make_pair(j1,i1) <= std::make_pair(m1,k1)))        \
                            {                                                             \
                                std::swap(m1, j1);                                        \
                                std::swap(k1, i1);                                        \
                            }                                                             \
                            arr(j, i, m, k) = arr(j1, i1, m1, k1);                        \
                            ++nDefined;                                                   \
                        }                                                                 \
         assert(nDefined == arr.length());                                                \
         return arr;                                                                      \
 } while (0);
 
     template <class Quadrature>
     ArrayND<double> oracleEpsCovarianceArray(
         const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad,
         const unsigned long nPoints, const unsigned maxdeg, const bool highOrder)
     {
         Private::validateNPoints("oracleEpsCovarianceArray", nPoints);
         npstat_EpsCovarianceArray_mainbody(oracleEpsCovariance, quad);
     }
 
     template <typename Numeric>
     ArrayND<double> sampleEpsCovarianceArray(
         const AbsClassicalOrthoPoly1D& poly, const Numeric* coords,
         const unsigned long nPoints, const unsigned maxdeg, const bool highOrder)
     {
         Private::validateNPoints("sampleEpsCovarianceArray", nPoints);
         npstat_EpsCovarianceArray_mainbody(sampleEpsCovariance, coords);
     }
 }