Index: trunk/tests/test_MathUtils.cc
===================================================================
--- trunk/tests/test_MathUtils.cc	(revision 867)
+++ trunk/tests/test_MathUtils.cc	(revision 868)
@@ -1,38 +1,138 @@
 #include <cmath>
+#include <vector>
+#include <utility>
+
 #include "UnitTest++.h"
 
 #include "test_utils.hh"
+
 #include "npstat/nm/MathUtils.hh"
+#include "npstat/nm/SimpleFunctors.hh"
 
 using namespace npstat;
 using namespace std;
 
 namespace {
+    class PolyFunctor : public npstat::Functor1<double,double>
+    {
+    public:
+        inline PolyFunctor(const double* coeffs, const unsigned degree)
+            : coeffs_(coeffs, coeffs+(degree+1U)) {}
+
+        inline double operator()(const double& x) const
+            {return polySeriesSum(&coeffs_[0], coeffs_.size()-1U, x);}
+
+    private:
+        std::vector<double> coeffs_;
+    };
+
     TEST(parabolicExtremum)
     {
         const double eps = 1.0e-11;
         const unsigned ntry = 100;
         double x = 0.0, y = 0.0;
 
         for (unsigned i=0; i<ntry; ++i)
         {
             double a = 0.0;
             while (fabs(a) < 0.1)
                 a = test_rng()*2.0 - 1.0;
             const double x0 = test_rng()*2.0 - 1.0;
             const double y0 = test_rng()*2.0 - 1.0;
 
             const bool isMinimum = a > 0.0;
             const double x1 = test_rng()*4.0 - 2.0;
             const double x2 = test_rng()*4.0 - 2.0;
             const double x3 = test_rng()*4.0 - 2.0;
             const double y1 = a*(x1 - x0)*(x1 - x0) + y0;
             const double y2 = a*(x2 - x0)*(x2 - x0) + y0;
             const double y3 = a*(x3 - x0)*(x3 - x0) + y0;
 
             CHECK_EQUAL(isMinimum, parabolicExtremum(x1,y1,x2,y2,x3,y3,&x,&y));
             CHECK_CLOSE(x0, x, eps);
             CHECK_CLOSE(y0, y, eps);
         }
     }
+
+    TEST(polyIntegralCoeffs)
+    {
+        const unsigned deg = 4;
+
+        double coeffs[deg+1U];
+        double intCoeffs[deg+2U];
+
+        for (unsigned i=0; i<10U; ++i)
+        {
+            for (unsigned j=0; j<=deg; ++j)
+                coeffs[j] = 2.0*test_rng() - 1.0;
+
+            PolyFunctor fcn(coeffs, deg);
+            polyIntegralCoeffs(coeffs, deg, intCoeffs);
+
+            for (unsigned j=0; j<20U; ++j)
+            {
+                const double x = 2.0*test_rng() - 1.0;
+                long double d;
+                polyAndDeriv(intCoeffs, deg+1U, x, 0, &d);
+                const long double ref = polySeriesSum(coeffs, deg, x);
+                CHECK_CLOSE(ref, d, 1.0e-12);
+
+                const double xmin = 2.0*test_rng() - 1.0;
+                const double xmax = 2.0*test_rng() - 1.0;
+                const double integ1 = polySeriesSum(intCoeffs, deg+1U, xmax) -
+                                      polySeriesSum(intCoeffs, deg+1U, xmin);
+                const double integ2 = simpson_integral(fcn, xmin, xmax);
+                CHECK_CLOSE(integ1, integ2, 1.0e-6);
+            }
+        }
+    }
+
+    TEST(chebyshevIntegralCoeffs)
+    {
+        const unsigned deg = 4;
+
+        double coeffs[deg+1U];
+        double intCoeffs[deg+2U];
+
+        double chebcoeffs[deg+1U];
+        double chebIntCoeffs[deg+2U];
+
+        for (unsigned itry=0; itry<10U; ++itry)
+        {
+            for (unsigned j=0; j<=deg; ++j)
+                coeffs[j] = 2.0*test_rng() - 1.0;
+
+            double xmin = 4.0*test_rng() - 2.0;
+            double xmax = 4.0*test_rng() - 2.0;
+            // double xmin = -1.0;
+            // double xmax = 1.0;
+            // if (xmin > xmax)
+            //    std::swap(xmin, xmax);
+
+            const PolyFunctor fcn(coeffs, deg);
+            chebyshevSeriesCoeffs(fcn, xmin, xmax, deg, chebcoeffs);
+            chebyshevIntegralCoeffs(chebcoeffs, deg, xmin, xmax, chebIntCoeffs);
+            polyIntegralCoeffs(coeffs, deg, intCoeffs);
+
+            for (unsigned j=0; j<20U; ++j)
+            {
+                const double x = xmin + test_rng()*(xmax - xmin);
+                const double ref = fcn(x);
+                CHECK_CLOSE(ref, polySeriesSum(coeffs, deg, x), 1.0e-14);
+
+                const double cheb = chebyshevSeriesSum(
+                    chebcoeffs, deg, xmin, xmax, x);
+                CHECK_CLOSE(ref, cheb, 1.0e-14);
+
+                const double refIntegral = polySeriesSum(intCoeffs, deg+1U, x) -
+                                           polySeriesSum(intCoeffs, deg+1U, xmin);
+
+                const double i0 = chebyshevSeriesSum(chebIntCoeffs, deg+1U, xmin, xmax, xmin);
+                CHECK_CLOSE(0.0, i0, 1.0e-12);
+
+                const double integ = chebyshevSeriesSum(chebIntCoeffs, deg+1U, xmin, xmax, x);
+                CHECK_CLOSE(refIntegral, integ, 1.0e-12);
+            }
+        }
+    }
 }
Index: trunk/tests/test_ContOrthoPoly1D.cc
===================================================================
--- trunk/tests/test_ContOrthoPoly1D.cc	(revision 867)
+++ trunk/tests/test_ContOrthoPoly1D.cc	(revision 868)
@@ -1,406 +1,406 @@
 #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"
 
 #include "npstat/stat/Distributions1D.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);
             }
     }
 
 #ifdef USE_LAPACK_QUAD
     TEST(ContOrthoPoly1D_orthonormalization_quad_STIELTJES)
     {
         const lapack_double eps = FLT128_EPSILON*100.0;
 
         const OrthoPolyMethod method = OPOLY_STIELTJES;
 
         const unsigned npoints = 64;
         const unsigned maxdeg = 10;
         const unsigned ntries = 5;
         const double xmin = -1.0;
         const double xmax = 1.0;
         const double range = xmax - xmin;
 
         std::vector<double> points(2U*npoints);
         std::vector<ContOrthoPoly1D::MeasurePoint> measure;
         measure.reserve(npoints);
 
         MersenneTwister rng;
 
         for (unsigned itry=0; itry<ntries; ++itry)
         {
             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);
 
             for (unsigned i=0; i<=maxdeg; ++i)
                 for (unsigned j=0; j<=maxdeg; ++j)
                 {
                     const lapack_double d = poly.empiricalKroneckerDelta(i, j);
                     CHECK_CLOSE(static_cast<lapack_double>(kdelta(i, j)), d, eps);
                 }
         }
     }
 
     TEST(ContOrthoPoly1D_orthonormalization_quad_LANCZOS)
     {
         const lapack_double eps = FLT128_EPSILON*100.0;
 
         const OrthoPolyMethod method = OPOLY_LANCZOS;
             
         const unsigned npoints = 64;
         const unsigned maxdeg = 10;
         const unsigned ntries = 5;
         const double xmin = -1.0;
         const double xmax = 1.0;
         const double range = xmax - xmin;
 
         std::vector<double> points(2U*npoints);
         std::vector<ContOrthoPoly1D::MeasurePoint> measure;
         measure.reserve(npoints);
 
         MersenneTwister rng;
 
         for (unsigned itry=0; itry<ntries; ++itry)
         {
             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);
 
             for (unsigned i=0; i<=maxdeg; ++i)
                 for (unsigned j=0; j<=maxdeg; ++j)
                 {
                     const lapack_double d = poly.empiricalKroneckerDelta(i, j);
                     CHECK_CLOSE(static_cast<lapack_double>(kdelta(i, j)), d, eps);
                 }
         }
     }
 #endif
 
     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, 1000*eps);
-                CHECK_CLOSE(fvalue, f2, 50*eps);
+                CHECK_CLOSE(fvalue, f2, 200*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_integrateEWPolyProduct)
     {
         const double eps = 1.0e-12;
 
         MersenneTwister rng;
         const unsigned npoints = 64;
         const unsigned maxdeg = 6;
         const unsigned ntries = 5;
         const unsigned nInteg = 64;
         std::vector<double> points(npoints);
         std::vector<double> aves(maxdeg + 1);
         Beta1D beta(0.0, 1.0, 2.0, 2.0);
         GaussLegendreQuadrature glq(nInteg);
  
         for (unsigned itry=0; itry<ntries; ++itry)
         {
             rng.run(&points[0], npoints, npoints);
             ContOrthoPoly1D poly(maxdeg, points);
 
             DensityFunctor1D df(beta);
             const Matrix<double>& mat = poly.extWeightProductAverages(
                 df, glq, maxdeg, 0.0, 1.0);
             poly.extWeightAverages(df, glq, maxdeg, 0.0, 1.0, &aves[0]);
 
             for (unsigned i=0; i<=maxdeg; ++i)
                 for (unsigned j=0; j<=maxdeg; ++j)
                 {
                     const double integ = poly.integrateEWPolyProduct(
                         DensityFunctor1D(beta), i, j, 0.0, 1.0, nInteg);
                     CHECK_CLOSE(integ, mat[i][j], eps);
                 }
 
             for (unsigned j=0; j<=maxdeg; ++j)
                 CHECK_CLOSE(aves[j], mat[0][j], 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);
         OrthoPoly1DDeg poly3(jac, 3);
 
         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);
             CHECK_CLOSE(jac.poly(3, x), poly3(x), eps);
         }
 
         // Test "sampleAverages" method
         const unsigned nrandom = 100;
         std::vector<double> r(nrandom);
         for (unsigned i=0; i<nrandom; ++i)
             r[i] = -1.0 + 2.0*test_rng();
         std::vector<double> result(maxdeg + 1);
         poly.sampleAverages(&r[0], r.size(), &result[0], maxdeg);
 
         Matrix<long double> mat(maxdeg + 1, maxdeg + 1, 0);
         std::vector<long double> acc(maxdeg + 1, 0.0L);
         std::vector<long double> tmp(maxdeg + 1);
         for (unsigned i=0; i<nrandom; ++i)
         {
             poly.allPolys(r[i], maxdeg, &tmp[0]);
             for (unsigned deg=0; deg<=maxdeg; ++deg)
             {
                 acc[deg] += tmp[deg];
                 for (unsigned k=0; k<=maxdeg; ++k)
                 {
                     const double prod = tmp[deg]*tmp[k];
                     mat[deg][k] += prod;
                 }
             }
         }
         for (unsigned deg=0; deg<=maxdeg; ++deg)
         {
             acc[deg] /= nrandom;
             CHECK_CLOSE(static_cast<double>(acc[deg]), result[deg], eps);
         }
 
         // Test "pairwiseSampleAverages" method
         mat /= nrandom;
         Matrix<double> averages = poly.sampleProductAverages(
             &r[0], r.size(), maxdeg);
         Matrix<double> refmat(mat);
         CHECK((averages - refmat).maxAbsValue() < eps);
     }
 
     TEST(ContOrthoPoly1D_jacobi_2)
     {
         typedef ContOrthoPoly1D::PreciseQuadrature Quadrature;
 
         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 = Quadrature::minimalExactRule(2*maxdeg+alpha+beta);
         Quadrature 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/swig/nm/MathUtils.i
===================================================================
--- trunk/npstat/swig/nm/MathUtils.i	(revision 867)
+++ trunk/npstat/swig/nm/MathUtils.i	(revision 868)
@@ -1,50 +1,56 @@
 %include stddecls.i
 %include numpy.i
 
 %include nm/SimpleFunctors.i
 
 %{
 #include "npstat/nm/MathUtils.hh"
 %}
 
 %apply (double* IN_ARRAY1, int DIM1) {
     (const double *c, const unsigned degreep1)
 };
 
 %ignore npstat::polySeriesSum;
 %rename(polySeriesSum) npstat::polySeriesSum2;
 
 %ignore npstat::polyAndDeriv;
 %rename(polyAndDeriv) npstat::polyAndDeriv2;
 
+%ignore npstat::polyIntegralCoeffs;
+%rename(polyIntegralCoeffs) npstat::polyIntegralCoeffs2;
+
 %ignore npstat::legendreSeriesSum;
 %rename(legendreSeriesSum) npstat::legendreSeriesSum2;
 
 %ignore npstat::hermiteSeriesSumProb;
 %rename(hermiteSeriesSumProb) npstat::hermiteSeriesSumProb2;
 
 %ignore npstat::hermiteSeriesSumPhys;
 %rename(hermiteSeriesSumPhys) npstat::hermiteSeriesSumPhys2;
 
 %ignore npstat::gegenbauerSeriesSum;
 %rename(gegenbauerSeriesSum) npstat::gegenbauerSeriesSum2;
 
 %ignore npstat::chebyshevSeriesSum;
 %rename(chebyshevSeriesSum) npstat::chebyshevSeriesSum2;
 
+%ignore npstat::chebyshevIntegralCoeffs;
+%rename(chebyshevIntegralCoeffs) npstat::chebyshevIntegralCoeffs2;
+
 %ignore npstat::solveQuadratic;
 %rename(solveQuadratic) npstat::solveQuadratic2;
 
 %ignore npstat::solveCubic;
 %rename(solveCubic) npstat::solveCubic2;
 
 %ignore npstat::parabolicExtremum;
 /* parabolicExtremum2 is in npstat/wrap/parabolicExtremum.hh */
 
 %include "npstat/nm/MathUtils.hh"
 
 namespace npstat {
     %template(chebyshevSeriesCoeffs)  chebyshevSeriesCoeffs2<Functor1RefHelper<double,double> >;
 }
 
 %clear (const double *c, const unsigned degreep1);
Index: trunk/npstat/nm/MathUtils.icc
===================================================================
--- trunk/npstat/nm/MathUtils.icc	(revision 867)
+++ trunk/npstat/nm/MathUtils.icc	(revision 868)
@@ -1,174 +1,234 @@
 #include <cmath>
 #include <cassert>
 #include <vector>
 
 namespace npstat {
     template<typename Numeric>
     long double polySeriesSum(const Numeric *a,
                               const unsigned degree, const long double x)
     {
         // Dimension of "a" should be at least degree+1
         assert(a);
         long double sum = 0.0L;
         for (int deg=degree; deg>=0; --deg)
         {
             sum *= x;
             sum += a[deg];
         }
         return sum;
     }
 
     template<typename Numeric>
     void polyAndDeriv(const Numeric *a, unsigned deg, const long double x,
                       long double *value, long double *deriv)
     {
         // Dimension of "a" should be at least deg+1
         assert(a);
         long double sum = 0.0L, der = 0.0L;
         for (; deg>=1; --deg)
         {
             sum *= x;
             der *= x;
             sum += a[deg];
             der += deg*a[deg];
         }
         if (value)
             *value = sum*x + a[0];
         if (deriv)
             *deriv = der;
     }
 
     template<typename Numeric>
+    void polyIntegralCoeffs(const Numeric* coeffs, const unsigned degree,
+                            Numeric* integralCoeffs)
+    {
+        assert(coeffs);
+        assert(integralCoeffs);
+        assert(coeffs != integralCoeffs);
+
+        integralCoeffs[0] = static_cast<Numeric>(0);
+        for (unsigned i=0; i<=degree; ++i)
+        {
+            const unsigned ip1 = i + 1U;
+            integralCoeffs[ip1] = coeffs[i]/ip1;
+        }
+    }
+
+    template<typename Numeric>
     long double legendreSeriesSum(const Numeric *a,
                                   const unsigned degree, const long double x)
     {
         assert(a);
         long double result = a[0], pminus2 = 1.0L, pminus1 = x;
         if (degree)
         {
             result += a[1]*x;
             for (unsigned i=2; i<=degree; ++i)
             {
                 const long double p = ((2*i-1U)*x*pminus1 - (i-1U)*pminus2)/i;
                 result += p*a[i];
                 pminus2 = pminus1;
                 pminus1 = p;
             }
         }
         return result;
     }
 
     template<typename Numeric>
     long double hermiteSeriesSumProb(const Numeric *a,
                                      const unsigned degree, const long double x)
     {
         assert(a);
         long double result = a[0], pminus2 = 1.0L, pminus1 = x;
         if (degree)
         {
             result += a[1]*x;
             for (unsigned i=2; i<=degree; ++i)
             {
                 const long double p = x*pminus1 - (i-1U)*pminus2;
                 result += p*a[i];
                 pminus2 = pminus1;
                 pminus1 = p;
             }
         }
         return result;
     }
 
     template<typename Numeric>
     long double hermiteSeriesSumPhys(const Numeric *a,
                                      const unsigned degree, const long double x)
     {
         assert(a);
         long double result = a[0], pminus2 = 1.0L, pminus1 = 2.0L*x;
         if (degree)
         {
             result += a[1]*2.0L*x;
             for (unsigned i=2; i<=degree; ++i)
             {
                 const long double p = 2.0L*(x*pminus1 - (i-1U)*pminus2);
                 result += p*a[i];
                 pminus2 = pminus1;
                 pminus1 = p;
             }
         }
         return result;
     }
 
     template<typename Numeric>
     long double gegenbauerSeriesSum(const Numeric *a, const unsigned degree,
                                     const long double lambda,
                                     const long double x)
     {
         assert(a);
         long double result = a[0], pminus2 = 1.0L, pminus1 = 2.0L*lambda*x;
         if (degree)
         {
             result += a[1]*pminus1;
             const long double lm = 2.0L*(lambda - 1.0L);
             for (unsigned i=2; i<=degree; ++i)
             {
                 const long double p = ((2*i+lm)*x*pminus1 - (i+lm)*pminus2)/i;
                 result += p*a[i];
                 pminus2 = pminus1;
                 pminus1 = p;
             }
         }
         return result;
     }
 
     template<typename Numeric>
-    long double chebyshevSeriesSum(const Numeric *a,
-                                   const unsigned degree, const long double x)
+    long double chebyshevSeriesSum(const Numeric *a, const unsigned degree,
+                                   const long double x)
     {
         assert(a);
         const long double twox = 2.0L*x;
 
         // Clenshaw recursion
         long double rp2 = 0.0L, rp1 = 0.0L, r = 0.0L;
         for (unsigned k = degree; k > 0U; --k)
         {
             r = twox*rp1 - rp2 + a[k];
             rp2 = rp1;
             rp1 = r;
         }
         return x*rp1 - rp2 + a[0];
     }
 
+    template<typename Numeric>
+    long double chebyshevSeriesSum(const Numeric *a, const unsigned degree,
+                                   const long double xmin, const long double xmax,
+                                   const long double x)
+    {
+        assert(xmin != xmax);
+        const long double xtrans = 2.0L*(x - xmin)/(xmax - xmin) - 1.0L;
+        return chebyshevSeriesSum(a, degree, xtrans);
+    }
+
     template<typename Functor, typename Numeric>
     void chebyshevSeriesCoeffs(const Functor& f, const long double xmin,
                                const long double xmax,
                                const unsigned degree, Numeric *coeffs)
     {
         const long double Pi = 3.141592653589793238462643383L;
         assert(coeffs);
         assert(xmin != xmax);
 
         // Simple direct formula is used for calculations
         // (faster calculations for large degrees can be done
         // by the cosine transform).
         const unsigned N = degree + 1U;
         const long double halfwidth = (xmax - xmin)/2.0L;
         std::vector<long double> fval(N);
 
         for (unsigned k=0; k<N; ++k)
         {
             const long double t = (cosl(Pi*(k + 0.5L)/N)+1.0L)*halfwidth+xmin;
             fval[k] = f(t);
         }
 
         const long double norm = 2.0L/N;
         for (unsigned j=0; j<N; ++j)
         {
             long double sum = 0.0L;
             for (unsigned k=0; k<N; ++k)
                 sum += fval[k]*cosl(Pi*j*(k + 0.5L)/N);
             if (!j)
                 sum /= 2.0L;
             coeffs[j] = static_cast<Numeric>(norm*sum);
         }
     }
+
+    template<typename Numeric>
+    void chebyshevIntegralCoeffs(const Numeric* coeffs, const unsigned degree,
+                                 const long double xmin, const long double xmax,
+                                 Numeric* integralCoeffs)
+    {
+        assert(coeffs);
+        assert(integralCoeffs);
+        assert(integralCoeffs != coeffs);
+        assert(xmin != xmax);
+
+        const long double norm = 0.25L*(xmax - xmin);
+        long double sum = 0.0L;
+        long double f = -1.0L;
+
+        integralCoeffs[1U] = norm*(2.0*coeffs[0] - coeffs[2]);
+        sum += integralCoeffs[1U];
+
+        for (unsigned j=2U; j<degree; ++j)
+        {
+            integralCoeffs[j] = norm*(coeffs[j-1] - coeffs[j+1])/j;
+            sum += f*integralCoeffs[j];
+            f *= -1.0L;
+        }
+
+        for (unsigned j=degree; j<degree+2; ++j)
+        {
+            integralCoeffs[j] = norm*(coeffs[j-1])/j;
+            sum += f*integralCoeffs[j];
+            f *= -1.0L;
+        }
+
+        integralCoeffs[0U] = sum;
+    }
 }
Index: trunk/npstat/nm/MathUtils.hh
===================================================================
--- trunk/npstat/nm/MathUtils.hh	(revision 867)
+++ trunk/npstat/nm/MathUtils.hh	(revision 868)
@@ -1,187 +1,236 @@
 #ifndef NPSTAT_MATHUTILS_HH_
 #define NPSTAT_MATHUTILS_HH_
 
 /*!
 // \file MathUtils.hh
 //
 // \brief Various simple mathematical utilities which did not end up inside
 //        dedicated headers
 //
 // Author: I. Volobouev
 //
 // March 2010
 */
 
 #include <utility>
 #include <vector>
 
 namespace npstat {
     /**
     // Solve the quadratic equation x*x + b*x + c == 0
     // in a numerically sound manner. Return the number of roots.
     */
     unsigned solveQuadratic(double b, double c,
                             double *x1, double *x2);
 
     /**
     // Find the real roots of the cubic:
     //   x**3 + p*x**2 + q*x + r = 0
     // The number of real roots is returned, and the roots are placed
     // into the "v3" array. Original code by Don Herbison-Evans (see
     // his article "Solving Quartics and Cubics for Graphics" in the
     // book "Graphics Gems V", page 3), with minimal adaptation for
     // this package by igv.
     */
     unsigned solveCubic(double p, double q, double r, double v3[3]);
 
     /**
     // Find an extremum of a parabola passing through the three given points.
     // The returned value is "true" for minimum and "false" for maximum.
     // std::invalid_argument will be thrown in case some of the x values
     // coincide or if the coordinates describe a straight line.
     */
     bool parabolicExtremum(double x0, double y0, double x1, double y1,
                            double x2, double y2, double* extremumCoordinate,
                            double* extremumValue);
 
     /**
     // Volume of the n-dimensional sphere of unit radius. Should be
     // multuplied by R^n to get the volume of the sphere with arbitrary radius.
     */
     double ndUnitSphereVolume(unsigned n);
 
     /**
     // Area of the sphere of unit radius embedded in the n-dimensional space.
     // Should be multuplied by R^(n-1) to get the area of the sphere with
     // arbitrary radius.
     */
     double ndUnitSphereArea(unsigned n);
 
     /**
     // Sum of polynomial series. The length of the
     // array of coefficients should be at least degree+1.
     // The highest degree coefficient is assumed to be
     // the last one in the "coeffs" array (0th degree
     // coefficient comes first).
     */
     template<typename Numeric>
     long double polySeriesSum(const Numeric *coeffs,
                               unsigned degree, long double x);
 
     /**
     // Sum and derivative of polynomial series. The length of the
     // array of coefficients should be at least degree+1.
     */
     template<typename Numeric>
     void polyAndDeriv(const Numeric *coeffs, unsigned degree,
                       long double x, long double *value, long double *deriv);
 
     /**
+    // Coefficients for the integral of polynomial series.
+    // The integration constant is set to 0. The length of the
+    // array of coefficients should be at least degree+1,
+    // and the length of the buffer for the integral coefficients
+    // should be at least degree+2.
+    */
+    template<typename Numeric>
+    void polyIntegralCoeffs(const Numeric* coeffs, unsigned degree,
+                            Numeric* integralCoeffs);
+
+    /**
     // Series using Legendre polynomials. 0th degree coefficient
     // comes first. Although any value of x can be specified, the result
     // is not going to be terribly meaningful in case |x| > 1.
     */
     template<typename Numeric>
     long double legendreSeriesSum(const Numeric *coeffs,
                                   unsigned degree, long double x);
 
     /**
     // Series for Hermite polynomials orthogonal with weight exp(-x*x/2).
     // These are sometimes called the "probabilists' Hermite polynomials".
     */
     template<typename Numeric>
     long double hermiteSeriesSumProb(const Numeric *coeffs,
                                      unsigned degree, long double x);
 
     /**
     // Series for Hermite polynomials orthogonal with weight exp(-x*x). These
     // are sometimes called the "physicists' Hermite polynomials". The weight
     // exp(-x*x) is also used in the "GaussHermiteQuadrature" class.
     */
     template<typename Numeric>
     long double hermiteSeriesSumPhys(const Numeric *coeffs,
                                      unsigned degree, long double x);
 
     /** Series for Gegenbauer polynomials. "lambda" is the parameter. */
     template<typename Numeric>
     long double gegenbauerSeriesSum(const Numeric *coeffs, unsigned degree,
                                     long double lambda, long double x);
 
     /** Series for the Chebyshev polynomials of the first kind */
     template<typename Numeric>
-    long double chebyshevSeriesSum(const Numeric *coeffs,
-                                   unsigned degree, long double x);
+    long double chebyshevSeriesSum(const Numeric *coeffs, unsigned degree,
+                                   long double x);
+
+    /**
+    // Series for the Chebyshev polynomials of the first kind
+    // on the given interval [xmin, xmax]
+    */
+    template<typename Numeric>
+    long double chebyshevSeriesSum(const Numeric *coeffs, unsigned degree,
+                                   long double xmin, long double xmax,
+                                   long double x);
 
     /**
     // Generate Chebyshev series coefficients for a given functor on the
     // interval [xmin, xmax]. The array of coefficients must be at least
     // degree+1 long. The functor will be given a long double argument.
     */
     template<typename Functor, typename Numeric>
     void chebyshevSeriesCoeffs(const Functor& f,
                                long double xmin, long double xmax,
                                unsigned degree, Numeric *coeffs);
 
+    /**
+    // Converts Chebyshev series coefficients into coefficients
+    // for the function integral. The 0's degree coefficient
+    // will be chosen in such a way that the integral is 0
+    // at xmin. The array of coefficients must be at least
+    // degree+1 long, and the buffer for the integral coefficients
+    // must be at least degree+2 long.
+    */
+    template<typename Numeric>
+    void chebyshevIntegralCoeffs(const Numeric* coeffs, unsigned degree,
+                                 long double xmin, long double xmax,
+                                 Numeric* integralCoeffs);
+
 #ifdef SWIG
     inline double polySeriesSum2(const double *c, const unsigned degreep1,
                                  const double x)
         {return polySeriesSum(c, degreep1-1U, x);}
 
     inline std::pair<double,double> polyAndDeriv2(
         const double *c, const unsigned degreep1,
         const double x)
     {
         long double p, d;
         polyAndDeriv(c, degreep1-1U, x, &p, &d);
         return std::pair<double,double>(p, d);
     }
 
+    inline std::vector<double> polyIntegralCoeffs2(const double *c, const unsigned degreep1)
+    {
+        std::vector<double> tmp(degreep1+1U);
+        polyIntegralCoeffs(c, degreep1-1U, &tmp[0]);
+        return tmp;
+    }
+
     inline double legendreSeriesSum2(const double *c, const unsigned degreep1,
                                      const double x)
         {return legendreSeriesSum(c, degreep1-1U, x);}
 
     inline double hermiteSeriesSumProb2(const double *c, const unsigned degreep1,
                                         const double x)
         {return hermiteSeriesSumProb(c, degreep1-1U, x);}
 
     inline double hermiteSeriesSumPhys2(const double *c, const unsigned degreep1,
                                         const double x)
         {return hermiteSeriesSumPhys(c, degreep1-1U, x);}
 
     inline double gegenbauerSeriesSum2(const double *c, const unsigned degreep1,
                                        const double lambda, const double x)
         {return gegenbauerSeriesSum(c, degreep1-1U, lambda, x);}
 
     inline double chebyshevSeriesSum2(const double *c, const unsigned degreep1,
+                                      const double xmin, const double xmax,
                                       const double x)
-        {return chebyshevSeriesSum(c, degreep1-1U, x);}
+        {return chebyshevSeriesSum(c, degreep1-1U, xmin, xmax, x);}
+
+    inline std::vector<double> chebyshevIntegralCoeffs2(const double *c, const unsigned degreep1,
+                                                        const double xmin, const double xmax)
+    {
+        std::vector<double> tmp(degreep1+1U);
+        chebyshevIntegralCoeffs(c, degreep1-1U, xmin, xmax, &tmp[0]);
+        return tmp;
+    }
 
     template<typename Functor>
     inline std::vector<double> chebyshevSeriesCoeffs2(const Functor& f,
         const double xmin, const double xmax, const unsigned degree)
     {
         std::vector<double> tmp(degree+1U);
         chebyshevSeriesCoeffs(f, xmin, xmax, degree, &tmp[0]);
         return tmp;
     }
 
     inline std::vector<double> solveQuadratic2(const double b, const double c)
     {
         double x[2];
         const unsigned nRoots = solveQuadratic(b, c, &x[0], &x[1]);
         return std::vector<double>(x, x+nRoots);
     }
 
     inline std::vector<double> solveCubic2(const double p, const double q,
                                            const double r)
     {
         double v3[3];
         const unsigned nRoots = solveCubic(p, q, r, v3);
         return std::vector<double>(v3, v3+nRoots);
     }
 #endif // SWIG
 }
 
 #include "npstat/nm/MathUtils.icc"
 
 #endif // NPSTAT_MATHUTILS_HH_