Index: trunk/npstat/stat/DensityOrthoPoly1D.hh
===================================================================
--- trunk/npstat/stat/DensityOrthoPoly1D.hh	(revision 698)
+++ trunk/npstat/stat/DensityOrthoPoly1D.hh	(revision 699)
@@ -1,101 +1,103 @@
 #ifndef NPSTAT_DENSITYORTHOPOLY1D_HH_
 #define NPSTAT_DENSITYORTHOPOLY1D_HH_
 
 /*!
 // \file DensityOrthoPoly1D.hh
 //
 // \brief Continuous orthonormal polynomial series in one dimension
 //        generated for a continuous density used as the weight
 //
 // Author: I. Volobouev
 //
 // July 2020
 */
 
 #include <vector>
 #include <cassert>
 #include <cmath>
 
 #include "npstat/nm/AbsClassicalOrthoPoly1D.hh"
 #include "npstat/nm/OrthoPolyMethod.hh"
 
 #include "npstat/stat/AbsDistribution1D.hh"
 
 namespace npstat {
     class DensityOrthoPoly1D : public AbsClassicalOrthoPoly1D
     {
     public:
         /**
         // "nIntegPoints" parameter is the number of integration points
         // for the Fejer quadrature. Note that Fejer quadrature with n
         // points can integrate exactly polynomials of degree n-1 and
         // that this quadrature will be used for calculating various
         // inner products in which the density is used as the weight.
         // Therefore, "nIntegPoints" should be at least
         // 2*maxDegree + 1 + C, where "C" is the polynomial degree
         // that can be used to approximate the density reasonably well.
         //
         // The support of the weight for the integration purpose will
         // be set to [distro.quantile(0.0), distro.quantile(1.0)]. If
         // this interval is very wide, the results may be meaningless
         // (monomials like x^n can exceed the floating point dynamic
         // range, etc).
         */
         DensityOrthoPoly1D(const AbsDistribution1D& distro,
                            unsigned maxDegree,
                            unsigned nIntegPoints,
                            OrthoPolyMethod m = OPOLY_STIELTJES);
 
+        DensityOrthoPoly1D(const DensityOrthoPoly1D&);
+        DensityOrthoPoly1D& operator=(const DensityOrthoPoly1D&);
+
         inline virtual ~DensityOrthoPoly1D() {}
 
         inline virtual DensityOrthoPoly1D* clone() const
             {return new DensityOrthoPoly1D(*this);}
 
         inline virtual long double weight(const long double x) const
             {return distro_->density(x);}
         inline virtual double xmin() const {return distro_->quantile(0.0);}
         inline virtual double xmax() const {return distro_->quantile(1.0);}
         inline virtual unsigned maxDegree() const {return maxdeg_;}
 
     protected:
         virtual std::pair<long double,long double>
         recurrenceCoeffs(unsigned deg) const;
 
     private:
         typedef std::pair<double,double> MeasurePt;
 
         struct Recurrence
         {
             inline Recurrence(const long double a,
                               const long double b)
                 : alpha(a), beta(b)
             {
                 assert(beta > 0.0L);
                 sqrbeta = sqrtl(beta);
             }
 
             long double alpha;
             long double beta;
             long double sqrbeta;
         };
 
+        // Disable default constructor
         DensityOrthoPoly1D();
-        DensityOrthoPoly1D(const DensityOrthoPoly1D&);
-        DensityOrthoPoly1D& operator=(const DensityOrthoPoly1D&);
 
         void calcRecurrenceCoeffs(const std::vector<MeasurePt>& measure,
                                   OrthoPolyMethod m);
         void calcRecurrenceStieltjes(const std::vector<MeasurePt>& measure);
         void calcRecurrenceLanczos(const std::vector<MeasurePt>& measure);
         std::pair<long double,long double> monicInnerProducts(
             const std::vector<MeasurePt>& measure, unsigned degree) const;
         long double monicpoly(unsigned degree, long double x) const;
 
         CPP11_auto_ptr<AbsDistribution1D> distro_;
         std::vector<Recurrence> rCoeffs_;
         long double wsum_;
         unsigned maxdeg_;
     };
 }
 
 #endif // NPSTAT_DENSITYORTHOPOLY1D_HH_
Index: trunk/npstat/stat/DensityOrthoPoly1D.cc
===================================================================
--- trunk/npstat/stat/DensityOrthoPoly1D.cc	(revision 698)
+++ trunk/npstat/stat/DensityOrthoPoly1D.cc	(revision 699)
@@ -1,160 +1,172 @@
 #include <stdexcept>
 
 #include "npstat/nm/FejerQuadrature.hh"
 #include "npstat/stat/DensityOrthoPoly1D.hh"
 
 namespace npstat {
     DensityOrthoPoly1D::DensityOrthoPoly1D(const AbsDistribution1D& distro,
                                            const unsigned maxDegree,
                                            const unsigned nIntegPoints,
                                            const OrthoPolyMethod m)
         : distro_(distro.clone()), wsum_(0.0L), maxdeg_(maxDegree)
     {
         if (!nIntegPoints) throw std::invalid_argument(
             "In npstat::DensityOrthoPoly1D constructor: the number of "
             "integration points must be positive");
         const FejerQuadrature quad(nIntegPoints);
         const DensityFunctor1D weight(distro);
         const std::vector<MeasurePt>& measure = quad.weightedIntegrationPoints(
             weight, distro_->quantile(0.0), distro_->quantile(1.0));
         assert(measure.size() == nIntegPoints);
         for (unsigned i=0; i<nIntegPoints; ++i)
             wsum_ += measure[i].second;
         calcRecurrenceCoeffs(measure, m);
     }
 
     void DensityOrthoPoly1D::calcRecurrenceCoeffs(
         const std::vector<MeasurePt>& measure,
         const OrthoPolyMethod m)
     {
         rCoeffs_.clear();
         rCoeffs_.reserve(maxdeg_ + 1);
         switch (m)
         {
         case OPOLY_STIELTJES:
             calcRecurrenceStieltjes(measure);
             break;
 
         case OPOLY_LANCZOS:
             calcRecurrenceLanczos(measure);
             break;
 
         default:
             throw std::runtime_error(
                 "In npstat::DensityOrthoPoly1D::calcRecurrenceCoeffs: "
                 "incomplete switch statement. This is a bug. Please report.");
         }
         assert(rCoeffs_.size() == maxdeg_ + 1);
     }
 
     void DensityOrthoPoly1D::calcRecurrenceStieltjes(const std::vector<MeasurePt>& measure)
     {
         long double prevSprod = 1.0L;
         for (unsigned deg=0; deg<=maxdeg_; ++deg)
         {
             const std::pair<long double,long double> ip = monicInnerProducts(measure, deg);
             rCoeffs_.push_back(Recurrence(ip.first/ip.second, ip.second/prevSprod));
             prevSprod = ip.second;
         }
     }
 
     void DensityOrthoPoly1D::calcRecurrenceLanczos(const std::vector<MeasurePt>& measure)
     {
         typedef long double Real;
 
         const unsigned mSize = measure.size();
         std::vector<Real> dmem(mSize*2U);
 
         Real* p0 = &dmem[0];
         for (unsigned i=0; i<mSize; ++i)
             p0[i] = measure[i].first;
 
         Real* p1 = &dmem[mSize];
         clearBuffer(p1, mSize);
         p1[0] = measure[0].second/wsum_;
 
         const unsigned mm1 = mSize - 1;
         for (unsigned n=0; n<mm1; ++n)
         {
             const unsigned np1 = n + 1U;
             Real xlam = measure[np1].first;
             Real pn = measure[np1].second/wsum_;
             Real gam = 1.0, sig = 0.0, t = 0.0;
 
             for (unsigned k=0; k<=np1; ++k)
             {
                 const Real rho = p1[k] + pn;
                 const Real tmp = gam*rho;
                 Real tsig = sig;
                 if (rho <= 0.0)
                 {
                     gam = 1.0;
                     sig = 0.0;
                 }
                 else
                 {
                     gam = p1[k]/rho;
                     sig = pn/rho;
                 }
                 const Real tk = sig*(p0[k] - xlam) - gam*t;
                 p0[k] -= tk - t;
                 t = tk;
                 if (sig < 0.0)
                     pn = tsig*p1[k];
                 else
                     pn = t*t/sig;
                 tsig = sig;
                 p1[k] = tmp;
             }
         }
 
         for (unsigned deg=0; deg<=maxdeg_; ++deg)
             rCoeffs_.push_back(Recurrence(p0[deg], p1[deg]));
     }
 
     std::pair<long double,long double>
     DensityOrthoPoly1D::monicInnerProducts(const std::vector<MeasurePt>& measure,
                                            const unsigned degree) const
     {
         long double sum = 0.0L, xsum = 0.0L;
         const unsigned mSize = measure.size();
         for (unsigned i = 0; i < mSize; ++i)
         {
             const long double x = measure[i].first;
             const long double p = monicpoly(degree, x);
             const long double pprod = p*p*measure[i].second;
             sum += pprod;
             xsum += x*pprod;
         }
         return std::pair<long double,long double>(xsum/wsum_, sum/wsum_);
     }
 
     long double DensityOrthoPoly1D::monicpoly(const unsigned degree,
                                               const long double x) const
     {
         long double polyk = 1.0L, polykm1 = 0.0L;
         for (unsigned k=0; k<degree; ++k)
         {
             const long double p = (x - rCoeffs_[k].alpha)*polyk -
                                    rCoeffs_[k].beta*polykm1;
             polykm1 = polyk;
             polyk = p;
         }
         return polyk;
     }
 
     std::pair<long double,long double>
     DensityOrthoPoly1D::recurrenceCoeffs(const unsigned deg) const
     {
         const Recurrence& r = rCoeffs_.at(deg);
         return std::pair<long double,long double>(r.alpha, r.sqrbeta);
     }
 
     DensityOrthoPoly1D::DensityOrthoPoly1D(const DensityOrthoPoly1D& r)
         : distro_(r.distro_->clone()),
           rCoeffs_(r.rCoeffs_),
           wsum_(r.wsum_),
           maxdeg_(r.maxdeg_)
     {
     }
+
+    DensityOrthoPoly1D& DensityOrthoPoly1D::operator=(const DensityOrthoPoly1D& r)
+    {
+        if (&r != this)
+        {
+            distro_ = CPP11_auto_ptr<AbsDistribution1D>(r.distro_->clone());
+            rCoeffs_ = r.rCoeffs_;
+            wsum_ = r.wsum_;
+            maxdeg_ = r.maxdeg_;
+        }
+        return *this;
+    }
 }