Index: trunk/npstat/nm/ContOrthoPoly1D.cc
===================================================================
--- trunk/npstat/nm/ContOrthoPoly1D.cc	(revision 909)
+++ trunk/npstat/nm/ContOrthoPoly1D.cc	(revision 910)
@@ -1,978 +1,979 @@
 #include "npstat/nm/ContOrthoPoly1D.hh"
 #include "npstat/nm/allocators.hh"
 #include "npstat/nm/StorablePolySeries1D.hh"
 #include "npstat/nm/Triple.hh"
 
 inline static int kdelta(const unsigned i, const unsigned j)
 {
     return i == j ? 1 : 0;
 }
 
 namespace npstat {
     const ContOrthoPoly1D::precise_type ContOrthoPoly1D::precise_zero = 0.0;
 
     void ContOrthoPoly1D::calcMeanXandWeightSums()
     {
         // Calculate wsum_, wsumsq_, and meanX_
         const unsigned long mSize = measure_.size();
         precise_type xsum = precise_zero;
         if (allWeightsEqual_)
         {
             const precise_type lw0 = measure_[0].second;
             wsum_ = mSize*lw0;
             wsumsq_ = mSize*lw0*lw0;
             for (unsigned long i = 0; i < mSize; ++i)
             {
                 const double x = measure_[i].first;
                 xsum += x;
                 if (x < minX_)
                     minX_ = x;
                 if (x > maxX_)
                     maxX_ = x;
             }
             meanX_ = xsum/mSize;
         }
         else
         {
             for (unsigned long i = 0; i < mSize; ++i)
             {
                 const double x = measure_[i].first;
                 const precise_type w = measure_[i].second;
                 wsum_ += w;
                 wsumsq_ += w*w;
                 xsum += w*x;
                 if (x < minX_)
                     minX_ = x;
                 if (x > maxX_)
                     maxX_ = x;
             }
             assert(wsum_ > precise_zero);
             meanX_ = xsum/wsum_;
         }
 
         // Shift the measure
         double shift = meanX_;
         for (unsigned long i = 0; i < mSize; ++i)
             measure_[i].first -= shift;
 
         // Try to improve the mean
         xsum = precise_zero;
         precise_type dx;
         if (allWeightsEqual_)
         {
             for (unsigned long i = 0; i < mSize; ++i)
                 xsum += measure_[i].first;
             dx = xsum/mSize;
         }
         else
         {
             for (unsigned long i = 0; i < mSize; ++i)
             {
                 const precise_type w = measure_[i].second;
                 xsum += w*measure_[i].first;
             }
             dx = xsum/wsum_;
         }
 
         // Shift the measure and the mean
         meanX_ += dx;
         shift = dx;
         for (unsigned long i = 0; i < mSize; ++i)
             measure_[i].first -= shift;
     }
 
     void ContOrthoPoly1D::calcRecurrenceCoeffs(const OrthoPolyMethod m)
     {
         rCoeffs_.clear();
         rCoeffs_.reserve(maxdeg_ + 1);
         switch (m)
         {
         case OPOLY_STIELTJES:
             calcRecurrenceStieltjes();
             break;
 
         case OPOLY_LANCZOS:
             calcRecurrenceLanczos();
             break;
 
         default:
             throw std::runtime_error(
                 "In npstat::ContOrthoPoly1D::calcRecurrenceCoeffs: "
                 "incomplete switch statement. This is a bug. Please report.");
         }
         assert(rCoeffs_.size() == maxdeg_ + 1);
     }
 
     double ContOrthoPoly1D::effectiveSampleSize() const
     {
         if (allWeightsEqual_)
             return measure_.size();
         else
             return wsum_*wsum_/wsumsq_;
     }
 
     ContOrthoPoly1D::precise_type ContOrthoPoly1D::monicpoly(
         const unsigned degree, const precise_type x) const
     {
         assert(degree <= maxdeg_);
 
         precise_type polyk = 1.0, polykm1 = precise_zero;
         for (unsigned k=0; k<degree; ++k)
         {
             const precise_type p = (x - rCoeffs_[k].alpha)*polyk -
                                    rCoeffs_[k].beta*polykm1;
             polykm1 = polyk;
             polyk = p;
         }
         return polyk;
     }
 
     ContOrthoPoly1D::precise_type ContOrthoPoly1D::normpoly(
         const unsigned degree, const precise_type x) const
     {
         assert(degree <= maxdeg_);
 
         precise_type polyk = 1.0, polykm1 = precise_zero;
         for (unsigned k=0; k<degree; ++k)
         {
             const precise_type p = ((x - rCoeffs_[k].alpha)*polyk -
                                    rCoeffs_[k].sqrbeta*polykm1)/rCoeffs_[k+1].sqrbeta;
             polykm1 = polyk;
             polyk = p;
         }
         return polyk;
     }
 
     double ContOrthoPoly1D::powerAverage(const unsigned deg,
                                          const unsigned power) const
     {
         switch (power)
         {
         case 0:
             return 1.0;
         case 1:
             return kdelta(deg, 0U);
         case 2:
             return 1.0;
         default:
             {
                 if (deg > maxdeg_) throw std::invalid_argument(
                     "In npstat::ContOrthoPoly1D::powerAverage: "
                     "degree argument is out of range");
                 const unsigned long mSize = measure_.size();
                 precise_type sum = precise_zero;
                 for (unsigned long i = 0; i < mSize; ++i)
                     sum += measure_[i].second*powl(normpoly(deg, measure_[i].first), power);
                 return sum/wsum_;
             }
         }
     }
 
     std::pair<ContOrthoPoly1D::precise_type,ContOrthoPoly1D::precise_type>
     ContOrthoPoly1D::twonormpoly(
         const unsigned deg1, const unsigned deg2, const precise_type x) const
     {
         precise_type p1 = precise_zero, p2 = precise_zero, polyk = 1.0, polykm1 = precise_zero;
         const unsigned degmax = std::max(deg1, deg2);
         for (unsigned k=0; k<degmax; ++k)
         {
             if (k == deg1)
                 p1 = polyk;
             if (k == deg2)
                 p2 = polyk;
             const precise_type p = ((x - rCoeffs_[k].alpha)*polyk -
                                    rCoeffs_[k].sqrbeta*polykm1)/rCoeffs_[k+1].sqrbeta;
             polykm1 = polyk;
             polyk = p;
         }
         if (deg1 == degmax)
             p1 = polyk;
         if (deg2 == degmax)
             p2 = polyk;
         return std::pair<precise_type,precise_type>(p1, p2);
     }
 
     ContOrthoPoly1D::precise_type ContOrthoPoly1D::normpolyprod(
         const unsigned* degrees, const unsigned nDeg,
         const unsigned degmax, const precise_type x) const
     {
         precise_type prod = 1.0, polyk = 1.0, polykm1 = precise_zero;
         for (unsigned k=0; k<degmax; ++k)
         {
             for (unsigned i=0; i<nDeg; ++i)
                 if (k == degrees[i])
                     prod *= polyk;
             const precise_type p = ((x - rCoeffs_[k].alpha)*polyk -
                                     rCoeffs_[k].sqrbeta*polykm1)/rCoeffs_[k+1].sqrbeta;
             polykm1 = polyk;
             polyk = p;
         }
         for (unsigned i=0; i<nDeg; ++i)
             if (degmax == degrees[i])
                 prod *= polyk;
         return prod;
     }
 
     std::pair<ContOrthoPoly1D::precise_type,ContOrthoPoly1D::precise_type>
     ContOrthoPoly1D::monicInnerProducts(const unsigned degree) const
     {
         precise_type sum = precise_zero, xsum = precise_zero;
         const unsigned long mSize = measure_.size();
         for (unsigned long i = 0; i < mSize; ++i)
         {
             const precise_type x = measure_[i].first;
             const precise_type p = monicpoly(degree, x);
             const precise_type pprod = p*p*measure_[i].second;
             sum += pprod;
             xsum += x*pprod;
         }
         return std::pair<precise_type,precise_type>(xsum/wsum_, sum/wsum_);
     }
 
     void ContOrthoPoly1D::weightCoeffs(double *coeffs,
                                        const unsigned maxdeg) const
     {
         if (maxdeg > maxdeg_) throw std::invalid_argument(
             "In npstat::ContOrthoPoly1D::weightsSquaredCoeffs: "
             "maximum degree out of range");
         assert(coeffs);
 
         std::vector<precise_type> scalarProducts(maxdeg+1U, precise_zero);
 
         const unsigned long mSize = measure_.size();
         for (unsigned long i = 0; i < mSize; ++i)
         {
             const double x = measure_[i].first;
             const double w = measure_[i].second;
             const precise_type f = w;
 
             precise_type polyk = 1.0, polykm1 = precise_zero;
             for (unsigned k=0; k<maxdeg; ++k)
             {
                 scalarProducts[k] += polyk*f*w;
                 const precise_type poly = ((x - rCoeffs_[k].alpha)*polyk -
                                           rCoeffs_[k].sqrbeta*polykm1)/rCoeffs_[k+1].sqrbeta;
                 polykm1 = polyk;
                 polyk = poly;
             }
             scalarProducts[maxdeg] += polyk*f*w;
         }
 
         for (unsigned deg=0; deg<=maxdeg; ++deg)
             coeffs[deg] = scalarProducts[deg]/wsum_;
     }
 
     void ContOrthoPoly1D::calcRecurrenceStieltjes()
     {
         precise_type prevSprod = 1.0;
         for (unsigned deg=0; deg<=maxdeg_; ++deg)
         {
             const std::pair<precise_type,precise_type> ip = monicInnerProducts(deg);
             rCoeffs_.push_back(Recurrence(ip.first/ip.second, ip.second/prevSprod));
             prevSprod = ip.second;
         }
     }
 
     void ContOrthoPoly1D::calcRecurrenceLanczos()
     {
         const unsigned long mSize = measure_.size();
         std::vector<precise_type> dmem(mSize*2UL);
 
         precise_type* p0 = &dmem[0];
         for (unsigned long i=0; i<mSize; ++i)
             p0[i] = measure_[i].first;
 
         precise_type* p1 = &dmem[mSize];
         clearBuffer(p1, mSize);
         p1[0] = measure_[0].second/wsum_;
 
         const unsigned long mm1 = mSize - 1;
         for (unsigned long n=0; n<mm1; ++n)
         {
             const unsigned long np1 = n + 1UL;
             precise_type xlam = measure_[np1].first;
             precise_type pn = measure_[np1].second/wsum_;
             precise_type gam = 1.0, sig = 0.0, t = 0.0;
 
             for (unsigned long k=0; k<=np1; ++k)
             {
                 const precise_type rho = p1[k] + pn;
                 const precise_type tmp = gam*rho;
                 precise_type tsig = sig;
                 if (rho <= 0.0)
                 {
                     gam = 1.0;
                     sig = 0.0;
                 }
                 else
                 {
                     gam = p1[k]/rho;
                     sig = pn/rho;
                 }
                 const precise_type tk = sig*(p0[k] - xlam) - gam*t;
                 p0[k] -= tk - t;
                 t = tk;
                 if (sig < 0.0)
                     pn = tsig*p1[k];
                 else
                     pn = t*t/sig;
                 tsig = sig;
                 p1[k] = tmp;
             }
         }
 
         for (unsigned deg=0; deg<=maxdeg_; ++deg)
             rCoeffs_.push_back(Recurrence(p0[deg], p1[deg]));
     }
 
     double ContOrthoPoly1D::poly(const unsigned deg, const double x) const
     {
         if (deg > maxdeg_) throw std::invalid_argument(
             "In npstat::ContOrthoPoly1D::poly: degree argument is out of range");
         return normpoly(deg, x - meanX_);
     }
 
     std::pair<double,double> ContOrthoPoly1D::polyPair(
         const unsigned deg1, const unsigned deg2, const double x) const
     {
         if (deg1 > maxdeg_ || deg2 > maxdeg_) throw std::invalid_argument(
             "In npstat::ContOrthoPoly1D::polyPair: degree argument is out of range");
         const std::pair<precise_type,precise_type>& ld =
             twonormpoly(deg1, deg2, x - meanX_);
         return std::pair<double,double>(ld.first, ld.second);
     }
 
     ContOrthoPoly1D::precise_type ContOrthoPoly1D::empiricalKroneckerDelta(
         const unsigned ideg1, const unsigned ideg2) const
     {
         if (ideg1 > maxdeg_ || ideg2 > maxdeg_) throw std::invalid_argument(
             "In npstat::ContOrthoPoly1D::empiricalKroneckerDelta: "
             "degree argument is out of range");
         precise_type sum = precise_zero;
         const unsigned long mSize = measure_.size();
         for (unsigned long i = 0; i < mSize; ++i)
         {
             const std::pair<precise_type, precise_type>& p =
                 twonormpoly(ideg1, ideg2, measure_[i].first);
             sum += measure_[i].second*p.first*p.second;
         }
         return sum/wsum_;
     }
 
     Matrix<ContOrthoPoly1D::precise_type>
     ContOrthoPoly1D::empiricalKroneckerMatrix(const unsigned maxdeg) const
     {
         if (maxdeg > maxdeg_) throw std::invalid_argument(
             "In npstat::ContOrthoPoly1D::empiricalKroneckerMatrix: "
             "matrix size is out of range");
         Matrix<ContOrthoPoly1D::precise_type> kd(maxdeg+1U, maxdeg+1U, 0);
         std::vector<ContOrthoPoly1D::precise_type> polyv(maxdeg+1U);
         ContOrthoPoly1D::precise_type* polyBuf = &polyv[0];
         const unsigned long mSize = measure_.size();
         for (unsigned long ipt = 0; ipt < mSize; ++ipt)
         {
             const double weight = measure_[ipt].second;
             getAllPolys(measure_[ipt].first, maxdeg, polyBuf);
             for (unsigned i=0; i<=maxdeg; ++i)
             {
                 const ContOrthoPoly1D::precise_type polyI = polyBuf[i];
                 ContOrthoPoly1D::precise_type* row = kd[i];
                 for (unsigned j=i; j<=maxdeg; ++j)
                     row[j] += weight*polyI*polyBuf[j];
             }
         }
         for (unsigned i=0; i<maxdeg; ++i)
             for (unsigned j=i+1U; j<=maxdeg; ++j)
                 kd[j][i] = kd[i][j];
-        return kd/wsum_;
+        kd /= wsum_;
+        return kd;
     }
 
     double ContOrthoPoly1D::jointPowerAverage(
         const unsigned ideg1, const unsigned power1,
         const unsigned ideg2, const unsigned power2) const
     {
         // Process various simple special cases first
         if (!power1)
             return powerAverage(ideg2, power2);
         if (!power2)
             return powerAverage(ideg1, power1);
         if (power1 == 1U && power2 == 1U)
             return kdelta(ideg1, ideg2);
 
         // General calculation
         if (ideg1 > maxdeg_ || ideg2 > maxdeg_) throw std::invalid_argument(
             "In npstat::ContOrthoPoly1D::jointPowerAverage: "
             "degree argument is out of range");
         precise_type sum = precise_zero;
         const unsigned long mSize = measure_.size();
         for (unsigned long i = 0; i < mSize; ++i)
         {
             const std::pair<precise_type, precise_type>& p =
                 twonormpoly(ideg1, ideg2, measure_[i].first);
             sum += measure_[i].second*powl(p.first,power1)*powl(p.second,power2);
         }
         return sum/wsum_;
     }
 
     double ContOrthoPoly1D::jointAverage(
         const unsigned* degrees, const unsigned nDegrees,
         const bool sorted) const
     {
         if (nDegrees)
         {
             assert(degrees);
 
             // See if we can avoid a direct calculation
             if (nDegrees == 1U)
                 return kdelta(degrees[0], 0U);
             if (nDegrees == 2U)
                 return kdelta(degrees[0], degrees[1]);
 
             // Check if we can remove leading zeros
             {
                 unsigned nonZeroDegs = nDegrees;
                 while (nonZeroDegs && !degrees[0])
                 {
                     ++degrees;
                     --nonZeroDegs;
                 }
                 if (nonZeroDegs < nDegrees)
                     return jointAverage(degrees, nonZeroDegs, sorted);
             }
 
             // Check if we can remove any other zeros
             if (!sorted)
             {
                 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(degBuf, nonZeroDegs, sorted);
                     }
                 }
             }
 
             unsigned degmax;
             if (sorted)
                 degmax = degrees[nDegrees-1U];
             else
                 degmax = *std::max_element(degrees, degrees+nDegrees);
             if (degmax > maxdeg_) throw std::invalid_argument(
                 "In npstat::ContOrthoPoly1D::jointAverage: "
                 "degree argument is out of range");
 
             precise_type sum = precise_zero;
             const unsigned long mSize = measure_.size();
             for (unsigned long i = 0; i < mSize; ++i)
                 sum += measure_[i].second*normpolyprod(
                     degrees, nDegrees, degmax, measure_[i].first);
             return sum/wsum_;
         }
         else
             return 1.0;
     }
 
     double ContOrthoPoly1D::empiricalKroneckerCovariance(
         const unsigned deg1, const unsigned deg2,
         const unsigned deg3, const unsigned deg4) const
     {
         double cov = 0.0;
         if (!((deg1 == 0 && deg2 == 0) || (deg3 == 0 && deg4 == 0)))
         {
             unsigned degs[4];
             degs[0] = deg1;
             degs[1] = deg2;
             degs[2] = deg3;
             degs[3] = deg4;
             cov = (jointAverage(degs, 4) - kdelta(deg1, deg2)*kdelta(deg3, deg4))/
                 effectiveSampleSize();
             if (std::min(deg1, deg2) == std::min(deg3, deg4) &&
                 std::max(deg1, deg2) == std::max(deg3, deg4))
                 if (cov < 0.0)
                     cov = 0.0;
         }
         return cov;
     }
 
     std::pair<double,double>
     ContOrthoPoly1D::recurrenceCoeffs(const unsigned deg) const
     {
         if (deg > maxdeg_) throw std::invalid_argument(
             "In npstat::ContOrthoPoly1D::recurrenceCoeffs: "
             "degree argument is out of range");
         const Recurrence& rc(rCoeffs_[deg]);
         return std::pair<double,double>(rc.alpha, rc.beta);
     }
 
     Matrix<ContOrthoPoly1D::precise_type>
     ContOrthoPoly1D::jacobiMatrix(const unsigned n) const
     {
         if (!n) throw std::invalid_argument(
             "In npstat::ContOrthoPoly1D::jacobiMatrix: "
             "can not build matrix of zero size");
         if (n > maxdeg_) throw std::invalid_argument(
             "In npstat::ContOrthoPoly1D::jacobiMatrix: "
             "matrix size is out of range");
         Matrix<precise_type> mat(n, n, 0);
         unsigned ip1 = 1;
         for (unsigned i=0; i<n; ++i, ++ip1)
         {
             mat[i][i] = rCoeffs_[i].alpha;
             if (i)
                 mat[i][i-1] = rCoeffs_[i].sqrbeta;
             if (ip1 < n)
                 mat[i][ip1] = rCoeffs_[ip1].sqrbeta;
         }
         return mat;
     }
 
     ContOrthoPoly1D::precise_type ContOrthoPoly1D::integratePoly(
         const unsigned deg1, const unsigned power,
         const double xmin, const double xmax) const
     {
         if (!deg1)
             return xmax - xmin;
         if (deg1 > maxdeg_) throw std::invalid_argument(
             "In npstat::ContOrthoPoly1D::integratePoly: "
             "degree argument is out of range");
         const unsigned nGood = PreciseQuadrature::minimalExactRule(deg1*power);
         if (!nGood) throw std::invalid_argument(
             "In npstat::ContOrthoPoly1D::integratePoly: "
             "product of poly degree and power is too large");
         PreciseQuadrature glq(nGood);
         PolyFcn fcn(*this, deg1, power);
         return glq.integrate(fcn, xmin - meanX_, xmax - meanX_);
     }
 
     ContOrthoPoly1D::precise_type ContOrthoPoly1D::integrateTripleProduct(
         const unsigned deg1, const unsigned deg2,
         const unsigned deg3, const double xmin,
         const double xmax) const
     {
         const unsigned maxdex = std::max(std::max(deg1, deg2), deg3);
         if (maxdex > maxdeg_) throw std::invalid_argument(
             "In npstat::ContOrthoPoly1D::integrateTripleProduct: "
             "degree argument is out of range");
         const unsigned sumdeg = deg1 + deg2 + deg3;
         const unsigned nGood = PreciseQuadrature::minimalExactRule(sumdeg);
         if (!nGood) throw std::invalid_argument(
             "In npstat::ContOrthoPoly1D::integrateTripleProduct: "
             "sum of the polynomial degrees is too large");
         PreciseQuadrature glq(nGood);
         TripleProdFcn fcn(*this, deg1, deg2, deg3);
         return glq.integrate(fcn, xmin - meanX_, xmax - meanX_);
     }
 
     Matrix<ContOrthoPoly1D::precise_type> ContOrthoPoly1D::tripleProductMatrix(
         const std::vector<std::pair<unsigned,unsigned> >& pairs,
         const unsigned maxdeg, const double xmin, const double xmax) const
     {
         typedef npstat::Triple<unsigned,unsigned,unsigned> Key3;
         typedef std::map<Key3,precise_type> TripleProductMap;
 
         const unsigned nPairs = pairs.size();
         if (!nPairs) throw std::invalid_argument(
             "In npstat::ContOrthoPoly1D::tripleProductMatrix: "
             "empty vector of degree specifications");
         const unsigned nPolys = maxdeg + 1U;
 
         npstat::Matrix<precise_type> design0(nPairs, nPolys);
         TripleProductMap tmap;
 
         // Avoid calculating the same integral more than once
         for (unsigned ipair=0; ipair<nPairs; ++ipair)
         {
             const unsigned n = pairs[ipair].first;
             const unsigned m = pairs[ipair].second;
             for (unsigned ipoly=0; ipoly<nPolys; ++ipoly)
             {
                 unsigned degs[3];
                 degs[0] = m;
                 degs[1] = n;
                 degs[2] = ipoly;
                 std::sort(degs, degs+3);
                 const Key3 key(degs[0], degs[1], degs[2]);
                 TripleProductMap::iterator it = tmap.find(key);
                 if (it == tmap.end())
                 {
                     const precise_type value = integrateTripleProduct(
                         n, m, ipoly, xmin, xmax);
                     std::pair<TripleProductMap::iterator,bool> where =
                         tmap.insert(std::make_pair(key, value));
                     assert(where.second);
                     it = where.first;
                 }
                 design0[ipair][ipoly] = it->second;
             }
         }
 
         return design0;
     }
 
     double ContOrthoPoly1D::cachedJointAverage(const unsigned deg1,
                                                const unsigned deg2,
                                                const unsigned deg3,
                                                const unsigned deg4) const
     {
         MemoKey key(deg1, deg2, deg3, deg4);
         return getCachedAverage(key);
     }
 
     double ContOrthoPoly1D::cachedJointAverage(const unsigned deg1,
                                                const unsigned deg2,
                                                const unsigned deg3,
                                                const unsigned deg4,
                                                const unsigned deg5,
                                                const unsigned deg6) const
     {
         MemoKey key(deg1, deg2, deg3, deg4, deg5, deg6);
         return getCachedAverage(key);
     }
 
     double ContOrthoPoly1D::cachedJointAverage(const unsigned deg1,
                                                const unsigned deg2,
                                                const unsigned deg3,
                                                const unsigned deg4,
                                                const unsigned deg5,
                                                const unsigned deg6,
                                                const unsigned deg7,
                                                const unsigned deg8) const
     {
         MemoKey key(deg1, deg2, deg3, deg4, deg5, deg6, deg7, deg8);
         return getCachedAverage(key);
     }
 
     double ContOrthoPoly1D::getCachedAverage(const MemoKey& key) const
     {
         double value;
         std::map<MemoKey,double>::const_iterator it = cachedAverages_.find(key);
         if (it == cachedAverages_.end())
         {
             value = jointAverage(key.degrees(), key.nDegrees(), true);
             cachedAverages_.insert(std::pair<MemoKey,double>(key, value));
         }
         else
             value = it->second;
         return value;
     }
 
     double ContOrthoPoly1D::cov4(const unsigned deg1, const unsigned deg2,
                                  const unsigned deg3, const unsigned deg4) const
     {
         const bool has0 = (deg1 == 0 && deg2 == 0) || (deg3 == 0 && deg4 == 0);
         double cov = 0.0;
         if (!has0)
         {
             const double N = effectiveSampleSize();
             cov = (cachedJointAverage(deg1, deg2, deg3, deg4) -
                    kdelta(deg1, deg2)*kdelta(deg3, deg4))/N;
             if (std::min(deg1, deg2) == std::min(deg3, deg4) &&
                 std::max(deg1, deg2) == std::max(deg3, deg4))
                 if (cov < 0.0)
                     cov = 0.0;
         }
         return cov;
     }
 
     double ContOrthoPoly1D::cov6(const unsigned a, const unsigned b,
                                  const unsigned c, const unsigned d,
                                  const unsigned e, const unsigned f) const
     {
         const bool has0 = (a == 0 && b == 0) || (c == 0 && d == 0) || (e == 0 && f == 0);
         double cov = 0.0;
         if (!has0)
         {
             double sum = cachedJointAverage(a, b, c, d, e, f);
             double add = 2.0;
 
             if (kdelta(a, b))
                 sum -= cachedJointAverage(c, d, e, f);
             else
                 add = 0.0;
             if (kdelta(c, d))
                 sum -= cachedJointAverage(a, b, e, f);
             else
                 add = 0.0;
             if (kdelta(e, f))
                 sum -= cachedJointAverage(a, b, c, d);
             else
                 add = 0.0;
 
             const double N = effectiveSampleSize();
             cov = (sum + add)/N/N;
         }
         return cov;
     }
 
     double ContOrthoPoly1D::slowCov8(const unsigned a, const unsigned b,
                                      const unsigned c, const unsigned d,
                                      const unsigned e, const unsigned f,
                                      const unsigned g, const unsigned h) const
     {
         const bool has0 = (a == 0 && b == 0) || (c == 0 && d == 0) ||
                           (e == 0 && f == 0) || (g == 0 && h == 0);
         double cov = 0.0;
         if (!has0)
         {
             const double pabcd = cachedJointAverage(a, b, c, d);
             const double pefgh = cachedJointAverage(e, f, g, h);
             const double pabef = cachedJointAverage(a, b, e, f);
             const double pcdgh = cachedJointAverage(c, d, g, h);
             const double pabgh = cachedJointAverage(a, b, g, h);
             const double pcdef = cachedJointAverage(c, d, e, f);
             const double pabcdef = cachedJointAverage(a, b, c, d, e, f);
             const double pabcdgh = cachedJointAverage(a, b, c, d, g, h);
             const double pabefgh = cachedJointAverage(a, b, e, f, g, h);
             const double pcdefgh = cachedJointAverage(c, d, e, f, g, h);
             const double pabcdefgh = cachedJointAverage(a, b, c, d, e, f, g, h);
 
             const double deltaprod = kdelta(a, b)*kdelta(c, d)*kdelta(e, f)*kdelta(g, h);
 
             const double tmp = kdelta(e,f)*kdelta(g,h)*pabcd +
                                kdelta(c,d)*kdelta(g,h)*pabef +
                                kdelta(c,d)*kdelta(e,f)*pabgh +
                                kdelta(a,b)*kdelta(c,d)*pefgh +
                                kdelta(a,b)*kdelta(e,f)*pcdgh +
                                kdelta(a,b)*kdelta(g,h)*pcdef;
 
             const double term1 = pabcd*pefgh + pabef*pcdgh +
                 pabgh*pcdef + 3.0*deltaprod - tmp;
 
             const double term2 = pabcdefgh - 6.0*deltaprod -
                 kdelta(a,b)*pcdefgh - kdelta(c,d)*pabefgh -
                 kdelta(e,f)*pabcdgh - kdelta(g,h)*pabcdef -
                 pabcd*pefgh - pabef*pcdgh - pabgh*pcdef + 2.0*tmp;
 
             const double nPoints = effectiveSampleSize();
             const double prod8 = (term1 + term2/nPoints)/nPoints/nPoints;
             cov = prod8 - cov4(a, b, c, d)*cov4(e, f, g, h);
         }
         return cov;
     }
 
     double ContOrthoPoly1D::covCov4(const unsigned a, const unsigned b,
                                     const unsigned c, const unsigned d,
                                     const unsigned e, const unsigned f,
                                     const unsigned g, const unsigned h) const
     {
         const bool has0 = (a == 0 && b == 0) || (c == 0 && d == 0) ||
                           (e == 0 && f == 0) || (g == 0 && h == 0);
         double cov = 0.0;
         if (!has0)
         {
             const double pabcdefgh = cachedJointAverage(a, b, c, d, e, f, g, h);
             const double pabcd = cachedJointAverage(a, b, c, d);
             const double pefgh = cachedJointAverage(e, f, g, h);
             const double N = effectiveSampleSize();
             cov = (pabcdefgh - pabcd*pefgh)/N/N/N;
         }
         return cov;
     }
 
     double ContOrthoPoly1D::cov8(const unsigned a, const unsigned b,
                                  const unsigned c, const unsigned d,
                                  const unsigned e, const unsigned f,
                                  const unsigned g, const unsigned h) const
     {
         const bool has0 = (a == 0 && b == 0) || (c == 0 && d == 0) ||
                           (e == 0 && f == 0) || (g == 0 && h == 0);
         double cov = 0.0;
         if (!has0)
         {
             const bool includeCubicPart = true;
 
             // First, calculate the O(N^-2) part
             const double pabef = cachedJointAverage(a, b, e, f);
             const double pcdgh = cachedJointAverage(c, d, g, h);
             const double pabgh = cachedJointAverage(a, b, g, h);
             const double pcdef = cachedJointAverage(c, d, e, f);
 
             const double deltaprod = kdelta(a, b)*kdelta(c, d)*kdelta(e, f)*kdelta(g, h);
             const double tmp2 = kdelta(c, d)*kdelta(g, h)*pabef +
                                 kdelta(c, d)*kdelta(e, f)*pabgh +
                                 kdelta(a, b)*kdelta(e, f)*pcdgh +
                                 kdelta(a, b)*kdelta(g, h)*pcdef;
 
             const double term2 = pabef*pcdgh + pabgh*pcdef + 2.0*deltaprod - tmp2;
             double term3 = 0.0;
 
             if (includeCubicPart)
             {
                 // Additional terms needed to calculate the O(N^-3) part
                 const double pabcdefgh = cachedJointAverage(a, b, c, d, e, f, g, h);
 
                 double sixsum = 0.0;
                 if (kdelta(a, b))
                     sixsum += cachedJointAverage(c, d, e, f, g, h);
                 if (kdelta(c, d))
                     sixsum += cachedJointAverage(a, b, e, f, g, h);
                 if (kdelta(e, f))
                     sixsum += cachedJointAverage(a, b, c, d, g, h);
                 if (kdelta(g, h))
                     sixsum += cachedJointAverage(a, b, c, d, e, f);
 
                 const double pabcd = cachedJointAverage(a, b, c, d);
                 const double pefgh = cachedJointAverage(e, f, g, h);
                 const double tmp3 = tmp2 + kdelta(e, f)*kdelta(g, h)*pabcd +
                                            kdelta(a, b)*kdelta(c, d)*pefgh;
 
                 term3 = pabcdefgh - 6.0*deltaprod - sixsum  -
                     (pabcd*pefgh + pabef*pcdgh + pabgh*pcdef) + 2.0*tmp3;
             }
 
             const double N = effectiveSampleSize();
             cov = (term2 + term3/N)/N/N;
 
             // const bool isVariance = ?;
             // if (isVariance)
             //     if (cov < 0.0)
             //         cov = 0.0;
         }
         return cov;
     }
 
     double ContOrthoPoly1D::epsExpectation(const unsigned m_in,
                                            const unsigned n_in,
                                            const bool highOrder) const
     {
         if (highOrder)
         {
             precise_type sum = precise_zero;
             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*cachedJointAverage(k, k, m, n);
                 }
                 if (m == n)
                     sum -= 1.0;
                 else
                     sum -= (cachedJointAverage(m, m, m, n) +
                             cachedJointAverage(m, n, n, n))/2.0;
             }
             return sum/effectiveSampleSize();
         }
         else
             return 0.0;
     }
 
     double ContOrthoPoly1D::epsCovariance(const unsigned m1_in,
                                           const unsigned n1_in,
                                           const unsigned m2_in,
                                           const unsigned n2_in,
                                           const bool highOrder) const
     {
         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);
 
             precise_type sum = precise_zero;
 
             // Process the -v_{m1,n1} term (i.e., the linear one) of eps_{m1,n1}
             sum += cov4(m2, n2, m1, n1);
             sum += cov6(m2, n2, n2, n2, m1, n1)/2.0;
             sum += cov6(m2, n2, m2, m2, m1, n1)/2.0;
             for (unsigned k=0; k<=n2; ++k)
             {
                 const double factor = k > m2 ? 1.0 : 2.0;
                 sum -= factor*cov6(k, m2, k, n2, m1, n1);
             }
 
             // Process the term -v_{m1,n1}/2 (v_{n1,n1} + v_{m1,m1}) of eps_{m1,n1}
             unsigned idx[2];
             idx[0] = n1;
             idx[1] = m1;
             for (unsigned ii=0; ii<2; ++ii)
             {
                 const unsigned mOrN = idx[ii];
                 sum += cov6(m1, n1, mOrN, mOrN, m2, n2)/2.0;
                 sum += cov8(m1, n1, mOrN, mOrN, m2, n2, n2, n2)/4.0;
                 sum += cov8(m1, n1, mOrN, mOrN, m2, n2, m2, m2)/4.0;
                 for (unsigned k=0; k<=n2; ++k)
                 {
                     const double factor = k > m2 ? 1.0 : 2.0;
                     sum -= factor/2.0*cov8(m1, n1, mOrN, mOrN, k, m2, k, n2);
                 }
             }
 
             // Process the sum in eps_{m1,n1}
             for (unsigned k1=0; k1<=n1; ++k1)
             {
                 const double f1 = k1 > m1 ? 1.0 : 2.0;
                 sum -= f1*cov6(k1, m1, k1, n1, m2, n2);
                 sum -= f1*cov8(k1, m1, k1, n1, m2, n2, n2, n2)/2.0;
                 sum -= f1*cov8(k1, m1, k1, n1, m2, n2, m2, m2)/2.0;
                 for (unsigned k=0; k<=n2; ++k)
                 {
                     const double factor = k > m2 ? 1.0 : 2.0;
                     sum += f1*factor*cov8(k1, m1, k1, n1, k, m2, k, n2);
                 }
             }
 
             return sum;
         }
         else
             return cov4(m2_in, n2_in, m1_in, n1_in);
     }
 
     Matrix<double> ContOrthoPoly1D::epsCovarianceMatrix(
         const std::vector<std::pair<unsigned,unsigned> >& pairs,
         const bool highOrder) const
     {
         if (pairs.empty()) throw std::invalid_argument(
             "In npstat::ContOrthoPoly1D::epsCovarianceMatrix: empty list of pairs");
 
         const unsigned nPairs = pairs.size();
         Matrix<double> cov(nPairs, nPairs);
         for (unsigned ipair1=0; ipair1<nPairs; ++ipair1)
         {
             const unsigned n = pairs[ipair1].first;
             const unsigned m = pairs[ipair1].second;
             for (unsigned ipair2=0; ipair2<=ipair1; ++ipair2)
             {
                 const double value = epsCovariance(
                     n, m, pairs[ipair2].first, pairs[ipair2].second, highOrder);
                 cov[ipair1][ipair2] = value;
                 if (ipair1 != ipair2)
                     cov[ipair2][ipair1] = value;
             }
         }
         return cov;
     }
 
     CPP11_auto_ptr<StorablePolySeries1D> ContOrthoPoly1D::makeStorablePolySeries(
         const double i_xmin, const double i_xmax,
         const double *coeffs, const unsigned maxdeg) const
     {
         const unsigned sz = rCoeffs_.size();
         std::vector<std::pair<precise_type,precise_type> > rc(sz);
         for (unsigned i=0; i<sz; ++i)
         {
             const Recurrence& r(rCoeffs_[i]);
             rc[i].first = r.alpha;
             rc[i].second = r.sqrbeta;
         }
         return CPP11_auto_ptr<StorablePolySeries1D>(
             new StorablePolySeries1D(rc, i_xmin, i_xmax, meanX_, coeffs, maxdeg));
     }
 }