Index: trunk/npstat/nm/Poly1D.cc
===================================================================
--- trunk/npstat/nm/Poly1D.cc	(revision 752)
+++ trunk/npstat/nm/Poly1D.cc	(revision 753)
@@ -1,433 +1,433 @@
 #include <cmath>
 #include <cfloat>
 #include <cassert>
 #include <stdexcept>
 #include <algorithm>
 #include <utility>
 
 #include "npstat/nm/Poly1D.hh"
 #include "npstat/nm/MathUtils.hh"
 
 static unsigned nSignChanges(const long double* seq, const unsigned len)
 {
     assert(len);
     assert(seq);
 
     unsigned nChanges = 0;
     int sign = 100;
     for (unsigned i=0; i<len; ++i)
         if (seq[i])
         {
             const int newsign = seq[i] > 0.0L ? 1 : -1;
             if (sign != 100)
             {
                 if (sign*newsign < 0)
                     ++nChanges;
             }
             sign = newsign;
         }
     return nChanges;
 }
 
 namespace npstat {
     Poly1D::Poly1D(const std::vector<double>& c)
         : coeffs_(c.begin(), c.end())
     {
         if (coeffs_.empty())
-            coeffs_.resize(1, 0.0);
+            coeffs_.resize(1, 0.0L);
         else
             truncateLeadingZeros();
     }
 
     Poly1D::Poly1D(const std::vector<long double>& c)
         : coeffs_(c)
     {
         if (coeffs_.empty())
-            coeffs_.resize(1, 0.0);
+            coeffs_.resize(1, 0.0L);
         else
             truncateLeadingZeros();
     }
 
     Poly1D::Poly1D(const unsigned degree, const long double degCoeff)
-        : coeffs_(degree+1U, 0.0)
+        : coeffs_(degree+1U, 0.0L)
     {
         coeffs_[degree] = degCoeff;
-        if (degree && degCoeff == 0.0)
+        if (degree && degCoeff == 0.0L)
             coeffs_.resize(1);
     }
 
     Poly1D::Poly1D(const long double* coeffs, const unsigned degree)
         : coeffs_(coeffs, coeffs + (degree+1U))
     {
         assert(coeffs);
         truncateLeadingZeros();
     }
 
     Poly1D::Poly1D(const double* coeffs, const unsigned degree)
         : coeffs_(coeffs, coeffs + (degree+1U))
     {
         assert(coeffs);
         truncateLeadingZeros();
     }
 
     long double Poly1D::operator()(const long double x) const
     {
         return polySeriesSum(&coeffs_[0], coeffs_.size()-1U, x);
     }
 
     Poly1D Poly1D::derivative() const
     {
         const unsigned myDeg = deg();
         if (myDeg)
         {
             Poly1D result(myDeg-1U, myDeg*coeffs_[myDeg]);
             for (unsigned i=0; i<myDeg; ++i)
                 result.coeffs_[i] = (i+1)*coeffs_[i+1];
             return result;
         }
         else
         {
             Poly1D dummy;
             return dummy;
         }
     }
 
     Poly1D Poly1D::integral(const long double c) const
     {
         const unsigned myDeg = deg();
         if (myDeg == 0U)
         {
             if (coeffs_[0])
             {
                 long double cr[2];
                 cr[0] = c;
                 cr[1] = coeffs_[0];
                 return Poly1D(cr, 1U);
             }
             else
                 return c;
         }
         else
         {
             Poly1D result(myDeg+1U, coeffs_[myDeg]/(myDeg+1U));
             for (unsigned i=0; i<myDeg; ++i)
                 result.coeffs_[i+1U] = coeffs_[i]/(i+1U);
             result.coeffs_[0] = c;
             return result;
         }
     }
 
     void Poly1D::truncate(const unsigned degree)
     {
         if (degree < deg())
         {
             coeffs_.resize(degree + 1U);
             truncateLeadingZeros();
         }
     }
 
     bool Poly1D::isClose(const Poly1D& r, const long double eps) const
     {
-        if (eps < 0.0) throw std::invalid_argument(
+        if (eps < 0.0L) throw std::invalid_argument(
             "In npstat::Poly1D::isClose: eps must not be negative");
         const unsigned maxDeg = std::max(deg(), r.deg());
         for (unsigned i=0; i<=maxDeg; ++i)
             if (std::abs(r[i] - (*this)[i]) > eps)
                 return false;
         return true;
     }
 
     Poly1D Poly1D::operator-() const
     {
         Poly1D poly(*this);
         const unsigned len = coeffs_.size();
         for (unsigned i=0; i<len; ++i)
-            poly.coeffs_[i] *= -1.0;
+            poly.coeffs_[i] *= -1.0L;
         return poly;
     }
 
     void Poly1D::truncateLeadingZeros()
     {
         const unsigned sz = coeffs_.size();
         if (sz > 1U)
         {
             unsigned firstNon0 = sz;
             for (unsigned i=sz-1U; i>0; --i)
                 if (coeffs_[i])
                 {
                     firstNon0 = i;
                     break;
                 }
             if (firstNon0 == sz)
                 coeffs_.resize(1);
             else if (firstNon0 < sz-1U)
                 coeffs_.resize(firstNon0+1U);
         }
     }
 
     Poly1D& Poly1D::operator*=(const Poly1D& r)
     {
         if (isNull())
             return *this;
 
         if (r.isNull())
         {
             coeffs_.resize(1);
             coeffs_[0] = 0.0L;
             return *this;
         }
 
         const unsigned mydeg = deg();
         const unsigned rdeg = r.deg();
         if (mydeg == 0U)
         {
             const long double c = coeffs_[0];
             coeffs_ = r.coeffs_;
             for (unsigned i=0; i<=rdeg; ++i)
                 coeffs_[i] *= c;
         }
         else if (rdeg == 0U)
         {
             const long double c = r.coeffs_[0];
             for (unsigned i=0; i<=mydeg; ++i)
                 coeffs_[i] *= c;
         }
         else
         {
             std::vector<long double> prod(mydeg + rdeg + 1U, 0.0L);
             for (unsigned i=0; i<=mydeg; ++i)
                 for (unsigned j=0; j<=rdeg; ++j)
                     prod[i + j] += coeffs_[i]*r.coeffs_[j];
             coeffs_.swap(prod);
         }
 
         return *this;
     }
 
     Poly1D Poly1D::operator*(const Poly1D& r) const
     {
         if (isNull() || r.isNull())
         {
             Poly1D dummy;
             return dummy;
         }
 
         const unsigned mydeg = deg();
         const unsigned rdeg = r.deg();
         if (mydeg == 0U)
         {
             Poly1D p(r);
             const long double c = coeffs_[0];
             for (unsigned i=0; i<=rdeg; ++i)
                 p.coeffs_[i] *= c;
             return p;
         }
         else if (rdeg == 0U)
         {
             Poly1D p(*this);
             const long double c = r.coeffs_[0];
             for (unsigned i=0; i<=mydeg; ++i)
                 p.coeffs_[i] *= c;
             return p;
         }
         else
         {
             std::vector<long double> prod(mydeg + rdeg + 1U, 0.0L);
             for (unsigned i=0; i<=mydeg; ++i)
                 for (unsigned j=0; j<=rdeg; ++j)
                     prod[i + j] += coeffs_[i]*r.coeffs_[j];
             return Poly1D(prod);
         }
     }
 
     Poly1D Poly1D::operator+(const Poly1D& r) const
     {
         const unsigned maxdeg = std::max(deg(), r.deg());
         Poly1D result(maxdeg, (*this)[maxdeg]+r[maxdeg]);
         for (unsigned i=0; i<maxdeg; ++i)
             result.coeffs_[i] = (*this)[i]+r[i];
         result.truncateLeadingZeros();
         return result;
     }
 
     Poly1D Poly1D::operator-(const Poly1D& r) const
     {
         const unsigned maxdeg = std::max(deg(), r.deg());
         Poly1D result(maxdeg, (*this)[maxdeg]-r[maxdeg]);
         for (unsigned i=0; i<maxdeg; ++i)
             result.coeffs_[i] = (*this)[i]-r[i];
         result.truncateLeadingZeros();
         return result;
     }
 
     Poly1D& Poly1D::operator+=(const Poly1D& r)
     {
         const unsigned maxdeg = std::max(deg(), r.deg());
         if (maxdeg > deg())
             coeffs_.resize(maxdeg + 1U);
         for (unsigned i=0; i<=maxdeg; ++i)
             coeffs_[i] += r[i];
         truncateLeadingZeros();
         return *this;
     }
 
     Poly1D& Poly1D::operator-=(const Poly1D& r)
     {
         const unsigned maxdeg = std::max(deg(), r.deg());
         if (maxdeg > deg())
             coeffs_.resize(maxdeg + 1U);
         for (unsigned i=0; i<=maxdeg; ++i)
             coeffs_[i] -= r[i];
         truncateLeadingZeros();
         return *this;
     }
 
     long double Poly1D::leadingCoefficient() const
     {
         const unsigned sz = coeffs_.size();
         assert(sz);
         const long double c = coeffs_.back();
         if (sz > 1U)
             assert(c);
         return c;
     }
 
     Poly1D Poly1D::operator/(const Poly1D& b) const
     {
         if (b.isNull()) throw std::invalid_argument(
             "In npstat::Poly1D::operator/: division by zero encountered");
 
         const unsigned d = b.deg();
         const long double c = b.leadingCoefficient();
         if (d)
         {
             Poly1D q;
             Poly1D r(*this);
             if (r.deg() >= d)
                 q.reserve(r.deg() - d);
             while (r.deg() >= d)
             {
                 const Poly1D s(r.deg()-d, r.leadingCoefficient()/c);
                 const Poly1D& sb = s*b;
                 q += s;
                 r -= sb;
                 r.truncate(sb.deg() - 1U);
             }
             return q;
         }
         else
         {
             Poly1D q(*this);
             const unsigned mydeg = deg();
             for (unsigned i=0; i<=mydeg; ++i)
                 q.coeffs_[i] /= c;
             return q;
         }
     }
 
     Poly1D Poly1D::operator%(const Poly1D& b) const
     {
         if (b.isNull()) throw std::invalid_argument(
             "In npstat::Poly1D::operator%: division by zero encountered");
 
         const unsigned d = b.deg();
         if (d)
         {
             Poly1D r(*this);
             const long double c = b.leadingCoefficient();
             while (r.deg() >= d)
             {
                 const Poly1D s(r.deg()-d, r.leadingCoefficient()/c);
                 const Poly1D& sb = s*b;
                 r -= sb;
                 r.truncate(sb.deg() - 1U);
             }
             return r;
         }
         else
         {
             Poly1D dummy;
             return dummy;
         }
     }
 
     unsigned Poly1D::nRoots(const long double a, const long double b) const
     {
         if (a >= b) throw std::invalid_argument(
             "In npstat::Poly1D::nRoots: invalid "
             "interval definition, must have a < b");
         const unsigned mydeg = deg();
         switch (mydeg)
         {
         case 0U:
             return 0U;
 
         case 1U:
         {
             const long double pa = (*this)(a);
             const long double pb = (*this)(b);
             const long double prod = pa*pb;
             if (prod < 0.0L)
                 return 1U;
             else if (prod > 0.0L)
                 return 0U;
             else if (pb == 0.0L)
                 return 1U;
             else
                 return 0U;
         }
 
         default:
         {
             // Construct the Sturm sequence
             std::vector<long double> buf(2U*(mydeg+1U));
             long double* avalues = &buf[0];
             long double* bvalues = avalues + (mydeg + 1U);
 
             Poly1D p_im1(*this);
             Poly1D p_i(derivative());
             Poly1D* ptr_im1 = &p_im1;
             Poly1D* ptr_i = &p_i;
             avalues[0] = p_im1(a);
             avalues[1] = p_i(a);
             bvalues[0] = p_im1(b);
             bvalues[1] = p_i(b);
             unsigned seqLength = 2U;
             for (;;)
             {
                 *ptr_im1 = -(*ptr_im1 % *ptr_i);
                 if (ptr_im1->isNull())
                     break;
                 assert(seqLength <= mydeg);
                 avalues[seqLength] = (*ptr_im1)(a);
                 bvalues[seqLength++] = (*ptr_im1)(b);
                 std::swap(ptr_im1, ptr_i);
             }
             const unsigned na = nSignChanges(avalues, seqLength);
             const unsigned nb = nSignChanges(bvalues, seqLength);
             assert(na >= nb);
             return na - nb;
         }
         }
     }
 
     Poly1D Poly1D::monicDeg0()
     {
         return 1.0L;
     }
 
     Poly1D Poly1D::monicDeg1(const long double b)
     {
         Poly1D mono(1, 1.0L);
         mono.coeffs_[0] = b;
         return mono;
     }
 
     Poly1D Poly1D::monicDeg2(const long double b, const long double c)
     {
         Poly1D mono(2, 1.0L);
         mono.coeffs_[0] = c;
         mono.coeffs_[1] = b;
         return mono;
     }
 }