Index: trunk/npstat/stat/DensityOrthoPoly1D.cc
===================================================================
--- trunk/npstat/stat/DensityOrthoPoly1D.cc	(revision 699)
+++ trunk/npstat/stat/DensityOrthoPoly1D.cc	(revision 700)
@@ -1,172 +1,175 @@
 #include <stdexcept>
 
 #include "npstat/nm/FejerQuadrature.hh"
 #include "npstat/stat/DensityOrthoPoly1D.hh"
 
 namespace npstat {
     DensityOrthoPoly1D::DensityOrthoPoly1D(const AbsDistribution1D& distro,
                                            const unsigned maxDegree,
                                            const unsigned nIntegPoints,
                                            const OrthoPolyMethod m)
         : distro_(distro.clone()), wsum_(0.0L), maxdeg_(maxDegree)
     {
         if (!nIntegPoints) throw std::invalid_argument(
             "In npstat::DensityOrthoPoly1D constructor: the number of "
             "integration points must be positive");
         const FejerQuadrature quad(nIntegPoints);
         const DensityFunctor1D weight(distro);
         const std::vector<MeasurePt>& measure = quad.weightedIntegrationPoints(
             weight, distro_->quantile(0.0), distro_->quantile(1.0));
         assert(measure.size() == nIntegPoints);
         for (unsigned i=0; i<nIntegPoints; ++i)
             wsum_ += measure[i].second;
+        if (wsum_ <= 0.0L) throw std::invalid_argument(
+            "In npstat::DensityOrthoPoly1D constructor: total weight "
+            "is not positive");
         calcRecurrenceCoeffs(measure, m);
     }
 
     void DensityOrthoPoly1D::calcRecurrenceCoeffs(
         const std::vector<MeasurePt>& measure,
         const OrthoPolyMethod m)
     {
         rCoeffs_.clear();
         rCoeffs_.reserve(maxdeg_ + 1);
         switch (m)
         {
         case OPOLY_STIELTJES:
             calcRecurrenceStieltjes(measure);
             break;
 
         case OPOLY_LANCZOS:
             calcRecurrenceLanczos(measure);
             break;
 
         default:
             throw std::runtime_error(
                 "In npstat::DensityOrthoPoly1D::calcRecurrenceCoeffs: "
                 "incomplete switch statement. This is a bug. Please report.");
         }
         assert(rCoeffs_.size() == maxdeg_ + 1);
     }
 
     void DensityOrthoPoly1D::calcRecurrenceStieltjes(const std::vector<MeasurePt>& measure)
     {
         long double prevSprod = 1.0L;
         for (unsigned deg=0; deg<=maxdeg_; ++deg)
         {
             const std::pair<long double,long double> ip = monicInnerProducts(measure, deg);
             rCoeffs_.push_back(Recurrence(ip.first/ip.second, ip.second/prevSprod));
             prevSprod = ip.second;
         }
     }
 
     void DensityOrthoPoly1D::calcRecurrenceLanczos(const std::vector<MeasurePt>& measure)
     {
         typedef long double Real;
 
         const unsigned mSize = measure.size();
         std::vector<Real> dmem(mSize*2U);
 
         Real* p0 = &dmem[0];
         for (unsigned i=0; i<mSize; ++i)
             p0[i] = measure[i].first;
 
         Real* p1 = &dmem[mSize];
         clearBuffer(p1, mSize);
         p1[0] = measure[0].second/wsum_;
 
         const unsigned mm1 = mSize - 1;
         for (unsigned n=0; n<mm1; ++n)
         {
             const unsigned np1 = n + 1U;
             Real xlam = measure[np1].first;
             Real pn = measure[np1].second/wsum_;
             Real gam = 1.0, sig = 0.0, t = 0.0;
 
             for (unsigned k=0; k<=np1; ++k)
             {
                 const Real rho = p1[k] + pn;
                 const Real tmp = gam*rho;
                 Real tsig = sig;
                 if (rho <= 0.0)
                 {
                     gam = 1.0;
                     sig = 0.0;
                 }
                 else
                 {
                     gam = p1[k]/rho;
                     sig = pn/rho;
                 }
                 const Real 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]));
     }
 
     std::pair<long double,long double>
     DensityOrthoPoly1D::monicInnerProducts(const std::vector<MeasurePt>& measure,
                                            const unsigned degree) const
     {
         long double sum = 0.0L, xsum = 0.0L;
         const unsigned mSize = measure.size();
         for (unsigned i = 0; i < mSize; ++i)
         {
             const long double x = measure[i].first;
             const long double p = monicpoly(degree, x);
             const long double pprod = p*p*measure[i].second;
             sum += pprod;
             xsum += x*pprod;
         }
         return std::pair<long double,long double>(xsum/wsum_, sum/wsum_);
     }
 
     long double DensityOrthoPoly1D::monicpoly(const unsigned degree,
                                               const long double x) const
     {
         long double polyk = 1.0L, polykm1 = 0.0L;
         for (unsigned k=0; k<degree; ++k)
         {
             const long double p = (x - rCoeffs_[k].alpha)*polyk -
                                    rCoeffs_[k].beta*polykm1;
             polykm1 = polyk;
             polyk = p;
         }
         return polyk;
     }
 
     std::pair<long double,long double>
     DensityOrthoPoly1D::recurrenceCoeffs(const unsigned deg) const
     {
         const Recurrence& r = rCoeffs_.at(deg);
         return std::pair<long double,long double>(r.alpha, r.sqrbeta);
     }
 
     DensityOrthoPoly1D::DensityOrthoPoly1D(const DensityOrthoPoly1D& r)
         : distro_(r.distro_->clone()),
           rCoeffs_(r.rCoeffs_),
           wsum_(r.wsum_),
           maxdeg_(r.maxdeg_)
     {
     }
 
     DensityOrthoPoly1D& DensityOrthoPoly1D::operator=(const DensityOrthoPoly1D& r)
     {
         if (&r != this)
         {
             distro_ = CPP11_auto_ptr<AbsDistribution1D>(r.distro_->clone());
             rCoeffs_ = r.rCoeffs_;
             wsum_ = r.wsum_;
             maxdeg_ = r.maxdeg_;
         }
         return *this;
     }
 }