Index: trunk/npstat/stat/SeriesCGF1D.cc
===================================================================
--- trunk/npstat/stat/SeriesCGF1D.cc	(revision 754)
+++ trunk/npstat/stat/SeriesCGF1D.cc	(revision 755)
@@ -1,92 +1,108 @@
 #include <cmath>
 #include <cassert>
 #include <stdexcept>
 
 #include "npstat/stat/SeriesCGF1D.hh"
 #include "npstat/rng/permutation.hh"
+#include "npstat/nm/Poly1D.hh"
 
 namespace npstat {
     SeriesCGF1D::SeriesCGF1D(const double* cumulants,
                              const unsigned nCumulants,
                              const double i_smin, const double i_smax)
         : cumulants_(cumulants, cumulants+nCumulants),
           smin_(i_smin),
           smax_(i_smax)
     {
         assert(cumulants);
         if (nCumulants < 2U) throw std::invalid_argument(
-            "In npstat::SeriesCGF1D constructor : not enough cumulants");
+            "In npstat::SeriesCGF1D constructor: not enough cumulants");
         if (cumulants[1] <= 0.0) throw std::invalid_argument(
-            "In npstat::SeriesCGF1D constructor : variance must be positive");
+            "In npstat::SeriesCGF1D constructor: variance must be positive");
         if (!(smin_ < 0.0)) throw std::invalid_argument(
-            "In npstat::SeriesCGF1D constructor : invalid support infimum");
+            "In npstat::SeriesCGF1D constructor: invalid support infimum");
         if (!(smax_ > 0.0)) throw std::invalid_argument(
-            "In npstat::SeriesCGF1D constructor : invalid support supremum");
+            "In npstat::SeriesCGF1D constructor: invalid support supremum");
     }
 
     bool SeriesCGF1D::isEqual(const AbsCGF1D& other) const
     {
         const SeriesCGF1D& r = static_cast<const SeriesCGF1D&>(other);
         return cumulants_ == r.cumulants_ && smin_ == r.smin_ && smax_ == r.smax_;
     }
 
     double SeriesCGF1D::derivative(const unsigned order, const double s) const
     {
         if (s == 0.0)
         {
             if (order)
             {
                 const unsigned idx = order - 1U;
                 if (idx < cumulants_.size())
                     return cumulants_[idx];
             }
         }
         else if (!(smin_ < s && s < smax_)) throw std::invalid_argument(
-            "In npstat::SeriesCGF1D::derivative : argument out of range");
+            "In npstat::SeriesCGF1D::derivative: argument out of range");
         else
         {
             const unsigned degree = cumulants_.size();
             if (order <= degree)
             {
                 const double* a = &cumulants_[0];
                 long double sum = 0.0L;
                 if (order)
                 {
                     const int shift = order-1U;
                     for (int deg=degree-order; deg>=0; --deg)
                     {
                         sum *= s;
                         sum += a[deg+shift]/ldfactorial(deg);
                     }
                 }
                 else
                 {
                     for (unsigned deg=degree; deg>0U; --deg)
                     {
                         sum += a[deg-1]/ldfactorial(deg);
                         sum *= s;
                     }
                 }
                 return sum;
             }
         }
         return 0.0;
     }
 
     SeriesCGF1D* SeriesCGF1D::shiftAndScale(const double mu,
                                             const double sigma) const
     {
         if (!(sigma > 0.0)) throw std::invalid_argument(
-            "In npstat::SeriesCGF1D::shiftAndScale : scale must be positive");
+            "In npstat::SeriesCGF1D::shiftAndScale: scale must be positive");
         std::vector<double> cum(cumulants_);
         const unsigned sz = cum.size();
         cum[0] -= mu;
         if (sigma != 1.0)
             for (unsigned i=0; i<sz; ++i)
                 cum[i] /= pow(sigma, i+1U);
 
         // Rethink this. It is not obvious how the
         // support should be modified.
         return new SeriesCGF1D(&cum[0], sz, smin_, smax_);
     }
+
+    bool SeriesCGF1D::isConvex(const double a, const double b) const
+    {
+        if (a >= b) throw std::invalid_argument(
+            "In npstat::SeriesCGF1D::sisConvex: invalid interval");
+        const unsigned degree = cumulants_.size();
+        std::vector<long double> coeffs(degree + 1U);
+        coeffs[0] = 0.0L;
+        for (unsigned deg=0; deg<degree; ++deg)
+            coeffs[deg+1] = cumulants_[deg]/ldfactorial(deg+1U);
+        const Poly1D cgf0(coeffs);
+        const Poly1D& cgf1(cgf0.derivative());
+        const Poly1D& cgf2(cgf1.derivative());
+        return cgf2((a + b)/2.0) > 0.0L && !cgf2.nRoots(a, b);
+    }
 }
Index: trunk/npstat/stat/SeriesCGF1D.hh
===================================================================
--- trunk/npstat/stat/SeriesCGF1D.hh	(revision 754)
+++ trunk/npstat/stat/SeriesCGF1D.hh	(revision 755)
@@ -1,59 +1,65 @@
 #ifndef NPSTAT_SERIESCGF1D_HH_
 #define NPSTAT_SERIESCGF1D_HH_
 
 /*!
 // \file SeriesCGF1D.hh
 //
 // \brief Univariate cumulant generating function defined by cumulants
 //
 // Author: I. Volobouev
 //
 // December 2019
 */
 
 #include <vector>
 
 #include "npstat/stat/AbsCGF1D.hh"
 
 namespace npstat {
     class SeriesCGF1D : public AbsCGF1D
     {
     public:
         /**
         // At least the first two cumulants (mean and variance)
         // must be provided in the "cumulants" array.
         // cumulants[0] must contain the mean, cumulants[1] the
         // variance, etc. "nCumulants" specifies the length of
         // the array (and also the highest cumulant order).
         // The variance must be positive (i.e., the degenerate
         // distribution of a single point mass is excluded from
         // consideration). The parameters "smin" and "smax"
         // should be chosen in such a way that the CGF remains
         // convex on the (smin, smax) interval.
         */
         SeriesCGF1D(const double* cumulants, unsigned nCumulants,
                     double smin, double smax);
 
         inline virtual SeriesCGF1D* clone() const
             {return new SeriesCGF1D(*this);}
 
         inline virtual ~SeriesCGF1D() {}
 
         virtual SeriesCGF1D* shiftAndScale(double mu, double sigma) const;
 
         inline virtual double smin() const {return smin_;}
         inline virtual double smax() const {return smax_;}
 
         virtual double derivative(unsigned order, double s) const;
 
+        /**
+        // Check whether the CGF is convex on the interval (a, b].
+        // This method ignores smin and smax.
+        */
+        virtual bool isConvex(double a, double b) const;
+
     protected:
         virtual bool isEqual(const AbsCGF1D& r) const;
 
     private:
         std::vector<double> cumulants_;
         double smin_;
         double smax_;
     };
 }
 
 #endif // NPSTAT_SERIESCGF1D_HH_