Pulled in the latest changes from the Eigen trunk

15 files changed, 1012 insertions, 567 deletions

diff --git a/unsupported/Eigen/CXX11/Core b/unsupported/Eigen/CXX11/Core
index bba3d578d..f6c3b49bb 100644
--- a/unsupported/Eigen/CXX11/Core
+++ b/unsupported/Eigen/CXX11/Core
@@ -34,6 +34,7 @@
 #if __cplusplus <= 199711L
 #include "src/Core/util/EmulateCXX11Meta.h"
 #else
+#include <array>
 #include "src/Core/util/CXX11Workarounds.h"
 #include "src/Core/util/CXX11Meta.h"
 #endif
diff --git a/unsupported/Eigen/MPRealSupport b/unsupported/Eigen/MPRealSupport
index d4b03647d..632de3854 100644
--- a/unsupported/Eigen/MPRealSupport
+++ b/unsupported/Eigen/MPRealSupport
@@ -27,6 +27,8 @@ namespace Eigen {
   * via the <a href="http://www.holoborodko.com/pavel/mpfr">MPFR C++</a>
   * library which itself is built upon <a href="http://www.mpfr.org/">MPFR</a>/<a href="http://gmplib.org/">GMP</a>.
   *
+  * \warning MPFR C++ is licensed under the GPL.
+  *
   * You can find a copy of MPFR C++ that is known to be compatible in the unsupported/test/mpreal folder.
   *
   * Here is an example:
@@ -86,9 +88,9 @@ int main()
     inline static Real epsilon  (const Real& x)                                         {    return mpfr::machine_epsilon(x); }
 
     inline static Real dummy_precision()   
-    { 
-        unsigned int weak_prec = ((mpfr::mpreal::get_default_prec()-1) * 90) / 100;
-        return mpfr::machine_epsilon(weak_prec);
+    {
+      mpfr_prec_t weak_prec = ((mpfr::mpreal::get_default_prec()-1) * 90) / 100;
+      return mpfr::machine_epsilon(weak_prec);
     }
   };
 
@@ -139,64 +141,53 @@ int main()
     public:
       typedef mpfr::mpreal ResScalar;
       enum {
-        nr = 2, // must be 2 for proper packing...
+        nr = 1,
         mr = 1,
-        WorkSpaceFactor = nr,
         LhsProgress = 1,
         RhsProgress = 1
       };
     };
 
-    template<typename Index, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
-    struct gebp_kernel<mpfr::mpreal,mpfr::mpreal,Index,mr,nr,ConjugateLhs,ConjugateRhs>
+    template<typename Index, bool ConjugateLhs, bool ConjugateRhs>
+    struct gebp_kernel<mpfr::mpreal,mpfr::mpreal,Index,1,1,ConjugateLhs,ConjugateRhs>
     {
       typedef mpfr::mpreal mpreal;
 
       EIGEN_DONT_INLINE
       void operator()(mpreal* res, Index resStride, const mpreal* blockA, const mpreal* blockB, Index rows, Index depth, Index cols, mpreal alpha,
-                      Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0, mpreal* /*unpackedB*/ = 0)
+                      Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0)
       {
-        mpreal acc1, acc2, tmp;
+        if(rows==0 || cols==0 || depth==0)
+          return;
+        
+        mpreal  acc1(0,mpfr_get_prec(blockA[0].mpfr_srcptr())),
+                tmp (0,mpfr_get_prec(blockA[0].mpfr_srcptr()));
         
         if(strideA==-1) strideA = depth;
         if(strideB==-1) strideB = depth;
 
-        for(Index j=0; j<cols; j+=nr)
+        for(Index i=0; i<rows; ++i)
         {
-          Index actual_nr = (std::min<Index>)(nr,cols-j);
-          mpreal *C1 = res + j*resStride;
-          mpreal *C2 = res + (j+1)*resStride;
-          for(Index i=0; i<rows; i++)
+          for(Index j=0; j<cols; ++j)
           {
-            mpreal *B = const_cast<mpreal*>(blockB) + j*strideB + offsetB*actual_nr;
-            mpreal *A = const_cast<mpreal*>(blockA) + i*strideA + offsetA;
+            mpreal *C1 = res + j*resStride;
+            
+            const mpreal *A = blockA + i*strideA + offsetA;
+            const mpreal *B = blockB + j*strideB + offsetB;
+            
             acc1 = 0;
-            acc2 = 0;
             for(Index k=0; k<depth; k++)
             {
-              mpfr_mul(tmp.mpfr_ptr(), A[k].mpfr_ptr(), B[0].mpfr_ptr(), mpreal::get_default_rnd());
+              mpfr_mul(tmp.mpfr_ptr(), A[k].mpfr_srcptr(), B[k].mpfr_srcptr(), mpreal::get_default_rnd());
               mpfr_add(acc1.mpfr_ptr(), acc1.mpfr_ptr(), tmp.mpfr_ptr(),  mpreal::get_default_rnd());
-              
-              if(actual_nr==2) {
-                mpfr_mul(tmp.mpfr_ptr(), A[k].mpfr_ptr(), B[1].mpfr_ptr(), mpreal::get_default_rnd());
-                mpfr_add(acc2.mpfr_ptr(), acc2.mpfr_ptr(), tmp.mpfr_ptr(),  mpreal::get_default_rnd());
-              }
-              
-              B+=actual_nr;
             }
             
-            mpfr_mul(acc1.mpfr_ptr(), acc1.mpfr_ptr(), alpha.mpfr_ptr(), mpreal::get_default_rnd());
-            mpfr_add(C1[i].mpfr_ptr(), C1[i].mpfr_ptr(), acc1.mpfr_ptr(),  mpreal::get_default_rnd());
-            
-            if(actual_nr==2) {
-              mpfr_mul(acc2.mpfr_ptr(), acc2.mpfr_ptr(), alpha.mpfr_ptr(), mpreal::get_default_rnd());
-              mpfr_add(C2[i].mpfr_ptr(), C2[i].mpfr_ptr(), acc2.mpfr_ptr(),  mpreal::get_default_rnd());
-            }
+            mpfr_mul(acc1.mpfr_ptr(), acc1.mpfr_srcptr(), alpha.mpfr_srcptr(), mpreal::get_default_rnd());
+            mpfr_add(C1[i].mpfr_ptr(), C1[i].mpfr_srcptr(), acc1.mpfr_srcptr(),  mpreal::get_default_rnd());
           }
         }
       }
     };
-
   } // end namespace internal
 }
 
diff --git a/unsupported/Eigen/SVD b/unsupported/Eigen/SVD
index 7cc059280..900a8aa60 100644
--- a/unsupported/Eigen/SVD
+++ b/unsupported/Eigen/SVD
@@ -29,10 +29,6 @@
 #include "../../Eigen/src/SVD/JacobiSVD_MKL.h"
 #endif
 
-#ifdef EIGEN2_SUPPORT
-#include "../../Eigen/src/Eigen2Support/SVD.h"
-#endif
-
 #include "../../Eigen/src/Core/util/ReenableStupidWarnings.h"
 
 #endif // EIGEN_SVD_MODULE_H
diff --git a/unsupported/Eigen/src/IterativeSolvers/GMRES.h b/unsupported/Eigen/src/IterativeSolvers/GMRES.h
index 073367506..c8c84069e 100644
--- a/unsupported/Eigen/src/IterativeSolvers/GMRES.h
+++ b/unsupported/Eigen/src/IterativeSolvers/GMRES.h
@@ -2,7 +2,7 @@
 // for linear algebra.
 //
 // Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
-// Copyright (C) 2012 Kolja Brix <brix@igpm.rwth-aaachen.de>
+// Copyright (C) 2012, 2014 Kolja Brix <brix@igpm.rwth-aaachen.de>
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
@@ -72,16 +72,20 @@ bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Precondition
 
 	VectorType p0 = rhs - mat*x;
 	VectorType r0 = precond.solve(p0);
-// 	RealScalar r0_sqnorm = r0.squaredNorm();
+ 
+	// is initial guess already good enough?
+	if(abs(r0.norm()) < tol) {
+		return true; 
+	}
 
 	VectorType w = VectorType::Zero(restart + 1);
 
-	FMatrixType H = FMatrixType::Zero(m, restart + 1);
+	FMatrixType H = FMatrixType::Zero(m, restart + 1); // Hessenberg matrix
 	VectorType tau = VectorType::Zero(restart + 1);
 	std::vector < JacobiRotation < Scalar > > G(restart);
 
 	// generate first Householder vector
-	VectorType e;
+	VectorType e(m-1);
 	RealScalar beta;
 	r0.makeHouseholder(e, tau.coeffRef(0), beta);
 	w(0)=(Scalar) beta;
diff --git a/unsupported/Eigen/src/IterativeSolvers/Scaling.h b/unsupported/Eigen/src/IterativeSolvers/Scaling.h
index 4fd439202..d113e6e90 100644
--- a/unsupported/Eigen/src/IterativeSolvers/Scaling.h
+++ b/unsupported/Eigen/src/IterativeSolvers/Scaling.h
@@ -9,6 +9,9 @@
 
 #ifndef EIGEN_ITERSCALING_H
 #define EIGEN_ITERSCALING_H
+
+namespace Eigen {
+
 /**
   * \ingroup IterativeSolvers_Module
   * \brief iterative scaling algorithm to equilibrate rows and column norms in matrices
@@ -41,8 +44,6 @@
   * 
   * \sa \ref IncompleteLUT 
   */
-namespace Eigen {
-using std::abs; 
 template<typename _MatrixType>
 class IterScaling
 {
@@ -71,6 +72,7 @@ class IterScaling
      */
     void compute (const MatrixType& mat)
     {
+      using std::abs;
       int m = mat.rows(); 
       int n = mat.cols();
       eigen_assert((m>0 && m == n) && "Please give a non - empty matrix");
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h
index 05df5a102..160120d03 100644
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h
@@ -176,8 +176,8 @@ void matrix_exp_pade17(const MatrixType &A, MatrixType &U, MatrixType &V)
   const MatrixType A4 = A2 * A2;
   const MatrixType A6 = A4 * A2;
   const MatrixType A8 = A4 * A4;
-  V = b[17] * m_tmp1 + b[15] * A6 + b[13] * A4 + b[11] * A2; // used for temporary storage
-  matrixType tmp = A8 * V;
+  V = b[17] * A8 + b[15] * A6 + b[13] * A4 + b[11] * A2; // used for temporary storage
+  MatrixType tmp = A8 * V;
   tmp += b[9] * A8 + b[7] * A6 + b[5] * A4 + b[3] * A2 
     + b[1] * MatrixType::Identity(A.rows(), A.cols());
   U.noalias() = A * tmp;
diff --git a/unsupported/Eigen/src/SparseExtra/MarketIO.h b/unsupported/Eigen/src/SparseExtra/MarketIO.h
index 1c40d3f7c..25ff4228d 100644
--- a/unsupported/Eigen/src/SparseExtra/MarketIO.h
+++ b/unsupported/Eigen/src/SparseExtra/MarketIO.h
@@ -133,6 +133,7 @@ template<typename SparseMatrixType>
 bool loadMarket(SparseMatrixType& mat, const std::string& filename)
 {
   typedef typename SparseMatrixType::Scalar Scalar;
+  typedef typename SparseMatrixType::Index Index;
   std::ifstream input(filename.c_str(),std::ios::in);
   if(!input)
     return false;
@@ -142,11 +143,11 @@ bool loadMarket(SparseMatrixType& mat, const std::string& filename)
   
   bool readsizes = false;
 
-  typedef Triplet<Scalar,int> T;
+  typedef Triplet<Scalar,Index> T;
   std::vector<T> elements;
   
-  int M(-1), N(-1), NNZ(-1);
-  int count = 0;
+  Index M(-1), N(-1), NNZ(-1);
+  Index count = 0;
   while(input.getline(buffer, maxBuffersize))
   {
     // skip comments   
@@ -169,7 +170,7 @@ bool loadMarket(SparseMatrixType& mat, const std::string& filename)
     }
     else
     { 
-      int i(-1), j(-1);
+      Index i(-1), j(-1);
       Scalar value; 
       if( internal::GetMarketLine(line, M, N, i, j, value) ) 
       {
diff --git a/unsupported/Eigen/src/Splines/Spline.h b/unsupported/Eigen/src/Splines/Spline.h
index 1b47992d6..c46f728bc 100644
--- a/unsupported/Eigen/src/Splines/Spline.h
+++ b/unsupported/Eigen/src/Splines/Spline.h
@@ -44,9 +44,15 @@ namespace Eigen
     
     /** \brief The data type used to store knot vectors. */
     typedef typename SplineTraits<Spline>::KnotVectorType KnotVectorType;
+
+    /** \brief The data type used to store parameter vectors. */
+    typedef typename SplineTraits<Spline>::ParameterVectorType ParameterVectorType;
     
     /** \brief The data type used to store non-zero basis functions. */
     typedef typename SplineTraits<Spline>::BasisVectorType BasisVectorType;
+
+    /** \brief The data type used to store the values of the basis function derivatives. */
+    typedef typename SplineTraits<Spline>::BasisDerivativeType BasisDerivativeType;
     
     /** \brief The data type representing the spline's control points. */
     typedef typename SplineTraits<Spline>::ControlPointVectorType ControlPointVectorType;
@@ -203,10 +209,25 @@ namespace Eigen
      **/
     static BasisVectorType BasisFunctions(Scalar u, DenseIndex degree, const KnotVectorType& knots);
 
+    /**
+     * \copydoc Spline::basisFunctionDerivatives
+     * \param degree The degree of the underlying spline
+     * \param knots The underlying spline's knot vector.
+     **/    
+    static BasisDerivativeType BasisFunctionDerivatives(
+      const Scalar u, const DenseIndex order, const DenseIndex degree, const KnotVectorType& knots);
 
   private:
     KnotVectorType m_knots; /*!< Knot vector. */
     ControlPointVectorType  m_ctrls; /*!< Control points. */
+
+    template <typename DerivativeType>
+    static void BasisFunctionDerivativesImpl(
+      const typename Spline<_Scalar, _Dim, _Degree>::Scalar u,
+      const DenseIndex order,
+      const DenseIndex p, 
+      const typename Spline<_Scalar, _Dim, _Degree>::KnotVectorType& U,
+      DerivativeType& N_);
   };
 
   template <typename _Scalar, int _Dim, int _Degree>
@@ -345,20 +366,24 @@ namespace Eigen
   }
 
   /* --------------------------------------------------------------------------------------------- */
-
-  template <typename SplineType, typename DerivativeType>
-  void basisFunctionDerivativesImpl(const SplineType& spline, typename SplineType::Scalar u, DenseIndex order, DerivativeType& N_)
+  
+  
+  template <typename _Scalar, int _Dim, int _Degree>
+  template <typename DerivativeType>
+  void Spline<_Scalar, _Dim, _Degree>::BasisFunctionDerivativesImpl(
+    const typename Spline<_Scalar, _Dim, _Degree>::Scalar u,
+    const DenseIndex order,
+    const DenseIndex p, 
+    const typename Spline<_Scalar, _Dim, _Degree>::KnotVectorType& U,
+    DerivativeType& N_)
   {
+    typedef Spline<_Scalar, _Dim, _Degree> SplineType;
     enum { Order = SplineTraits<SplineType>::OrderAtCompileTime };
 
     typedef typename SplineTraits<SplineType>::Scalar Scalar;
     typedef typename SplineTraits<SplineType>::BasisVectorType BasisVectorType;
-    typedef typename SplineTraits<SplineType>::KnotVectorType KnotVectorType;
-
-    const KnotVectorType& U = spline.knots();
-
-    const DenseIndex p = spline.degree();
-    const DenseIndex span = spline.span(u);
+  
+    const DenseIndex span = SplineType::Span(u, p, U);
 
     const DenseIndex n = (std::min)(p, order);
 
@@ -455,8 +480,8 @@ namespace Eigen
   typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::BasisDerivativeType
     Spline<_Scalar, _Dim, _Degree>::basisFunctionDerivatives(Scalar u, DenseIndex order) const
   {
-    typename SplineTraits< Spline >::BasisDerivativeType der;
-    basisFunctionDerivativesImpl(*this, u, order, der);
+    typename SplineTraits<Spline<_Scalar, _Dim, _Degree> >::BasisDerivativeType der;
+    BasisFunctionDerivativesImpl(u, order, degree(), knots(), der);
     return der;
   }
 
@@ -465,8 +490,21 @@ namespace Eigen
   typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::BasisDerivativeType
     Spline<_Scalar, _Dim, _Degree>::basisFunctionDerivatives(Scalar u, DenseIndex order) const
   {
-    typename SplineTraits< Spline, DerivativeOrder >::BasisDerivativeType der;
-    basisFunctionDerivativesImpl(*this, u, order, der);
+    typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::BasisDerivativeType der;
+    BasisFunctionDerivativesImpl(u, order, degree(), knots(), der);
+    return der;
+  }
+
+  template <typename _Scalar, int _Dim, int _Degree>
+  typename SplineTraits<Spline<_Scalar, _Dim, _Degree> >::BasisDerivativeType
+  Spline<_Scalar, _Dim, _Degree>::BasisFunctionDerivatives(
+    const typename Spline<_Scalar, _Dim, _Degree>::Scalar u,
+    const DenseIndex order,
+    const DenseIndex degree,
+    const typename Spline<_Scalar, _Dim, _Degree>::KnotVectorType& knots)
+  {
+    typename SplineTraits<Spline>::BasisDerivativeType der;
+    BasisFunctionDerivativesImpl(u, order, degree, knots, der);
     return der;
   }
 }
diff --git a/unsupported/Eigen/src/Splines/SplineFitting.h b/unsupported/Eigen/src/Splines/SplineFitting.h
index 0265d532c..d3c245fa9 100644
--- a/unsupported/Eigen/src/Splines/SplineFitting.h
+++ b/unsupported/Eigen/src/Splines/SplineFitting.h
@@ -10,10 +10,14 @@
 #ifndef EIGEN_SPLINE_FITTING_H
 #define EIGEN_SPLINE_FITTING_H
 
+#include <algorithm>
+#include <functional>
 #include <numeric>
+#include <vector>
 
 #include "SplineFwd.h"
 
+#include <Eigen/LU>
 #include <Eigen/QR>
 
 namespace Eigen
@@ -50,6 +54,129 @@ namespace Eigen
   }
 
   /**
+   * \brief Computes knot averages when derivative constraints are present.
+   * Note that this is a technical interpretation of the referenced article
+   * since the algorithm contained therein is incorrect as written.
+   * \ingroup Splines_Module
+   *
+   * \param[in] parameters The parameters at which the interpolation B-Spline
+   *            will intersect the given interpolation points. The parameters
+   *            are assumed to be a non-decreasing sequence.
+   * \param[in] degree The degree of the interpolating B-Spline. This must be
+   *            greater than zero.
+   * \param[in] derivativeIndices The indices corresponding to parameters at
+   *            which there are derivative constraints. The indices are assumed
+   *            to be a non-decreasing sequence.
+   * \param[out] knots The calculated knot vector. These will be returned as a
+   *             non-decreasing sequence
+   *
+   * \sa Les A. Piegl, Khairan Rajab, Volha Smarodzinana. 2008.
+   * Curve interpolation with directional constraints for engineering design. 
+   * Engineering with Computers
+   **/
+  template <typename KnotVectorType, typename ParameterVectorType, typename IndexArray>
+  void KnotAveragingWithDerivatives(const ParameterVectorType& parameters,
+                                    const unsigned int degree,
+                                    const IndexArray& derivativeIndices,
+                                    KnotVectorType& knots)
+  {
+    typedef typename ParameterVectorType::Scalar Scalar;
+
+    DenseIndex numParameters = parameters.size();
+    DenseIndex numDerivatives = derivativeIndices.size();
+
+    if (numDerivatives < 1)
+    {
+      KnotAveraging(parameters, degree, knots);
+      return;
+    }
+
+    DenseIndex startIndex;
+    DenseIndex endIndex;
+  
+    DenseIndex numInternalDerivatives = numDerivatives;
+    
+    if (derivativeIndices[0] == 0)
+    {
+      startIndex = 0;
+      --numInternalDerivatives;
+    }
+    else
+    {
+      startIndex = 1;
+    }
+    if (derivativeIndices[numDerivatives - 1] == numParameters - 1)
+    {
+      endIndex = numParameters - degree;
+      --numInternalDerivatives;
+    }
+    else
+    {
+      endIndex = numParameters - degree - 1;
+    }
+
+    // There are (endIndex - startIndex + 1) knots obtained from the averaging
+    // and 2 for the first and last parameters.
+    DenseIndex numAverageKnots = endIndex - startIndex + 3;
+    KnotVectorType averageKnots(numAverageKnots);
+    averageKnots[0] = parameters[0];
+
+    int newKnotIndex = 0;
+    for (DenseIndex i = startIndex; i <= endIndex; ++i)
+      averageKnots[++newKnotIndex] = parameters.segment(i, degree).mean();
+    averageKnots[++newKnotIndex] = parameters[numParameters - 1];
+
+    newKnotIndex = -1;
+  
+    ParameterVectorType temporaryParameters(numParameters + 1);
+    KnotVectorType derivativeKnots(numInternalDerivatives);
+    for (unsigned int i = 0; i < numAverageKnots - 1; ++i)
+    {
+      temporaryParameters[0] = averageKnots[i];
+      ParameterVectorType parameterIndices(numParameters);
+      int temporaryParameterIndex = 1;
+      for (int j = 0; j < numParameters; ++j)
+      {
+        Scalar parameter = parameters[j];
+        if (parameter >= averageKnots[i] && parameter < averageKnots[i + 1])
+        {
+          parameterIndices[temporaryParameterIndex] = j;
+          temporaryParameters[temporaryParameterIndex++] = parameter;
+        }
+      }
+      temporaryParameters[temporaryParameterIndex] = averageKnots[i + 1];
+
+      for (int j = 0; j <= temporaryParameterIndex - 2; ++j)
+      {
+        for (DenseIndex k = 0; k < derivativeIndices.size(); ++k)
+        {
+          if (parameterIndices[j + 1] == derivativeIndices[k]
+              && parameterIndices[j + 1] != 0
+              && parameterIndices[j + 1] != numParameters - 1)
+          {
+            derivativeKnots[++newKnotIndex] = temporaryParameters.segment(j, 3).mean();
+            break;
+          }
+        }
+      }
+    }
+    
+    KnotVectorType temporaryKnots(averageKnots.size() + derivativeKnots.size());
+
+    std::merge(averageKnots.data(), averageKnots.data() + averageKnots.size(),
+               derivativeKnots.data(), derivativeKnots.data() + derivativeKnots.size(),
+               temporaryKnots.data());
+
+    // Number of control points (one for each point and derivative) plus spline order.
+    DenseIndex numKnots = numParameters + numDerivatives + degree + 1;
+    knots.resize(numKnots);
+
+    knots.head(degree).fill(temporaryKnots[0]);
+    knots.tail(degree).fill(temporaryKnots.template tail<1>()[0]);
+    knots.segment(degree, temporaryKnots.size()) = temporaryKnots;
+  }
+
+  /**
    * \brief Computes chord length parameters which are required for spline interpolation.
    * \ingroup Splines_Module
    *
@@ -86,6 +213,7 @@ namespace Eigen
   struct SplineFitting
   {
     typedef typename SplineType::KnotVectorType KnotVectorType;
+    typedef typename SplineType::ParameterVectorType ParameterVectorType;
 
     /**
      * \brief Fits an interpolating Spline to the given data points.
@@ -109,6 +237,52 @@ namespace Eigen
      **/
     template <typename PointArrayType>
     static SplineType Interpolate(const PointArrayType& pts, DenseIndex degree, const KnotVectorType& knot_parameters);
+
+    /**
+     * \brief Fits an interpolating spline to the given data points and
+     * derivatives.
+     * 
+     * \param points The points for which an interpolating spline will be computed.
+     * \param derivatives The desired derivatives of the interpolating spline at interpolation
+     *                    points.
+     * \param derivativeIndices An array indicating which point each derivative belongs to. This
+     *                          must be the same size as @a derivatives.
+     * \param degree The degree of the interpolating spline.
+     *
+     * \returns A spline interpolating @a points with @a derivatives at those points.
+     *
+     * \sa Les A. Piegl, Khairan Rajab, Volha Smarodzinana. 2008.
+     * Curve interpolation with directional constraints for engineering design. 
+     * Engineering with Computers
+     **/
+    template <typename PointArrayType, typename IndexArray>
+    static SplineType InterpolateWithDerivatives(const PointArrayType& points,
+                                                 const PointArrayType& derivatives,
+                                                 const IndexArray& derivativeIndices,
+                                                 const unsigned int degree);
+
+    /**
+     * \brief Fits an interpolating spline to the given data points and derivatives.
+     * 
+     * \param points The points for which an interpolating spline will be computed.
+     * \param derivatives The desired derivatives of the interpolating spline at interpolation points.
+     * \param derivativeIndices An array indicating which point each derivative belongs to. This
+     *                          must be the same size as @a derivatives.
+     * \param degree The degree of the interpolating spline.
+     * \param parameters The parameters corresponding to the interpolation points.
+     *
+     * \returns A spline interpolating @a points with @a derivatives at those points.
+     *
+     * \sa Les A. Piegl, Khairan Rajab, Volha Smarodzinana. 2008.
+     * Curve interpolation with directional constraints for engineering design. 
+     * Engineering with Computers
+     */
+    template <typename PointArrayType, typename IndexArray>
+    static SplineType InterpolateWithDerivatives(const PointArrayType& points,
+                                                 const PointArrayType& derivatives,
+                                                 const IndexArray& derivativeIndices,
+                                                 const unsigned int degree,
+                                                 const ParameterVectorType& parameters);
   };
 
   template <typename SplineType>
@@ -151,6 +325,106 @@ namespace Eigen
     ChordLengths(pts, chord_lengths);
     return Interpolate(pts, degree, chord_lengths);
   }
+  
+  template <typename SplineType>
+  template <typename PointArrayType, typename IndexArray>
+  SplineType 
+  SplineFitting<SplineType>::InterpolateWithDerivatives(const PointArrayType& points,
+                                                        const PointArrayType& derivatives,
+                                                        const IndexArray& derivativeIndices,
+                                                        const unsigned int degree,
+                                                        const ParameterVectorType& parameters)
+  {
+    typedef typename SplineType::KnotVectorType::Scalar Scalar;      
+    typedef typename SplineType::ControlPointVectorType ControlPointVectorType;
+
+    typedef Matrix<Scalar, Dynamic, Dynamic> MatrixType;
+
+    const DenseIndex n = points.cols() + derivatives.cols();
+    
+    KnotVectorType knots;
+
+    KnotAveragingWithDerivatives(parameters, degree, derivativeIndices, knots);
+    
+    // fill matrix
+    MatrixType A = MatrixType::Zero(n, n);
+
+    // Use these dimensions for quicker populating, then transpose for solving.
+    MatrixType b(points.rows(), n);
+
+    DenseIndex startRow;
+    DenseIndex derivativeStart;
+
+    // End derivatives.
+    if (derivativeIndices[0] == 0)
+    {
+      A.template block<1, 2>(1, 0) << -1, 1;
+      
+      Scalar y = (knots(degree + 1) - knots(0)) / degree;
+      b.col(1) = y*derivatives.col(0);
+      
+      startRow = 2;
+      derivativeStart = 1;
+    }
+    else
+    {
+      startRow = 1;
+      derivativeStart = 0;
+    }
+    if (derivativeIndices[derivatives.cols() - 1] == points.cols() - 1)
+    {
+      A.template block<1, 2>(n - 2, n - 2) << -1, 1;
+
+      Scalar y = (knots(knots.size() - 1) - knots(knots.size() - (degree + 2))) / degree;
+      b.col(b.cols() - 2) = y*derivatives.col(derivatives.cols() - 1);
+    }
+    
+    DenseIndex row = startRow;
+    DenseIndex derivativeIndex = derivativeStart;
+    for (DenseIndex i = 1; i < parameters.size() - 1; ++i)
+    {
+      const DenseIndex span = SplineType::Span(parameters[i], degree, knots);
+
+      if (derivativeIndices[derivativeIndex] == i)
+      {
+        A.block(row, span - degree, 2, degree + 1)
+          = SplineType::BasisFunctionDerivatives(parameters[i], 1, degree, knots);
+
+        b.col(row++) = points.col(i);
+        b.col(row++) = derivatives.col(derivativeIndex++);
+      }
+      else
+      {
+        A.row(row++).segment(span - degree, degree + 1)
+          = SplineType::BasisFunctions(parameters[i], degree, knots);
+      }
+    }
+    b.col(0) = points.col(0);
+    b.col(b.cols() - 1) = points.col(points.cols() - 1);
+    A(0,0) = 1;
+    A(n - 1, n - 1) = 1;
+    
+    // Solve
+    FullPivLU<MatrixType> lu(A);
+    ControlPointVectorType controlPoints = lu.solve(MatrixType(b.transpose())).transpose();
+
+    SplineType spline(knots, controlPoints);
+    
+    return spline;
+  }
+  
+  template <typename SplineType>
+  template <typename PointArrayType, typename IndexArray>
+  SplineType
+  SplineFitting<SplineType>::InterpolateWithDerivatives(const PointArrayType& points,
+                                                        const PointArrayType& derivatives,
+                                                        const IndexArray& derivativeIndices,
+                                                        const unsigned int degree)
+  {
+    ParameterVectorType parameters;
+    ChordLengths(points, parameters);
+    return InterpolateWithDerivatives(points, derivatives, derivativeIndices, degree, parameters);
+  }
 }
 
 #endif // EIGEN_SPLINE_FITTING_H
diff --git a/unsupported/Eigen/src/Splines/SplineFwd.h b/unsupported/Eigen/src/Splines/SplineFwd.h
index 9ea23a9a1..0a95fbf3e 100644
--- a/unsupported/Eigen/src/Splines/SplineFwd.h
+++ b/unsupported/Eigen/src/Splines/SplineFwd.h
@@ -48,6 +48,9 @@ namespace Eigen
       
       /** \brief The data type used to store knot vectors. */
       typedef Array<Scalar,1,Dynamic> KnotVectorType;
+
+      /** \brief The data type used to store parameter vectors. */
+      typedef Array<Scalar,1,Dynamic> ParameterVectorType;
       
       /** \brief The data type representing the spline's control points. */
       typedef Array<Scalar,Dimension,Dynamic> ControlPointVectorType;
diff --git a/unsupported/doc/Overview.dox b/unsupported/doc/Overview.dox
index d048377df..45464a545 100644
--- a/unsupported/doc/Overview.dox
+++ b/unsupported/doc/Overview.dox
@@ -1,14 +1,15 @@
+/// \brief Namespace containing all symbols from the %Eigen library.
 namespace Eigen {
 
-/** \mainpage Eigen's unsupported modules
+/** \mainpage %Eigen's unsupported modules
 
-This is the API documentation for Eigen's unsupported modules.
+This is the API documentation for %Eigen's unsupported modules.
 
 These modules are contributions from various users. They are provided "as is", without any support.
 
 Click on the \e Modules tab at the top of this page to get a list of all unsupported modules.
 
-Don't miss the <a href="..//index.html">official Eigen documentation</a>.
+Don't miss the <a href="../index.html">official Eigen documentation</a>.
 
 */
 
@@ -18,8 +19,10 @@ Don't miss the <a href="..//index.html">official Eigen documentation</a>.
 
 The unsupported modules are contributions from various users. They are
 provided "as is", without any support. Nevertheless, some of them are
-subject to be included in Eigen in the future.
+subject to be included in %Eigen in the future.
 
 */
 
+/// \internal \brief Namespace containing low-level routines from the %Eigen library.
+namespace internal {}
 }
diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt
index 7458128fb..406564673 100644
--- a/unsupported/test/CMakeLists.txt
+++ b/unsupported/test/CMakeLists.txt
@@ -100,7 +100,7 @@ if(EIGEN_TEST_CXX11)
   ei_add_test(cxx11_meta "-std=c++0x")
   ei_add_test(cxx11_tensor_simple "-std=c++0x")
   ei_add_test(cxx11_tensor_symmetry "-std=c++0x")
-  ei_add_test(cxx11_tensor_assign "-std=c++0x")
+#  ei_add_test(cxx11_tensor_assign "-std=c++0x")
   ei_add_test(cxx11_tensor_comparisons "-std=c++0x")
   ei_add_test(cxx11_tensor_contraction "-std=c++0x")
   ei_add_test(cxx11_tensor_convolution "-std=c++0x")
@@ -108,7 +108,7 @@ if(EIGEN_TEST_CXX11)
 #  ei_add_test(cxx11_tensor_fixed_size "-std=c++0x")
   ei_add_test(cxx11_tensor_lvalue "-std=c++0x")
   ei_add_test(cxx11_tensor_map "-std=c++0x")
-  ei_add_test(cxx11_tensor_morphing "-std=c++0x")
+#  ei_add_test(cxx11_tensor_morphing "-std=c++0x")
 #  ei_add_test(cxx11_tensor_device  "-std=c++0x")
 #  ei_add_test(cxx11_tensor_fixed_size "-std=c++0x")
 #  ei_add_test(cxx11_tensor_thread_pool "-std=c++0x")
diff --git a/unsupported/test/mpreal/mpreal.h b/unsupported/test/mpreal/mpreal.h
index 70d37ab71..dddda7dd9 100644
--- a/unsupported/test/mpreal/mpreal.h
+++ b/unsupported/test/mpreal/mpreal.h
@@ -1,59 +1,47 @@
 /*

-    Multi-precision floating point number class for C++. 

+    MPFR C++: Multi-precision floating point number class for C++. 

     Based on MPFR library:    http://mpfr.org

 

     Project homepage:    http://www.holoborodko.com/pavel/mpfr

     Contact e-mail:      pavel@holoborodko.com

 

-    Copyright (c) 2008-2013 Pavel Holoborodko

+    Copyright (c) 2008-2014 Pavel Holoborodko

 

     Contributors:

     Dmitriy Gubanov, Konstantin Holoborodko, Brian Gladman, 

     Helmut Jarausch, Fokko Beekhof, Ulrich Mutze, Heinz van Saanen, 

     Pere Constans, Peter van Hoof, Gael Guennebaud, Tsai Chia Cheng, 

-    Alexei Zubanov, Jauhien Piatlicki, Victor Berger, John Westwood.

+    Alexei Zubanov, Jauhien Piatlicki, Victor Berger, John Westwood,

+    Petr Aleksandrov, Orion Poplawski, Charles Karney.

 

-    ****************************************************************************

-    This library is free software; you can redistribute it and/or

-    modify it under the terms of the GNU Lesser General Public

-    License as published by the Free Software Foundation; either

-    version 2.1 of the License, or (at your option) any later version.

-

-    This library is distributed in the hope that it will be useful,

-    but WITHOUT ANY WARRANTY; without even the implied warranty of

-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

-    Lesser General Public License for more details.

+    Licensing:

+    (A) MPFR C++ is under GNU General Public License ("GPL").

+    

+    (B) Non-free licenses may also be purchased from the author, for users who 

+        do not want their programs protected by the GPL.

 

-    You should have received a copy of the GNU Lesser General Public

-    License along with this library; if not, write to the Free Software

-    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA

+        The non-free licenses are for users that wish to use MPFR C++ in 

+        their products but are unwilling to release their software 

+        under the GPL (which would require them to release source code 

+        and allow free redistribution).

 

-    ****************************************************************************

-    Redistribution and use in source and binary forms, with or without

-    modification, are permitted provided that the following conditions

-    are met:

-    

-    1. Redistributions of source code must retain the above copyright

-    notice, this list of conditions and the following disclaimer.

-    

-    2. Redistributions in binary form must reproduce the above copyright

-    notice, this list of conditions and the following disclaimer in the

-    documentation and/or other materials provided with the distribution.

+        Such users can purchase an unlimited-use license from the author.

+        Contact us for more details.

     

-    3. The name of the author may not be used to endorse or promote products

-    derived from this software without specific prior written permission.

-

-    THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND

-    ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE

-    IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE

-    ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE

-    FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL

-    DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS

-    OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)

-    HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT

-    LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY

-    OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF

-    SUCH DAMAGE.

+    GNU General Public License ("GPL") copyright permissions statement:

+    **************************************************************************

+    This program is free software: you can redistribute it and/or modify

+    it under the terms of the GNU General Public License as published by

+    the Free Software Foundation, either version 3 of the License, or

+    (at your option) any later version.

+

+    This program is distributed in the hope that it will be useful,

+    but WITHOUT ANY WARRANTY; without even the implied warranty of

+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

+    GNU General Public License for more details.

+

+    You should have received a copy of the GNU General Public License

+    along with this program.  If not, see <http://www.gnu.org/licenses/>.

 */

 

 #ifndef __MPREAL_H__

@@ -71,6 +59,15 @@
 // Options

 #define MPREAL_HAVE_INT64_SUPPORT               // Enable int64_t support if possible. Available only for MSVC 2010 & GCC. 

 #define MPREAL_HAVE_MSVC_DEBUGVIEW              // Enable Debugger Visualizer for "Debug" builds in MSVC.

+#define MPREAL_HAVE_DYNAMIC_STD_NUMERIC_LIMITS  // Enable extended std::numeric_limits<mpfr::mpreal> specialization.

+                                                // Meaning that "digits", "round_style" and similar members are defined as functions, not constants.

+                                                // See std::numeric_limits<mpfr::mpreal> at the end of the file for more information.

+

+// Library version

+#define MPREAL_VERSION_MAJOR 3

+#define MPREAL_VERSION_MINOR 5

+#define MPREAL_VERSION_PATCHLEVEL 9

+#define MPREAL_VERSION_STRING "3.5.9"

 

 // Detect compiler using signatures from http://predef.sourceforge.net/

 #if defined(__GNUC__) && defined(__INTEL_COMPILER)

@@ -83,17 +80,29 @@
     #define IsInf(x) std::isinf(x)              // GNU C/C++ (and/or other compilers), just hope for C99 conformance

 #endif

 

-// Detect support for r-value references (move semantic). Borrowed from Eigen.

 // A Clang feature extension to determine compiler features.

-// We use it to determine 'cxx_rvalue_references'

 #ifndef __has_feature

-# define __has_feature(x) 0

+    #define __has_feature(x) 0

 #endif

 

+// Detect support for r-value references (move semantic). Borrowed from Eigen.

 #if (__has_feature(cxx_rvalue_references) || \

-    defined(__GXX_EXPERIMENTAL_CXX0X__) || \

-    (defined(_MSC_VER) && _MSC_VER >= 1600))

-#define MPREAL_HAVE_MOVE_SUPPORT

+       defined(__GXX_EXPERIMENTAL_CXX0X__) || __cplusplus >= 201103L || \

+      (defined(_MSC_VER) && _MSC_VER >= 1600))

+

+    #define MPREAL_HAVE_MOVE_SUPPORT

+

+    // Use fields in mpfr_t structure to check if it was initialized / set dummy initialization 

+    #define mpfr_is_initialized(x)      (0 != (x)->_mpfr_d)

+    #define mpfr_set_uninitialized(x)   ((x)->_mpfr_d = 0 )

+#endif

+

+// Detect support for explicit converters. 

+#if (__has_feature(cxx_explicit_conversions) || \

+       defined(__GXX_EXPERIMENTAL_CXX0X__) || __cplusplus >= 201103L || \

+      (defined(_MSC_VER) && _MSC_VER >= 1800))

+

+    #define MPREAL_HAVE_EXPLICIT_CONVERTERS

 #endif

 

 // Detect available 64-bit capabilities

@@ -112,7 +121,7 @@
         #endif

 

     #elif defined (__GNUC__) && defined(__linux__)

-        #if defined(__amd64__) || defined(__amd64) || defined(__x86_64__) || defined(__x86_64) || defined(__ia64) || defined(__itanium__) || defined(_M_IA64)

+        #if defined(__amd64__) || defined(__amd64) || defined(__x86_64__) || defined(__x86_64) || defined(__ia64) || defined(__itanium__) || defined(_M_IA64) || defined (__PPC64__)

             #undef MPREAL_HAVE_INT64_SUPPORT    // Remove all shaman dances for x64 builds since

             #undef MPFR_USE_INTMAX_T            // GCC already supports x64 as of "long int" is 64-bit integer, nothing left to do

         #else

@@ -175,7 +184,7 @@ public:
     mpreal(const int u,               mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd());

     

     // Construct mpreal from mpfr_t structure.

-    // shared = true allows to avoid deep copy, so that mpreal and 'u' shared same data & pointers.    

+    // shared = true allows to avoid deep copy, so that mpreal and 'u' share the same data & pointers.    

     mpreal(const mpfr_t  u, bool shared = false);   

 

 #if defined (MPREAL_HAVE_INT64_SUPPORT)

@@ -315,14 +324,34 @@ public:
     friend bool operator == (const mpreal& a, const double b);

 

     // Type Conversion operators

+    bool            toBool      (mp_rnd_t mode = GMP_RNDZ)    const;

     long            toLong      (mp_rnd_t mode = GMP_RNDZ)    const;

     unsigned long   toULong     (mp_rnd_t mode = GMP_RNDZ)    const;

+    float           toFloat     (mp_rnd_t mode = GMP_RNDN)    const;

     double          toDouble    (mp_rnd_t mode = GMP_RNDN)    const;

     long double     toLDouble   (mp_rnd_t mode = GMP_RNDN)    const;

 

+#if defined (MPREAL_HAVE_EXPLICIT_CONVERTERS)

+    explicit operator bool               () const { return toBool();       }

+    explicit operator int                () const { return toLong();       }

+    explicit operator long               () const { return toLong();       }

+    explicit operator long long          () const { return toLong();       }

+    explicit operator unsigned           () const { return toULong();      }

+    explicit operator unsigned long      () const { return toULong();      }

+    explicit operator unsigned long long () const { return toULong();      }

+    explicit operator float              () const { return toFloat();      }

+    explicit operator double             () const { return toDouble();     }

+    explicit operator long double        () const { return toLDouble();    }

+#endif

+

 #if defined (MPREAL_HAVE_INT64_SUPPORT)

     int64_t         toInt64     (mp_rnd_t mode = GMP_RNDZ)    const;

     uint64_t        toUInt64    (mp_rnd_t mode = GMP_RNDZ)    const;

+

+    #if defined (MPREAL_HAVE_EXPLICIT_CONVERTERS)

+    explicit operator int64_t   () const { return toInt64();      }

+    explicit operator uint64_t  () const { return toUInt64();     }

+    #endif

 #endif

 

     // Get raw pointers so that mpreal can be directly used in raw mpfr_* functions

@@ -331,122 +360,124 @@ public:
     ::mpfr_srcptr mpfr_srcptr() const;

 

     // Convert mpreal to string with n significant digits in base b

-    // n = 0 -> convert with the maximum available digits 

-    std::string        toString(int n = 0, int b = 10, mp_rnd_t mode = mpreal::get_default_rnd()) const;

+    // n = -1 -> convert with the maximum available digits

+    std::string toString(int n = -1, int b = 10, mp_rnd_t mode = mpreal::get_default_rnd()) const;

 

 #if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))

-    std::string        toString(const std::string& format) const;

+    std::string toString(const std::string& format) const;

 #endif

 

+    std::ostream& output(std::ostream& os) const;

+

     // Math Functions

-    friend const mpreal sqr (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal sqrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal sqrt(const unsigned long int v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal cbrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal root(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal pow (const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal pow (const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal pow (const mpreal& a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal pow (const mpreal& a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal pow (const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal pow (const unsigned long int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal fabs(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-

-    friend const mpreal abs(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

+    friend const mpreal sqr (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal sqrt(const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal sqrt(const unsigned long int v, mp_rnd_t rnd_mode);

+    friend const mpreal cbrt(const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal root(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode);

+    friend const mpreal pow (const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode);

+    friend const mpreal pow (const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode);

+    friend const mpreal pow (const mpreal& a, const unsigned long int b, mp_rnd_t rnd_mode);

+    friend const mpreal pow (const mpreal& a, const long int b, mp_rnd_t rnd_mode);

+    friend const mpreal pow (const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode);

+    friend const mpreal pow (const unsigned long int a, const unsigned long int b, mp_rnd_t rnd_mode);

+    friend const mpreal fabs(const mpreal& v, mp_rnd_t rnd_mode);

+

+    friend const mpreal abs(const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode);

+    friend inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode);

+    friend inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode);

+    friend inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode);

+    friend inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode);

     friend int cmpabs(const mpreal& a,const mpreal& b);

     

-    friend const mpreal log  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal log2 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal log10(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal exp  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 

-    friend const mpreal exp2 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal exp10(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal log1p(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal expm1(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-

-    friend const mpreal cos(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal sin(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal tan(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal sec(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal csc(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal cot(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-

-    friend const mpreal acos  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal asin  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal atan  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal acot  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal asec  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal acsc  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-

-    friend const mpreal cosh  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal sinh  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal tanh  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal sech  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal csch  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal coth  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal acosh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal asinh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal atanh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal acoth (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal asech (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal acsch (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-

-    friend const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-

-    friend const mpreal fac_ui (unsigned long int v,  mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal eint   (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-

-    friend const mpreal gamma    (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal lngamma  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal lgamma   (const mpreal& v, int *signp = 0, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal zeta     (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal erf      (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal erfc     (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal besselj0 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 

-    friend const mpreal besselj1 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 

-    friend const mpreal besseljn (long n, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal bessely0 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal bessely1 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal besselyn (long n, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 

-    friend const mpreal fma      (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal fms      (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal agm      (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal sum      (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

+    friend const mpreal log  (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal log2 (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal log10(const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal exp  (const mpreal& v, mp_rnd_t rnd_mode); 

+    friend const mpreal exp2 (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal exp10(const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal log1p(const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal expm1(const mpreal& v, mp_rnd_t rnd_mode);

+

+    friend const mpreal cos(const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal sin(const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal tan(const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal sec(const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal csc(const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal cot(const mpreal& v, mp_rnd_t rnd_mode);

+    friend int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode);

+

+    friend const mpreal acos  (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal asin  (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal atan  (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode);

+    friend const mpreal acot  (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal asec  (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal acsc  (const mpreal& v, mp_rnd_t rnd_mode);

+

+    friend const mpreal cosh  (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal sinh  (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal tanh  (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal sech  (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal csch  (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal coth  (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal acosh (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal asinh (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal atanh (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal acoth (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal asech (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal acsch (const mpreal& v, mp_rnd_t rnd_mode);

+

+    friend const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode);

+

+    friend const mpreal fac_ui (unsigned long int v,  mp_prec_t prec, mp_rnd_t rnd_mode);

+    friend const mpreal eint   (const mpreal& v, mp_rnd_t rnd_mode);

+

+    friend const mpreal gamma    (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal lngamma  (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal lgamma   (const mpreal& v, int *signp, mp_rnd_t rnd_mode);

+    friend const mpreal zeta     (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal erf      (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal erfc     (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal besselj0 (const mpreal& v, mp_rnd_t rnd_mode); 

+    friend const mpreal besselj1 (const mpreal& v, mp_rnd_t rnd_mode); 

+    friend const mpreal besseljn (long n, const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal bessely0 (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal bessely1 (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal besselyn (long n, const mpreal& v, mp_rnd_t rnd_mode); 

+    friend const mpreal fma      (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode);

+    friend const mpreal fms      (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode);

+    friend const mpreal agm      (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode);

+    friend const mpreal sum      (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode);

     friend int sgn(const mpreal& v); // returns -1 or +1

 

 // MPFR 2.4.0 Specifics

 #if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))

-    friend int          sinh_cosh   (mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal li2         (const mpreal& v,                       mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal fmod        (const mpreal& x, const mpreal& y,      mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal rec_sqrt    (const mpreal& v,                       mp_rnd_t rnd_mode = mpreal::get_default_rnd());

+    friend int          sinh_cosh   (mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal li2         (const mpreal& v,                       mp_rnd_t rnd_mode);

+    friend const mpreal fmod        (const mpreal& x, const mpreal& y,      mp_rnd_t rnd_mode);

+    friend const mpreal rec_sqrt    (const mpreal& v,                       mp_rnd_t rnd_mode);

 

     // MATLAB's semantic equivalents

-    friend const mpreal rem (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); // Remainder after division

-    friend const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); // Modulus after division

+    friend const mpreal rem (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); // Remainder after division

+    friend const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); // Modulus after division

 #endif

 

 // MPFR 3.0.0 Specifics

 #if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))

-    friend const mpreal digamma (const mpreal& v,        mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal ai      (const mpreal& v,        mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal urandom (gmp_randstate_t& state, mp_rnd_t rnd_mode = mpreal::get_default_rnd());     // use gmp_randinit_default() to init state, gmp_randclear() to clear

-    friend const mpreal grandom (gmp_randstate_t& state, mp_rnd_t rnd_mode = mpreal::get_default_rnd());     // use gmp_randinit_default() to init state, gmp_randclear() to clear

-    friend const mpreal grandom (unsigned int seed = 0);

+    friend const mpreal digamma (const mpreal& v,        mp_rnd_t rnd_mode);

+    friend const mpreal ai      (const mpreal& v,        mp_rnd_t rnd_mode);

+    friend const mpreal urandom (gmp_randstate_t& state, mp_rnd_t rnd_mode);     // use gmp_randinit_default() to init state, gmp_randclear() to clear

+    friend const mpreal grandom (gmp_randstate_t& state, mp_rnd_t rnd_mode);     // use gmp_randinit_default() to init state, gmp_randclear() to clear

+    friend const mpreal grandom (unsigned int seed);

 #endif

     

     // Uniformly distributed random number generation in [0,1] using

     // Mersenne-Twister algorithm by default.

     // Use parameter to setup seed, e.g.: random((unsigned)time(NULL))

     // Check urandom() for more precise control.

-    friend const mpreal random(unsigned int seed = 0);

+    friend const mpreal random(unsigned int seed);

 

     // Exponent and mantissa manipulation

     friend const mpreal frexp(const mpreal& v, mp_exp_t* exp);    

@@ -458,31 +489,31 @@ public:
 

     // Constants

     // don't forget to call mpfr_free_cache() for every thread where you are using const-functions

-    friend const mpreal const_log2      (mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal const_pi        (mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal const_euler     (mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal const_catalan   (mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t rnd_mode = mpreal::get_default_rnd());

+    friend const mpreal const_log2      (mp_prec_t prec, mp_rnd_t rnd_mode);

+    friend const mpreal const_pi        (mp_prec_t prec, mp_rnd_t rnd_mode);

+    friend const mpreal const_euler     (mp_prec_t prec, mp_rnd_t rnd_mode);

+    friend const mpreal const_catalan   (mp_prec_t prec, mp_rnd_t rnd_mode);

 

     // returns +inf iff sign>=0 otherwise -inf

-    friend const mpreal const_infinity(int sign = 1, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t rnd_mode = mpreal::get_default_rnd());

+    friend const mpreal const_infinity(int sign, mp_prec_t prec);

 

     // Output/ Input

     friend std::ostream& operator<<(std::ostream& os, const mpreal& v);

     friend std::istream& operator>>(std::istream& is, mpreal& v);

 

     // Integer Related Functions

-    friend const mpreal rint (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

+    friend const mpreal rint (const mpreal& v, mp_rnd_t rnd_mode);

     friend const mpreal ceil (const mpreal& v);

     friend const mpreal floor(const mpreal& v);

     friend const mpreal round(const mpreal& v);

     friend const mpreal trunc(const mpreal& v);

-    friend const mpreal rint_ceil   (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal rint_floor  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal rint_round  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal rint_trunc  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal frac        (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal remainder   (         const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal remquo      (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

+    friend const mpreal rint_ceil   (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal rint_floor  (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal rint_round  (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal rint_trunc  (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal frac        (const mpreal& v, mp_rnd_t rnd_mode);

+    friend const mpreal remainder   (         const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode);

+    friend const mpreal remquo      (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode);

     

     // Miscellaneous Functions

     friend const mpreal nexttoward (const mpreal& x, const mpreal& y);

@@ -549,8 +580,8 @@ public:
     // Efficient swapping of two mpreal values - needed for std algorithms

     friend void swap(mpreal& x, mpreal& y);

     

-    friend const mpreal fmax(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-    friend const mpreal fmin(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

+    friend const mpreal fmax(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode);

+    friend const mpreal fmin(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode);

 

 private:

     // Human friendly Debug Preview in Visual Studio.

@@ -579,16 +610,16 @@ public:
 // Default constructor: creates mp number and initializes it to 0.

 inline mpreal::mpreal() 

 { 

-    mpfr_init2(mp,mpreal::get_default_prec()); 

-    mpfr_set_ui(mp,0,mpreal::get_default_rnd());

+    mpfr_init2 (mpfr_ptr(), mpreal::get_default_prec()); 

+    mpfr_set_ui(mpfr_ptr(), 0, mpreal::get_default_rnd());

 

     MPREAL_MSVC_DEBUGVIEW_CODE;

 }

 

 inline mpreal::mpreal(const mpreal& u) 

 {

-    mpfr_init2(mp,mpfr_get_prec(u.mp));

-    mpfr_set(mp,u.mp,mpreal::get_default_rnd());

+    mpfr_init2(mpfr_ptr(),mpfr_get_prec(u.mpfr_srcptr()));

+    mpfr_set  (mpfr_ptr(),u.mpfr_srcptr(),mpreal::get_default_rnd());

 

     MPREAL_MSVC_DEBUGVIEW_CODE;

 }

@@ -596,8 +627,7 @@ inline mpreal::mpreal(const mpreal& u)
 #ifdef MPREAL_HAVE_MOVE_SUPPORT

 inline mpreal::mpreal(mpreal&& other)

 {

-    mpfr_ptr()->_mpfr_d = 0;                // make sure "other" holds no pinter to actual data 

-

+    mpfr_set_uninitialized(mpfr_ptr());     // make sure "other" holds no pinter to actual data 

     mpfr_swap(mpfr_ptr(), other.mpfr_ptr());

 

     MPREAL_MSVC_DEBUGVIEW_CODE;

@@ -616,12 +646,12 @@ inline mpreal::mpreal(const mpfr_t  u, bool shared)
 {

     if(shared)

     {

-        std::memcpy(mp, u, sizeof(mpfr_t));

+        std::memcpy(mpfr_ptr(), u, sizeof(mpfr_t));

     }

     else

     {

-        mpfr_init2(mp, mpfr_get_prec(u));

-        mpfr_set  (mp, u, mpreal::get_default_rnd());

+        mpfr_init2(mpfr_ptr(), mpfr_get_prec(u));

+        mpfr_set  (mpfr_ptr(), u, mpreal::get_default_rnd());

     }

 

     MPREAL_MSVC_DEBUGVIEW_CODE;

@@ -629,40 +659,40 @@ inline mpreal::mpreal(const mpfr_t  u, bool shared)
 

 inline mpreal::mpreal(const mpf_t u)

 {

-    mpfr_init2(mp,(mp_prec_t) mpf_get_prec(u)); // (gmp: mp_bitcnt_t) unsigned long -> long (mpfr: mp_prec_t)

-    mpfr_set_f(mp,u,mpreal::get_default_rnd());

+    mpfr_init2(mpfr_ptr(),(mp_prec_t) mpf_get_prec(u)); // (gmp: mp_bitcnt_t) unsigned long -> long (mpfr: mp_prec_t)

+    mpfr_set_f(mpfr_ptr(),u,mpreal::get_default_rnd());

 

     MPREAL_MSVC_DEBUGVIEW_CODE;

 }

 

 inline mpreal::mpreal(const mpz_t u, mp_prec_t prec, mp_rnd_t mode)

 {

-    mpfr_init2(mp,prec);

-    mpfr_set_z(mp,u,mode);

+    mpfr_init2(mpfr_ptr(), prec);

+    mpfr_set_z(mpfr_ptr(), u, mode);

 

     MPREAL_MSVC_DEBUGVIEW_CODE;

 }

 

 inline mpreal::mpreal(const mpq_t u, mp_prec_t prec, mp_rnd_t mode)

 {

-    mpfr_init2(mp,prec);

-    mpfr_set_q(mp,u,mode);

+    mpfr_init2(mpfr_ptr(), prec);

+    mpfr_set_q(mpfr_ptr(), u, mode);

 

     MPREAL_MSVC_DEBUGVIEW_CODE;

 }

 

 inline mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode)

 {

-     mpfr_init2(mp, prec);

+     mpfr_init2(mpfr_ptr(), prec);

 

 #if (MPREAL_DOUBLE_BITS_OVERFLOW > -1)

-	if(fits_in_bits(u, MPREAL_DOUBLE_BITS_OVERFLOW))

-	{

-		mpfr_set_d(mp, u, mode);

-	}else

-		throw conversion_overflow();

+  if(fits_in_bits(u, MPREAL_DOUBLE_BITS_OVERFLOW))

+  {

+    mpfr_set_d(mpfr_ptr(), u, mode);

+  }else

+    throw conversion_overflow();

 #else

-	mpfr_set_d(mp, u, mode);

+  mpfr_set_d(mpfr_ptr(), u, mode);

 #endif

 

     MPREAL_MSVC_DEBUGVIEW_CODE;

@@ -670,40 +700,40 @@ inline mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode)
 

 inline mpreal::mpreal(const long double u, mp_prec_t prec, mp_rnd_t mode)

 { 

-    mpfr_init2(mp,prec);

-    mpfr_set_ld(mp,u,mode);

+    mpfr_init2 (mpfr_ptr(), prec);

+    mpfr_set_ld(mpfr_ptr(), u, mode);

 

     MPREAL_MSVC_DEBUGVIEW_CODE;

 }

 

 inline mpreal::mpreal(const unsigned long int u, mp_prec_t prec, mp_rnd_t mode)

 { 

-    mpfr_init2(mp,prec);

-    mpfr_set_ui(mp,u,mode);

+    mpfr_init2 (mpfr_ptr(), prec);

+    mpfr_set_ui(mpfr_ptr(), u, mode);

 

     MPREAL_MSVC_DEBUGVIEW_CODE;

 }

 

 inline mpreal::mpreal(const unsigned int u, mp_prec_t prec, mp_rnd_t mode)

 { 

-    mpfr_init2(mp,prec);

-    mpfr_set_ui(mp,u,mode);

+    mpfr_init2 (mpfr_ptr(), prec);

+    mpfr_set_ui(mpfr_ptr(), u, mode);

 

     MPREAL_MSVC_DEBUGVIEW_CODE;

 }

 

 inline mpreal::mpreal(const long int u, mp_prec_t prec, mp_rnd_t mode)

 { 

-    mpfr_init2(mp,prec);

-    mpfr_set_si(mp,u,mode);

+    mpfr_init2 (mpfr_ptr(), prec);

+    mpfr_set_si(mpfr_ptr(), u, mode);

 

     MPREAL_MSVC_DEBUGVIEW_CODE;

 }

 

 inline mpreal::mpreal(const int u, mp_prec_t prec, mp_rnd_t mode)

 { 

-    mpfr_init2(mp,prec);

-    mpfr_set_si(mp,u,mode);

+    mpfr_init2 (mpfr_ptr(), prec);

+    mpfr_set_si(mpfr_ptr(), u, mode);

 

     MPREAL_MSVC_DEBUGVIEW_CODE;

 }

@@ -711,16 +741,16 @@ inline mpreal::mpreal(const int u, mp_prec_t prec, mp_rnd_t mode)
 #if defined (MPREAL_HAVE_INT64_SUPPORT)

 inline mpreal::mpreal(const uint64_t u, mp_prec_t prec, mp_rnd_t mode)

 {

-    mpfr_init2(mp,prec);

-    mpfr_set_uj(mp, u, mode); 

+    mpfr_init2 (mpfr_ptr(), prec);

+    mpfr_set_uj(mpfr_ptr(), u, mode); 

 

     MPREAL_MSVC_DEBUGVIEW_CODE;

 }

 

 inline mpreal::mpreal(const int64_t u, mp_prec_t prec, mp_rnd_t mode)

 {

-    mpfr_init2(mp,prec);

-    mpfr_set_sj(mp, u, mode); 

+    mpfr_init2 (mpfr_ptr(), prec);

+    mpfr_set_sj(mpfr_ptr(), u, mode); 

 

     MPREAL_MSVC_DEBUGVIEW_CODE;

 }

@@ -728,24 +758,26 @@ inline mpreal::mpreal(const int64_t u, mp_prec_t prec, mp_rnd_t mode)
 

 inline mpreal::mpreal(const char* s, mp_prec_t prec, int base, mp_rnd_t mode)

 {

-    mpfr_init2(mp, prec);

-    mpfr_set_str(mp, s, base, mode); 

+    mpfr_init2  (mpfr_ptr(), prec);

+    mpfr_set_str(mpfr_ptr(), s, base, mode); 

 

     MPREAL_MSVC_DEBUGVIEW_CODE;

 }

 

 inline mpreal::mpreal(const std::string& s, mp_prec_t prec, int base, mp_rnd_t mode)

 {

-    mpfr_init2(mp, prec);

-    mpfr_set_str(mp, s.c_str(), base, mode); 

+    mpfr_init2  (mpfr_ptr(), prec);

+    mpfr_set_str(mpfr_ptr(), s.c_str(), base, mode); 

 

     MPREAL_MSVC_DEBUGVIEW_CODE;

 }

 

 inline void mpreal::clear(::mpfr_ptr x)

 {

-    // Use undocumented field in mpfr_t structure to check if it was initialized 

-    if(0 != x->_mpfr_d) mpfr_clear(x);

+#ifdef MPREAL_HAVE_MOVE_SUPPORT

+    if(mpfr_is_initialized(x)) 

+#endif

+    mpfr_clear(x);

 }

 

 inline mpreal::~mpreal() 

@@ -821,6 +853,9 @@ const mpreal sqrt(const int v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
 const mpreal sqrt(const long double v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

 const mpreal sqrt(const double v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

 

+// abs

+inline const mpreal abs(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd());

+

 //////////////////////////////////////////////////////////////////////////

 // pow

 const mpreal pow(const mpreal& a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

@@ -873,6 +908,11 @@ const mpreal pow(const double a, const unsigned int b, mp_rnd_t rnd_mode = mprea
 const mpreal pow(const double a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

 const mpreal pow(const double a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

 

+inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

+inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

+inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

+inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

+

 //////////////////////////////////////////////////////////////////////////

 // Estimate machine epsilon for the given precision

 // Returns smallest eps such that 1.0 + eps != 1.0

@@ -920,15 +960,15 @@ inline mpreal& mpreal::operator=(const mpreal& v)
 {

     if (this != &v)

     {

-		mp_prec_t tp = mpfr_get_prec(mp);

-		mp_prec_t vp = mpfr_get_prec(v.mp);

+    mp_prec_t tp = mpfr_get_prec(  mpfr_srcptr());

+    mp_prec_t vp = mpfr_get_prec(v.mpfr_srcptr());

 

-		if(tp != vp){

-			clear(mpfr_ptr());

-			mpfr_init2(mpfr_ptr(), vp);

-		}

+    if(tp != vp){

+      clear(mpfr_ptr());

+      mpfr_init2(mpfr_ptr(), vp);

+    }

 

-        mpfr_set(mp, v.mp, mpreal::get_default_rnd());

+        mpfr_set(mpfr_ptr(), v.mpfr_srcptr(), mpreal::get_default_rnd());

 

         MPREAL_MSVC_DEBUGVIEW_CODE;

     }

@@ -937,7 +977,7 @@ inline mpreal& mpreal::operator=(const mpreal& v)
 

 inline mpreal& mpreal::operator=(const mpf_t v)

 {

-    mpfr_set_f(mp, v, mpreal::get_default_rnd());

+    mpfr_set_f(mpfr_ptr(), v, mpreal::get_default_rnd());

     

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

@@ -945,7 +985,7 @@ inline mpreal& mpreal::operator=(const mpf_t v)
 

 inline mpreal& mpreal::operator=(const mpz_t v)

 {

-    mpfr_set_z(mp, v, mpreal::get_default_rnd());

+    mpfr_set_z(mpfr_ptr(), v, mpreal::get_default_rnd());

     

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

@@ -953,7 +993,7 @@ inline mpreal& mpreal::operator=(const mpz_t v)
 

 inline mpreal& mpreal::operator=(const mpq_t v)

 {

-    mpfr_set_q(mp, v, mpreal::get_default_rnd());

+    mpfr_set_q(mpfr_ptr(), v, mpreal::get_default_rnd());

 

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

@@ -961,7 +1001,7 @@ inline mpreal& mpreal::operator=(const mpq_t v)
 

 inline mpreal& mpreal::operator=(const long double v)        

 {    

-    mpfr_set_ld(mp, v, mpreal::get_default_rnd());

+    mpfr_set_ld(mpfr_ptr(), v, mpreal::get_default_rnd());

 

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

@@ -970,22 +1010,22 @@ inline mpreal& mpreal::operator=(const long double v)
 inline mpreal& mpreal::operator=(const double v)                

 {   

 #if (MPREAL_DOUBLE_BITS_OVERFLOW > -1)

-	if(fits_in_bits(v, MPREAL_DOUBLE_BITS_OVERFLOW))

-	{

-		mpfr_set_d(mp,v,mpreal::get_default_rnd());

-	}else

-		throw conversion_overflow();

+  if(fits_in_bits(v, MPREAL_DOUBLE_BITS_OVERFLOW))

+  {

+    mpfr_set_d(mpfr_ptr(),v,mpreal::get_default_rnd());

+  }else

+    throw conversion_overflow();

 #else

-	mpfr_set_d(mp,v,mpreal::get_default_rnd());

+  mpfr_set_d(mpfr_ptr(),v,mpreal::get_default_rnd());

 #endif

 

-	MPREAL_MSVC_DEBUGVIEW_CODE;

+  MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator=(const unsigned long int v)    

 {    

-    mpfr_set_ui(mp, v, mpreal::get_default_rnd());    

+    mpfr_set_ui(mpfr_ptr(), v, mpreal::get_default_rnd());    

 

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

@@ -993,7 +1033,7 @@ inline mpreal& mpreal::operator=(const unsigned long int v)
 

 inline mpreal& mpreal::operator=(const unsigned int v)        

 {    

-    mpfr_set_ui(mp, v, mpreal::get_default_rnd());    

+    mpfr_set_ui(mpfr_ptr(), v, mpreal::get_default_rnd());    

 

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

@@ -1001,7 +1041,7 @@ inline mpreal& mpreal::operator=(const unsigned int v)
 

 inline mpreal& mpreal::operator=(const long int v)            

 {    

-    mpfr_set_si(mp, v, mpreal::get_default_rnd());    

+    mpfr_set_si(mpfr_ptr(), v, mpreal::get_default_rnd());    

 

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

@@ -1009,7 +1049,7 @@ inline mpreal& mpreal::operator=(const long int v)
 

 inline mpreal& mpreal::operator=(const int v)

 {    

-    mpfr_set_si(mp, v, mpreal::get_default_rnd());    

+    mpfr_set_si(mpfr_ptr(), v, mpreal::get_default_rnd());    

 

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

@@ -1026,11 +1066,11 @@ inline mpreal& mpreal::operator=(const char* s)
 

     mpfr_t t;

 

-    mpfr_init2(t, mpfr_get_prec(mp));

+    mpfr_init2(t, mpfr_get_prec(mpfr_srcptr()));

 

     if(0 == mpfr_set_str(t, s, 10, mpreal::get_default_rnd()))

     {

-        mpfr_set(mp, t, mpreal::get_default_rnd()); 

+        mpfr_set(mpfr_ptr(), t, mpreal::get_default_rnd()); 

         MPREAL_MSVC_DEBUGVIEW_CODE;

     }

 

@@ -1049,11 +1089,11 @@ inline mpreal& mpreal::operator=(const std::string& s)
 

     mpfr_t t;

 

-    mpfr_init2(t, mpfr_get_prec(mp));

+    mpfr_init2(t, mpfr_get_prec(mpfr_srcptr()));

 

     if(0 == mpfr_set_str(t, s.c_str(), 10, mpreal::get_default_rnd()))

     {

-        mpfr_set(mp, t, mpreal::get_default_rnd()); 

+        mpfr_set(mpfr_ptr(), t, mpreal::get_default_rnd()); 

         MPREAL_MSVC_DEBUGVIEW_CODE;

     }

 

@@ -1066,7 +1106,7 @@ inline mpreal& mpreal::operator=(const std::string& s)
 // + Addition

 inline mpreal& mpreal::operator+=(const mpreal& v)

 {

-    mpfr_add(mp,mp,v.mp,mpreal::get_default_rnd());

+    mpfr_add(mpfr_ptr(), mpfr_srcptr(), v.mpfr_srcptr(), mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

@@ -1080,14 +1120,14 @@ inline mpreal& mpreal::operator+=(const mpf_t u)
 

 inline mpreal& mpreal::operator+=(const mpz_t u)

 {

-    mpfr_add_z(mp,mp,u,mpreal::get_default_rnd());

+    mpfr_add_z(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator+=(const mpq_t u)

 {

-    mpfr_add_q(mp,mp,u,mpreal::get_default_rnd());

+    mpfr_add_q(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

@@ -1102,7 +1142,7 @@ inline mpreal& mpreal::operator+= (const long double u)
 inline mpreal& mpreal::operator+= (const double u)

 {

 #if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))

-    mpfr_add_d(mp,mp,u,mpreal::get_default_rnd());

+    mpfr_add_d(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());

 #else

     *this += mpreal(u);

 #endif

@@ -1113,28 +1153,28 @@ inline mpreal& mpreal::operator+= (const double u)
 

 inline mpreal& mpreal::operator+=(const unsigned long int u)

 {

-    mpfr_add_ui(mp,mp,u,mpreal::get_default_rnd());

+    mpfr_add_ui(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator+=(const unsigned int u)

 {

-    mpfr_add_ui(mp,mp,u,mpreal::get_default_rnd());

+    mpfr_add_ui(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator+=(const long int u)

 {

-    mpfr_add_si(mp,mp,u,mpreal::get_default_rnd());

+    mpfr_add_si(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator+=(const int u)

 {

-    mpfr_add_si(mp,mp,u,mpreal::get_default_rnd());

+    mpfr_add_si(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

@@ -1154,9 +1194,9 @@ inline const mpreal mpreal::operator+()const    {    return mpreal(*this); }
 

 inline const mpreal operator+(const mpreal& a, const mpreal& b)

 {

-	mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));

-	mpfr_add(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());

-	return c;

+  mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));

+  mpfr_add(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());

+  return c;

 }

 

 inline mpreal& mpreal::operator++() 

@@ -1187,21 +1227,21 @@ inline const mpreal mpreal::operator-- (int)
 // - Subtraction

 inline mpreal& mpreal::operator-=(const mpreal& v)

 {

-    mpfr_sub(mp,mp,v.mp,mpreal::get_default_rnd());

+    mpfr_sub(mpfr_ptr(),mpfr_srcptr(),v.mpfr_srcptr(),mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator-=(const mpz_t v)

 {

-    mpfr_sub_z(mp,mp,v,mpreal::get_default_rnd());

+    mpfr_sub_z(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator-=(const mpq_t v)

 {

-    mpfr_sub_q(mp,mp,v,mpreal::get_default_rnd());

+    mpfr_sub_q(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

@@ -1216,7 +1256,7 @@ inline mpreal& mpreal::operator-=(const long double v)
 inline mpreal& mpreal::operator-=(const double v)

 {

 #if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))

-    mpfr_sub_d(mp,mp,v,mpreal::get_default_rnd());

+    mpfr_sub_d(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

 #else

     *this -= mpreal(v);    

 #endif

@@ -1227,28 +1267,28 @@ inline mpreal& mpreal::operator-=(const double v)
 

 inline mpreal& mpreal::operator-=(const unsigned long int v)

 {

-    mpfr_sub_ui(mp,mp,v,mpreal::get_default_rnd());

+    mpfr_sub_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator-=(const unsigned int v)

 {

-    mpfr_sub_ui(mp,mp,v,mpreal::get_default_rnd());

+    mpfr_sub_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator-=(const long int v)

 {

-    mpfr_sub_si(mp,mp,v,mpreal::get_default_rnd());

+    mpfr_sub_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator-=(const int v)

 {

-    mpfr_sub_si(mp,mp,v,mpreal::get_default_rnd());

+    mpfr_sub_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

@@ -1256,15 +1296,15 @@ inline mpreal& mpreal::operator-=(const int v)
 inline const mpreal mpreal::operator-()const

 {

     mpreal u(*this);

-    mpfr_neg(u.mp,u.mp,mpreal::get_default_rnd());

+    mpfr_neg(u.mpfr_ptr(),u.mpfr_srcptr(),mpreal::get_default_rnd());

     return u;

 }

 

 inline const mpreal operator-(const mpreal& a, const mpreal& b)

 {

-    mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));

-    mpfr_sub(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());

-    return c;

+  mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));

+  mpfr_sub(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());

+  return c;

 }

 

 inline const mpreal operator-(const double  b, const mpreal& a)

@@ -1312,21 +1352,21 @@ inline const mpreal operator-(const int b, const mpreal& a)
 // * Multiplication

 inline mpreal& mpreal::operator*= (const mpreal& v)

 {

-    mpfr_mul(mp,mp,v.mp,mpreal::get_default_rnd());

+    mpfr_mul(mpfr_ptr(),mpfr_srcptr(),v.mpfr_srcptr(),mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator*=(const mpz_t v)

 {

-    mpfr_mul_z(mp,mp,v,mpreal::get_default_rnd());

+    mpfr_mul_z(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator*=(const mpq_t v)

 {

-    mpfr_mul_q(mp,mp,v,mpreal::get_default_rnd());

+    mpfr_mul_q(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

@@ -1341,7 +1381,7 @@ inline mpreal& mpreal::operator*=(const long double v)
 inline mpreal& mpreal::operator*=(const double v)

 {

 #if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))

-    mpfr_mul_d(mp,mp,v,mpreal::get_default_rnd());

+    mpfr_mul_d(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

 #else

     *this *= mpreal(v);    

 #endif

@@ -1351,58 +1391,58 @@ inline mpreal& mpreal::operator*=(const double v)
 

 inline mpreal& mpreal::operator*=(const unsigned long int v)

 {

-    mpfr_mul_ui(mp,mp,v,mpreal::get_default_rnd());

+    mpfr_mul_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator*=(const unsigned int v)

 {

-    mpfr_mul_ui(mp,mp,v,mpreal::get_default_rnd());

+    mpfr_mul_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator*=(const long int v)

 {

-    mpfr_mul_si(mp,mp,v,mpreal::get_default_rnd());

+    mpfr_mul_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator*=(const int v)

 {

-    mpfr_mul_si(mp,mp,v,mpreal::get_default_rnd());

+    mpfr_mul_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline const mpreal operator*(const mpreal& a, const mpreal& b)

 {

-	mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));

-	mpfr_mul(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());

-	return c;

+  mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));

+  mpfr_mul(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());

+  return c;

 }

 

 //////////////////////////////////////////////////////////////////////////

 // / Division

 inline mpreal& mpreal::operator/=(const mpreal& v)

 {

-    mpfr_div(mp,mp,v.mp,mpreal::get_default_rnd());

+    mpfr_div(mpfr_ptr(),mpfr_srcptr(),v.mpfr_srcptr(),mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator/=(const mpz_t v)

 {

-    mpfr_div_z(mp,mp,v,mpreal::get_default_rnd());

+    mpfr_div_z(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator/=(const mpq_t v)

 {

-    mpfr_div_q(mp,mp,v,mpreal::get_default_rnd());

+    mpfr_div_q(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

@@ -1417,7 +1457,7 @@ inline mpreal& mpreal::operator/=(const long double v)
 inline mpreal& mpreal::operator/=(const double v)

 {

 #if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))

-    mpfr_div_d(mp,mp,v,mpreal::get_default_rnd());

+    mpfr_div_d(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

 #else

     *this /= mpreal(v);    

 #endif

@@ -1427,72 +1467,72 @@ inline mpreal& mpreal::operator/=(const double v)
 

 inline mpreal& mpreal::operator/=(const unsigned long int v)

 {

-    mpfr_div_ui(mp,mp,v,mpreal::get_default_rnd());

+    mpfr_div_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator/=(const unsigned int v)

 {

-    mpfr_div_ui(mp,mp,v,mpreal::get_default_rnd());

+    mpfr_div_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator/=(const long int v)

 {

-    mpfr_div_si(mp,mp,v,mpreal::get_default_rnd());

+    mpfr_div_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator/=(const int v)

 {

-    mpfr_div_si(mp,mp,v,mpreal::get_default_rnd());

+    mpfr_div_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline const mpreal operator/(const mpreal& a, const mpreal& b)

 {

-	mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));

-	mpfr_div(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());

-	return c;

+  mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_srcptr()), mpfr_get_prec(b.mpfr_srcptr())));

+  mpfr_div(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());

+  return c;

 }

 

 inline const mpreal operator/(const unsigned long int b, const mpreal& a)

 {

-    mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));

+    mpreal x(0, mpfr_get_prec(a.mpfr_srcptr()));

     mpfr_ui_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());

     return x;

 }

 

 inline const mpreal operator/(const unsigned int b, const mpreal& a)

 {

-    mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));

+    mpreal x(0, mpfr_get_prec(a.mpfr_srcptr()));

     mpfr_ui_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());

     return x;

 }

 

 inline const mpreal operator/(const long int b, const mpreal& a)

 {

-    mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));

-    mpfr_si_div(x.mpfr_ptr(), b, a.mpfr_srcptr(),mpreal::get_default_rnd());

+    mpreal x(0, mpfr_get_prec(a.mpfr_srcptr()));

+    mpfr_si_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());

     return x;

 }

 

 inline const mpreal operator/(const int b, const mpreal& a)

 {

-    mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));

-    mpfr_si_div(x.mpfr_ptr(), b, a.mpfr_srcptr(),mpreal::get_default_rnd());

+    mpreal x(0, mpfr_get_prec(a.mpfr_srcptr()));

+    mpfr_si_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());

     return x;

 }

 

 inline const mpreal operator/(const double  b, const mpreal& a)

 {

 #if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))

-    mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));

-    mpfr_d_div(x.mpfr_ptr(), b, a.mpfr_srcptr(),mpreal::get_default_rnd());

+    mpreal x(0, mpfr_get_prec(a.mpfr_srcptr()));

+    mpfr_d_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());

     return x;

 #else

     mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));

@@ -1505,56 +1545,56 @@ inline const mpreal operator/(const double  b, const mpreal& a)
 // Shifts operators - Multiplication/Division by power of 2

 inline mpreal& mpreal::operator<<=(const unsigned long int u)

 {

-    mpfr_mul_2ui(mp,mp,u,mpreal::get_default_rnd());

+    mpfr_mul_2ui(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator<<=(const unsigned int u)

 {

-    mpfr_mul_2ui(mp,mp,static_cast<unsigned long int>(u),mpreal::get_default_rnd());

+    mpfr_mul_2ui(mpfr_ptr(),mpfr_srcptr(),static_cast<unsigned long int>(u),mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator<<=(const long int u)

 {

-    mpfr_mul_2si(mp,mp,u,mpreal::get_default_rnd());

+    mpfr_mul_2si(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator<<=(const int u)

 {

-    mpfr_mul_2si(mp,mp,static_cast<long int>(u),mpreal::get_default_rnd());

+    mpfr_mul_2si(mpfr_ptr(),mpfr_srcptr(),static_cast<long int>(u),mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator>>=(const unsigned long int u)

 {

-    mpfr_div_2ui(mp,mp,u,mpreal::get_default_rnd());

+    mpfr_div_2ui(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator>>=(const unsigned int u)

 {

-    mpfr_div_2ui(mp,mp,static_cast<unsigned long int>(u),mpreal::get_default_rnd());

+    mpfr_div_2ui(mpfr_ptr(),mpfr_srcptr(),static_cast<unsigned long int>(u),mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator>>=(const long int u)

 {

-    mpfr_div_2si(mp,mp,u,mpreal::get_default_rnd());

+    mpfr_div_2si(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::operator>>=(const int u)

 {

-    mpfr_div_2si(mp,mp,static_cast<long int>(u),mpreal::get_default_rnd());

+    mpfr_div_2si(mpfr_ptr(),mpfr_srcptr(),static_cast<long int>(u),mpreal::get_default_rnd());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

@@ -1603,7 +1643,7 @@ inline const mpreal operator>>(const mpreal& v, const int k)
 inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode)

 {

     mpreal x(v);

-    mpfr_mul_2ui(x.mp,v.mp,k,rnd_mode);

+    mpfr_mul_2ui(x.mpfr_ptr(),v.mpfr_srcptr(),k,rnd_mode);

     return x;

 }

 

@@ -1611,61 +1651,63 @@ inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_m
 inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode)

 {

     mpreal x(v);

-    mpfr_mul_2si(x.mp,v.mp,k,rnd_mode);

+    mpfr_mul_2si(x.mpfr_ptr(),v.mpfr_srcptr(),k,rnd_mode);

     return x;

 }

 

 inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode)

 {

     mpreal x(v);

-    mpfr_div_2ui(x.mp,v.mp,k,rnd_mode);

+    mpfr_div_2ui(x.mpfr_ptr(),v.mpfr_srcptr(),k,rnd_mode);

     return x;

 }

 

 inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode)

 {

     mpreal x(v);

-    mpfr_div_2si(x.mp,v.mp,k,rnd_mode);

+    mpfr_div_2si(x.mpfr_ptr(),v.mpfr_srcptr(),k,rnd_mode);

     return x;

 }

 

 //////////////////////////////////////////////////////////////////////////

 //Boolean operators

-inline bool operator >  (const mpreal& a, const mpreal& b){    return (mpfr_greater_p(a.mp,b.mp)      !=0);    }

-inline bool operator >= (const mpreal& a, const mpreal& b){    return (mpfr_greaterequal_p(a.mp,b.mp) !=0);    }

-inline bool operator <  (const mpreal& a, const mpreal& b){    return (mpfr_less_p(a.mp,b.mp)         !=0);    }

-inline bool operator <= (const mpreal& a, const mpreal& b){    return (mpfr_lessequal_p(a.mp,b.mp)    !=0);    }

-inline bool operator == (const mpreal& a, const mpreal& b){    return (mpfr_equal_p(a.mp,b.mp)        !=0);    }

-inline bool operator != (const mpreal& a, const mpreal& b){    return (mpfr_lessgreater_p(a.mp,b.mp)  !=0);    }

-

-inline bool operator == (const mpreal& a, const unsigned long int b ){    return (mpfr_cmp_ui(a.mp,b) == 0);    }

-inline bool operator == (const mpreal& a, const unsigned int b      ){    return (mpfr_cmp_ui(a.mp,b) == 0);    }

-inline bool operator == (const mpreal& a, const long int b          ){    return (mpfr_cmp_si(a.mp,b) == 0);    }

-inline bool operator == (const mpreal& a, const int b               ){    return (mpfr_cmp_si(a.mp,b) == 0);    }

-inline bool operator == (const mpreal& a, const long double b       ){    return (mpfr_cmp_ld(a.mp,b) == 0);    }

-inline bool operator == (const mpreal& a, const double b            ){    return (mpfr_cmp_d(a.mp,b)  == 0);    }

-

-

-inline bool isnan    (const mpreal& v){    return (mpfr_nan_p(v.mp)     != 0);    }

-inline bool isinf    (const mpreal& v){    return (mpfr_inf_p(v.mp)     != 0);    }

-inline bool isfinite (const mpreal& v){    return (mpfr_number_p(v.mp)  != 0);    }

-inline bool iszero   (const mpreal& v){    return (mpfr_zero_p(v.mp)    != 0);    }

-inline bool isint    (const mpreal& v){    return (mpfr_integer_p(v.mp) != 0);    }

+inline bool operator >  (const mpreal& a, const mpreal& b){    return (mpfr_greater_p       (a.mpfr_srcptr(),b.mpfr_srcptr()) !=0 );    }

+inline bool operator >= (const mpreal& a, const mpreal& b){    return (mpfr_greaterequal_p  (a.mpfr_srcptr(),b.mpfr_srcptr()) !=0 );    }

+inline bool operator <  (const mpreal& a, const mpreal& b){    return (mpfr_less_p          (a.mpfr_srcptr(),b.mpfr_srcptr()) !=0 );    }

+inline bool operator <= (const mpreal& a, const mpreal& b){    return (mpfr_lessequal_p     (a.mpfr_srcptr(),b.mpfr_srcptr()) !=0 );    }

+inline bool operator == (const mpreal& a, const mpreal& b){    return (mpfr_equal_p         (a.mpfr_srcptr(),b.mpfr_srcptr()) !=0 );    }

+inline bool operator != (const mpreal& a, const mpreal& b){    return (mpfr_lessgreater_p   (a.mpfr_srcptr(),b.mpfr_srcptr()) !=0 );    }

+

+inline bool operator == (const mpreal& a, const unsigned long int b ){    return (mpfr_cmp_ui(a.mpfr_srcptr(),b) == 0 );    }

+inline bool operator == (const mpreal& a, const unsigned int b      ){    return (mpfr_cmp_ui(a.mpfr_srcptr(),b) == 0 );    }

+inline bool operator == (const mpreal& a, const long int b          ){    return (mpfr_cmp_si(a.mpfr_srcptr(),b) == 0 );    }

+inline bool operator == (const mpreal& a, const int b               ){    return (mpfr_cmp_si(a.mpfr_srcptr(),b) == 0 );    }

+inline bool operator == (const mpreal& a, const long double b       ){    return (mpfr_cmp_ld(a.mpfr_srcptr(),b) == 0 );    }

+inline bool operator == (const mpreal& a, const double b            ){    return (mpfr_cmp_d (a.mpfr_srcptr(),b) == 0 );    }

+

+

+inline bool isnan    (const mpreal& op){    return (mpfr_nan_p    (op.mpfr_srcptr()) != 0 );    }

+inline bool isinf    (const mpreal& op){    return (mpfr_inf_p    (op.mpfr_srcptr()) != 0 );    }

+inline bool isfinite (const mpreal& op){    return (mpfr_number_p (op.mpfr_srcptr()) != 0 );    }

+inline bool iszero   (const mpreal& op){    return (mpfr_zero_p   (op.mpfr_srcptr()) != 0 );    }

+inline bool isint    (const mpreal& op){    return (mpfr_integer_p(op.mpfr_srcptr()) != 0 );    }

 

 #if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))

-inline bool isregular(const mpreal& v){    return (mpfr_regular_p(v.mp));}

+inline bool isregular(const mpreal& op){    return (mpfr_regular_p(op.mpfr_srcptr()));}

 #endif 

 

 //////////////////////////////////////////////////////////////////////////

 // Type Converters

-inline long             mpreal::toLong   (mp_rnd_t mode)  const    {    return mpfr_get_si(mp, mode);    }

-inline unsigned long    mpreal::toULong  (mp_rnd_t mode)  const    {    return mpfr_get_ui(mp, mode);    }

-inline double           mpreal::toDouble (mp_rnd_t mode)  const    {    return mpfr_get_d (mp, mode);    }

-inline long double      mpreal::toLDouble(mp_rnd_t mode)  const    {    return mpfr_get_ld(mp, mode);    }

+inline bool             mpreal::toBool (mp_rnd_t /*mode*/) const   {    return  mpfr_zero_p (mpfr_srcptr()) == 0;     }

+inline long             mpreal::toLong   (mp_rnd_t mode)  const    {    return  mpfr_get_si (mpfr_srcptr(), mode);    }

+inline unsigned long    mpreal::toULong  (mp_rnd_t mode)  const    {    return  mpfr_get_ui (mpfr_srcptr(), mode);    }

+inline float            mpreal::toFloat  (mp_rnd_t mode)  const    {    return  mpfr_get_flt(mpfr_srcptr(), mode);    }

+inline double           mpreal::toDouble (mp_rnd_t mode)  const    {    return  mpfr_get_d  (mpfr_srcptr(), mode);    }

+inline long double      mpreal::toLDouble(mp_rnd_t mode)  const    {    return  mpfr_get_ld (mpfr_srcptr(), mode);    }

 

 #if defined (MPREAL_HAVE_INT64_SUPPORT)

-inline int64_t      mpreal::toInt64 (mp_rnd_t mode)    const{    return mpfr_get_sj(mp, mode);    }

-inline uint64_t     mpreal::toUInt64(mp_rnd_t mode)    const{    return mpfr_get_uj(mp, mode);    }

+inline int64_t      mpreal::toInt64 (mp_rnd_t mode)    const{    return mpfr_get_sj(mpfr_srcptr(), mode);    }

+inline uint64_t     mpreal::toUInt64(mp_rnd_t mode)    const{    return mpfr_get_uj(mpfr_srcptr(), mode);    }

 #endif

 

 inline ::mpfr_ptr     mpreal::mpfr_ptr()             { return mp; }

@@ -1689,7 +1731,7 @@ inline std::string mpreal::toString(const std::string& format) const
 

     if( !format.empty() )

     {

-        if(!(mpfr_asprintf(&s,format.c_str(),mp) < 0))

+        if(!(mpfr_asprintf(&s, format.c_str(), mpfr_srcptr()) < 0))

         {

             out = std::string(s);

 

@@ -1704,20 +1746,19 @@ inline std::string mpreal::toString(const std::string& format) const
 

 inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const

 {

+    // TODO: Add extended format specification (f, e, rounding mode) as it done in output operator

     (void)b;

     (void)mode;

 

 #if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))

 

-    // Use MPFR native function for output

-    char format[128];

-    int digits;

+    std::ostringstream format;

 

-    digits = n > 0 ? n : bits2digits(mpfr_get_prec(mp));

-

-    sprintf(format,"%%.%dRNg",digits);        // Default format

+    int digits = (n >= 0) ? n : bits2digits(mpfr_get_prec(mpfr_srcptr()));

+    

+    format << "%." << digits << "RNg";

 

-    return toString(std::string(format));

+    return toString(format.str());

 

 #else

 

@@ -1735,8 +1776,8 @@ inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const
     if(mpfr_zero_p(mp)) return "0";

     if(mpfr_nan_p(mp))  return "NaN";

 

-    s  = mpfr_get_str(NULL,&exp,b,0,mp,mode);

-    ns = mpfr_get_str(NULL,&exp,b,n,mp,mode);

+    s  = mpfr_get_str(NULL, &exp, b, 0, mp, mode);

+    ns = mpfr_get_str(NULL, &exp, b, (std::max)(0,n), mp, mode);

 

     if(s!=NULL && ns!=NULL)

     {

@@ -1821,17 +1862,43 @@ inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const
 

 //////////////////////////////////////////////////////////////////////////

 // I/O

+inline std::ostream& mpreal::output(std::ostream& os) const 

+{

+    std::ostringstream format;

+    const std::ios::fmtflags flags = os.flags();

+

+    format << ((flags & std::ios::showpos) ? "%+" : "%");

+    if (os.precision() >= 0)

+        format << '.' << os.precision() << "R*"

+               << ((flags & std::ios::floatfield) == std::ios::fixed ? 'f' :

+                   (flags & std::ios::floatfield) == std::ios::scientific ? 'e' :

+                   'g');

+    else

+        format << "R*e";

+

+    char *s = NULL;

+    if(!(mpfr_asprintf(&s, format.str().c_str(),

+                        mpfr::mpreal::get_default_rnd(),

+                        mpfr_srcptr())

+        < 0))

+    {

+        os << std::string(s);

+        mpfr_free_str(s);

+    }

+    return os;

+}

+

 inline std::ostream& operator<<(std::ostream& os, const mpreal& v)

 {

-    return os<<v.toString(static_cast<int>(os.precision()));

+    return v.output(os);

 }

 

 inline std::istream& operator>>(std::istream &is, mpreal& v)

 {

-    // ToDo, use cout::hexfloat and other flags to setup base

+    // TODO: use cout::hexfloat and other flags to setup base

     std::string tmp;

     is >> tmp;

-    mpfr_set_str(v.mp, tmp.c_str(), 10, mpreal::get_default_rnd());

+    mpfr_set_str(v.mpfr_ptr(), tmp.c_str(), 10, mpreal::get_default_rnd());

     return is;

 }

 

@@ -1844,53 +1911,53 @@ inline mp_prec_t digits2bits(int d)
 {

     const double LOG2_10 = 3.3219280948873624;

 

-    return (mp_prec_t) std::ceil( d * LOG2_10 );

+    return mp_prec_t(std::ceil( d * LOG2_10 ));

 }

 

 inline int bits2digits(mp_prec_t b)

 {

     const double LOG10_2 = 0.30102999566398119;

 

-    return (int) std::floor( b * LOG10_2 );

+    return int(std::floor( b * LOG10_2 ));

 }

 

 //////////////////////////////////////////////////////////////////////////

 // Set/Get number properties

-inline int sgn(const mpreal& v)

+inline int sgn(const mpreal& op)

 {

-    int r = mpfr_signbit(v.mp);

-    return (r>0?-1:1);

+    int r = mpfr_signbit(op.mpfr_srcptr());

+    return (r > 0? -1 : 1);

 }

 

 inline mpreal& mpreal::setSign(int sign, mp_rnd_t RoundingMode)

 {

-    mpfr_setsign(mp,mp,(sign < 0 ? 1 : 0),RoundingMode);

+    mpfr_setsign(mpfr_ptr(), mpfr_srcptr(), (sign < 0 ? 1 : 0), RoundingMode);

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline int mpreal::getPrecision() const

 {

-    return mpfr_get_prec(mp);

+    return int(mpfr_get_prec(mpfr_srcptr()));

 }

 

 inline mpreal& mpreal::setPrecision(int Precision, mp_rnd_t RoundingMode)

 {

-    mpfr_prec_round(mp, Precision, RoundingMode);

+    mpfr_prec_round(mpfr_ptr(), Precision, RoundingMode);

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

 

 inline mpreal& mpreal::setInf(int sign) 

 { 

-    mpfr_set_inf(mp,sign);

+    mpfr_set_inf(mpfr_ptr(), sign);

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }    

 

 inline mpreal& mpreal::setNan() 

 {

-    mpfr_set_nan(mp);

+    mpfr_set_nan(mpfr_ptr());

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return *this;

 }

@@ -1899,9 +1966,9 @@ inline mpreal&    mpreal::setZero(int sign)
 {

 

 #if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))

-    mpfr_set_zero(mp, sign);

+    mpfr_set_zero(mpfr_ptr(), sign);

 #else

-    mpfr_set_si(mp, 0, (mpfr_get_default_rounding_mode)());

+    mpfr_set_si(mpfr_ptr(), 0, (mpfr_get_default_rounding_mode)());

     setSign(sign);

 #endif 

 

@@ -1911,23 +1978,23 @@ inline mpreal&    mpreal::setZero(int sign)
 

 inline mp_prec_t mpreal::get_prec() const

 {

-    return mpfr_get_prec(mp);

+    return mpfr_get_prec(mpfr_srcptr());

 }

 

 inline void mpreal::set_prec(mp_prec_t prec, mp_rnd_t rnd_mode)

 {

-    mpfr_prec_round(mp,prec,rnd_mode);

+    mpfr_prec_round(mpfr_ptr(),prec,rnd_mode);

     MPREAL_MSVC_DEBUGVIEW_CODE;

 }

 

 inline mp_exp_t mpreal::get_exp ()

 {

-    return mpfr_get_exp(mp);

+    return mpfr_get_exp(mpfr_srcptr());

 }

 

 inline int mpreal::set_exp (mp_exp_t e)

 {

-    int x = mpfr_set_exp(mp, e);

+    int x = mpfr_set_exp(mpfr_ptr(), e);

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return x;

 }

@@ -1945,7 +2012,7 @@ inline const mpreal ldexp(const mpreal& v, mp_exp_t exp)
     mpreal x(v);

 

     // rounding is not important since we just increasing the exponent

-    mpfr_mul_2si(x.mp,x.mp,exp,mpreal::get_default_rnd()); 

+    mpfr_mul_2si(x.mpfr_ptr(), x.mpfr_srcptr(), exp, mpreal::get_default_rnd()); 

     return x;

 }

 

@@ -1960,9 +2027,9 @@ inline mpreal machine_epsilon(const mpreal& x)
     /* the smallest eps such that x + eps != x */

     if( x < 0)

     {

-        return nextabove(-x)+x;

+        return nextabove(-x) + x;

     }else{

-        return nextabove(x)-x;

+        return nextabove( x) - x;

     }

 }

 

@@ -1997,22 +2064,22 @@ inline bool isEqualFuzzy(const mpreal& a, const mpreal& b)
 

 inline const mpreal modf(const mpreal& v, mpreal& n)

 {

-    mpreal frac(v);

+    mpreal f(v);

 

     // rounding is not important since we are using the same number

-    mpfr_frac(frac.mp,frac.mp,mpreal::get_default_rnd());    

-    mpfr_trunc(n.mp,v.mp);

-    return frac;

+    mpfr_frac (f.mpfr_ptr(),f.mpfr_srcptr(),mpreal::get_default_rnd());    

+    mpfr_trunc(n.mpfr_ptr(),v.mpfr_srcptr());

+    return f;

 }

 

 inline int mpreal::check_range (int t, mp_rnd_t rnd_mode)

 {

-    return mpfr_check_range(mp,t,rnd_mode);

+    return mpfr_check_range(mpfr_ptr(),t,rnd_mode);

 }

 

 inline int mpreal::subnormalize (int t,mp_rnd_t rnd_mode)

 {

-    int r = mpfr_subnormalize(mp,t,rnd_mode);

+    int r = mpfr_subnormalize(mpfr_ptr(),t,rnd_mode);

     MPREAL_MSVC_DEBUGVIEW_CODE;

     return r;

 }

@@ -2065,8 +2132,11 @@ inline mp_exp_t mpreal::get_emax_max (void)
         mpfr_##f(y.mpfr_ptr(), x.mpfr_srcptr(), r);           \

         return y; 

 

-inline const mpreal sqr  (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(sqr );    }

-inline const mpreal sqrt (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(sqrt);    }

+inline const mpreal sqr  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd())

+{   MPREAL_UNARY_MATH_FUNCTION_BODY(sqr );    }

+

+inline const mpreal sqrt (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd())

+{   MPREAL_UNARY_MATH_FUNCTION_BODY(sqrt);    }

 

 inline const mpreal sqrt(const unsigned long int x, mp_rnd_t r)

 {

@@ -2092,14 +2162,14 @@ inline const mpreal sqrt(const int v, mp_rnd_t rnd_mode)
     else        return mpreal().setNan(); // NaN

 }

 

-inline const mpreal root(const mpreal& x, unsigned long int k, mp_rnd_t r)

+inline const mpreal root(const mpreal& x, unsigned long int k, mp_rnd_t r = mpreal::get_default_rnd())

 {

     mpreal y(0, mpfr_get_prec(x.mpfr_srcptr())); 

     mpfr_root(y.mpfr_ptr(), x.mpfr_srcptr(), k, r);  

     return y; 

 }

 

-inline const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t r)

+inline const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t r = mpreal::get_default_rnd())

 {

     mpreal y(0, mpfr_get_prec(a.mpfr_srcptr()));

     mpfr_dim(y.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), r);

@@ -2111,7 +2181,7 @@ inline int cmpabs(const mpreal& a,const mpreal& b)
     return mpfr_cmpabs(a.mpfr_ptr(), b.mpfr_srcptr());

 }

 

-inline int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode)

+inline int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

 {

     return mpfr_sin_cos(s.mpfr_ptr(), c.mpfr_ptr(), v.mpfr_srcptr(), rnd_mode);

 }

@@ -2119,84 +2189,85 @@ inline int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode)
 inline const mpreal sqrt  (const long double v, mp_rnd_t rnd_mode)    {   return sqrt(mpreal(v),rnd_mode);    }

 inline const mpreal sqrt  (const double v, mp_rnd_t rnd_mode)         {   return sqrt(mpreal(v),rnd_mode);    }

 

-inline const mpreal cbrt  (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(cbrt );    }

-inline const mpreal fabs  (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(abs  );    }

-inline const mpreal abs   (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(abs  );    }

-inline const mpreal log   (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(log  );    }

-inline const mpreal log2  (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(log2 );    }

-inline const mpreal log10 (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(log10);    }

-inline const mpreal exp   (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(exp  );    }

-inline const mpreal exp2  (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(exp2 );    }

-inline const mpreal exp10 (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(exp10);    }

-inline const mpreal cos   (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(cos  );    }

-inline const mpreal sin   (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(sin  );    }

-inline const mpreal tan   (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(tan  );    }

-inline const mpreal sec   (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(sec  );    }

-inline const mpreal csc   (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(csc  );    }

-inline const mpreal cot   (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(cot  );    }

-inline const mpreal acos  (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(acos );     }

-inline const mpreal asin  (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(asin );     }

-inline const mpreal atan  (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(atan );     }

-

-inline const mpreal acot  (const mpreal& v, mp_rnd_t r) {   return atan (1/v, r);                      }

-inline const mpreal asec  (const mpreal& v, mp_rnd_t r) {   return acos (1/v, r);                      }

-inline const mpreal acsc  (const mpreal& v, mp_rnd_t r) {   return asin (1/v, r);                      }

-inline const mpreal acoth (const mpreal& v, mp_rnd_t r) {   return atanh(1/v, r);                      }

-inline const mpreal asech (const mpreal& v, mp_rnd_t r) {   return acosh(1/v, r);                      }

-inline const mpreal acsch (const mpreal& v, mp_rnd_t r) {   return asinh(1/v, r);                      }

-

-inline const mpreal cosh  (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(cosh );    }

-inline const mpreal sinh  (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(sinh );    }

-inline const mpreal tanh  (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(tanh );    }

-inline const mpreal sech  (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(sech );    }

-inline const mpreal csch  (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(csch );    }

-inline const mpreal coth  (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(coth );    }

-inline const mpreal acosh (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(acosh);    }

-inline const mpreal asinh (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(asinh);    }

-inline const mpreal atanh (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(atanh);    }

-

-inline const mpreal log1p   (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(log1p  );    }

-inline const mpreal expm1   (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(expm1  );    }

-inline const mpreal eint    (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(eint   );    }

-inline const mpreal gamma   (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(gamma  );    }

-inline const mpreal lngamma (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(lngamma);    }

-inline const mpreal zeta    (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(zeta   );    }

-inline const mpreal erf     (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(erf    );    }

-inline const mpreal erfc    (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(erfc   );    }

-inline const mpreal besselj0(const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(j0     );    }

-inline const mpreal besselj1(const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(j1     );    }

-inline const mpreal bessely0(const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(y0     );    }

-inline const mpreal bessely1(const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(y1     );    }

-

-inline const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode)

+inline const mpreal cbrt  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(cbrt );    }

+inline const mpreal fabs  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(abs  );    }

+inline const mpreal abs   (const mpreal& x, mp_rnd_t r)                             {   MPREAL_UNARY_MATH_FUNCTION_BODY(abs  );    }

+inline const mpreal log   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(log  );    }

+inline const mpreal log2  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(log2 );    }

+inline const mpreal log10 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(log10);    }

+inline const mpreal exp   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(exp  );    }

+inline const mpreal exp2  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(exp2 );    }

+inline const mpreal exp10 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(exp10);    }

+inline const mpreal cos   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(cos  );    }

+inline const mpreal sin   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(sin  );    }

+inline const mpreal tan   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(tan  );    }

+inline const mpreal sec   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(sec  );    }

+inline const mpreal csc   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(csc  );    }

+inline const mpreal cot   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(cot  );    }

+inline const mpreal acos  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(acos );    }

+inline const mpreal asin  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(asin );    }

+inline const mpreal atan  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(atan );    }

+

+inline const mpreal acot  (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) {   return atan (1/v, r);                      }

+inline const mpreal asec  (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) {   return acos (1/v, r);                      }

+inline const mpreal acsc  (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) {   return asin (1/v, r);                      }

+inline const mpreal acoth (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) {   return atanh(1/v, r);                      }

+inline const mpreal asech (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) {   return acosh(1/v, r);                      }

+inline const mpreal acsch (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) {   return asinh(1/v, r);                      }

+

+inline const mpreal cosh  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(cosh );    }

+inline const mpreal sinh  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(sinh );    }

+inline const mpreal tanh  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(tanh );    }

+inline const mpreal sech  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(sech );    }

+inline const mpreal csch  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(csch );    }

+inline const mpreal coth  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(coth );    }

+inline const mpreal acosh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(acosh);    }

+inline const mpreal asinh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(asinh);    }

+inline const mpreal atanh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(atanh);    }

+

+inline const mpreal log1p   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(log1p  );    }

+inline const mpreal expm1   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(expm1  );    }

+inline const mpreal eint    (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(eint   );    }

+inline const mpreal gamma   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(gamma  );    }

+inline const mpreal lngamma (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(lngamma);    }

+inline const mpreal zeta    (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(zeta   );    }

+inline const mpreal erf     (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(erf    );    }

+inline const mpreal erfc    (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(erfc   );    }

+inline const mpreal besselj0(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(j0     );    }

+inline const mpreal besselj1(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(j1     );    }

+inline const mpreal bessely0(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(y0     );    }

+inline const mpreal bessely1(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(y1     );    }

+

+inline const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

 {

     mpreal a(0,(std::max)(y.getPrecision(), x.getPrecision()));

     mpfr_atan2(a.mpfr_ptr(), y.mpfr_srcptr(), x.mpfr_srcptr(), rnd_mode);

     return a;

 }

 

-inline const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)

+inline const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

 {

     mpreal a(0,(std::max)(y.getPrecision(), x.getPrecision()));

     mpfr_hypot(a.mpfr_ptr(), x.mpfr_srcptr(), y.mpfr_srcptr(), rnd_mode);

     return a;

 }

 

-inline const mpreal remainder (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)

+inline const mpreal remainder (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

 {    

     mpreal a(0,(std::max)(y.getPrecision(), x.getPrecision()));

     mpfr_remainder(a.mpfr_ptr(), x.mpfr_srcptr(), y.mpfr_srcptr(), rnd_mode);

     return a;

 }

 

-inline const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)

+inline const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

 {

     mpreal a(0,(std::max)(y.getPrecision(), x.getPrecision()));

     mpfr_remquo(a.mpfr_ptr(),q, x.mpfr_srcptr(), y.mpfr_srcptr(), rnd_mode);

     return a;

 }

 

-inline const mpreal fac_ui (unsigned long int v, mp_prec_t prec, mp_rnd_t rnd_mode)

+inline const mpreal fac_ui (unsigned long int v, mp_prec_t prec     = mpreal::get_default_prec(),

+                                           mp_rnd_t  rnd_mode = mpreal::get_default_rnd())

 {

     mpreal x(0, prec);

     mpfr_fac_ui(x.mpfr_ptr(),v,rnd_mode);

@@ -2204,33 +2275,33 @@ inline const mpreal fac_ui (unsigned long int v, mp_prec_t prec, mp_rnd_t rnd_mo
 }

 

 

-inline const mpreal lgamma (const mpreal& v, int *signp, mp_rnd_t rnd_mode)

+inline const mpreal lgamma (const mpreal& v, int *signp = 0, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

 {

     mpreal x(v);

     int tsignp;

 

-    if(signp)   mpfr_lgamma(x.mp,  signp,v.mp,rnd_mode);

-    else        mpfr_lgamma(x.mp,&tsignp,v.mp,rnd_mode);

+    if(signp)   mpfr_lgamma(x.mpfr_ptr(),  signp,v.mpfr_srcptr(),rnd_mode);

+    else        mpfr_lgamma(x.mpfr_ptr(),&tsignp,v.mpfr_srcptr(),rnd_mode);

 

     return x;

 }

 

 

-inline const mpreal besseljn (long n, const mpreal& x, mp_rnd_t r)

+inline const mpreal besseljn (long n, const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd())

 {

-    mpreal  y(0, mpfr_get_prec(x.mpfr_srcptr()));

+    mpreal  y(0, x.getPrecision());

     mpfr_jn(y.mpfr_ptr(), n, x.mpfr_srcptr(), r);

     return y;

 }

 

-inline const mpreal besselyn (long n, const mpreal& x, mp_rnd_t r)

+inline const mpreal besselyn (long n, const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd())

 {

-    mpreal  y(0, mpfr_get_prec(x.mpfr_srcptr()));

+    mpreal  y(0, x.getPrecision());

     mpfr_yn(y.mpfr_ptr(), n, x.mpfr_srcptr(), r);

     return y;

 }

 

-inline const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode)

+inline const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

 {

     mpreal a;

     mp_prec_t p1, p2, p3;

@@ -2245,7 +2316,7 @@ inline const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, m
     return a;

 }

 

-inline const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode)

+inline const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

 {

     mpreal a;

     mp_prec_t p1, p2, p3;

@@ -2260,7 +2331,7 @@ inline const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, m
     return a;

 }

 

-inline const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode)

+inline const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

 {

     mpreal a;

     mp_prec_t p1, p2;

@@ -2275,7 +2346,7 @@ inline const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode)
     return a;

 }

 

-inline const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode)

+inline const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

 {

     mpreal x;

     mpfr_ptr* t;

@@ -2292,23 +2363,23 @@ inline const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_m
 // MPFR 2.4.0 Specifics

 #if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))

 

-inline int sinh_cosh(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode)

+inline int sinh_cosh(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

 {

     return mpfr_sinh_cosh(s.mp,c.mp,v.mp,rnd_mode);

 }

 

-inline const mpreal li2 (const mpreal& x, mp_rnd_t r) 

+inline const mpreal li2 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) 

 {   

     MPREAL_UNARY_MATH_FUNCTION_BODY(li2);    

 }

 

-inline const mpreal rem (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)

+inline const mpreal rem (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

 {

     /*  R = rem(X,Y) if Y != 0, returns X - n * Y where n = trunc(X/Y). */

     return fmod(x, y, rnd_mode);

 }

 

-inline const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)

+inline const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

 {

     (void)rnd_mode;

     

@@ -2333,7 +2404,7 @@ inline const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
     return m;

 }

 

-inline const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)

+inline const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

 {

     mpreal a;

     mp_prec_t yp, xp;

@@ -2348,7 +2419,7 @@ inline const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
     return a;

 }

 

-inline const mpreal rec_sqrt(const mpreal& v, mp_rnd_t rnd_mode)

+inline const mpreal rec_sqrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

 {

     mpreal x(v);

     mpfr_rec_sqrt(x.mp,v.mp,rnd_mode);

@@ -2359,41 +2430,41 @@ inline const mpreal rec_sqrt(const mpreal& v, mp_rnd_t rnd_mode)
 //////////////////////////////////////////////////////////////////////////

 // MPFR 3.0.0 Specifics

 #if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))

-inline const mpreal digamma (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(digamma);     }

-inline const mpreal ai      (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(ai);          }

+inline const mpreal digamma (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(digamma);     }

+inline const mpreal ai      (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(ai);          }

 #endif // MPFR 3.0.0 Specifics

 

 //////////////////////////////////////////////////////////////////////////

 // Constants

-inline const mpreal const_log2 (mp_prec_t p, mp_rnd_t r)

+inline const mpreal const_log2 (mp_prec_t p = mpreal::get_default_prec(), mp_rnd_t r = mpreal::get_default_rnd())

 {

     mpreal x(0, p);

     mpfr_const_log2(x.mpfr_ptr(), r);

     return x;

 }

 

-inline const mpreal const_pi (mp_prec_t p, mp_rnd_t r)

+inline const mpreal const_pi (mp_prec_t p = mpreal::get_default_prec(), mp_rnd_t r = mpreal::get_default_rnd())

 {

     mpreal x(0, p);

     mpfr_const_pi(x.mpfr_ptr(), r);

     return x;

 }

 

-inline const mpreal const_euler (mp_prec_t p, mp_rnd_t r)

+inline const mpreal const_euler (mp_prec_t p = mpreal::get_default_prec(), mp_rnd_t r = mpreal::get_default_rnd())

 {

     mpreal x(0, p);

     mpfr_const_euler(x.mpfr_ptr(), r);

     return x;

 }

 

-inline const mpreal const_catalan (mp_prec_t p, mp_rnd_t r)

+inline const mpreal const_catalan (mp_prec_t p = mpreal::get_default_prec(), mp_rnd_t r = mpreal::get_default_rnd())

 {

     mpreal x(0, p);

     mpfr_const_catalan(x.mpfr_ptr(), r);

     return x;

 }

 

-inline const mpreal const_infinity (int sign, mp_prec_t p, mp_rnd_t /*r*/)

+inline const mpreal const_infinity (int sign = 1, mp_prec_t p = mpreal::get_default_prec())

 {

     mpreal x(0, p);

     mpfr_set_inf(x.mpfr_ptr(), sign);

@@ -2430,12 +2501,12 @@ inline const mpreal trunc(const mpreal& v)
     return x;

 }

 

-inline const mpreal rint       (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(rint      );     }

-inline const mpreal rint_ceil  (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(rint_ceil );     }

-inline const mpreal rint_floor (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(rint_floor);     }

-inline const mpreal rint_round (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(rint_round);     }

-inline const mpreal rint_trunc (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(rint_trunc);     }

-inline const mpreal frac       (const mpreal& x, mp_rnd_t r) {   MPREAL_UNARY_MATH_FUNCTION_BODY(frac      );     }

+inline const mpreal rint       (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(rint      );     }

+inline const mpreal rint_ceil  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(rint_ceil );     }

+inline const mpreal rint_floor (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(rint_floor);     }

+inline const mpreal rint_round (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(rint_round);     }

+inline const mpreal rint_trunc (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(rint_trunc);     }

+inline const mpreal frac       (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(frac      );     }

 

 //////////////////////////////////////////////////////////////////////////

 // Miscellaneous Functions

@@ -2443,14 +2514,14 @@ inline void         swap (mpreal& a, mpreal& b)            {    mpfr_swap(a.mp,b
 inline const mpreal (max)(const mpreal& x, const mpreal& y){    return (x>y?x:y);       }

 inline const mpreal (min)(const mpreal& x, const mpreal& y){    return (x<y?x:y);       }

 

-inline const mpreal fmax(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)

+inline const mpreal fmax(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

 {

     mpreal a;

     mpfr_max(a.mp,x.mp,y.mp,rnd_mode);

     return a;

 }

 

-inline const mpreal fmin(const mpreal& x, const mpreal& y,  mp_rnd_t rnd_mode)

+inline const mpreal fmin(const mpreal& x, const mpreal& y,  mp_rnd_t rnd_mode = mpreal::get_default_rnd())

 {

     mpreal a;

     mpfr_min(a.mp,x.mp,y.mp,rnd_mode);

@@ -2485,16 +2556,16 @@ inline const mpreal urandomb (gmp_randstate_t& state)
     return x;

 }

 

-#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))

+#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,1,0))

 // use gmp_randinit_default() to init state, gmp_randclear() to clear

-inline const mpreal urandom (gmp_randstate_t& state, mp_rnd_t rnd_mode)

+inline const mpreal urandom (gmp_randstate_t& state, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

 {

     mpreal x;

     mpfr_urandom(x.mp,state,rnd_mode);

     return x;

 }

 

-inline const mpreal grandom (gmp_randstate_t& state, mp_rnd_t rnd_mode)

+inline const mpreal grandom (gmp_randstate_t& state, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

 {

     mpreal x;

     mpfr_grandom(x.mp, NULL, state, rnd_mode);

@@ -2516,7 +2587,7 @@ inline const mpreal random2 (mp_size_t size, mp_exp_t exp)
 // a = random(seed); <- initialization & first random number generation

 // a = random();     <- next random numbers generation

 // seed != 0

-inline const mpreal random(unsigned int seed)

+inline const mpreal random(unsigned int seed = 0)

 {

 

 #if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))

@@ -2541,7 +2612,7 @@ inline const mpreal random(unsigned int seed)
 }

 

 #if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))

-inline const mpreal grandom(unsigned int seed)

+inline const mpreal grandom(unsigned int seed = 0)

 {

     static gmp_randstate_t state;

     static bool isFirstTime = true;

@@ -2578,21 +2649,21 @@ inline bool mpreal::fits_in_bits(double x, int n)
     return IsInf(x) || (std::modf ( std::ldexp ( std::frexp ( x, &i ), n ), &t ) == 0.0);

 }

 

-inline const mpreal pow(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode)

+inline const mpreal pow(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

 {

     mpreal x(a);

     mpfr_pow(x.mp,x.mp,b.mp,rnd_mode);

     return x;

 }

 

-inline const mpreal pow(const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode)

+inline const mpreal pow(const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

 {

     mpreal x(a);

     mpfr_pow_z(x.mp,x.mp,b,rnd_mode);

     return x;

 }

 

-inline const mpreal pow(const mpreal& a, const unsigned long int b, mp_rnd_t rnd_mode)

+inline const mpreal pow(const mpreal& a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

 {

     mpreal x(a);

     mpfr_pow_ui(x.mp,x.mp,b,rnd_mode);

@@ -2604,7 +2675,7 @@ inline const mpreal pow(const mpreal& a, const unsigned int b, mp_rnd_t rnd_mode
     return pow(a,static_cast<unsigned long int>(b),rnd_mode);

 }

 

-inline const mpreal pow(const mpreal& a, const long int b, mp_rnd_t rnd_mode)

+inline const mpreal pow(const mpreal& a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

 {

     mpreal x(a);

     mpfr_pow_si(x.mp,x.mp,b,rnd_mode);

@@ -2626,7 +2697,7 @@ inline const mpreal pow(const mpreal& a, const double b, mp_rnd_t rnd_mode)
     return pow(a,mpreal(b),rnd_mode);

 }

 

-inline const mpreal pow(const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode)

+inline const mpreal pow(const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

 {

     mpreal x(a);

     mpfr_ui_pow(x.mp,a,b.mp,rnd_mode);

@@ -2641,13 +2712,13 @@ inline const mpreal pow(const unsigned int a, const mpreal& b, mp_rnd_t rnd_mode
 inline const mpreal pow(const long int a, const mpreal& b, mp_rnd_t rnd_mode)

 {

     if (a>=0)     return pow(static_cast<unsigned long int>(a),b,rnd_mode);

-    else        return pow(mpreal(a),b,rnd_mode);

+    else          return pow(mpreal(a),b,rnd_mode);

 }

 

 inline const mpreal pow(const int a, const mpreal& b, mp_rnd_t rnd_mode)

 {

     if (a>=0)     return pow(static_cast<unsigned long int>(a),b,rnd_mode);

-    else        return pow(mpreal(a),b,rnd_mode);

+    else          return pow(mpreal(a),b,rnd_mode);

 }

 

 inline const mpreal pow(const long double a, const mpreal& b, mp_rnd_t rnd_mode)

@@ -2676,13 +2747,13 @@ inline const mpreal pow(const unsigned long int a, const unsigned int b, mp_rnd_
 inline const mpreal pow(const unsigned long int a, const long int b, mp_rnd_t rnd_mode)

 {

     if(b>0)    return pow(a,static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui

-    else    return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow

+    else       return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow

 }

 

 inline const mpreal pow(const unsigned long int a, const int b, mp_rnd_t rnd_mode)

 {

     if(b>0)    return pow(a,static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui

-    else    return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow

+    else       return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow

 }

 

 inline const mpreal pow(const unsigned long int a, const long double b, mp_rnd_t rnd_mode)

@@ -2879,7 +2950,7 @@ inline const mpreal pow(const double a, const int b, mp_rnd_t rnd_mode)
 // Non-throwing swap C++ idiom: http://en.wikibooks.org/wiki/More_C%2B%2B_Idioms/Non-throwing_swap

 namespace std

 {

-	// only allowed to extend namespace std with specializations

+  // we are allowed to extend namespace std with specializations only

     template <>

     inline void swap(mpfr::mpreal& x, mpfr::mpreal& y) 

     { 

@@ -2907,20 +2978,6 @@ namespace std
         static const bool tinyness_before   = true;

 

         static const float_denorm_style has_denorm  = denorm_absent;

-        

-        inline static float_round_style round_style()

-        {

-            mp_rnd_t r = mpfr::mpreal::get_default_rnd();

-

-            switch (r)

-            {

-                case GMP_RNDN: return round_to_nearest;

-                case GMP_RNDZ: return round_toward_zero; 

-                case GMP_RNDU: return round_toward_infinity; 

-                case GMP_RNDD: return round_toward_neg_infinity; 

-                default: return round_indeterminate;

-            }

-        }

 

         inline static mpfr::mpreal (min)    (mp_prec_t precision = mpfr::mpreal::get_default_prec()) {  return  mpfr::minval(precision);  }

         inline static mpfr::mpreal (max)    (mp_prec_t precision = mpfr::mpreal::get_default_prec()) {  return  mpfr::maxval(precision);  }

@@ -2928,7 +2985,7 @@ namespace std
 

         // Returns smallest eps such that 1 + eps != 1 (classic machine epsilon)

         inline static mpfr::mpreal epsilon(mp_prec_t precision = mpfr::mpreal::get_default_prec()) {  return  mpfr::machine_epsilon(precision); }

-		

+    

         // Returns smallest eps such that x + eps != x (relative machine epsilon)

         inline static mpfr::mpreal epsilon(const mpfr::mpreal& x) {  return mpfr::machine_epsilon(x);  }

 

@@ -2951,10 +3008,30 @@ namespace std
         MPREAL_PERMISSIVE_EXPR static const int min_exponent10 = (int) (MPFR_EMIN_DEFAULT * 0.3010299956639811); 

         MPREAL_PERMISSIVE_EXPR static const int max_exponent10 = (int) (MPFR_EMAX_DEFAULT * 0.3010299956639811); 

 

-        // Should be constant according to standard, but 'digits' depends on precision in MPFR

+#ifdef MPREAL_HAVE_DYNAMIC_STD_NUMERIC_LIMITS

+

+        // Following members should be constant according to standard, but they can be variable in MPFR

+        // So we define them as functions here. 

+        //

+        // This is preferable way for std::numeric_limits<mpfr::mpreal> specialization.

+        // But it is incompatible with standard std::numeric_limits and might not work with other libraries, e.g. boost. 

+        // See below for compatible implementation. 

+        inline static float_round_style round_style()

+        {

+            mp_rnd_t r = mpfr::mpreal::get_default_rnd();

+

+            switch (r)

+            {

+            case GMP_RNDN: return round_to_nearest;

+            case GMP_RNDZ: return round_toward_zero; 

+            case GMP_RNDU: return round_toward_infinity; 

+            case GMP_RNDD: return round_toward_neg_infinity; 

+            default: return round_indeterminate;

+            }

+        }

 

-        inline static int digits()                        {    return mpfr::mpreal::get_default_prec();    }

-        inline static int digits(const mpfr::mpreal& x)   {    return x.getPrecision();                    }

+        inline static int digits()                        {    return int(mpfr::mpreal::get_default_prec());    }

+        inline static int digits(const mpfr::mpreal& x)   {    return x.getPrecision();                         }

 

         inline static int digits10(mp_prec_t precision = mpfr::mpreal::get_default_prec())

         {

@@ -2970,6 +3047,25 @@ namespace std
         {

             return digits10(precision);

         }

+#else

+        // Digits and round_style are NOT constants when it comes to mpreal.

+        // If possible, please use functions digits() and round_style() defined above.

+        //

+        // These (default) values are preserved for compatibility with existing libraries, e.g. boost.

+        // Change them accordingly to your application. 

+        //

+        // For example, if you use 256 bits of precision uniformly in your program, then:

+        // digits       = 256

+        // digits10     = 77 

+        // max_digits10 = 78

+        // 

+        // Approximate formula for decimal digits is: digits10 = floor(log10(2) * digits). See bits2digits() for more details.

+

+        static const std::float_round_style round_style = round_to_nearest;

+        static const int digits       = 53;

+        static const int digits10     = 15;

+        static const int max_digits10 = 16;

+#endif

     };

 

 }

diff --git a/unsupported/test/mpreal_support.cpp b/unsupported/test/mpreal_support.cpp
index bc00382be..1aa9e786a 100644
--- a/unsupported/test/mpreal_support.cpp
+++ b/unsupported/test/mpreal_support.cpp
@@ -32,12 +32,11 @@ void test_mpreal_support()
     VERIFY_IS_APPROX(A.array().abs2().sqrt(), A.array().abs());
     VERIFY_IS_APPROX(A.array().sin(),         sin(A.array()));
     VERIFY_IS_APPROX(A.array().cos(),         cos(A.array()));
-    
 
     // Cholesky
     X = S.selfadjointView<Lower>().llt().solve(B);
     VERIFY_IS_APPROX((S.selfadjointView<Lower>()*X).eval(),B);
-
+    
     // partial LU
     X = A.lu().solve(B);
     VERIFY_IS_APPROX((A*X).eval(),B);
diff --git a/unsupported/test/splines.cpp b/unsupported/test/splines.cpp
index a7eb3e0c4..97665af96 100644
--- a/unsupported/test/splines.cpp
+++ b/unsupported/test/splines.cpp
@@ -13,23 +13,23 @@
 
 namespace Eigen {
   
-// lets do some explicit instantiations and thus
-// force the compilation of all spline functions...
-template class Spline<double, 2, Dynamic>;
-template class Spline<double, 3, Dynamic>;
+  // lets do some explicit instantiations and thus
+  // force the compilation of all spline functions...
+  template class Spline<double, 2, Dynamic>;
+  template class Spline<double, 3, Dynamic>;
 
-template class Spline<double, 2, 2>;
-template class Spline<double, 2, 3>;
-template class Spline<double, 2, 4>;
-template class Spline<double, 2, 5>;
+  template class Spline<double, 2, 2>;
+  template class Spline<double, 2, 3>;
+  template class Spline<double, 2, 4>;
+  template class Spline<double, 2, 5>;
 
-template class Spline<float, 2, Dynamic>;
-template class Spline<float, 3, Dynamic>;
+  template class Spline<float, 2, Dynamic>;
+  template class Spline<float, 3, Dynamic>;
 
-template class Spline<float, 3, 2>;
-template class Spline<float, 3, 3>;
-template class Spline<float, 3, 4>;
-template class Spline<float, 3, 5>;
+  template class Spline<float, 3, 2>;
+  template class Spline<float, 3, 3>;
+  template class Spline<float, 3, 4>;
+  template class Spline<float, 3, 5>;
 
 }
 
@@ -234,11 +234,48 @@ void check_global_interpolation2d()
   }
 }
 
+void check_global_interpolation_with_derivatives2d()
+{
+  typedef Spline2d::PointType PointType;
+  typedef Spline2d::KnotVectorType KnotVectorType;
+
+  const unsigned int numPoints = 100;
+  const unsigned int dimension = 2;
+  const unsigned int degree = 3;
+
+  ArrayXXd points = ArrayXXd::Random(dimension, numPoints);
+
+  KnotVectorType knots;
+  Eigen::ChordLengths(points, knots);
+
+  ArrayXXd derivatives = ArrayXXd::Random(dimension, numPoints);
+  VectorXd derivativeIndices(numPoints);
+
+  for (Eigen::DenseIndex i = 0; i < numPoints; ++i)
+      derivativeIndices(i) = static_cast<double>(i);
+
+  const Spline2d spline = SplineFitting<Spline2d>::InterpolateWithDerivatives(
+    points, derivatives, derivativeIndices, degree);  
+    
+  for (Eigen::DenseIndex i = 0; i < points.cols(); ++i)
+  {
+    PointType point = spline(knots(i));
+    PointType referencePoint = points.col(i);
+    VERIFY_IS_APPROX(point, referencePoint);
+    PointType derivative = spline.derivatives(knots(i), 1).col(1);
+    PointType referenceDerivative = derivatives.col(i);
+    VERIFY_IS_APPROX(derivative, referenceDerivative);
+  }
+}
 
 void test_splines()
 {
-  CALL_SUBTEST( eval_spline3d() );
-  CALL_SUBTEST( eval_spline3d_onbrks() );
-  CALL_SUBTEST( eval_closed_spline2d() );
-  CALL_SUBTEST( check_global_interpolation2d() );
+  for (int i = 0; i < g_repeat; ++i)
+  {
+    CALL_SUBTEST( eval_spline3d() );
+    CALL_SUBTEST( eval_spline3d_onbrks() );
+    CALL_SUBTEST( eval_closed_spline2d() );
+    CALL_SUBTEST( check_global_interpolation2d() );
+    CALL_SUBTEST( check_global_interpolation_with_derivatives2d() );
+  }
 }