Index: trunk/tests/test_ContOrthoPoly1D.cc
===================================================================
--- trunk/tests/test_ContOrthoPoly1D.cc	(revision 609)
+++ trunk/tests/test_ContOrthoPoly1D.cc	(revision 610)
@@ -1,244 +1,244 @@
 #include "UnitTest++.h"
 #include "test_utils.hh"
 
 #include "npstat/nm/ContOrthoPoly1D.hh"
 #include "npstat/nm/ClassicalOrthoPolys1D.hh"
 #include "npstat/nm/FejerQuadrature.hh"
 #include "npstat/nm/StorablePolySeries1D.hh"
 
 #include "npstat/rng/MersenneTwister.hh"
 
 using namespace npstat;
 using namespace std;
 
 inline static int kdelta(const unsigned i, const unsigned j)
 {
     return i == j ? 1 : 0;
 }
 
 namespace {
     TEST(ContOrthoPoly1D_orthonormalization)
     {
         const double eps = 1.0e-15;
 
         const OrthoPolyMethod method[] = {OPOLY_STIELTJES, OPOLY_LANCZOS};
         const unsigned nMethods = sizeof(method)/sizeof(method[0]);
             
         const unsigned npoints = 64;
         const unsigned maxdeg = 10;
         const unsigned ntries = 10;
 
         std::vector<double> points(2U*npoints);
         std::vector<ContOrthoPoly1D::MeasurePoint> measure;
         measure.reserve(npoints);
         std::vector<double> coeffs(maxdeg + 1);
 
         MersenneTwister rng;
 
         for (unsigned imeth=0; imeth<nMethods; ++imeth)
             for (unsigned itry=0; itry<ntries; ++itry)
             {
                 const double xmin = rng()*5.0 - 6.0;
                 const double xmax = rng()*2.0 + 1.0;
                 const double range = xmax - xmin;
                 rng.run(&points[0], 2U*npoints, 2U*npoints);
 
                 measure.clear();
                 for (unsigned i=0; i<npoints; ++i)
                 {
                     ContOrthoPoly1D::MeasurePoint p(points[2*i]*range + xmin, points[2*i+1]);
                     measure.push_back(p);
                 }
                 ContOrthoPoly1D poly(maxdeg, measure, method[imeth]);
 
                 for (unsigned i=0; i<=maxdeg; ++i)
                     for (unsigned j=0; j<=maxdeg; ++j)
                     {
                         const double d = poly.empiricalKroneckerDelta(i, j);
                         CHECK_CLOSE(static_cast<double>(kdelta(i, j)), d, eps);
                     }
 
                 measure.clear();
                 for (unsigned i=0; i<npoints; ++i)
                 {
                     ContOrthoPoly1D::MeasurePoint p(points[2*i]*range + xmin, 1.0);
                     measure.push_back(p);
                 }
                 ContOrthoPoly1D poly2(maxdeg, measure, method[imeth]);
 
                 poly2.weightCoeffs(&coeffs[0], maxdeg);
                 for (unsigned i=0; i<=maxdeg; ++i)
                     CHECK_CLOSE(static_cast<double>(kdelta(i, 0)), coeffs[i], eps);
             }
     }
 
     TEST(ContOrthoPoly1D_weightCoeffs)
     {
         const double eps = 1.0e-7;
 
         const unsigned maxdeg = 10;
         const unsigned ntries = 10;
         const unsigned npoints = maxdeg + 1U;
 
         std::vector<double> points(2U*npoints);
         std::vector<ContOrthoPoly1D::MeasurePoint> measure;
         measure.reserve(npoints);
         std::vector<double> coeffs(maxdeg + 1);
 
         MersenneTwister rng;
 
         for (unsigned itry=0; itry<ntries; ++itry)
         {
             const double xmin = rng()*5.0 - 6.0;
             const double xmax = rng()*2.0 + 1.0;
             const double range = xmax - xmin;
             rng.run(&points[0], 2U*npoints, 2U*npoints);
 
             measure.clear();
             for (unsigned i=0; i<npoints; ++i)
             {
                 ContOrthoPoly1D::MeasurePoint p(points[2*i]*range + xmin, points[2*i+1]);
                 measure.push_back(p);
             }
             ContOrthoPoly1D poly(maxdeg, measure, OPOLY_LANCZOS);
             poly.weightCoeffs(&coeffs[0], maxdeg);
             CPP11_auto_ptr<StorablePolySeries1D> poly2(
                 poly.makeStorablePolySeries(xmin, xmax));
 
             for (unsigned deg=0; deg<=maxdeg; ++deg)
                 for (unsigned ipow=0; ipow<4; ++ipow)
                 {
                     const double i1 = poly.integratePoly(deg, ipow, xmin, xmax);
                     const double i2 = poly2->integratePoly(deg, ipow);
                     if (fabs(i2) > 1.0e-10)
                         CHECK_CLOSE(1.0, i1/i2, 1.0e-10);
                 }
 
             for (unsigned i=0; i<npoints; ++i)
             {
                 const double x = measure[i].first;
                 const double w = measure[i].second;
                 const double fvalue = poly.series(&coeffs[0], maxdeg, x);
                 const double f2 = poly2->series(&coeffs[0], maxdeg, x);
-                CHECK_CLOSE(w, fvalue, 10*eps);
+                CHECK_CLOSE(w, fvalue, 100*eps);
                 CHECK_CLOSE(fvalue, f2, eps);
             }
         }
     }
 
     TEST(ContOrthoPoly1D_cov8)
     {
         const double eps = 1.0e-12;
 
         MersenneTwister rng;
         const unsigned npoints = 64;
         const unsigned maxdeg = 6;
         const unsigned ntries = 5;
         const unsigned degtries = 100;
         std::vector<double> points(npoints);
 
         for (unsigned itry=0; itry<ntries; ++itry)
         {
             rng.run(&points[0], npoints, npoints);
             ContOrthoPoly1D poly(maxdeg, points);
 
             for (unsigned ideg=0; ideg<degtries; ++ideg)
             {
                 unsigned d[8];
                 bool has0 = true;
                 while (has0)
                 {
                     for (unsigned i=0; i<8; ++i)
                         d[i] = (maxdeg + 0.99)*rng();
                     has0 = (d[0] == 0 && d[1] == 0) ||
                            (d[2] == 0 && d[3] == 0) ||
                            (d[4] == 0 && d[5] == 0) ||
                            (d[6] == 0 && d[7] == 0);
                 }
                 const double cov8 = poly.cov8(d[0],d[1],d[2],d[3],d[4],d[5],d[6],d[7]);
                 const double check = poly.slowCov8(d[0],d[1],d[2],d[3],d[4],d[5],d[6],d[7]);
                 CHECK_CLOSE(check, cov8, eps);
             }
         }
     }
 
     TEST(ContOrthoPoly1D_epsCovariance)
     {
         const double eps = 1.0e-15;
         const unsigned npoints = 64;
         const unsigned maxdeg = 6;
         const unsigned ntries = 5;
 
         MersenneTwister rng;
         std::vector<double> points(npoints);
 
         for (unsigned itry=0; itry<ntries; ++itry)
         {
             rng.run(&points[0], npoints, npoints);
             ContOrthoPoly1D poly(maxdeg, points);
 
             for (unsigned i=0; i<=maxdeg; ++i)
                 for (unsigned j=0; j<=i; ++j)
                 {
                     CHECK(poly.epsExpectation(i, j, false) == 0.0);
 
                     for (unsigned k=0; k<=maxdeg; ++k)
                         for (unsigned m=0; m<=k; ++m)
                         {
                             const double eps1 = poly.empiricalKroneckerCovariance(i, j, k, m);
                             const double eps2 = poly.epsCovariance(i, j, k, m, false);
                             CHECK_CLOSE(eps1, eps2, eps);
                         }
                 }
         }
     }
 
     // Test that we can reproduce Jacobi polynomials using ContOrthoPoly1D
     TEST(ContOrthoPoly1D_jacobi_1)
     {
         const double eps = 1.0e-13;
 
         const unsigned alpha = 3;
         const unsigned beta = 4;
         const unsigned maxdeg = 5;
         const unsigned npoints = 64;
 
         JacobiOrthoPoly1D jac(alpha, beta);
         OrthoPoly1DWeight weight(jac);
 
         const unsigned npt = FejerQuadrature::minimalExactRule(2*maxdeg+alpha+beta);
         FejerQuadrature fej(npt);
         ContOrthoPoly1D poly(maxdeg, fej.weightedIntegrationPoints(weight, -1.0, 1.0), OPOLY_LANCZOS);
 
         for (unsigned i=0; i<npoints; ++i)
         {
             const double x = test_rng()*2.0 - 1.0;
             for (unsigned deg=0; deg<=maxdeg; ++deg)
                 CHECK_CLOSE(jac.poly(deg, x), poly.poly(deg, x), eps);
         }
     }
 
     TEST(ContOrthoPoly1D_jacobi_2)
     {
         const double eps = 1.0e-13;
 
         const unsigned alpha = 3;
         const unsigned beta = 4;
         const unsigned maxdeg = 5;
         const unsigned npoints = 64;
 
         JacobiOrthoPoly1D jac(alpha, beta);
         OrthoPoly1DWeight weight(jac);
 
         const unsigned npt = GaussLegendreQuadrature::minimalExactRule(2*maxdeg+alpha+beta);
         GaussLegendreQuadrature glq(npt);
         ContOrthoPoly1D poly(maxdeg, glq.weightedIntegrationPoints(weight, -1.0, 1.0), OPOLY_LANCZOS);
 
         for (unsigned i=0; i<npoints; ++i)
         {
             const double x = test_rng()*2.0 - 1.0;
             for (unsigned deg=0; deg<=maxdeg; ++deg)
                 CHECK_CLOSE(jac.poly(deg, x), poly.poly(deg, x), eps);
         }
     }
 }
Index: trunk/npstat/stat/EdgeworthSeries1D.cc
===================================================================
--- trunk/npstat/stat/EdgeworthSeries1D.cc	(revision 609)
+++ trunk/npstat/stat/EdgeworthSeries1D.cc	(revision 610)
@@ -1,610 +1,644 @@
 #include <cmath>
 
 #include "npstat/nm/MathUtils.hh"
 #include "npstat/nm/SpecialFunctions.hh"
 #include "npstat/nm/findRootUsingBisections.hh"
 #include "npstat/nm/GaussHermiteQuadrature.hh"
 
 #include "npstat/stat/EdgeworthSeries1D.hh"
 #include "npstat/stat/distributionReadError.hh"
 
 #include "geners/binaryIO.hh"
 #include "geners/vectorIO.hh"
 
-#define MAXCOEFFS 36
+#define EDGEWORTH_MAXORDER 14U
+#define MAXCOEFFS (EDGEWORTH_MAXORDER*3U)
 #define LSQR2 1.41421356237309504880169L
 #define SQRTWOPIL 2.50662827463100050241577L
 
 static long double ldgexceedance(const long double x)
 {
     return erfcl(x/LSQR2)/2.0L;
 }
 
 static long double ldgcdf(const long double x)
 {
     return erfcl(-x/LSQR2)/2.0L;
 }
 
 namespace npstat {
     EdgeworthSeries1D::EdgeworthSeries1D(
         const std::vector<double>& cumulants,
         const EdgeworthSeriesMethod m, const unsigned i_order,
         const bool i_slrMode)
         : cumulants_(cumulants), m_(m), order_(i_order), slrMode_(i_slrMode)
     {
         unsigned maxOrder = 4U;
         if (m == EDGEWORTH_CLASSICAL && !slrMode_)
-            maxOrder = 12U;
+            maxOrder = EDGEWORTH_MAXORDER;
         if (order_ > maxOrder) throw std::invalid_argument(
             "In npstat::EdgeworthSeries1D constructor: "
             "order parameter out of range");
         if (order_)
         {
             const unsigned nMin = minNCumulants();
             if (cumulants_.size() < nMin)
                 throw std::invalid_argument(
                     "In npstat::EdgeworthSeries1D constructor: "
                     "not enough cumulants provided");
             if (nMin > 1U)
                 if (cumulants_[1] <= 0.0)
                     throw std::invalid_argument(
                         "In npstat::EdgeworthSeries1D constructor: "
                         "variance must be positive");
         }
     }
 
     double EdgeworthSeries1D::edgMean() const
     {
         if (minNCumulants() && m_ == EDGEWORTH_CLASSICAL)
             return cumulants_.at(0);
         else
             return 0.0;
     }
 
     double EdgeworthSeries1D::edgStdev() const
     {
         if (minNCumulants() > 1)
             return normalizedSigma();
         else
             return 1.0;
     }
 
     unsigned EdgeworthSeries1D::minNCumulants() const
     {
         if (slrMode_)
             return order_;
         else
             return order_ ? order_ + 2U : 0U;
     }
 
     bool EdgeworthSeries1D::isEqual(const AbsDistribution1D& o) const
     {
         const EdgeworthSeries1D& r = static_cast<const EdgeworthSeries1D&>(o);
         return cumulants_ == r.cumulants_ && m_ == r.m_ &&
                order_ == r.order_ && slrMode_ == r.slrMode_;
     }
 
     bool EdgeworthSeries1D::write(std::ostream& os) const
     {
         const unsigned m = static_cast<unsigned>(m_);
         const unsigned char slr = slrMode_;
         gs::write_pod(os, m);
         gs::write_pod(os, order_);
         gs::write_pod(os, slr);
         return gs::write_item(os, cumulants_, false) && !os.fail();
     }
 
     EdgeworthSeries1D* EdgeworthSeries1D::read(
         const gs::ClassId& id, std::istream& in)
     {
         static const gs::ClassId myClassId(
             gs::ClassId::makeId<EdgeworthSeries1D>());
         myClassId.ensureSameId(id);
 
         unsigned m, order;
         unsigned char slr;
         gs::read_pod(in, &m);
         gs::read_pod(in, &order);
         gs::read_pod(in, &slr);
 
         std::vector<double> cumulants;
         gs::restore_item(in, &cumulants, false);
 
         if (in.fail())
         {
             distributionReadError(in, classname());
             return 0;
         }
         else
             return new EdgeworthSeries1D(
                 cumulants, static_cast<EdgeworthSeriesMethod>(m), order, slr);
     }
 
     double EdgeworthSeries1D::densityFactor(const double xin) const
     {
         long double coeffs[MAXCOEFFS+1];
         unsigned maxdeg = 0;
 
         if (order_)
         {
             if (slrMode_)
                 hermiteCoeffsSLR(coeffs+1, &maxdeg);
             else
                 hermiteCoeffsNormal(coeffs+1, &maxdeg);
             ++maxdeg;
         }
         coeffs[0] = 0.0L;
         return 1.0 + hermiteSeriesSumProb(coeffs, maxdeg, normalizedCoord(xin))/normalizedSigma();
     }
 
     double EdgeworthSeries1D::normalizedSigma() const
     {
         double sigma = 1.0;
         if (order_ > 1U && m_ == EDGEWORTH_CLASSICAL)
             sigma = sqrt(cumulants_.at(1));
         return sigma;
     }
 
     double EdgeworthSeries1D::density(const double xin) const
     {
         const double x = normalizedCoord(xin);
         return exp(-x*x/2)/SQRTWOPIL/normalizedSigma()*densityFactor(xin);
     }
 
     double EdgeworthSeries1D::cdf(const double xin) const
     {
         const double z = normalizedCoord(xin);
         return ldgcdf(z) - expl(-z*z/2.0)/SQRTWOPIL/normalizedSigma()*cdfFactor(xin);
     }
 
     double EdgeworthSeries1D::exceedance(const double xin) const
     {
         const double z = normalizedCoord(xin);
         return ldgexceedance(z) + expl(-z*z/2.0)/SQRTWOPIL/normalizedSigma()*cdfFactor(xin);
     }
 
     double EdgeworthSeries1D::quantile(const double r1) const
     {
         if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error(
             "In npstat::EdgeworthSeries1D::quantile: "
             "cdf argument outside of [0, 1] interval");
         if (r1 == 0.0 || r1 == 1.0)
             return inverseGaussCdf(r1);
         const double x0 = inverseGaussCdf(0.0);
         const double x1 = inverseGaussCdf(1.0);
         CdfFunctor1D fcn(*this);
         double q = 0.0;
         if (!findRootUsingBisections(fcn, r1, x0, x1, 1.0e-15, &q))
             throw std::runtime_error("In npstat::EdgeworthSeries1D::quantile: "
                                      "root finding failed");
         return q;
     }
 
     double EdgeworthSeries1D::inverseExceedance(const double r1) const
     {
         if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error(
             "In npstat::EdgeworthSeries1D::inverseExceedance: "
             "exceedance argument outside of [0, 1] interval");
         if (r1 == 0.0 || r1 == 1.0)
             return inverseGaussCdf(1.0 - r1);
         const double x0 = inverseGaussCdf(0.0);
         const double x1 = inverseGaussCdf(1.0);
         ExceedanceFunctor1D fcn(*this);
         double q = 0.0;
         if (!findRootUsingBisections(fcn, r1, x0, x1, 1.0e-15, &q))
             throw std::runtime_error("In npstat::EdgeworthSeries1D::inverseExceedance: "
                                      "root finding failed");
         return q;
     }
 
     double EdgeworthSeries1D::normalizedCoord(const double xin) const
     {
         double x = xin;
         if (order_)
         {
             switch (m_)
             {
             case EDGEWORTH_SEVERINI:
                 break;
 
             case EDGEWORTH_CLASSICAL:
                 {
                     double sigma = 1.0;
                     if (order_ > 1U)
                         sigma = sqrt(cumulants_.at(1));
                     x = (xin - cumulants_.at(0))/sigma;
                 }
                 break;
 
             default:
                 assert(!"Incomplete swich statement. This is a bug. Please report.");
             }
         }
         return x;
     }
 
     double EdgeworthSeries1D::cdfFactor(const double xin) const
     {
         if (order_)
         {
             long double coeffs[MAXCOEFFS];
             unsigned maxdeg;
             if (slrMode_)
                 hermiteCoeffsSLR(coeffs, &maxdeg);
             else
                 hermiteCoeffsNormal(coeffs, &maxdeg);
             return hermiteSeriesSumProb(coeffs, maxdeg, normalizedCoord(xin));
         }
         else
             return 0.0;
     }
 
     double EdgeworthSeries1D::empiricalCentralMoment(const unsigned k) const
     {
         switch (k)
         {
         case 0U:
             return 1.0;
 
         case 1U:
             return 0.0;
 
         default:
             if (k > 200U)
                 throw std::invalid_argument("In npstat::EdgeworthSeries1D::empiricalCentralMoment: "
                                             "moment order is too high");
             GaussHermiteQuadrature quad(128);
             EdgeworthSeriesMomentFcn fcn(*this, k);
             return quad.integrateProb(edgMean(), edgStdev(), fcn);
         }
     }
 
     void EdgeworthSeries1D::hermiteCoeffsSLR(long double* coeffs, unsigned* maxdeg) const
     {
         for (unsigned i=0; i<MAXCOEFFS; ++i)
             coeffs[i] = 0.0L;
         *maxdeg = 0U;
 
         assert(order_);
         const unsigned nk = cumulants_.size();
         assert(nk >= order_);
         *maxdeg = order_ - 1U;
         assert(MAXCOEFFS > *maxdeg);
         const double k1 = m_ == EDGEWORTH_SEVERINI ? cumulants_[0] : 0.0;
         coeffs[0] = k1;
 
         if (order_ > 1U)
         {
             const double k2 = m_ == EDGEWORTH_SEVERINI ? cumulants_[1] : 1.0;
             const double k2m1 = k2 - 1.0;
             coeffs[1] = (k1*k1 + k2m1)/2.0;
             
             if (order_ > 2U)
             {
                 const double k3 = cumulants_[2];
                 coeffs[2] = (k1*(k1*k1 + 3*k2m1) + k3)/6;
                 
                 if (order_ > 3U)
                 {
                     const double k4 = cumulants_[3];
                     const double k1sq = k1*k1;
                     coeffs[3] = (k1sq*k1sq + 6*k1sq*k2m1 + 3*k2m1*k2m1 + 4*k1*k3 + k4)/24;
                     
                     if (order_ > 4U)
                         // We should not be here
                         assert(0);
                 }
 
                 if (m_ == EDGEWORTH_CLASSICAL)
                 {
                     const double sigma = sqrt(cumulants_[1]);
                     for (unsigned i=2; i<=*maxdeg; ++i)
                         coeffs[i] /= powl(sigma, i);
                 }
             }
         }
     }
 
     void EdgeworthSeries1D::hermiteCoeffsClassical(long double* coeffs, unsigned* maxdeg) const
     {
         for (unsigned i=0; i<MAXCOEFFS; ++i)
             coeffs[i] = 0.0L;
         *maxdeg = 0U;
 
         assert(order_);
         const unsigned nk = cumulants_.size();
         assert(nk > 2U);
 
         const double kappa3 = cumulants_[2];
 
         *maxdeg = 2U;
-        assert(MAXCOEFFS > *maxdeg);
         coeffs[2] = kappa3/6.0;
 
         if (order_ >= 2U)
         {
             *maxdeg = 5U;
-            assert(MAXCOEFFS > *maxdeg);
 
             const double kappa4 = cumulants_[3];
             const double kappa3_2 = kappa3*kappa3;
             coeffs[3] += kappa4/24;
             coeffs[5] += kappa3_2/72;
 
             if (order_ >= 3U)
             {
                 *maxdeg = 8U;
-                assert(MAXCOEFFS > *maxdeg);
 
                 const double kappa5 = cumulants_[4];
                 const double kappa3_3 = kappa3_2*kappa3;
                 coeffs[4] += kappa5/120;
                 coeffs[6] += (kappa3*kappa4)/144;
                 coeffs[8] += kappa3_3/1296;
 
                 if (order_ >= 4U)
                 {
                     *maxdeg = 11U;
-                    assert(MAXCOEFFS > *maxdeg);
 
                     const double kappa6 = cumulants_[5];
                     const double kappa3_4 = kappa3_3*kappa3;
                     const double kappa4_2 = kappa4*kappa4;
                     coeffs[5] += kappa6/720;
                     coeffs[7] += kappa4_2/1152 + (kappa3*kappa5)/720;
                     coeffs[9] += (kappa3_2*kappa4)/1728;
                     coeffs[11] += kappa3_4/31104;
 
                     if (order_ >= 5U)
                     {
                         *maxdeg = 14U;
-                        assert(MAXCOEFFS > *maxdeg);
 
                         const double kappa7 = cumulants_[6];
                         const double kappa3_5 = kappa3_4*kappa3;
                         coeffs[6] += kappa7/5040;
                         coeffs[8] += (kappa4*kappa5)/2880 + (kappa3*kappa6)/4320;
                         coeffs[10] += (kappa3*kappa4_2)/6912 + (kappa3_2*kappa5)/8640;
                         coeffs[12] += (kappa3_3*kappa4)/31104;
                         coeffs[14] += kappa3_5/933120;
 
                         if (order_ >= 6U)
                         {
                             *maxdeg = 17U;
-                            assert(MAXCOEFFS > *maxdeg);
 
                             const double kappa8 = cumulants_[7];
                             const double kappa5_2 = kappa5*kappa5;
                             const double kappa4_3 = kappa4_2*kappa4;
                             const double kappa3_6 = kappa3_5*kappa3;
                             coeffs[7] += kappa8/40320;
                             coeffs[9] += kappa5_2/28800 + (kappa4*kappa6)/17280 + (kappa3*kappa7)/30240;
                             coeffs[11] += kappa4_3/82944 + (kappa3*kappa4*kappa5)/17280 + (kappa3_2*kappa6)/51840;
                             coeffs[13] += (kappa3_2*kappa4_2)/82944 + (kappa3_3*kappa5)/155520;
                             coeffs[15] += (kappa3_4*kappa4)/746496;
                             coeffs[17] += kappa3_6/33592320;
 
                             if (order_ >= 7U)
                             {
                                 *maxdeg = 20U;
-                                assert(MAXCOEFFS > *maxdeg);
 
                                 const double kappa9 = cumulants_[8];
                                 const double kappa3_7 = kappa3_6*kappa3;
                                 coeffs[8] += kappa9/362880.0;
                                 coeffs[10] += (kappa5*kappa6)/86400 + (kappa4*kappa7)/120960 + (kappa3*kappa8)/241920.0;
                                 coeffs[12] += (kappa4_2*kappa5)/138240 + (kappa3*kappa5_2)/172800 + (kappa3*kappa4*kappa6)/103680 + (kappa3_2*kappa7)/362880.0;
                                 coeffs[14] += (kappa3*kappa4_3)/497664 + (kappa3_2*kappa4*kappa5)/207360 + (kappa3_3*kappa6)/933120.0;
                                 coeffs[16] += (kappa3_3*kappa4_2)/1492992 + (kappa3_4*kappa5)/3732480.0;
                                 coeffs[18] += (kappa3_5*kappa4)/22394880.0;
                                 coeffs[20] += kappa3_7/1410877440.0;
 
                                 if (order_ >= 8U)
                                 {
                                     *maxdeg = 23U;
-                                    assert(MAXCOEFFS > *maxdeg);
 
                                     const double kappa10 = cumulants_[9];
                                     const double kappa6_2 = kappa6*kappa6;
                                     const double kappa3_8 = kappa3_7*kappa3;
                                     const double kappa4_4 = kappa4_3*kappa4;
                                     coeffs[9] += kappa10/3628800.0;
                                     coeffs[11] += kappa6_2/1036800 + (kappa5*kappa7)/604800 + (kappa4*kappa8)/967680 + (kappa3*kappa9)/2177280.0;
                                     coeffs[13] += (kappa4*kappa5_2)/691200 + (kappa4_2*kappa6)/829440 + (kappa3*kappa5*kappa6)/518400 + (kappa3*kappa4*kappa7)/725760 + (kappa3_2*kappa8)/2903040.0;
                                     coeffs[15] += kappa4_4/7962624 + (kappa3*kappa4_2*kappa5)/829440 + (kappa3_2*kappa5_2)/2073600 + (kappa3_2*kappa4*kappa6)/1244160 + (kappa3_3*kappa7)/6531840.0;
                                     coeffs[17] += (kappa3_2*kappa4_3)/5971968 + (kappa3_3*kappa4*kappa5)/3732480 + (kappa3_4*kappa6)/22394880.0;
                                     coeffs[19] += (kappa3_4*kappa4_2)/35831808 + (kappa3_5*kappa5)/111974400.0;
                                     coeffs[21] += (kappa3_6*kappa4)/806215680.0;
                                     coeffs[23] += kappa3_8/67722117120.0;
 
                                     if (order_ >= 9U)
                                     {
                                         *maxdeg = 26U;
-                                        assert(MAXCOEFFS > *maxdeg);
 
                                         const double kappa11 = cumulants_[10];
                                         const double kappa5_3 = kappa5_2*kappa5;
                                         const double kappa3_9 = kappa3_8*kappa3;
                                         coeffs[10] += kappa11/39916800.0;
                                         coeffs[12] += (kappa6*kappa7)/3628800.0 + (kappa5*kappa8)/4838400.0 + (kappa4*kappa9)/8709120.0 + (kappa3*kappa10)/21772800.0;
                                         coeffs[14] += kappa5_3/10368000.0 + (kappa4*kappa5*kappa6)/2073600.0 + (kappa3*kappa6_2)/6220800.0 + (kappa4_2*kappa7)/5806080.0 + (kappa3*kappa5*kappa7)/3628800.0 + (kappa3*kappa4*kappa8)/5806080.0 + (kappa3_2*kappa9)/26127360.0;
                                         coeffs[16] += (kappa4_3*kappa5)/9953280.0 + (kappa3*kappa4*kappa5_2)/4147200.0 + (kappa3*kappa4_2*kappa6)/4976640.0 + (kappa3_2*kappa5*kappa6)/6220800.0 + (kappa3_2*kappa4*kappa7)/8709120.0 + (kappa3_3*kappa8)/52254720.0;
                                         coeffs[18] += (kappa3*kappa4_4)/47775744 + (kappa3_2*kappa4_2*kappa5)/9953280.0 + (kappa3_3*kappa5_2)/37324800.0 + (kappa3_3*kappa4*kappa6)/22394880.0 + (kappa3_4*kappa7)/156764160.0;
                                         coeffs[20] += (kappa3_3*kappa4_3)/107495424 + (kappa3_4*kappa4*kappa5)/89579520.0 + (kappa3_5*kappa6)/671846400.0;
                                         coeffs[22] += (kappa3_5*kappa4_2)/1074954240.0 + (kappa3_6*kappa5)/4031078400.0;
                                         coeffs[24] += (kappa3_7*kappa4)/33861058560.0;
                                         coeffs[26] += kappa3_9/3656994324480.0;
 
                                         if (order_ >= 10U)
                                         {
                                             *maxdeg = 29U;
-                                            assert(MAXCOEFFS > *maxdeg);
 
                                             const double kappa12 = cumulants_[11];
                                             const double kappa7_2 = kappa7*kappa7;
                                             const double kappa3_10 = kappa3_9*kappa3;
                                             const double kappa4_5 = kappa4_4*kappa4;
                                             coeffs[11] += kappa12/479001600.0;
                                             coeffs[13] += kappa7_2/50803200.0 + (kappa6*kappa8)/29030400.0 + (kappa5*kappa9)/43545600.0 + (kappa4*kappa10)/87091200.0 + (kappa3*kappa11)/239500800.0;
                                             coeffs[15] += (kappa5_2*kappa6)/20736000.0 + (kappa4*kappa6_2)/24883200.0 + (kappa4*kappa5*kappa7)/14515200.0 + (kappa3*kappa6*kappa7)/21772800.0 + (kappa4_2*kappa8)/46448640.0 + (kappa3*kappa5*kappa8)/29030400.0 + (kappa3*kappa4*kappa9)/52254720.0 + (kappa3_2*kappa10)/261273600.0;
                                             coeffs[17] += (kappa4_2*kappa5_2)/33177600.0 + (kappa3*kappa5_3)/62208000.0 + (kappa4_3*kappa6)/59719680.0 + (kappa3*kappa4*kappa5*kappa6)/12441600.0 + (kappa3_2*kappa6_2)/74649600.0 + (kappa3*kappa4_2*kappa7)/34836480.0 + (kappa3_2*kappa5*kappa7)/43545600.0 + (kappa3_2*kappa4*kappa8)/69672960.0 + (kappa3_3*kappa9)/470292480.0;
                                             coeffs[19] += kappa4_5/955514880.0 + (kappa3*kappa4_3*kappa5)/59719680.0 + (kappa3_2*kappa4*kappa5_2)/49766400.0 + (kappa3_2*kappa4_2*kappa6)/59719680.0 + (kappa3_3*kappa5*kappa6)/111974400.0 + (kappa3_3*kappa4*kappa7)/156764160.0 + (kappa3_4*kappa8)/1254113280.0;
                                             coeffs[21] += (kappa3_2*kappa4_4)/573308928 + (kappa3_3*kappa4_2*kappa5)/179159040.0 + (kappa3_4*kappa5_2)/895795200.0 + (kappa3_4*kappa4*kappa6)/537477120.0 + (kappa3_5*kappa7)/4702924800.0;
                                             coeffs[23] += (kappa3_4*kappa4_3)/2579890176.0 + (kappa3_5*kappa4*kappa5)/2687385600.0 + (kappa3_6*kappa6)/24186470400.0;
                                             coeffs[25] += (kappa3_6*kappa4_2)/38698352640.0 + (kappa3_7*kappa5)/169305292800.0;
                                             coeffs[27] += (kappa3_8*kappa4)/1625330810880.0;
                                             coeffs[29] += kappa3_10/219419659468800.0;
 
                                             if (order_ >= 11U)
                                             {
                                                 *maxdeg = 32U;
-                                                assert(MAXCOEFFS > *maxdeg);
 
                                                 const double kappa13 = cumulants_[12];
                                                 const double kappa3_11 = kappa3_10*kappa3;
-                                                coeffs[12] += kappa13/6227020800.0;
-                                                coeffs[14] += (kappa7*kappa8)/203212800.0 + (kappa6*kappa9)/261273600.0 + (kappa5*kappa10)/435456000.0 + (kappa4*kappa11)/958003200.0 + (kappa3*kappa12)/2874009600.0;
-                                                coeffs[16] += (kappa5*kappa6_2)/124416000.0 + (kappa5_2*kappa7)/145152000.0 + (kappa4*kappa6*kappa7)/87091200.0 + (kappa3*kappa7_2)/304819200.0 + (kappa4*kappa5*kappa8)/116121600.0 + (kappa3*kappa6*kappa8)/174182400.0 + (kappa4_2*kappa9)/418037760.0 + (kappa3*kappa5*kappa9)/261273600.0 + (kappa3*kappa4*kappa10)/522547200.0 + (kappa3_2*kappa11)/2874009600.0;
-                                                coeffs[18] += (kappa4*kappa5_3)/248832000.0 + (kappa4_2*kappa5*kappa6)/99532800.0 + (kappa3*kappa5_2*kappa6)/124416000.0 + (kappa3*kappa4*kappa6_2)/149299200.0 + (kappa4_3*kappa7)/418037760.0 + (kappa3*kappa4*kappa5*kappa7)/87091200.0 + (kappa3_2*kappa6*kappa7)/261273600.0 + (kappa3*kappa4_2*kappa8)/278691840.0 + (kappa3_2*kappa5*kappa8)/348364800.0 + (kappa3_2*kappa4*kappa9)/627056640.0 + (kappa3_3*kappa10)/4702924800.0;
-                                                coeffs[20] += (kappa4_4*kappa5)/955514880.0 + (kappa3*kappa4_2*kappa5_2)/199065600.0 + (kappa3_2*kappa5_3)/746496000.0 + (kappa3*kappa4_3*kappa6)/358318080.0 + (kappa3_2*kappa4*kappa5*kappa6)/149299200.0 + (kappa3_3*kappa6_2)/1343692800.0 + (kappa3_2*kappa4_2*kappa7)/418037760.0 + (kappa3_3*kappa5*kappa7)/783820800.0 + (kappa3_3*kappa4*kappa8)/1254113280.0 + (kappa3_4*kappa9)/11287019520.0;
-                                                coeffs[22] += (kappa3*kappa4_5)/5733089280.0 + (kappa3_2*kappa4_3*kappa5)/716636160.0 + (kappa3_3*kappa4*kappa5_2)/895795200.0 + (kappa3_3*kappa4_2*kappa6)/1074954240.0 + (kappa3_4*kappa5*kappa6)/2687385600.0 + (kappa3_4*kappa4*kappa7)/3762339840.0 + (kappa3_5*kappa8)/37623398400.0;
-                                                coeffs[24] += (kappa3_3*kappa4_4)/10319560704 + (kappa3_4*kappa4_2*kappa5)/4299816960.0 + (kappa3_5*kappa5_2)/26873856000.0 + (kappa3_5*kappa4*kappa6)/16124313600.0 + (kappa3_6*kappa7)/169305292800.0;
-                                                coeffs[26] += (kappa3_5*kappa4_3)/77396705280.0 + (kappa3_6*kappa4*kappa5)/96745881600.0 + (kappa3_7*kappa6)/1015831756800.0;
-                                                coeffs[28] += (kappa3_7*kappa4_2)/1625330810880.0 + (kappa3_8*kappa5)/8126654054400.0;
-                                                coeffs[30] += (kappa3_9*kappa4)/87767863787520.0;
-                                                coeffs[32] += kappa3_11/14481697524940800.0;
+                                                coeffs[12] += kappa13/6227020800.0L;
+                                                coeffs[14] += (kappa7*kappa8)/203212800.0L + (kappa6*kappa9)/261273600.0L + (kappa5*kappa10)/435456000.0L + (kappa4*kappa11)/958003200.0L + (kappa3*kappa12)/2874009600.0L;
+                                                coeffs[16] += (kappa5*kappa6_2)/124416000.0L + (kappa5_2*kappa7)/145152000.0L + (kappa4*kappa6*kappa7)/87091200.0L + (kappa3*kappa7_2)/304819200.0L + (kappa4*kappa5*kappa8)/116121600.0L + (kappa3*kappa6*kappa8)/174182400.0L + (kappa4_2*kappa9)/418037760.0L + (kappa3*kappa5*kappa9)/261273600.0L + (kappa3*kappa4*kappa10)/522547200.0L + (kappa3_2*kappa11)/2874009600.0L;
+                                                coeffs[18] += (kappa4*kappa5_3)/248832000.0L + (kappa4_2*kappa5*kappa6)/99532800.0L + (kappa3*kappa5_2*kappa6)/124416000.0L + (kappa3*kappa4*kappa6_2)/149299200.0L + (kappa4_3*kappa7)/418037760.0L + (kappa3*kappa4*kappa5*kappa7)/87091200.0L + (kappa3_2*kappa6*kappa7)/261273600.0L + (kappa3*kappa4_2*kappa8)/278691840.0L + (kappa3_2*kappa5*kappa8)/348364800.0L + (kappa3_2*kappa4*kappa9)/627056640.0L + (kappa3_3*kappa10)/4702924800.0L;
+                                                coeffs[20] += (kappa4_4*kappa5)/955514880.0L + (kappa3*kappa4_2*kappa5_2)/199065600.0L + (kappa3_2*kappa5_3)/746496000.0L + (kappa3*kappa4_3*kappa6)/358318080.0L + (kappa3_2*kappa4*kappa5*kappa6)/149299200.0L + (kappa3_3*kappa6_2)/1343692800.0L + (kappa3_2*kappa4_2*kappa7)/418037760.0L + (kappa3_3*kappa5*kappa7)/783820800.0L + (kappa3_3*kappa4*kappa8)/1254113280.0L + (kappa3_4*kappa9)/11287019520.0L;
+                                                coeffs[22] += (kappa3*kappa4_5)/5733089280.0L + (kappa3_2*kappa4_3*kappa5)/716636160.0L + (kappa3_3*kappa4*kappa5_2)/895795200.0L + (kappa3_3*kappa4_2*kappa6)/1074954240.0L + (kappa3_4*kappa5*kappa6)/2687385600.0L + (kappa3_4*kappa4*kappa7)/3762339840.0L + (kappa3_5*kappa8)/37623398400.0L;
+                                                coeffs[24] += (kappa3_3*kappa4_4)/10319560704 + (kappa3_4*kappa4_2*kappa5)/4299816960.0L + (kappa3_5*kappa5_2)/26873856000.0L + (kappa3_5*kappa4*kappa6)/16124313600.0L + (kappa3_6*kappa7)/169305292800.0L;
+                                                coeffs[26] += (kappa3_5*kappa4_3)/77396705280.0L + (kappa3_6*kappa4*kappa5)/96745881600.0L + (kappa3_7*kappa6)/1015831756800.0L;
+                                                coeffs[28] += (kappa3_7*kappa4_2)/1625330810880.0L + (kappa3_8*kappa5)/8126654054400.0L;
+                                                coeffs[30] += (kappa3_9*kappa4)/87767863787520.0L;
+                                                coeffs[32] += kappa3_11/14481697524940800.0L;
 
                                                 if (order_ >= 12U)
                                                 {
                                                     *maxdeg = 35U;
-                                                    assert(MAXCOEFFS > *maxdeg);
 
                                                     const double kappa14 = cumulants_[13];
                                                     const double kappa8_2 = kappa8*kappa8;
                                                     const double kappa6_3 = kappa6_2*kappa6;
                                                     const double kappa3_12 = kappa3_11*kappa3;
                                                     const double kappa5_4 = kappa5_3*kappa5;
                                                     const double kappa4_6 = kappa4_5*kappa4;
-                                                    coeffs[13] += kappa14/87178291200.0;
-                                                    coeffs[15] += kappa8_2/3251404800.0 + (kappa7*kappa9)/1828915200.0 + (kappa6*kappa10)/2612736000.0 + (kappa5*kappa11)/4790016000.0 + (kappa4*kappa12)/11496038400.0 + (kappa3*kappa13)/37362124800.0;
-                                                    coeffs[17] += kappa6_3/2239488000.0 + (kappa5*kappa6*kappa7)/435456000.0 + (kappa4*kappa7_2)/1219276800.0 + (kappa5_2*kappa8)/1161216000.0 + (kappa4*kappa6*kappa8)/696729600.0 + (kappa3*kappa7*kappa8)/1219276800.0 + (kappa4*kappa5*kappa9)/1045094400.0 + (kappa3*kappa6*kappa9)/1567641600.0 + (kappa4_2*kappa10)/4180377600.0 + (kappa3*kappa5*kappa10)/2612736000.0 + (kappa3*kappa4*kappa11)/5748019200.0 + (kappa3_2*kappa12)/34488115200.0;
-                                                    coeffs[19] += kappa5_4/4976640000.0 + (kappa4*kappa5_2*kappa6)/497664000.0 + (kappa4_2*kappa6_2)/1194393600.0 + (kappa3*kappa5*kappa6_2)/746496000.0 + (kappa4_2*kappa5*kappa7)/696729600.0 + (kappa3*kappa5_2*kappa7)/870912000.0 + (kappa3*kappa4*kappa6*kappa7)/522547200.0 + (kappa3_2*kappa7_2)/3657830400.0 + (kappa4_3*kappa8)/3344302080.0 + (kappa3*kappa4*kappa5*kappa8)/696729600.0 + (kappa3_2*kappa6*kappa8)/2090188800.0 + (kappa3*kappa4_2*kappa9)/2508226560.0 + (kappa3_2*kappa5*kappa9)/3135283200.0 + (kappa3_2*kappa4*kappa10)/6270566400.0 + (kappa3_3*kappa11)/51732172800.0;
-                                                    coeffs[21] += (kappa4_3*kappa5_2)/2388787200.0 + (kappa3*kappa4*kappa5_3)/1492992000.0 + (kappa4_4*kappa6)/5733089280.0 + (kappa3*kappa4_2*kappa5*kappa6)/597196800.0 + (kappa3_2*kappa5_2*kappa6)/1492992000.0 + (kappa3_2*kappa4*kappa6_2)/1791590400.0 + (kappa3*kappa4_3*kappa7)/2508226560.0 + (kappa3_2*kappa4*kappa5*kappa7)/1045094400.0 + (kappa3_3*kappa6*kappa7)/4702924800.0 + (kappa3_2*kappa4_2*kappa8)/3344302080.0 + (kappa3_3*kappa5*kappa8)/6270566400.0 + (kappa3_3*kappa4*kappa9)/11287019520.0 + (kappa3_4*kappa10)/112870195200.0;
-                                                    coeffs[23] += kappa4_6/137594142720.0 + (kappa3*kappa4_4*kappa5)/5733089280.0 + (kappa3_2*kappa4_2*kappa5_2)/2388787200.0 + (kappa3_3*kappa5_3)/13436928000.0 + (kappa3_2*kappa4_3*kappa6)/4299816960.0 + (kappa3_3*kappa4*kappa5*kappa6)/2687385600.0 + (kappa3_4*kappa6_2)/32248627200.0 + (kappa3_3*kappa4_2*kappa7)/7524679680.0 + (kappa3_4*kappa5*kappa7)/18811699200.0 + (kappa3_4*kappa4*kappa8)/30098718720.0 + (kappa3_5*kappa9)/338610585600.0;
-                                                    coeffs[25] += (kappa3_2*kappa4_5)/68797071360.0 + (kappa3_3*kappa4_3*kappa5)/12899450880.0 + (kappa3_4*kappa4*kappa5_2)/21499084800.0 + (kappa3_4*kappa4_2*kappa6)/25798901760.0 + (kappa3_5*kappa5*kappa6)/80621568000.0 + (kappa3_5*kappa4*kappa7)/112870195200.0 + (kappa3_6*kappa8)/1354442342400.0;
-                                                    coeffs[27] += (kappa3_4*kappa4_4)/247669456896 + (kappa3_5*kappa4_2*kappa5)/128994508800.0 + (kappa3_6*kappa5_2)/967458816000.0 + (kappa3_6*kappa4*kappa6)/580475289600.0 + (kappa3_7*kappa7)/7110822297600.0;
-                                                    coeffs[29] += (kappa3_6*kappa4_3)/2786281390080.0 + (kappa3_7*kappa4*kappa5)/4063327027200.0 + (kappa3_8*kappa6)/48759924326400.0;
-                                                    coeffs[31] += (kappa3_8*kappa4_2)/78015878922240.0 + (kappa3_9*kappa5)/438839318937600.0;
-                                                    coeffs[33] += (kappa3_10*kappa4)/5266071827251200.0;
-                                                    coeffs[35] += kappa3_12/1042682221795737600.0;
+                                                    coeffs[13] += kappa14/87178291200.0L;
+                                                    coeffs[15] += kappa8_2/3251404800.0L + (kappa7*kappa9)/1828915200.0L + (kappa6*kappa10)/2612736000.0L + (kappa5*kappa11)/4790016000.0L + (kappa4*kappa12)/11496038400.0L + (kappa3*kappa13)/37362124800.0L;
+                                                    coeffs[17] += kappa6_3/2239488000.0L + (kappa5*kappa6*kappa7)/435456000.0L + (kappa4*kappa7_2)/1219276800.0L + (kappa5_2*kappa8)/1161216000.0L + (kappa4*kappa6*kappa8)/696729600.0L + (kappa3*kappa7*kappa8)/1219276800.0L + (kappa4*kappa5*kappa9)/1045094400.0L + (kappa3*kappa6*kappa9)/1567641600.0L + (kappa4_2*kappa10)/4180377600.0L + (kappa3*kappa5*kappa10)/2612736000.0L + (kappa3*kappa4*kappa11)/5748019200.0L + (kappa3_2*kappa12)/34488115200.0L;
+                                                    coeffs[19] += kappa5_4/4976640000.0L + (kappa4*kappa5_2*kappa6)/497664000.0L + (kappa4_2*kappa6_2)/1194393600.0L + (kappa3*kappa5*kappa6_2)/746496000.0L + (kappa4_2*kappa5*kappa7)/696729600.0L + (kappa3*kappa5_2*kappa7)/870912000.0L + (kappa3*kappa4*kappa6*kappa7)/522547200.0L + (kappa3_2*kappa7_2)/3657830400.0L + (kappa4_3*kappa8)/3344302080.0L + (kappa3*kappa4*kappa5*kappa8)/696729600.0L + (kappa3_2*kappa6*kappa8)/2090188800.0L + (kappa3*kappa4_2*kappa9)/2508226560.0L + (kappa3_2*kappa5*kappa9)/3135283200.0L + (kappa3_2*kappa4*kappa10)/6270566400.0L + (kappa3_3*kappa11)/51732172800.0L;
+                                                    coeffs[21] += (kappa4_3*kappa5_2)/2388787200.0L + (kappa3*kappa4*kappa5_3)/1492992000.0L + (kappa4_4*kappa6)/5733089280.0L + (kappa3*kappa4_2*kappa5*kappa6)/597196800.0L + (kappa3_2*kappa5_2*kappa6)/1492992000.0L + (kappa3_2*kappa4*kappa6_2)/1791590400.0L + (kappa3*kappa4_3*kappa7)/2508226560.0L + (kappa3_2*kappa4*kappa5*kappa7)/1045094400.0L + (kappa3_3*kappa6*kappa7)/4702924800.0L + (kappa3_2*kappa4_2*kappa8)/3344302080.0L + (kappa3_3*kappa5*kappa8)/6270566400.0L + (kappa3_3*kappa4*kappa9)/11287019520.0L + (kappa3_4*kappa10)/112870195200.0L;
+                                                    coeffs[23] += kappa4_6/137594142720.0L + (kappa3*kappa4_4*kappa5)/5733089280.0L + (kappa3_2*kappa4_2*kappa5_2)/2388787200.0L + (kappa3_3*kappa5_3)/13436928000.0L + (kappa3_2*kappa4_3*kappa6)/4299816960.0L + (kappa3_3*kappa4*kappa5*kappa6)/2687385600.0L + (kappa3_4*kappa6_2)/32248627200.0L + (kappa3_3*kappa4_2*kappa7)/7524679680.0L + (kappa3_4*kappa5*kappa7)/18811699200.0L + (kappa3_4*kappa4*kappa8)/30098718720.0L + (kappa3_5*kappa9)/338610585600.0L;
+                                                    coeffs[25] += (kappa3_2*kappa4_5)/68797071360.0L + (kappa3_3*kappa4_3*kappa5)/12899450880.0L + (kappa3_4*kappa4*kappa5_2)/21499084800.0L + (kappa3_4*kappa4_2*kappa6)/25798901760.0L + (kappa3_5*kappa5*kappa6)/80621568000.0L + (kappa3_5*kappa4*kappa7)/112870195200.0L + (kappa3_6*kappa8)/1354442342400.0L;
+                                                    coeffs[27] += (kappa3_4*kappa4_4)/247669456896 + (kappa3_5*kappa4_2*kappa5)/128994508800.0L + (kappa3_6*kappa5_2)/967458816000.0L + (kappa3_6*kappa4*kappa6)/580475289600.0L + (kappa3_7*kappa7)/7110822297600.0L;
+                                                    coeffs[29] += (kappa3_6*kappa4_3)/2786281390080.0L + (kappa3_7*kappa4*kappa5)/4063327027200.0L + (kappa3_8*kappa6)/48759924326400.0L;
+                                                    coeffs[31] += (kappa3_8*kappa4_2)/78015878922240.0L + (kappa3_9*kappa5)/438839318937600.0L;
+                                                    coeffs[33] += (kappa3_10*kappa4)/5266071827251200.0L;
+                                                    coeffs[35] += kappa3_12/1042682221795737600.0L;
 
                                                     if (order_ >= 13U)
-                                                        assert(0);
+                                                    {
+                                                        *maxdeg = 38U;
+
+                                                        const double kappa15 = cumulants_[14];
+                                                        const double kappa4_7 = kappa4_3*kappa4_4;
+                                                        const double kappa3_13 = kappa3_6*kappa3_7;
+                                                        const double kappa3_14 = kappa3_7*kappa3_7;
+                                                        const double kappa9_2 = kappa9*kappa9;
+                                                        coeffs[14] += kappa15/1307674368000.0L;
+                                                        coeffs[16] += (kappa8*kappa9)/14631321600.0L + (kappa7*kappa10)/18289152000.0L + (kappa6*kappa11)/28740096000.0L + (kappa5*kappa12)/57480192000.0L + (kappa4*kappa13)/149448499200.0L + (kappa3*kappa14)/523069747200.0L;
+                                                        coeffs[18] += (kappa6_2*kappa7)/5225472000.0L + (kappa5*kappa7_2)/6096384000.0L + (kappa5*kappa6*kappa8)/3483648000.0L + (kappa4*kappa7*kappa8)/4877107200.0L + (kappa3*kappa8_2)/19508428800.0L + (kappa5_2*kappa9)/10450944000.0L + (kappa4*kappa6*kappa9)/6270566400.0L + (kappa3*kappa7*kappa9)/10973491200.0L + (kappa4*kappa5*kappa10)/10450944000.0L + (kappa3*kappa6*kappa10)/15676416000.0L + (kappa4_2*kappa11)/45984153600.0L + (kappa3*kappa5*kappa11)/28740096000.0L + (kappa3*kappa4*kappa12)/68976230400.0L + (kappa3_2*kappa13)/448345497600.0L;
+                                                        coeffs[20] += (kappa5_3*kappa6)/7464960000.0L + (kappa4*kappa5*kappa6_2)/2985984000.0L + (kappa3*kappa6_3)/13436928000.0L + (kappa4*kappa5_2*kappa7)/3483648000.0L + (kappa4_2*kappa6*kappa7)/4180377600.0L + (kappa3*kappa5*kappa6*kappa7)/2612736000.0L + (kappa3*kappa4*kappa7_2)/7315660800.0L + (kappa4_2*kappa5*kappa8)/5573836800.0L + (kappa3*kappa5_2*kappa8)/6967296000.0L + (kappa3*kappa4*kappa6*kappa8)/4180377600.0L + (kappa3_2*kappa7*kappa8)/14631321600.0L + (kappa4_3*kappa9)/30098718720.0L + (kappa3*kappa4*kappa5*kappa9)/6270566400.0L + (kappa3_2*kappa6*kappa9)/18811699200.0L + (kappa3*kappa4_2*kappa10)/25082265600.0L + (kappa3_2*kappa5*kappa10)/31352832000.0L + (kappa3_2*kappa4*kappa11)/68976230400.0L + (kappa3_3*kappa12)/620786073600.0L;
+                                                        coeffs[22] += (kappa4_2*kappa5_3)/11943936000.0L + (kappa3*kappa5_4)/29859840000.0L + (kappa4_3*kappa5*kappa6)/7166361600.0L + (kappa3*kappa4*kappa5_2*kappa6)/2985984000.0L + (kappa3*kappa4_2*kappa6_2)/7166361600.0L + (kappa3_2*kappa5*kappa6_2)/8957952000.0L + (kappa4_4*kappa7)/40131624960.0L + (kappa3*kappa4_2*kappa5*kappa7)/4180377600.0L + (kappa3_2*kappa5_2*kappa7)/10450944000.0L + (kappa3_2*kappa4*kappa6*kappa7)/6270566400.0L + (kappa3_3*kappa7_2)/65840947200.0L + (kappa3*kappa4_3*kappa8)/20065812480.0L + (kappa3_2*kappa4*kappa5*kappa8)/8360755200.0L + (kappa3_3*kappa6*kappa8)/37623398400.0L + (kappa3_2*kappa4_2*kappa9)/30098718720.0L + (kappa3_3*kappa5*kappa9)/56435097600.0L + (kappa3_3*kappa4*kappa10)/112870195200.0L + (kappa3_4*kappa11)/1241572147200.0L;
+                                                        coeffs[24] += (kappa4_5*kappa5)/114661785600.0L + (kappa3*kappa4_3*kappa5_2)/14332723200.0L + (kappa3_2*kappa4*kappa5_3)/17915904000.0L + (kappa3*kappa4_4*kappa6)/34398535680.0L + (kappa3_2*kappa4_2*kappa5*kappa6)/7166361600.0L + (kappa3_3*kappa5_2*kappa6)/26873856000.0L + (kappa3_3*kappa4*kappa6_2)/32248627200.0L + (kappa3_2*kappa4_3*kappa7)/30098718720.0L + (kappa3_3*kappa4*kappa5*kappa7)/18811699200.0L + (kappa3_4*kappa6*kappa7)/112870195200.0L + (kappa3_3*kappa4_2*kappa8)/60197437440.0L + (kappa3_4*kappa5*kappa8)/150493593600.0L + (kappa3_4*kappa4*kappa9)/270888468480.0L + (kappa3_5*kappa10)/3386105856000.0L;
+                                                        coeffs[26] += (kappa3*kappa4_6)/825564856320.0L + (kappa3_2*kappa4_4*kappa5)/68797071360.0L + (kappa3_3*kappa4_2*kappa5_2)/42998169600.0L + (kappa3_4*kappa5_3)/322486272000.0L + (kappa3_3*kappa4_3*kappa6)/77396705280.0L + (kappa3_4*kappa4*kappa5*kappa6)/64497254400.0L + (kappa3_5*kappa6_2)/967458816000.0L + (kappa3_4*kappa4_2*kappa7)/180592312320.0L + (kappa3_5*kappa5*kappa7)/564350976000.0L + (kappa3_5*kappa4*kappa8)/902961561600.0L + (kappa3_6*kappa9)/12189981081600.0L;
+                                                        coeffs[28] += (kappa3_3*kappa4_5)/1238347284480.0L + (kappa3_4*kappa4_3*kappa5)/309586821120.0L + (kappa3_5*kappa4*kappa5_2)/644972544000.0L + (kappa3_5*kappa4_2*kappa6)/773967052800.0L + (kappa3_6*kappa5*kappa6)/2902376448000.0L + (kappa3_6*kappa4*kappa7)/4063327027200.0L + (kappa3_7*kappa8)/56886578380800.0L;
+                                                        coeffs[30] += (kappa3_5*kappa4_4)/7430083706880.0L + (kappa3_6*kappa4_2*kappa5)/4643802316800.0L + (kappa3_7*kappa5_2)/40633270272000.0L + (kappa3_7*kappa4*kappa6)/24379962163200.0L + (kappa3_8*kappa7)/341319470284800.0L;
+                                                        coeffs[32] += (kappa3_7*kappa4_3)/117023818383360.0L + (kappa3_8*kappa4*kappa5)/195039697305600.0L + (kappa3_9*kappa6)/2633035913625600.0L;
+                                                        coeffs[34] += (kappa3_9*kappa4_2)/4212857461800960.0L + (kappa3_10*kappa5)/26330359136256000.0L;
+                                                        coeffs[36] += (kappa3_11*kappa4)/347560740598579200.0L;
+                                                        coeffs[38] += kappa3_13/81329213300067532800.0L;
+
+                                                        if (order_ >= 14U)
+                                                        {
+                                                            *maxdeg = 41U;
+
+                                                            const double kappa16 = cumulants_[15];
+                                                            coeffs[15] += kappa16/20922789888000.0L;
+                                                            coeffs[17] += kappa9_2/263363788800.0L + (kappa8*kappa10)/146313216000.0L + (kappa7*kappa11)/201180672000.0L + (kappa6*kappa12)/344881152000.0L + (kappa5*kappa13)/747242496000.0L + (kappa4*kappa14)/2092278988800.0L + (kappa3*kappa15)/7846046208000.0L;
+                                                            coeffs[19] += (kappa6*kappa7_2)/36578304000.0L + (kappa6_2*kappa8)/41803776000.0L + (kappa5*kappa7*kappa8)/24385536000.0L + (kappa4*kappa8_2)/78033715200.0L + (kappa5*kappa6*kappa9)/31352832000.0L + (kappa4*kappa7*kappa9)/43893964800.0L + (kappa3*kappa8*kappa9)/87787929600.0L + (kappa5_2*kappa10)/104509440000.0L + (kappa4*kappa6*kappa10)/62705664000.0L + (kappa3*kappa7*kappa10)/109734912000.0L + (kappa4*kappa5*kappa11)/114960384000.0L + (kappa3*kappa6*kappa11)/172440576000.0L + (kappa4_2*kappa12)/551809843200.0L + (kappa3*kappa5*kappa12)/344881152000.0L + (kappa3*kappa4*kappa13)/896690995200.0L + (kappa3_2*kappa14)/6276836966400.0L;
+                                                            coeffs[21] += (kappa5_2*kappa6_2)/29859840000.0L + (kappa4*kappa6_3)/53747712000.0L + (kappa5_3*kappa7)/52254720000.0L + (kappa4*kappa5*kappa6*kappa7)/10450944000.0L + (kappa3*kappa6_2*kappa7)/31352832000.0L + (kappa4_2*kappa7_2)/58525286400.0L + (kappa3*kappa5*kappa7_2)/36578304000.0L + (kappa4*kappa5_2*kappa8)/27869184000.0L + (kappa4_2*kappa6*kappa8)/33443020800.0L + (kappa3*kappa5*kappa6*kappa8)/20901888000.0L + (kappa3*kappa4*kappa7*kappa8)/29262643200.0L + (kappa3_2*kappa8_2)/234101145600.0L + (kappa4_2*kappa5*kappa9)/50164531200.0L + (kappa3*kappa5_2*kappa9)/62705664000.0L + (kappa3*kappa4*kappa6*kappa9)/37623398400.0L + (kappa3_2*kappa7*kappa9)/131681894400.0L + (kappa4_3*kappa10)/300987187200.0L + (kappa3*kappa4*kappa5*kappa10)/62705664000.0L + (kappa3_2*kappa6*kappa10)/188116992000.0L + (kappa3*kappa4_2*kappa11)/275904921600.0L + (kappa3_2*kappa5*kappa11)/344881152000.0L + (kappa3_2*kappa4*kappa12)/827714764800.0L + (kappa3_3*kappa13)/8070218956800.0L;
+                                                            coeffs[23] += (kappa4*kappa5_4)/119439360000.0L + (kappa4_2*kappa5_2*kappa6)/23887872000.0L + (kappa3*kappa5_3*kappa6)/44789760000.0L + (kappa4_3*kappa6_2)/85996339200.0L + (kappa3*kappa4*kappa5*kappa6_2)/17915904000.0L + (kappa3_2*kappa6_3)/161243136000.0L + (kappa4_3*kappa5*kappa7)/50164531200.0L + (kappa3*kappa4*kappa5_2*kappa7)/20901888000.0L + (kappa3*kappa4_2*kappa6*kappa7)/25082265600.0L + (kappa3_2*kappa5*kappa6*kappa7)/31352832000.0L + (kappa3_2*kappa4*kappa7_2)/87787929600.0L + (kappa4_4*kappa8)/321052999680.0L + (kappa3*kappa4_2*kappa5*kappa8)/33443020800.0L + (kappa3_2*kappa5_2*kappa8)/83607552000.0L + (kappa3_2*kappa4*kappa6*kappa8)/50164531200.0L + (kappa3_3*kappa7*kappa8)/263363788800.0L + (kappa3*kappa4_3*kappa9)/180592312320.0L + (kappa3_2*kappa4*kappa5*kappa9)/75246796800.0L + (kappa3_3*kappa6*kappa9)/338610585600.0L + (kappa3_2*kappa4_2*kappa10)/300987187200.0L + (kappa3_3*kappa5*kappa10)/564350976000.0L + (kappa3_3*kappa4*kappa11)/1241572147200.0L + (kappa3_4*kappa12)/14898865766400.0L;
+                                                            coeffs[25] += (kappa4_4*kappa5_2)/229323571200.0L + (kappa3*kappa4_2*kappa5_3)/71663616000.0L + (kappa3_2*kappa5_4)/358318080000.0L + (kappa4_5*kappa6)/687970713600.0L + (kappa3*kappa4_3*kappa5*kappa6)/42998169600.0L + (kappa3_2*kappa4*kappa5_2*kappa6)/35831808000.0L + (kappa3_2*kappa4_2*kappa6_2)/85996339200.0L + (kappa3_3*kappa5*kappa6_2)/161243136000.0L + (kappa3*kappa4_4*kappa7)/240789749760.0L + (kappa3_2*kappa4_2*kappa5*kappa7)/50164531200.0L + (kappa3_3*kappa5_2*kappa7)/188116992000.0L + (kappa3_3*kappa4*kappa6*kappa7)/112870195200.0L + (kappa3_4*kappa7_2)/1580182732800.0L + (kappa3_2*kappa4_3*kappa8)/240789749760.0L + (kappa3_3*kappa4*kappa5*kappa8)/150493593600.0L + (kappa3_4*kappa6*kappa8)/902961561600.0L + (kappa3_3*kappa4_2*kappa9)/541776936960.0L + (kappa3_4*kappa5*kappa9)/1354442342400.0L + (kappa3_4*kappa4*kappa10)/2708884684800.0L + (kappa3_5*kappa11)/37247164416000.0L;
+                                                            coeffs[27] += kappa4_7/23115815976960.0L + (kappa3*kappa4_5*kappa5)/687970713600.0L + (kappa3_2*kappa4_3*kappa5_2)/171992678400.0L + (kappa3_3*kappa4*kappa5_3)/322486272000.0L + (kappa3_2*kappa4_4*kappa6)/412782428160.0L + (kappa3_3*kappa4_2*kappa5*kappa6)/128994508800.0L + (kappa3_4*kappa5_2*kappa6)/644972544000.0L + (kappa3_4*kappa4*kappa6_2)/773967052800.0L + (kappa3_3*kappa4_3*kappa7)/541776936960.0L + (kappa3_4*kappa4*kappa5*kappa7)/451480780800.0L + (kappa3_5*kappa6*kappa7)/3386105856000.0L + (kappa3_4*kappa4_2*kappa8)/1444738498560.0L + (kappa3_5*kappa5*kappa8)/4514807808000.0L + (kappa3_5*kappa4*kappa9)/8126654054400.0L + (kappa3_6*kappa10)/121899810816000.0L;
+                                                            coeffs[29] += (kappa3_2*kappa4_6)/9906778275840.0L + (kappa3_3*kappa4_4*kappa5)/1238347284480.0L + (kappa3_4*kappa4_2*kappa5_2)/1031956070400.0L + (kappa3_5*kappa5_3)/9674588160000.0L + (kappa3_4*kappa4_3*kappa6)/1857520926720.0L + (kappa3_5*kappa4*kappa5*kappa6)/1934917632000.0L + (kappa3_6*kappa6_2)/34828517376000.0L + (kappa3_5*kappa4_2*kappa7)/5417769369600.0L + (kappa3_6*kappa5*kappa7)/20316635136000.0L + (kappa3_6*kappa4*kappa8)/32506616217600.0L + (kappa3_7*kappa9)/511979205427200.0L;
+                                                            coeffs[31] += (kappa3_4*kappa4_5)/29720334827520.0L + (kappa3_5*kappa4_3*kappa5)/9287604633600.0L + (kappa3_6*kappa4*kappa5_2)/23219011584000.0L + (kappa3_6*kappa4_2*kappa6)/27862813900800.0L + (kappa3_7*kappa5*kappa6)/121899810816000.0L + (kappa3_7*kappa4*kappa7)/170659735142400.0L + (kappa3_8*kappa8)/2730555762278400.0L;
+                                                            coeffs[33] += (kappa3_6*kappa4_4)/267483013447680.0L + (kappa3_7*kappa4_2*kappa5)/195039697305600.0L + (kappa3_8*kappa5_2)/1950396973056000.0L + (kappa3_8*kappa4*kappa6)/1170238183833600.0L + (kappa3_9*kappa7)/18431251395379200.0L;
+                                                            coeffs[35] += (kappa3_8*kappa4_3)/5617143282401280.0L + (kappa3_9*kappa4*kappa5)/10532143654502400.0L + (kappa3_10*kappa6)/157982154817536000.0L;
+                                                            coeffs[37] += (kappa3_10*kappa4_2)/252771447708057600.0L + (kappa3_11*kappa5)/1737803702992896000.0L;
+                                                            coeffs[39] += (kappa3_12*kappa4)/25024373323097702400.0L;
+                                                            coeffs[41] += kappa3_14/6831653917205672755200.0L;
+
+                                                            if (order_ > EDGEWORTH_MAXORDER)
+                                                                assert(0);
+                                                        }
+                                                    }
                                                 }
                                             }
                                         }
                                     }
                                 }
                             }
                         }
                     }
                 }
             }
         }
 
         if (order_ > 1U)
         {
             const double sigma = sqrt(cumulants_[1]);
             if (sigma != 1.0)
                 for (unsigned i=2; i<=*maxdeg; ++i)
                     coeffs[i] /= powl(sigma, i);
         }
     }
 
     void EdgeworthSeries1D::hermiteCoeffsNormal(long double* coeffs, unsigned* maxdeg) const
     {
         if (m_ == EDGEWORTH_CLASSICAL)
         {
             hermiteCoeffsClassical(coeffs, maxdeg);
             return;
         }
 
         for (unsigned i=0; i<MAXCOEFFS; ++i)
             coeffs[i] = 0.0L;
         *maxdeg = 0U;
 
         assert(order_);
         const unsigned nk = cumulants_.size();
         assert(nk > 2U);
 
         const double k1 = cumulants_[0];
         const double k2 = cumulants_[1];
         const double k3 = cumulants_[2];
 
         *maxdeg = 2U;
         assert(MAXCOEFFS > *maxdeg);
         coeffs[0] = k1;
         coeffs[2] = k3/6.0;
 
         if (order_ > 1U)
         {
             assert(nk > 3U);
             const double k4 = cumulants_[3];
             const double k2m1 = k2 - 1.0;
 
             *maxdeg = 5;
             assert(MAXCOEFFS > *maxdeg);
             coeffs[1] += (k1*k1 + k2m1)/2.0;
             coeffs[3] += (4*k1*k3 + k4)/24.0;
             coeffs[5] += k3*k3/72.0;
 
             if (order_ > 2U)
             {
                 assert(nk > 4U);
                 const double k5 = cumulants_[4];
 
                 *maxdeg = 8;
                 assert(MAXCOEFFS > *maxdeg);
                 coeffs[2] += k1*(k1*k1 + 3*k2m1)/6;
                 coeffs[4] += (10*k1*k1*k3 + 10*k2m1*k3 + 5*k1*k4 + k5)/120;
                 coeffs[6] += (k3*(2*k1*k3 + k4))/144;
                 coeffs[8] += k3*k3*k3/1296;
 
                 if (order_ > 3U)
                 {
                     assert(nk > 5U);
                     const double k6 = cumulants_[5];
                     const double k1sq = k1*k1;
                     const double k3sq = k3*k3;
 
                     *maxdeg = 11;
                     assert(MAXCOEFFS > *maxdeg);
                     coeffs[3] += (k1sq*k1sq + 6*k1sq*k2m1 + 3*k2m1*k2m1)/24;
                     coeffs[5] += (20*k1*k1sq*k3 + 60*k1*k2m1*k3 + 15*k1sq*k4 + 15*k2m1*k4 + 6*k1*k5 + k6)/720;
                     coeffs[7] += (40*k1sq*k3sq + 40*k2m1*k3sq + 40*k1*k3*k4 + 5*k4*k4 + 8*k3*k5)/5760;
                     coeffs[9] += (k3sq*(4*k1*k3 + 3*k4))/5184;
                     coeffs[11] += k3sq*k3sq/31104;
 
                     if (order_ > 4U)
                         assert(0);
                 }
             }
         }
     }
 
     double EdgeworthSeriesMomentFcn::operator()(const double& x) const
     {
         const double fact = series_.densityFactor(x);
         switch (degree_)
         {
         case 0:
             return fact;
 
         case 1:
             return x*fact;
 
         default:
             return pow(x - mean_, degree_)*fact;
         }
     }
 }
Index: trunk/NEWS
===================================================================
--- trunk/NEWS	(revision 609)
+++ trunk/NEWS	(revision 610)
@@ -1,737 +1,737 @@
 Version 4.12.0 - development
 
-* Increased the maximum order of "classical" Edgeworth expansions to 12.
+* Increased the maximum order of "classical" Edgeworth expansions to 14.
 
 * Increased the maximum order of conversions between 1-d central moments
   and cumulants to 20.
 
 
 Version 4.11.0 - July 22 2019, by I. Volobouev
 
 * Added support for cumulant calculations for various Wald statistics
   in the Poisson process model.
 
 * Added functions convertCumulantsToCentralMoments and
   convertCentralMomentsToCumulants (header file cumulantConversion.hh).
 
 * Added function "cumulantUncertainties".
 
 
 Version 4.10.0 - July 11 2019, by I. Volobouev
 
 * Added SemiInfGaussianQuadrature class.
 
 * Added functions arrayMoment, arrayMoments, and arrayCentralMoments.
 
 * Added enum EdgeworthSeriesMethod and class EdgeworthSeries1D.
 
 * Added DeltaMixture1D class.
 
 * Added enum LikelihoodStatisticType.
 
 * Added functions "mixtureModelCumulants" and "poissonProcessCumulants"
   in the header likelihoodStatisticCumulants.hh.
 
 
 Version 4.9.0 - Dec 18 2018, by I. Volobouev
 
 * Added "integratePoly" and "jointIntegral" methods to the
   AbsClassicalOrthoPoly1D class.
 
 * Added utility functions "truncatedInverseSqrt" and "matrixIndexPairs".
 
 * Added a number of functions to "orthoPoly1DVProducts.hh".
 
 * Added classes ChebyshevOrthoPoly1st and ChebyshevOrthoPoly2nd inheriting
   from AbsClassicalOrthoPoly1D.
 
 * Added class HermiteProbOrthoPoly1D.
 
 * Added FejerQuadrature class.
 
 * Added classe IsoscelesTriangle1D and Logistic1D.
 
 * Added classes AbsKDE1DKernel and KDE1DHOSymbetaKernel.
 
 * Added static function "optimalDegreeHart" to OSDE1D class.
 
 
 Version 4.8.0 - Jul 9 2018, by I. Volobouev
 
 * Added ShiftedLegendreOrthoPoly1D class.
 
 * Added Lanczos method to generate recurrence coefficients for the
   ContOrthoPoly1D class.
 
 * Added npstat/stat/orthoPoly1DVProducts.hh file with various utilities
   for statistical analyis of chaos polynomials.
 
 
 Version 4.7.0 - Jan 13 2018, by I. Volobouev
 
 * Added "UniPareto1D" distribution (uniform with Pareto tail to the right).
 
 * More coordinate/weight cases for the GaussLegendreQuadrature class.
 
 * Added ContOrthoPoly1D class -- continuous orthogonal polynomials with
   discrete measure.
 
 * Added functions "linearLeastSquares" and "tdSymEigen" to the Matrix class.
 
 * Added OSDE1D class.
 
 * Added classes LocationScaleFamily1D and SinhAsinhTransform1D.
 
 * Added new functors (CdfFunctor1D, etc) as AbsDistribution1D helpers.
 
 * Small fix in StatAccumulatorArr.cc.
 
 
 Version 4.6.0 - Jan 23 2017, by I. Volobouev
 
 * Updated 1-d LOrPE cross validation code (classes AbsBandwidthCV,
   BandwidthCVLeastSquares1D, BandwidthCVPseudoLogli1D) for use with
   weighted samples in the case the sample itself is available at the
   point cross validation is run.
 
 
 Version 4.5.0 - Aug 01 2016, by I. Volobouev
 
 * A small fix in OrthoPolyND.icc (switched from cycling over unsigned to
   unsigned long in the scalar product function).
 
 * Implemented Gauss-Hermite quadrature with Gaussian density weight.
 
 * Changed the meaning of Quadratic1D and LogQuadratic1D parameters to be
   consistent with Legendre polynomial coefficients on [-1, 1] (new
   parameters are 1/2 of old).
 
 * Added class MinuitUnbinnedFitFcn1D (to interfaces).
 
 * Added function "findRootNewtonRaphson".
 
 * Added "statUncertainties" header with various functions.
 
 
 Version 4.4.0 - May 9 2016, by I. Volobouev
 
 * Added "timestamp" function.
 
 * Improved implementation of BinnedDensity1D::unscaledQuantile function.
   Certain problems caused by round-off errors are now fixed.
 
 * Added the capability to use the single closest parameter cells (thus
   disabling actual interpolation between parameter values, for speed)
   to "GridInterpolatedDistribution".
 
 
 Version 4.3.0 - March 19 2016, by I. Volobouev
 
 * Put additional assert statements inside OrthoPolyND.icc to prevent
   phony "array subscript is above array bounds" messages in g++ 4.9.2.
 
 * Improved CmdLine.hh.
 
 * Additional methods in CopulaInterpolationND and GridInterpolatedDistribution.
 
 * Added function "volumeDensityFromBinnedRadial".
 
 * Added convenience method "globalFilter" to the OrthoPoly1D class.
 
 * Initiated the transition of the Python API from Python 2 to Python 3.
 
 
 Version 4.2.0 - October 29 2015, by I. Volobouev
 
 * Added interpolation methods for the marginals to classes
   "CopulaInterpolationND" and "GridInterpolatedDistribution".
 
 * Removed assert on underflow in the "igamc" function. Now
   in such cases this function simply returns 0.0.
 
 
 Version 4.1.0 - July 27 2015, by I. Volobouev
 
 * Made a few additional methods virtual in AbsNtuple.
 
 * Declared methods "columnIndices" of AbsNtuple const (as they should be).
 
 * Added function "weightedCopulaHisto" to build copulas for samples of
   weighted points.
 
 * Added function "weightedVariableBandwidthSmooth1D" to use variable
   bandwidth smoothing with weighted histograms.
 
 * Added "AbsWeightedMarginalSmoother" interface class for smoothing samples
   of weighted points. Modified classes ConstantBandwidthSmoother1D,
   JohnsonKDESmoother, and LOrPEMarginalSmoother to support this interface.
 
 * Added class "VariableBandwidthSmoother1D" which implements both
   AbsMarginalSmoother and AbsWeightedMarginalSmoother interfaces.
 
 * Implemented cross-validation for weighted samples.
 
 * Implemented "buildInterpolatedCompositeDistroND" for generic construction
   of multivariate transfer functions.
 
 * Implemented "buildInterpolatedDistro1DNP" for generic construction
   of univariate transfer functions.
 
 
 Version 4.0.1 - June 17 2015, by I. Volobouev
 
 * Added "dump_qmc" example executable.
 
 
 Version 4.0.0 - June 10 2015, by I. Volobouev
 
 * Complete overhaul of 1-d filter-building code. Addition of new boundary
   methods is a lot easier now. The user API for choosing a LOrPE boundary
   method is encapsulated in the new "BoundaryHandling" class.
 
 * Implemented a number of new filter builders with different boundary
   treatments.
 
 * Updated the "LocalPolyFilter1D" class so that it holds the local
   bandwidth factors derived by the filter builders.
 
 
 Version 3.8.0 - June 1 2015, by I. Volobouev
 
 * Implemented class ConstSqFilter1DBuilder (declared in the header file
   npstat/stat/Filter1DBuilders.hh). The "BoundaryMethod" enum has been
   extended accordingly. Other files using this enum have been updated.
 
 * Implemented class FoldingSqFilter1DBuilder. Similar to ConstSqFilter1DBuilder
   but it also folds the kernel in addition to stretching it.
 
 * Added virtual destructors to a number of classes.
 
 * Added "externalMemArrayND" with various signatures to allow the use of
   ArrayND with memory not managed by ArrayND.
 
 * Added move constructors and move assignment operators to ArrayND and
   Matrix classes.
 
 
 Version 3.7.0 - May 10 2015, by I. Volobouev
 
 * Better numerical derivative calculation in InterpolatedDistribution1D.cc.
 
 * Added class "LocalMultiFilter1D" for fast generation of filters which
   correspond to each orthogonal polynomial separately.
 
 * Added a function calculating the area of n-dimensional sphere.
 
 * Added I/O capabilities to the RadialProfileND class.
 
 * Added class "LogRatioTransform1D".
 
 * Added utility function "multiFill1DHistoWithCDFWeights" (header file
   histoUtils.hh).
 
 * Avoiding underflow of the incomplete gamma in "convertToSphericalRandom".
 
 
 Version 3.6.0 - April 6 2015, by I. Volobouev
 
 * Fixed Marsaglia's code calculating the Anderson-Darling statistics
   (it was breaking down for large values of z).
 
 * Added high-level driver function "simpleVariableBandwidthSmooth1D" to
   automatically build the pilot estimate for "variableBandwidthSmooth1D".
 
 * Swithched to log of sigma as Minuit parameter in "minuitFitJohnsonCurves"
   instead of sigma (linear sigma would sometimes break the fit when Minuit
   would come up with 0 or negative trial value for it).
 
 * Extended "MinuitDensityFitFcn1D" class so that it could be used to fit
   non-uniformly binned histograms.
 
 * Adapted "minuitFitJohnsonCurves" function so that it could be used with
   histograms templated upon DualHistoAxis.
 
 * Added functions "fillArrayCentersPreservingAreas" and
   "canFillArrayCentersPreservingAreas".
 
 * Implemented an interface to monotonous coordinate transformations with
   the intent of using them in constructing densities. Implemented a number
   of transforms.
 
 * Implemented class "TransformedDistribution1D".
 
 * Added class "VerticallyInterpolatedDistro1D1P".
 
 * Added utility function "fill1DHistoWithCDFWeights".
 
 
 Version 3.5.0 - February 21 2015, by I. Volobouev
 
 * Added "symPDEigenInv" method to the Matrix class.
 
 * Added "variableCount" method to unfolding bandwidth scanner classes.
 
 * Increased the tolerance parameters in JohnsonSu::initialize and in
   JohnsonFit constructor.
 
 * Bug fix in function "fillHistoFromText".
 
 
 Version 3.4.4 - January 13 2015, by I. Volobouev
 
 * Corrected handling of some "assert" statements so that the code compiles
   correctly with the -DNDEBUG option.
 
 
 Version 3.4.3 - January 5 2015, by I. Volobouev
 
 * Implemented class MirroredGauss1D.
 
 * Added method "getOracleData" to class UnfoldingBandwidthScanner1D.
 
 * Bug fix in FoldingFilter1DBuilder::makeOrthoPoly.
 
 
 Version 3.4.2 - December 15 2014, by I. Volobouev
 
 * Implemented InterpolatedDistro1D1P class.
 
 
 Version 3.4.1 - November 07 2014, by I. Volobouev
 
 * Implemented "divideTransforms" function for deconvolutions.
 
 * Implemented the Moyal distribution.
 
 * Added "fillHistoFromText" utility function.
 
 * Added "apply_lorpe_1d" example.
 
 
 Version 3.4.0 - October 01 2014, by I. Volobouev
 
 * Implemented Hadamard product and Hadamard ratio for matrices.
 
 * Bug fix in the "solve_linear_system" lapack interface function.
 
 
 Version 3.3.1 - August 08 2014, by I. Volobouev
 
 * Terminate iterative refinement of the unfolding error propagation
   matrix early in case the solution does not seem to improve.
 
 
 Version 3.3.0 - August 05 2014, by I. Volobouev
 
 * Added correction factors to the determination of the number of fitted
   parameters in various unfolding procedures.
 
 
 Version 3.2.0 - July 25 2014, by I. Volobouev
 
 * Added "gaussianResponseMatrix" function for non-uniform binning.
 
 * Added Pareto distribution.
 
 * Implemented EMS unfolding with sparse matrices.
 
 * Added methods "getObservedShape" and "getUnfoldedShape" to the AbsUnfoldND
   class.
 
 * Bug fix in the assignment operator of ProductDistributionND class. Made
   class ProductDistributionND persistent.
 
 * Bug fix in the error propagation for unfolding, in the code which takes
   into account the extra normalization constant.
 
 * Added "productResponseMatrix" function to assist in making sparse response
   matrices.
 
 * Bug fix in the factory constructor of the Cauchy1D class.
 
 * Completed implementation of the "RatioOfNormals" class.
 
 
 Version 3.1.0 - June 29 2014, by I. Volobouev
 
 * Improved (again) random number generator for the 1-d Gaussian distribution.
   Something about expectation values of normalized Hermite polynomials over
   random numbers made by this generator is still not fully understood. The
   standard deviation of these expectations decreases with the polynomial
   order (while it should stay constant). It is possible that the numbers
   of points used are simply insufficient to sample the tails correctly.
 
 * Implemented smoothed expectation-maximization (a.k.a. D'Agostini) unfolding
   for 1-d distributions in classes SmoothedEMUnfold1D and MultiscaleEMUnfold1D.
   In certain usage scenarios, MultiscaleEMUnfold1D can be more efficient than
   SmoothedEMUnfold1D.
 
 * Implemented smoothed expectation-maximization unfolding for multivariate
   distributions in a class SmoothedEMUnfoldND.
 
 * Added class "UnfoldingBandwidthScanner1D" to study 1-d unfolding behavior
   as a function of filter bandwidth.
 
 * Added class "UnfoldingBandwidthScannerND" to study multivariate unfolding
   behavior as a function of provided bandwidth values.
 
 * Added DummyLocalPolyFilter1D class useful when a filter is needed which
   does not smooth anything.
 
 * Added function "poissonLogLikelihood" (header file npstat/stat/arrayStats.hh).
 
 * Added function "pooledDiscreteTabulated1D" (header file
   npstat/stat/DiscreteDistributions1D.hh).
 
 * Implemented class UGaussConvolution1D (convolution of uniform distribution
   with a Gaussian).
 
 * Implemented gamma distribution (class Gamma1D).
 
 * Defined interface for comparing binned distributions, AbsBinnedComparison1D.
 
 * Implemented several comparison classes for comparing binned distributions:
   PearsonsChiSquared, BinnedKSTest1D, BinnedADTest1D. Class BinnedKSTest1D
   pulled in dependence on the "kstest" package.
 
 * Made classes LocalPolyFilter1D, LocalPolyFilterND, and SequentialPolyFilterND
   persistent.
 
 * Added code generating dense filter matrices to LocalPolyFilterND and
   SequentialPolyFilterND (as needed for unfolding).
 
 * Made class MemoizingSymbetaFilterProvider persistent.
 
 * Implemented function goldenSectionSearchOnAGrid (header file 
   npstat/nm/goldenSectionSearch.hh).
 
 * Implemented function parabolicExtremum (header npstat/nm/MathUtils.hh).
 
 * Added class DistributionMix1D (header npstat/stat/DistributionMix1D.hh).
 
 * Added interface to solving A*X = B, with matrices X and B, to the Matrix
   class (method "solveLinearSystems").
 
 * Added "reshape" methods to the ArrayND class.
 
 * Added "gaussianResponseMatrix" function.
 
 * Added a section on unfolding to the documentation.
 
 * Added "ems_unfold_1d" example program.
 
 
 Version 3.0.0 - March 14 2014, by I. Volobouev
 
 * Added interface to the LAPACK SVD routines.
 
 * Added function "lorpeMise1D" to calculate MISE for arbitrary distributions.
 
 * Added FoldingFilter1DBuilder class.
 
 * Changed interfaces for several high-level functions to use
   FoldingFilter1DBuilder. The major version number got bumped up.
 
 * Split DensityScan1D.hh away from AbsDistribution1D.hh.
 
 
 Version 2.7.0 - March 10 2014, by I. Volobouev
 
 * Added code to optimize operations with diagonal matrices.
 
 * Added discretizedDistance.hh file for simple L1 and L2 distance calculations
   with numeric arrays.
 
 * Added base class for future unfolding code.
 
 * The "reset" method of the Matrix class was renamed into "uninitialize"
   in order to be consistent with ArrayND.
 
 * Added function "multinomialCovariance1D".
 
 * Added "polyTimesWeight" method to the OrthoPoly1D class.
 
 * Added methods "TtimesThis" and "timesT" to the Matrix class. These
   methods are more efficient than transpose followed by multiplication.
 
 
 Version 2.6.0 - January 30 2014, by I. Volobouev
 
 * Added function "lorpeBackgroundCVDensity1D" which linearizes calculation
   of the cross validation likelihood in semiparametric fits. Argument
   "linearizeCrossValidation" was added to MinuitSemiparametricFitFcn1D
   constructor, "lorpeBackground1D" function, etc.
 
 * Added the ability to build filters with center point removed to classes
   WeightTableFilter1DBuilder and StretchingFilter1DBuilder. The function
   "symbetaLOrPEFilter1D" now has an appropriate switch.
 
 * Added "removeRowAndColumn" method to the Matrix class.
 
 * Added CircularBuffer class.
 
 * Added various plugin bandwidth functions which work with non-integer
   polynomial degrees.
 
 * Switched to the Legendre polynomial basis for calculating all 1-d
   orthogonal polynomials (instead of monomial basis).
 
 * Added MemoizingSymbetaFilterProvider class.
 
 * Added "operator+=" method to the MultivariateSumAccumulator class.
 
 * Simplified implementation of the PolyFilterCollection1D class.
   File PolyFilterCollection1D.icc is removed.
 
 * Added "RatioOfNormals" 1-d distribution function. Only the density is
   currently implemented but not the CDF.
 
 * Added ExpMapper1d class.
 
 
 Version 2.5.0 - October 15 2013, by I. Volobouev
 
 * Added "getFilterMatrix" method to the LocalPolyFilter1D class.
 
 * Added "genEigen" method to the Matrix class (for determination of
   eigenvalues and eigenvectors of general real matrices).
 
 * Refactored the LAPACK interface so that interface functions to floats
   are automatically generated from interface functions to doubles. See
   the comment at the end of the "lapack_interface.icc" file for the shell
   commands to do this.
 
 
 Version 2.4.0 - October 6 2013, by I. Volobouev
 
 * Added functions "lorpeBackground1D", "lorpeBgCVPseudoLogli1D", and
   "lorpeBgLogli1D".
 
 * Added minuit interface classes "MinuitLOrPEBgCVFcn1D" and
   "MinuitSemiparametricFitFcn1D".
 
 * Added "ScalableDensityConstructor1D" class for use with Minuit interface
   functions.
 
 * Added classes AbsSymbetaFilterProvider and SymbetaPolyCollection1D.
 
 
 Version 2.3.0 - October 1 2013, by I. Volobouev
 
 * Allowed point dimensionality to be larger than the histogram dimensionality
   in the "empiricalCopulaHisto" function.
 
 * Added "keepAllFilters" method to AbsFilter1DBuilder and all derived classes.
 
 * Implemented exclusion regions for WeightTableFilter1DBuilder and
   StretchingFilter1DBuilder.
 
 * "symbetaLOrPEFilter1D" function (in the header LocalPolyFilter1D.hh)
   is updated to take the exclusion mask argument.
 
 * Added "continuousDegreeTaper" function which can do something meaningful
   with the continuous LOrPE degree parameter.
 
 
 Version 2.2.0 - June 30 2013, by I. Volobouev
 
 * Added classes DiscreteBernsteinPoly1D and BernsteinFilter1DBuilder.
 
 * Added classes DiscreteBeta1D and BetaFilter1DBuilder.
 
 * Added BifurcatedGauss1D class to model Gaussian-like distributions with
   different sigmas on the left and right sides.
 
 * Added virtual destructors to the classes declared in the Filter1DBuilders.hh
   header.
 
 * Added a method to the Matrix template to calculate Frobenius norm.
 
 * Added methods to the Matrix template to calculate row and column sums.
 
 * Added "directSum" method to the Matrix template.
 
 * Added constructor from a subrange of another matrix to the Matrix template.
 
 * Added code to the LocalPolyFilter1D class that generates a doubly
   stochastic filter out of an arbitrary filter.
 
 * Added "npstat/nm/definiteIntegrals.hh" header and corresponding .cc file
   for various infrequently used integrals.
 
 * Added "betaKernelsBandwidth" function.
 
 
 Version 2.1.0 - June 20 2013, by I. Volobouev
 
 * Fixed couple problems which showed up in the robust regression code
   due to compiler update.
 
 * Fixed CensoredQuantileRegressionOnKDTree::process method (needed this->
   dereference for some member).
 
 
 Version 2.0.0 - June 15 2013, by I. Volobouev
 
 * Updated to use "Geners" version 1.3.0. A few interfaces were changed
   (API for the string archives was removed because Geners own string archive
   facilities are now adequate) so the major version number was bumped up.
 
 
 Version 1.6.0 - June 12 2013, by I. Volobouev
 
 * Updated some documentation.
 
 * Updated fitCompositeJohnson.icc to use simplified histogram constructors.
 
 * Bug fix in the "minuitLocalQuantileRegression1D" function.
 
 * Changed the "quantileBinFromCdf" function to use unsigned long type for
   array indices.
 
 * Added "weightedLocalQuantileRegression1D" function (namespace npsi) for
   local regression with single predictor on weighted points.
 
 
 Version 1.5.0 - May 23 2013, by I. Volobouev
 
 * Added interfaces to LAPACK routines DSYEVD, DSYEVR, and corresponding
   single precision versions.
 
 * Added the "symPSDefEffectiveRank" method to the Matrix class for
   calculating effective ranks of symmetric positive semidefinite matrices.
 
 * Added converting constructor and assignment operator to the Matrix class.
 
 * Run the Gram-Schmidt procedure twice when orthogonal polynomials are
   derived in order to improve orthogonality.
 
 
 Version 1.4.0 - May 20 2013, by I. Volobouev
 
 * Added the "append" method to the AbsNtuple class.
 
 
 Version 1.3.0 - May 10 2013, by I. Volobouev
 
 * Added the code for Hermite polynomial series.
 
 * Improved random number generator for the 1-d Gaussian distribution.
 
 * Added a framework for discrete 1-d distributions as well as two
   concrete distribution classes (Poisson1D, DiscreteTabulated1D).
 
 * Added functions "readCompressedStringArchiveExt" and
   "writeCompressedStringArchiveExt" which can read/write either compressed
   or uncompressed string archives, distinguished by file extension.
 
 
 Version 1.2.1 - March 22 2013, by I. Volobouev
 
 * Improved CmdLine.hh in the "examples/C++" directory.
 
 * Added class QuantileTable1D.
 
 * Added classes LeftCensoredDistribution and RightCensoredDistribution.
 
 
 Version 1.2.0 - March 13 2013, by I. Volobouev
 
 * Added convenience "fill" methods to work with the ntuples which have
   small number of columns (up to 10).
 
 * Fixed a bug in AbsRandomGenerator for univariate generators making
   multivariate points.
 
 * Added LOrPEMarginalSmoother class.
 
 
 Version 1.1.1 - March 11 2013, by I. Volobouev
 
 * Added utility function "symbetaLOrPEFilter1D" which creates 1-d LOrPE
   filters using kernels from the symmetric beta family (and the Gaussian).
 
 * Added high level driver function "lorpeSmooth1D".
 
 * Allowed variables with zero variances for calculation of correlation
   coefficients in "MultivariateSumsqAccumulator". Such variables will
   have zero correlation coefficients with all other variables.
 
 * Added rebinning constructor to the HistoND class.
 
 
 Version 1.1.0 - March 8 2013, by I. Volobouev
 
 * Changed NUHistoAxis::fltBinNumber method to produce correct results
   with interpolation degree 0. It is not yet obvious which method would
   work best for higher interpolation degrees.
 
 * Added functions for converting between StringArchive and python bytearray.
   They have been placed in a new header: wrap/stringArchiveToBinary.hh.
 
 * Added methods "exportMemSlice" and "importMemSlice" to ArrayND. These
   methods allow for filling array slices from unstructured memory buffers
   and for exporting array slices to such memory buffers.
 
 * Added "simpleColumnNames" function (header file AbsNtuple.hh) to generate
   trivial column names when ntuple column names are not important.
 
 * Added functions "neymanPearsonWindow1D" and "signalToBgMaximum1D".
   They are declared in a new header npstat/neymanPearsonWindow1D.hh.
 
 
 Version 1.0.5 - December 17 2012, by I. Volobouev
 
 * Flush string archives before writing them out in stringArchiveIO.cc.
 
 * Added class TruncatedDistribution1D.
 
 
 Version 1.0.4 - November 14 2012, by I. Volobouev
 
 * Added utilities for reading/writing Geners string archives to files.
 
 * Added BinSummary class.
 
 * Doxygen documentation improved. Every header file in stat, nm, rng,
   and interfaces now has a brief description.
 
 
 Version 1.0.3 - September 27 2012, by I. Volobouev
 
 * Fixed some bugs related to moving StorableMultivariateFunctor code
   from "nm" to "stat".
 
 
 Version 1.0.2 - August 6 2012, by I. Volobouev
 
 * Added converting copy constructor to the "LinInterpolatedTableND" class.
 
 * Added StorableMultivariateFunctor class (together with the corresponding
   reader class).
 
 * Added StorableInterpolationFunctor class which inherits from the above
   and can be used with interpolation tables.
 
 * Added StorableHistoNDFunctor class which inherits from
   StorableMultivariateFunctor and can be used to interpolate histogram bins.
 
 * Added "transpose" method to HistoND class.
 
 * Created DualAxis class.
 
 * Created DualHistoAxis class.
 
 * Added conversion functions between histogram and grid axes.
 
 * Added "mergeTwoHistos" function for smooth merging of two histograms.
 
 * Added "ProductSymmetricBetaNDCdf" functor to be used as weight in
   merging histograms.
 
 * Added CoordinateSelector class.
 
 
 Version 1.0.1 - June 29 2012, by I. Volobouev
 
 * Implemented class LinInterpolatedTableND with related supporting code.