Index: trunk/npstat/stat/AbsCGF1D.cc
===================================================================
--- trunk/npstat/stat/AbsCGF1D.cc	(revision 686)
+++ trunk/npstat/stat/AbsCGF1D.cc	(revision 687)
@@ -1,164 +1,164 @@
 #include <cassert>
 #include <limits>
 #include <cmath>
 
 #include "npstat/stat/AbsCGF1D.hh"
 #include "npstat/nm/findRootUsingBisections.hh"
 
 namespace {
     class AbsCGF1DDeriv : public npstat::Functor1<double, double>
     {
     public:
         inline AbsCGF1DDeriv(const npstat::AbsCGF1D& cgf)
             : cgf_(cgf) {}
         inline double operator()(const double& s) const
             {return cgf_.derivative(1U, s);}
 
     private:
         const npstat::AbsCGF1D& cgf_;
     };
 }
 
 namespace npstat {
     double AbsCGF1D::saddlepoint(const double x) const
     {
         static const double factor = sqrt(2.0);
 
         const double infimum = smin();
         assert(infimum < 0.0);
         const double supremum = smax();
         assert(supremum > 0.0);
 
         // The first derivative is monotonously increasing
         const double d0 = derivative(1U, 0.0);
         if (x == d0)
             return 0.0;
 
         double xmin = 0.0, xmax = 0.0;
         double dmin = d0, dmax = d0;
         const double eps = std::numeric_limits<double>::epsilon();
         if (x > d0)
         {
             // Need to find useful xmax
             xmax = 1.0;
             while (xmax < supremum)
             {
                 dmax = derivative(1U, xmax);
                 if (x <= dmax)
                     break;
                 if (xmax >= std::numeric_limits<double>::max()/1.001/factor)
                     break;
                 xmax *= factor;
             }
             if (x > dmax)
             {
                 xmax = supremum/(1.0 + eps);
                 dmax = derivative(1U, xmax);
             }
             if (!(x <= dmax)) throw std::invalid_argument(
                 "Failed to solve the saddlepoint equation "
                 "within the CGF support");
         }
         else
         {
             // Need to find useful xmin
             xmin = -1.0;
             while (xmin > infimum)
             {
                 dmin = derivative(1U, xmin);
                 if (x >= dmin)
                     break;
                 if (xmin < -std::numeric_limits<double>::max()/1.001/factor)
                     break;
                 xmin *= factor;
             }
             if (x < dmin)
             {
                 xmin = infimum/(1.0 + eps);
                 dmin = derivative(1U, xmin);
             }
             if (!(x >= dmin)) throw std::invalid_argument(
                 "Failed to solve the saddlepoint equation "
                 "within the CGF support");
         }
 
         assert(x >= dmin);
         assert(x <= dmax);
 
         if (x == dmin)
             return xmin;
         if (x == dmax)
             return xmax;
 
         double root;
         AbsCGF1DDeriv deriv(*this);
         const bool status = findRootUsingBisections(deriv, x, xmin, xmax,
                                                     2.0*eps, &root);
         assert(status);
         return root;
     }
 
     double AbsCGF1D::standardizedCumulant(const unsigned order,
                                           const double s) const
     {
         if (order)
         {
             if (!(smin() < s && s < smax())) throw std::invalid_argument(
                 "In npstat::AbsCGF1D::standardizedCumulant: "
                 "argument outside of CGF support");
             const double var = derivative(2U, s);
             assert(var > 0.0);
             return derivative(order, s)/pow(var, order/2.0);
         }
         return 0.0;
     }
 
     double AbsCGF1D::saddlepointDensityApprox(
         const unsigned order, const double x) const
     {
         if (!order) throw std::invalid_argument(
             "In npstat::AbsCGF1D::saddlepointDensityApproximation: "
             "approximation order must be positive");
         if (order > 4U) throw std::invalid_argument(
             "In npstat::AbsCGF1D::saddlepointDensityApproximation: "
             "unsupported approximation order (too high)");
         const double shat = saddlepoint(x);
         const double Khat = derivative(0U, shat);
         const double varhat = derivative(2U, shat);
         assert(varhat > 0.0);
-        double factor = 1.0;
+        long double factor = 1.0L;
         if (order > 1U)
         {
             const double k3 = standardizedCumulant(3U, shat);
             const double k4 = standardizedCumulant(4U, shat);
             const double k3_2 = k3*k3;
             factor += (k4/8 - 5*k3_2/24);
 
             if (order > 2U)
             {
                 const double k5 = standardizedCumulant(5U, shat);
                 const double k6 = standardizedCumulant(6U, shat);
                 const double k3_4 = k3_2*k3_2;
                 const double k4_2 = k4*k4;
                 factor += ((385*k3_4)/1152 - (35*k3_2*k4)/64 +
                            (35*k4_2)/384 + (7*k3*k5)/48 - k6/48);
 
                 if (order > 3U)
                 {
                     const double k7 = standardizedCumulant(7U, shat);
                     const double k8 = standardizedCumulant(8U, shat);
                     const double k3_3 = k3*k3_2;
                     const double k3_6 = k3_3*k3_3;
                     const double k4_3 = k4*k4_2;
                     const double k5_2 = k5*k5;
                     factor += ((-85085*k3_6)/82944 + (25025*k3_4*k4)/9216 - 
                                (5005*k3_2*k4_2)/3072 + (385*k4_3)/3072 - 
                                (1001*k3_3*k5)/1152 + (77*k3*k4*k5)/128 - 
                                (21*k5_2)/640 + (77*k3_2*k6)/384 - 
                                (7*k4*k6)/128 - (k3*k7)/32 + k8/384);
                 }
             }
         }
-        return factor/sqrt(2.0*M_PI*varhat)*exp(Khat - shat*x);
+        return static_cast<double>(factor)/sqrt(2.0*M_PI*varhat)*exp(Khat - shat*x);
     }
 }