Index: trunk/npstat/stat/PolynomialDistro1D.cc
===================================================================
--- trunk/npstat/stat/PolynomialDistro1D.cc	(revision 887)
+++ trunk/npstat/stat/PolynomialDistro1D.cc	(revision 888)
@@ -1,187 +1,209 @@
 #include <cassert>
 #include <stdexcept>
 #include <sstream>
 #include <climits>
 #include <cfloat>
-#include <algorithm>
 
 #include "npstat/nm/findRootUsingBisections.hh"
 #include "npstat/nm/MathUtils.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)
     {
         setupCoeffs(coeffs, maxdeg);
         checkPositive(nCheck);
     }
 
     PolynomialDistro1D::PolynomialDistro1D(
         const double location, const double scale,
         const std::vector<double>& coeffs,
         const unsigned nCheck)
         : AbsScalableDistribution1D(location, scale)
     {
         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];
         for (unsigned i=maxdeg; i; --i)
         {
             if (allCoeffs_[i] == 0.0)
                 allCoeffs_.pop_back();
             else
                 break;
         }
 
         const unsigned deg = allCoeffs_.size() - 1;
         integCoeffs_.resize(deg + 2U);
         poly_.integralCoeffs(&allCoeffs_[0], deg, &integCoeffs_[0]);
     }
 
     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;
-        unsigned imin = 0;
-        double dmin = DBL_MAX;
-        for (unsigned i=0; i<nCheck; ++i)
+        const unsigned deg = allCoeffs_.size() - 1;
+        if (deg > 1U)
         {
-            double x = step*i;
-            if (i == nM1)
-                x = 1.0;
-            const double d = unscaledDensity(x);
-            if (d < dmin)
+            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;
+            unsigned imin = 0;
+            double dmin = DBL_MAX;
+            for (unsigned i=0; i<nCheck; ++i)
             {
-                imin = i;
-                dmin = d;
+                double x = step*i;
+                if (i == nM1)
+                    x = 1.0;
+                const double d = unscaledDensity(x);
+                if (d < dmin)
+                {
+                    imin = i;
+                    dmin = d;
+                }
+                if (d <= 0.0)
+                {
+                    std::ostringstream os;
+                    os.precision(17);
+                    os << "In npstat::PolynomialDistro1D::checkPositive: "
+                       << "density is not positive at x = " << x;
+                    throw std::invalid_argument(os.str());
+                }
             }
-            if (d <= 0.0)
+            if (imin > 0U && imin < nM1)
             {
-                std::ostringstream os;
-                os.precision(17);
-                os << "In npstat::PolynomialDistro1D::checkPositive: "
-                   << "density is not positive at x = " << x;
-                throw std::invalid_argument(os.str());
+                const double x = step*imin;
+                const double xleft = step*(imin-1U);
+                double xright = step*(imin+1U);
+                if (xright > 1.0)
+                    xright = 1.0;
+                double xmin, ymin;
+                const bool minFound = parabolicExtremum(
+                    xleft, unscaledDensity(xleft),
+                    x, dmin,
+                    xright, unscaledDensity(xright),
+                    &xmin, &ymin);
+                assert(minFound);
+                if (xmin < 0.0)
+                    xmin = 0.0;
+                else if (xmin > 1.0)
+                    xmin = 1.0;
+                const double d = unscaledDensity(xmin);
+                if (d <= 0.0)
+                {
+                    std::ostringstream os;
+                    os.precision(17);
+                    os << "In npstat::PolynomialDistro1D::checkPositive: "
+                       << "density is not positive at x = " << xmin;
+                    throw std::invalid_argument(os.str());
+                }
             }
         }
-        if (imin > 0U && imin < nM1 && allCoeffs_.size() > 2U)
+        else
         {
-            const double x = step*imin;
-            const double xleft = step*(imin-1U);
-            double xright = step*(imin+1U);
-            if (xright > 1.0)
-                xright = 1.0;
-            double xmin, ymin;
-            const bool minFound = parabolicExtremum(
-                xleft, unscaledDensity(xleft),
-                x, dmin,
-                xright, unscaledDensity(xright),
-                &xmin, &ymin);
-            assert(minFound);
-            const double d = unscaledDensity(xmin);
-            if (d <= 0.0)
+            const double xset[2] = {0.0, 1.0};
+            for (unsigned i=0; i<2U; ++i)
             {
-                std::ostringstream os;
-                os.precision(17);
-                os << "In npstat::PolynomialDistro1D::checkPositive: "
-                   << "density is not positive at x = " << xmin;
-                throw std::invalid_argument(os.str());
+                const double x = xset[i];
+                if (unscaledDensity(x) < 0.0)
+                {
+                    std::ostringstream os;
+                    os << "In npstat::PolynomialDistro1D::checkPositive: "
+                       << "density is negative 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 poly_.series(&integCoeffs_[0], integCoeffs_.size()-1, 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 1.0 - unscaledCdf(x);
     }
 
     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;
     }
 }