Index: trunk/npstat/stat/OrthoPolyGOFTest.icc
===================================================================
--- trunk/npstat/stat/OrthoPolyGOFTest.icc	(revision 732)
+++ trunk/npstat/stat/OrthoPolyGOFTest.icc	(revision 733)
@@ -1,144 +1,140 @@
 #include <cassert>
 #include <algorithm>
 #include <stdexcept>
 #include <cmath>
 
 #include "npstat/nm/truncatedInverseSqrt.hh"
 #include "npstat/nm/sumOfSquares.hh"
 
 namespace npstat {
     template <class Quadrature>
     OrthoPolyGOFTest::OrthoPolyGOFTest(
         const Quadrature& q,
         const AbsClassicalOrthoPoly1D& poly,
         const unsigned i_maxdeg,
         const double* importanceWeights,
         const unsigned nWeights)
         : poly_(poly.clone()),
           importanceWeights_(i_maxdeg+1U, 1.0),
-          buf_((i_maxdeg+1U)*(i_maxdeg+2U)/2U),
+          buf_((i_maxdeg+1U)*(i_maxdeg+2U)/2U - 1U),
           statBuf_(2U*i_maxdeg),
           maxdeg_(i_maxdeg)
     {
         if (nWeights)
         {
             if (nWeights != i_maxdeg + 1U) throw std::invalid_argument(
                 "In npstat::OrthoPolyGOFTest constructor: "
                 "inconsistent number of importance weights");
             assert(importanceWeights);
             for (unsigned i=0; i<nWeights; ++i)
             {
                 if (importanceWeights[i] <= 0.0) throw std::invalid_argument(
                     "In npstat::OrthoPolyGOFTest constructor: "
                     "all importance weights must be positive");
                 importanceWeights_[i] = importanceWeights[i];
             }
         }
         calculateNormalizingMatrix(q, &invsqr_);
     }
 
     template <class Quadrature>
     void OrthoPolyGOFTest::calculateNormalizingMatrix(
         const Quadrature& q, Matrix<lapack_double>* result) const
     {
         assert(result);
 
         const unsigned maxdegP1 = maxdeg_ + 1U;
-        const unsigned nRows = maxdegP1*(maxdegP1+1U)/2U;
+        const unsigned nRows = maxdegP1*(maxdegP1+1U)/2U - 1U;
         Matrix<lapack_double> covmat(nRows, nRows);
 
         unsigned iRow = 0;
         for (unsigned i=0; i<=maxdeg_; ++i)
         {
             const long double prod_i = importanceWeights_[i];
-            for (unsigned k=i; k<=maxdeg_; ++k, ++iRow)
+            for (unsigned k=i ? i : 1U; k<=maxdeg_; ++k, ++iRow)
             {
                 assert(iRow < nRows);
                 const double k_del_ik = k == i ? 1.0 : 0.0;
                 const long double prod_k = prod_i*importanceWeights_[k];
 
                 unsigned iCol = 0;
                 for (unsigned j=0; j<=maxdeg_; ++j)
                 {
                     const long double prod_j = prod_k*importanceWeights_[j];
-                    for (unsigned m=j; m<=maxdeg_; ++m, ++iCol)
+                    for (unsigned m=j ? j : 1U; m<=maxdeg_; ++m, ++iCol)
                     {
                         if (iCol >= iRow)
                         {
                             assert(iCol < nRows);
                             const long double prod_m = prod_j*importanceWeights_[m];
-                            double cov = 0.0;
-                            if (!((i == 0 && k == 0) || (j == 0 && m == 0)))
-                            {
-                                const double k_del_jm = j == m ? 1.0 : 0.0;
-                                cov = poly_->jointAverage(q, i, k, j, m) -
-                                    k_del_ik*k_del_jm;
-                                if (i == j && k == m && cov < 0.0)
-                                    cov = 0.0;
-                            }
+                            const double k_del_jm = j == m ? 1.0 : 0.0;
+                            double cov = poly_->jointAverage(q, i, k, j, m) -
+                                k_del_ik*k_del_jm;
+                            if (i == j && k == m && cov < 0.0)
+                                cov = 0.0;
                             covmat[iRow][iCol] = cov*prod_m;
                         }
                     }
                 }
                 assert(iCol == nRows);
             }
         }
         assert(iRow == nRows);
 
         iRow = 0;
         for (unsigned i=0; i<=maxdeg_; ++i)
         {
-            for (unsigned k=i; k<=maxdeg_; ++k, ++iRow)
+            for (unsigned k=i ? i : 1U; k<=maxdeg_; ++k, ++iRow)
             {
                 unsigned iCol = 0;
                 for (unsigned j=0; j<=maxdeg_; ++j)
                 {
-                    for (unsigned m=j; m<=maxdeg_; ++m, ++iCol)
+                    for (unsigned m=j ? j : 1U; m<=maxdeg_; ++m, ++iCol)
                         if (iCol < iRow)
                             covmat[iRow][iCol] = covmat[iCol][iRow];
                 }       
             }
         }
 
         truncatedInverseSqrt(covmat, 2U*maxdeg_, 0.0, result);
         assert(result->nRows() == 2U*maxdeg_);
     }
 
     template<typename Numeric>
     void OrthoPolyGOFTest::normalizedDeviations(
         const Numeric* data, const unsigned long lenData,
         double* deviations, const unsigned lenDeviations) const
     {
         const unsigned twomaxdeg = 2U*maxdeg_;
 
         assert(lenData);
         assert(data);
         assert(lenDeviations >= twomaxdeg);
         assert(deviations);
 
         const Matrix<double>& unnorm = poly_-> sampleProductAverages(
             data, lenData, maxdeg_);
         double* buf = const_cast<double*>(&buf_[0]);
 
         unsigned iRow = 0;
         for (unsigned i=0; i<=maxdeg_; ++i)
-            for (unsigned k=i; k<=maxdeg_; ++k, ++iRow)
+            for (unsigned k=i ? i : 1U; k<=maxdeg_; ++k, ++iRow)
                 buf[iRow] = unnorm[i][k]*importanceWeights_[i]*importanceWeights_[k];
 
         assert(iRow == buf_.size());
         invsqr_.timesVector(buf, iRow, deviations, twomaxdeg);
 
         const double sqrlen = std::sqrt(static_cast<double>(lenData));
         for (unsigned i=0; i<twomaxdeg; ++i)
             deviations[i] *= sqrlen;
     }
 
     template<typename Numeric>
     double OrthoPolyGOFTest::testStatistic(
         const Numeric* data, const unsigned long lenData) const
     {
         double* statbuf = const_cast<double*>(&statBuf_[0]);
         normalizedDeviations(data, lenData, statbuf, statBuf_.size());
         return sumOfSquares(statbuf, statBuf_.size());
     }
 }
Index: trunk/npstat/stat/PolynomialDistro1D.cc
===================================================================
--- trunk/npstat/stat/PolynomialDistro1D.cc	(revision 732)
+++ trunk/npstat/stat/PolynomialDistro1D.cc	(revision 733)
@@ -1,143 +1,143 @@
 #include <cassert>
 #include <stdexcept>
 #include <sstream>
 #include <climits>
 
 #include "npstat/nm/findRootUsingBisections.hh"
 
 #include "npstat/stat/PolynomialDistro1D.hh"
 #include "npstat/stat/distributionReadError.hh"
 
 namespace npstat {
     PolynomialDistro1D::PolynomialDistro1D(
         const double location, const double scale,
         const double* coeffs, const unsigned maxdeg,
         const unsigned nCheck)
         : AbsScalableDistribution1D(location, scale),
           glq_(GaussLegendreQuadrature::minimalExactRule(maxdeg))
     {
         setupCoeffs(coeffs, maxdeg);
         checkPositive(nCheck);
     }
 
     PolynomialDistro1D::PolynomialDistro1D(
         const double location, const double scale,
         const std::vector<double>& coeffs,
         const unsigned nCheck)
         : AbsScalableDistribution1D(location, scale),
           glq_(GaussLegendreQuadrature::minimalExactRule(coeffs.size()))
     {
         setupCoeffs(coeffs.empty() ? (const double*)0 : &coeffs[0],
                     coeffs.size());
         checkPositive(nCheck);
     }
 
     void PolynomialDistro1D::setupCoeffs(const double* coeffs,
                                          const unsigned maxdeg)
     {
         if (maxdeg)
             assert(coeffs);
         allCoeffs_.resize(maxdeg + 1U);
         allCoeffs_[0] = 1.0;
         for (unsigned i=0; i<maxdeg; ++i)
             allCoeffs_[i+1U] = coeffs[i];
     }
 
     void PolynomialDistro1D::checkPositive(const unsigned nCheck) const
     {
         if (nCheck < 2U) throw std::invalid_argument(
             "In npstat::PolynomialDistro1D::checkPositive: "
             "insufficient number of density check points");
         const unsigned nM1 = nCheck - 1U;
         const double step = 1.0/nM1;
         for (unsigned i=0; i<nCheck; ++i)
         {
             double x = step*i;
             if (i == nM1)
                 x = 1.0;
             if (unscaledDensity(x) <= 0.0)
             {
                 std::ostringstream os;
                 os.precision(17);
                 os << "In npstat::PolynomialDistro1D::checkPositive: "
-                    << "density is not positive at x = " << x;
+                   << "density is not positive at x = " << x;
                 throw std::invalid_argument(os.str());
             }
         }
     }
 
     double PolynomialDistro1D::unscaledDensity(const double x) const
     {
         if (x < 0.0 || x > 1.0)
             return 0.0;
         else
             return poly_.series(&allCoeffs_[0], allCoeffs_.size()-1, x);
     }
 
     double PolynomialDistro1D::unscaledCdf(const double x) const
     {
         if (x <= 0.0)
             return 0.0;
         else if (x >= 1.0)
             return 1.0;
         else
             return glq_.integrate(UnscaledDensityFunctor1D(*this), 0.0, x);
     }
 
     double PolynomialDistro1D::unscaledQuantile(const double r1) const
     {
         if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error(
             "In npstat::PolynomialDistro1D::unscaledQuantile: "
             "cdf argument outside of [0, 1] interval");
         if (r1 == 0.0)
             return 0.0;
         if (r1 == 1.0)
             return 1.0;
         double q;
         const bool status = findRootUsingBisections(
             UnscaledCdfFunctor1D(*this), r1, 0.0, 1.0, 4.0*DBL_EPSILON, &q);
         assert(status);
         return q;
     }
 
     double PolynomialDistro1D::unscaledExceedance(const double x) const
     {
         if (x <= 0.0)
             return 1.0;
         else if (x >= 1.0)
             return 0.0;
         else
             return glq_.integrate(UnscaledDensityFunctor1D(*this), x, 1.0);
     }
 
     bool PolynomialDistro1D::isEqual(const AbsDistribution1D& other) const
     {
         const PolynomialDistro1D& r = static_cast<const PolynomialDistro1D&>(other);
         return AbsScalableDistribution1D::isEqual(r) && allCoeffs_ == r.allCoeffs_;
     }
 
     bool PolynomialDistro1D::write(std::ostream& os) const
     {
         AbsScalableDistribution1D::write(os);
         gs::write_pod_vector(os, allCoeffs_);
         return !os.fail();
     }
 
     PolynomialDistro1D* PolynomialDistro1D::read(const gs::ClassId& id, std::istream& in)
     {
         static const gs::ClassId current(gs::ClassId::makeId<PolynomialDistro1D>());
         current.ensureSameId(id);
 
         double location, scale;
         if (AbsScalableDistribution1D::read(in, &location, &scale))
         {
             std::vector<double> table;
             gs::read_pod_vector(in, &table);
             if (!in.fail() && table.size())
                 return new PolynomialDistro1D(location, scale,
                                               &table[1], table.size()-1, 2U);
         }
         distributionReadError(in, classname());
         return 0;
     }
 }