Merge "Rebase Eigen to 3.2.4."

92 files changed, 1533 insertions, 839 deletions

diff --git a/Android.mk b/Android.mk
index 4078ec1b8..272cd956e 100644
--- a/Android.mk
+++ b/Android.mk
@@ -11,6 +11,3 @@
 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 # See the License for the specific language governing permissions and
 # limitations under the License.
-
-# This empty file is here solely for the purpose of optimizing the Android build
-# Please keep it there, and empty, thanks.
diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h
index c52b7d1a6..02ab93880 100644
--- a/Eigen/src/Cholesky/LDLT.h
+++ b/Eigen/src/Cholesky/LDLT.h
@@ -442,6 +442,7 @@ LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const MatrixType& a)
   m_transpositions.resize(size);
   m_isInitialized = false;
   m_temporary.resize(size);
+  m_sign = internal::ZeroSign;
 
   internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, m_sign);
 
@@ -502,7 +503,6 @@ struct solve_retval<LDLT<_MatrixType,_UpLo>, Rhs>
     using std::abs;
     using std::max;
     typedef typename LDLTType::MatrixType MatrixType;
-    typedef typename LDLTType::Scalar Scalar;
     typedef typename LDLTType::RealScalar RealScalar;
     const typename Diagonal<const MatrixType>::RealReturnType vectorD(dec().vectorD());
     // In some previous versions, tolerance was set to the max of 1/highest and the maximal diagonal entry * epsilon
diff --git a/Eigen/src/Core/ArrayWrapper.h b/Eigen/src/Core/ArrayWrapper.h
index a791bc358..b4641e2a0 100644
--- a/Eigen/src/Core/ArrayWrapper.h
+++ b/Eigen/src/Core/ArrayWrapper.h
@@ -29,6 +29,11 @@ struct traits<ArrayWrapper<ExpressionType> >
   : public traits<typename remove_all<typename ExpressionType::Nested>::type >
 {
   typedef ArrayXpr XprKind;
+  // Let's remove NestByRefBit
+  enum {
+    Flags0 = traits<typename remove_all<typename ExpressionType::Nested>::type >::Flags,
+    Flags = Flags0 & ~NestByRefBit
+  };
 };
 }
 
@@ -149,6 +154,11 @@ struct traits<MatrixWrapper<ExpressionType> >
  : public traits<typename remove_all<typename ExpressionType::Nested>::type >
 {
   typedef MatrixXpr XprKind;
+  // Let's remove NestByRefBit
+  enum {
+    Flags0 = traits<typename remove_all<typename ExpressionType::Nested>::type >::Flags,
+    Flags = Flags0 & ~NestByRefBit
+  };
 };
 }
 
diff --git a/Eigen/src/Core/DenseBase.h b/Eigen/src/Core/DenseBase.h
index c5800f6c8..04862f374 100644
--- a/Eigen/src/Core/DenseBase.h
+++ b/Eigen/src/Core/DenseBase.h
@@ -462,8 +462,10 @@ template<typename Derived> class DenseBase
     template<int p> RealScalar lpNorm() const;
 
     template<int RowFactor, int ColFactor>
-    const Replicate<Derived,RowFactor,ColFactor> replicate() const;
-    const Replicate<Derived,Dynamic,Dynamic> replicate(Index rowFacor,Index colFactor) const;
+    inline const Replicate<Derived,RowFactor,ColFactor> replicate() const;
+    
+    typedef Replicate<Derived,Dynamic,Dynamic> ReplicateReturnType;
+    inline const ReplicateReturnType replicate(Index rowFacor,Index colFactor) const;
 
     typedef Reverse<Derived, BothDirections> ReverseReturnType;
     typedef const Reverse<const Derived, BothDirections> ConstReverseReturnType;
diff --git a/Eigen/src/Core/Diagonal.h b/Eigen/src/Core/Diagonal.h
index aab8007b3..68cf6d4b0 100644
--- a/Eigen/src/Core/Diagonal.h
+++ b/Eigen/src/Core/Diagonal.h
@@ -190,18 +190,18 @@ MatrixBase<Derived>::diagonal() const
   *
   * \sa MatrixBase::diagonal(), class Diagonal */
 template<typename Derived>
-inline typename MatrixBase<Derived>::template DiagonalIndexReturnType<DynamicIndex>::Type
+inline typename MatrixBase<Derived>::DiagonalDynamicIndexReturnType
 MatrixBase<Derived>::diagonal(Index index)
 {
-  return typename DiagonalIndexReturnType<DynamicIndex>::Type(derived(), index);
+  return DiagonalDynamicIndexReturnType(derived(), index);
 }
 
 /** This is the const version of diagonal(Index). */
 template<typename Derived>
-inline typename MatrixBase<Derived>::template ConstDiagonalIndexReturnType<DynamicIndex>::Type
+inline typename MatrixBase<Derived>::ConstDiagonalDynamicIndexReturnType
 MatrixBase<Derived>::diagonal(Index index) const
 {
-  return typename ConstDiagonalIndexReturnType<DynamicIndex>::Type(derived(), index);
+  return ConstDiagonalDynamicIndexReturnType(derived(), index);
 }
 
 /** \returns an expression of the \a DiagIndex-th sub or super diagonal of the matrix \c *this
diff --git a/Eigen/src/Core/GeneralProduct.h b/Eigen/src/Core/GeneralProduct.h
index 2a59d9464..9e805a80f 100644
--- a/Eigen/src/Core/GeneralProduct.h
+++ b/Eigen/src/Core/GeneralProduct.h
@@ -232,7 +232,7 @@ EIGEN_DONT_INLINE void outer_product_selector_run(const ProductType& prod, Dest&
   // FIXME not very good if rhs is real and lhs complex while alpha is real too
   const Index cols = dest.cols();
   for (Index j=0; j<cols; ++j)
-    func(dest.col(j), prod.rhs().coeff(j) * prod.lhs());
+    func(dest.col(j), prod.rhs().coeff(0,j) * prod.lhs());
 }
 
 // Row major
@@ -243,7 +243,7 @@ EIGEN_DONT_INLINE void outer_product_selector_run(const ProductType& prod, Dest&
   // FIXME not very good if lhs is real and rhs complex while alpha is real too
   const Index rows = dest.rows();
   for (Index i=0; i<rows; ++i)
-    func(dest.row(i), prod.lhs().coeff(i) * prod.rhs());
+    func(dest.row(i), prod.lhs().coeff(i,0) * prod.rhs());
 }
 
 template<typename Lhs, typename Rhs>
diff --git a/Eigen/src/Core/MapBase.h b/Eigen/src/Core/MapBase.h
index ab50c9b81..cebed2bb6 100644
--- a/Eigen/src/Core/MapBase.h
+++ b/Eigen/src/Core/MapBase.h
@@ -168,6 +168,7 @@ template<typename Derived> class MapBase<Derived, ReadOnlyAccessors>
 template<typename Derived> class MapBase<Derived, WriteAccessors>
   : public MapBase<Derived, ReadOnlyAccessors>
 {
+    typedef MapBase<Derived, ReadOnlyAccessors> ReadOnlyMapBase;
   public:
 
     typedef MapBase<Derived, ReadOnlyAccessors> Base;
@@ -230,11 +231,13 @@ template<typename Derived> class MapBase<Derived, WriteAccessors>
 
     Derived& operator=(const MapBase& other)
     {
-      Base::Base::operator=(other);
+      ReadOnlyMapBase::Base::operator=(other);
       return derived();
     }
 
-    using Base::Base::operator=;
+    // In theory we could simply refer to Base:Base::operator=, but MSVC does not like Base::Base,
+    // see bugs 821 and 920.
+    using ReadOnlyMapBase::Base::operator=;
 };
 
 #undef EIGEN_STATIC_ASSERT_INDEX_BASED_ACCESS
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index 344b38f2f..cc3279746 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -215,7 +215,7 @@ template<typename Derived> class MatrixBase
 
     typedef Diagonal<Derived> DiagonalReturnType;
     DiagonalReturnType diagonal();
-	typedef typename internal::add_const<Diagonal<const Derived> >::type ConstDiagonalReturnType;
+    typedef typename internal::add_const<Diagonal<const Derived> >::type ConstDiagonalReturnType;
     ConstDiagonalReturnType diagonal() const;
 
     template<int Index> struct DiagonalIndexReturnType { typedef Diagonal<Derived,Index> Type; };
@@ -223,16 +223,12 @@ template<typename Derived> class MatrixBase
 
     template<int Index> typename DiagonalIndexReturnType<Index>::Type diagonal();
     template<int Index> typename ConstDiagonalIndexReturnType<Index>::Type diagonal() const;
+    
+    typedef Diagonal<Derived,DynamicIndex> DiagonalDynamicIndexReturnType;
+    typedef typename internal::add_const<Diagonal<const Derived,DynamicIndex> >::type ConstDiagonalDynamicIndexReturnType;
 
-    // Note: The "MatrixBase::" prefixes are added to help MSVC9 to match these declarations with the later implementations.
-    // On the other hand they confuse MSVC8...
-    #if (defined _MSC_VER) && (_MSC_VER >= 1500) // 2008 or later
-    typename MatrixBase::template DiagonalIndexReturnType<DynamicIndex>::Type diagonal(Index index);
-    typename MatrixBase::template ConstDiagonalIndexReturnType<DynamicIndex>::Type diagonal(Index index) const;
-    #else
-    typename DiagonalIndexReturnType<DynamicIndex>::Type diagonal(Index index);
-    typename ConstDiagonalIndexReturnType<DynamicIndex>::Type diagonal(Index index) const;
-    #endif
+    DiagonalDynamicIndexReturnType diagonal(Index index);
+    ConstDiagonalDynamicIndexReturnType diagonal(Index index) const;
 
     #ifdef EIGEN2_SUPPORT
     template<unsigned int Mode> typename internal::eigen2_part_return_type<Derived, Mode>::type part();
diff --git a/Eigen/src/Core/PermutationMatrix.h b/Eigen/src/Core/PermutationMatrix.h
index 1297b8413..f26f3e5cc 100644
--- a/Eigen/src/Core/PermutationMatrix.h
+++ b/Eigen/src/Core/PermutationMatrix.h
@@ -555,7 +555,10 @@ struct permut_matrix_product_retval
       const Index n = Side==OnTheLeft ? rows() : cols();
       // FIXME we need an is_same for expression that is not sensitive to constness. For instance
       // is_same_xpr<Block<const Matrix>, Block<Matrix> >::value should be true.
-      if(is_same<MatrixTypeNestedCleaned,Dest>::value && extract_data(dst) == extract_data(m_matrix))
+      if(    is_same<MatrixTypeNestedCleaned,Dest>::value
+          && blas_traits<MatrixTypeNestedCleaned>::HasUsableDirectAccess
+          && blas_traits<Dest>::HasUsableDirectAccess
+          && extract_data(dst) == extract_data(m_matrix))
       {
         // apply the permutation inplace
         Matrix<bool,PermutationType::RowsAtCompileTime,1,0,PermutationType::MaxRowsAtCompileTime> mask(m_permutation.size());
diff --git a/Eigen/src/Core/ProductBase.h b/Eigen/src/Core/ProductBase.h
index a494b5f87..cf74470a9 100644
--- a/Eigen/src/Core/ProductBase.h
+++ b/Eigen/src/Core/ProductBase.h
@@ -85,7 +85,14 @@ class ProductBase : public MatrixBase<Derived>
 
   public:
 
+#ifndef EIGEN_NO_MALLOC
+    typedef typename Base::PlainObject BasePlainObject;
+    typedef Matrix<Scalar,RowsAtCompileTime==1?1:Dynamic,ColsAtCompileTime==1?1:Dynamic,BasePlainObject::Options> DynPlainObject;
+    typedef typename internal::conditional<(BasePlainObject::SizeAtCompileTime==Dynamic) || (BasePlainObject::SizeAtCompileTime*int(sizeof(Scalar)) < int(EIGEN_STACK_ALLOCATION_LIMIT)),
+                                           BasePlainObject, DynPlainObject>::type PlainObject;
+#else
     typedef typename Base::PlainObject PlainObject;
+#endif
 
     ProductBase(const Lhs& a_lhs, const Rhs& a_rhs)
       : m_lhs(a_lhs), m_rhs(a_rhs)
@@ -180,7 +187,12 @@ namespace internal {
 template<typename Lhs, typename Rhs, int Mode, int N, typename PlainObject>
 struct nested<GeneralProduct<Lhs,Rhs,Mode>, N, PlainObject>
 {
-  typedef PlainObject const& type;
+  typedef typename GeneralProduct<Lhs,Rhs,Mode>::PlainObject const& type;
+};
+template<typename Lhs, typename Rhs, int Mode, int N, typename PlainObject>
+struct nested<const GeneralProduct<Lhs,Rhs,Mode>, N, PlainObject>
+{
+  typedef typename GeneralProduct<Lhs,Rhs,Mode>::PlainObject const& type;
 };
 }
 
diff --git a/Eigen/src/Core/Ref.h b/Eigen/src/Core/Ref.h
index cd6d949c4..df87b1af4 100644
--- a/Eigen/src/Core/Ref.h
+++ b/Eigen/src/Core/Ref.h
@@ -188,6 +188,8 @@ template<typename PlainObjectType, int Options, typename StrideType> class Ref
   : public RefBase<Ref<PlainObjectType, Options, StrideType> >
 {
     typedef internal::traits<Ref> Traits;
+    template<typename Derived>
+    inline Ref(const PlainObjectBase<Derived>& expr);
   public:
 
     typedef RefBase<Ref> Base;
@@ -196,20 +198,21 @@ template<typename PlainObjectType, int Options, typename StrideType> class Ref
 
     #ifndef EIGEN_PARSED_BY_DOXYGEN
     template<typename Derived>
-    inline Ref(PlainObjectBase<Derived>& expr,
-               typename internal::enable_if<bool(Traits::template match<Derived>::MatchAtCompileTime),Derived>::type* = 0)
+    inline Ref(PlainObjectBase<Derived>& expr)
     {
-      Base::construct(expr);
+      EIGEN_STATIC_ASSERT(static_cast<bool>(Traits::template match<Derived>::MatchAtCompileTime), STORAGE_LAYOUT_DOES_NOT_MATCH);
+      Base::construct(expr.derived());
     }
     template<typename Derived>
-    inline Ref(const DenseBase<Derived>& expr,
-               typename internal::enable_if<bool(internal::is_lvalue<Derived>::value&&bool(Traits::template match<Derived>::MatchAtCompileTime)),Derived>::type* = 0,
-               int = Derived::ThisConstantIsPrivateInPlainObjectBase)
+    inline Ref(const DenseBase<Derived>& expr)
     #else
     template<typename Derived>
     inline Ref(DenseBase<Derived>& expr)
     #endif
     {
+      EIGEN_STATIC_ASSERT(static_cast<bool>(internal::is_lvalue<Derived>::value), THIS_EXPRESSION_IS_NOT_A_LVALUE__IT_IS_READ_ONLY);
+      EIGEN_STATIC_ASSERT(static_cast<bool>(Traits::template match<Derived>::MatchAtCompileTime), STORAGE_LAYOUT_DOES_NOT_MATCH);
+      enum { THIS_EXPRESSION_IS_NOT_A_LVALUE__IT_IS_READ_ONLY = Derived::ThisConstantIsPrivateInPlainObjectBase};
       Base::construct(expr.const_cast_derived());
     }
 
diff --git a/Eigen/src/Core/Replicate.h b/Eigen/src/Core/Replicate.h
index dde86a834..ac4537c14 100644
--- a/Eigen/src/Core/Replicate.h
+++ b/Eigen/src/Core/Replicate.h
@@ -135,7 +135,7 @@ template<typename MatrixType,int RowFactor,int ColFactor> class Replicate
   */
 template<typename Derived>
 template<int RowFactor, int ColFactor>
-inline const Replicate<Derived,RowFactor,ColFactor>
+const Replicate<Derived,RowFactor,ColFactor>
 DenseBase<Derived>::replicate() const
 {
   return Replicate<Derived,RowFactor,ColFactor>(derived());
@@ -150,7 +150,7 @@ DenseBase<Derived>::replicate() const
   * \sa VectorwiseOp::replicate(), DenseBase::replicate<int,int>(), class Replicate
   */
 template<typename Derived>
-inline const Replicate<Derived,Dynamic,Dynamic>
+const typename DenseBase<Derived>::ReplicateReturnType
 DenseBase<Derived>::replicate(Index rowFactor,Index colFactor) const
 {
   return Replicate<Derived,Dynamic,Dynamic>(derived(),rowFactor,colFactor);
diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h
index 845ae1aec..4d65392c6 100644
--- a/Eigen/src/Core/TriangularMatrix.h
+++ b/Eigen/src/Core/TriangularMatrix.h
@@ -380,19 +380,19 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView
     EIGEN_STRONG_INLINE TriangularView& operator=(const ProductBase<ProductDerived, Lhs,Rhs>& other)
     {
       setZero();
-      return assignProduct(other,1);
+      return assignProduct(other.derived(),1);
     }
     
     template<typename ProductDerived, typename Lhs, typename Rhs>
     EIGEN_STRONG_INLINE TriangularView& operator+=(const ProductBase<ProductDerived, Lhs,Rhs>& other)
     {
-      return assignProduct(other,1);
+      return assignProduct(other.derived(),1);
     }
     
     template<typename ProductDerived, typename Lhs, typename Rhs>
     EIGEN_STRONG_INLINE TriangularView& operator-=(const ProductBase<ProductDerived, Lhs,Rhs>& other)
     {
-      return assignProduct(other,-1);
+      return assignProduct(other.derived(),-1);
     }
     
     
@@ -400,25 +400,34 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView
     EIGEN_STRONG_INLINE TriangularView& operator=(const ScaledProduct<ProductDerived>& other)
     {
       setZero();
-      return assignProduct(other,other.alpha());
+      return assignProduct(other.derived(),other.alpha());
     }
     
     template<typename ProductDerived>
     EIGEN_STRONG_INLINE TriangularView& operator+=(const ScaledProduct<ProductDerived>& other)
     {
-      return assignProduct(other,other.alpha());
+      return assignProduct(other.derived(),other.alpha());
     }
     
     template<typename ProductDerived>
     EIGEN_STRONG_INLINE TriangularView& operator-=(const ScaledProduct<ProductDerived>& other)
     {
-      return assignProduct(other,-other.alpha());
+      return assignProduct(other.derived(),-other.alpha());
     }
     
   protected:
     
     template<typename ProductDerived, typename Lhs, typename Rhs>
     EIGEN_STRONG_INLINE TriangularView& assignProduct(const ProductBase<ProductDerived, Lhs,Rhs>& prod, const Scalar& alpha);
+    
+    template<int Mode, bool LhsIsTriangular,
+         typename Lhs, bool LhsIsVector,
+         typename Rhs, bool RhsIsVector>
+    EIGEN_STRONG_INLINE TriangularView& assignProduct(const TriangularProduct<Mode, LhsIsTriangular, Lhs, LhsIsVector, Rhs, RhsIsVector>& prod, const Scalar& alpha)
+    {
+      lazyAssign(alpha*prod.eval());
+      return *this;
+    }
 
     MatrixTypeNested m_matrix;
 };
diff --git a/Eigen/src/Core/arch/NEON/Complex.h b/Eigen/src/Core/arch/NEON/Complex.h
index f183d31de..8d9255eef 100644
--- a/Eigen/src/Core/arch/NEON/Complex.h
+++ b/Eigen/src/Core/arch/NEON/Complex.h
@@ -110,7 +110,7 @@ template<> EIGEN_STRONG_INLINE Packet2cf ploaddup<Packet2cf>(const std::complex<
 template<> EIGEN_STRONG_INLINE void pstore <std::complex<float> >(std::complex<float> *   to, const Packet2cf& from) { EIGEN_DEBUG_ALIGNED_STORE pstore((float*)to, from.v); }
 template<> EIGEN_STRONG_INLINE void pstoreu<std::complex<float> >(std::complex<float> *   to, const Packet2cf& from) { EIGEN_DEBUG_UNALIGNED_STORE pstoreu((float*)to, from.v); }
 
-template<> EIGEN_STRONG_INLINE void prefetch<std::complex<float> >(const std::complex<float> *   addr) { __pld((float *)addr); }
+template<> EIGEN_STRONG_INLINE void prefetch<std::complex<float> >(const std::complex<float> *   addr) { EIGEN_ARM_PREFETCH((float *)addr); }
 
 template<> EIGEN_STRONG_INLINE std::complex<float>  pfirst<Packet2cf>(const Packet2cf& a)
 {
diff --git a/Eigen/src/Core/arch/NEON/PacketMath.h b/Eigen/src/Core/arch/NEON/PacketMath.h
index 163bac215..94dfab330 100644
--- a/Eigen/src/Core/arch/NEON/PacketMath.h
+++ b/Eigen/src/Core/arch/NEON/PacketMath.h
@@ -48,9 +48,18 @@ typedef uint32x4_t  Packet4ui;
   #define EIGEN_INIT_NEON_PACKET2(X, Y)       {X, Y}
   #define EIGEN_INIT_NEON_PACKET4(X, Y, Z, W) {X, Y, Z, W}
 #endif
-    
-#ifndef __pld
-#define __pld(x) asm volatile ( "   pld [%[addr]]\n" :: [addr] "r" (x) : "cc" );
+
+// arm64 does have the pld instruction. If available, let's trust the __builtin_prefetch built-in function
+// which available on LLVM and GCC (at least)
+#if EIGEN_HAS_BUILTIN(__builtin_prefetch) || defined(__GNUC__)
+  #define EIGEN_ARM_PREFETCH(ADDR) __builtin_prefetch(ADDR);
+#elif defined __pld
+  #define EIGEN_ARM_PREFETCH(ADDR) __pld(ADDR)
+#elif !defined(__aarch64__)
+  #define EIGEN_ARM_PREFETCH(ADDR) __asm__ __volatile__ ( "   pld [%[addr]]\n" :: [addr] "r" (ADDR) : "cc" );
+#else
+  // by default no explicit prefetching
+  #define EIGEN_ARM_PREFETCH(ADDR)
 #endif
 
 template<> struct packet_traits<float>  : default_packet_traits
@@ -209,8 +218,8 @@ template<> EIGEN_STRONG_INLINE void pstore<int>(int*       to, const Packet4i& f
 template<> EIGEN_STRONG_INLINE void pstoreu<float>(float*  to, const Packet4f& from) { EIGEN_DEBUG_UNALIGNED_STORE vst1q_f32(to, from); }
 template<> EIGEN_STRONG_INLINE void pstoreu<int>(int*      to, const Packet4i& from) { EIGEN_DEBUG_UNALIGNED_STORE vst1q_s32(to, from); }
 
-template<> EIGEN_STRONG_INLINE void prefetch<float>(const float* addr) { __pld(addr); }
-template<> EIGEN_STRONG_INLINE void prefetch<int>(const int*     addr) { __pld(addr); }
+template<> EIGEN_STRONG_INLINE void prefetch<float>(const float* addr) { EIGEN_ARM_PREFETCH(addr); }
+template<> EIGEN_STRONG_INLINE void prefetch<int>(const int*     addr) { EIGEN_ARM_PREFETCH(addr); }
 
 // FIXME only store the 2 first elements ?
 template<> EIGEN_STRONG_INLINE float  pfirst<Packet4f>(const Packet4f& a) { float EIGEN_ALIGN16 x[4]; vst1q_f32(x, a); return x[0]; }
diff --git a/Eigen/src/Core/arch/SSE/MathFunctions.h b/Eigen/src/Core/arch/SSE/MathFunctions.h
index 99cbd0d95..d16f30bb0 100644
--- a/Eigen/src/Core/arch/SSE/MathFunctions.h
+++ b/Eigen/src/Core/arch/SSE/MathFunctions.h
@@ -52,7 +52,7 @@ Packet4f plog<Packet4f>(const Packet4f& _x)
 
   Packet4i emm0;
 
-  Packet4f invalid_mask = _mm_cmplt_ps(x, _mm_setzero_ps());
+  Packet4f invalid_mask = _mm_cmpnge_ps(x, _mm_setzero_ps()); // not greater equal is true if x is NaN
   Packet4f iszero_mask = _mm_cmpeq_ps(x, _mm_setzero_ps());
 
   x = pmax(x, p4f_min_norm_pos);  /* cut off denormalized stuff */
@@ -166,7 +166,7 @@ Packet4f pexp<Packet4f>(const Packet4f& _x)
   emm0 = _mm_cvttps_epi32(fx);
   emm0 = _mm_add_epi32(emm0, p4i_0x7f);
   emm0 = _mm_slli_epi32(emm0, 23);
-  return pmul(y, _mm_castsi128_ps(emm0));
+  return pmax(pmul(y, Packet4f(_mm_castsi128_ps(emm0))), _x);
 }
 template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
 Packet2d pexp<Packet2d>(const Packet2d& _x)
@@ -239,7 +239,7 @@ Packet2d pexp<Packet2d>(const Packet2d& _x)
   emm0 = _mm_add_epi32(emm0, p4i_1023_0);
   emm0 = _mm_slli_epi32(emm0, 20);
   emm0 = _mm_shuffle_epi32(emm0, _MM_SHUFFLE(1,2,0,3));
-  return pmul(x, _mm_castsi128_pd(emm0));
+  return pmax(pmul(x, Packet2d(_mm_castsi128_pd(emm0))), _x);
 }
 
 /* evaluation of 4 sines at onces, using SSE2 intrinsics.
diff --git a/Eigen/src/Core/products/CoeffBasedProduct.h b/Eigen/src/Core/products/CoeffBasedProduct.h
index c06a0df1c..421f925e1 100644
--- a/Eigen/src/Core/products/CoeffBasedProduct.h
+++ b/Eigen/src/Core/products/CoeffBasedProduct.h
@@ -90,6 +90,7 @@ struct traits<CoeffBasedProduct<LhsNested,RhsNested,NestingFlags> >
             | (SameType && (CanVectorizeLhs || CanVectorizeRhs) ? PacketAccessBit : 0),
 
       CoeffReadCost = InnerSize == Dynamic ? Dynamic
+                    : InnerSize == 0 ? 0
                     : InnerSize * (NumTraits<Scalar>::MulCost + LhsCoeffReadCost + RhsCoeffReadCost)
                       + (InnerSize - 1) * NumTraits<Scalar>::AddCost,
 
@@ -133,7 +134,7 @@ class CoeffBasedProduct
     };
 
     typedef internal::product_coeff_impl<CanVectorizeInner ? InnerVectorizedTraversal : DefaultTraversal,
-                                   Unroll ? InnerSize-1 : Dynamic,
+                                   Unroll ? (InnerSize==0 ? 0 : InnerSize-1) : Dynamic,
                                    _LhsNested, _RhsNested, Scalar> ScalarCoeffImpl;
 
     typedef CoeffBasedProduct<LhsNested,RhsNested,NestByRefBit> LazyCoeffBasedProductType;
@@ -184,7 +185,7 @@ class CoeffBasedProduct
     {
       PacketScalar res;
       internal::product_packet_impl<Flags&RowMajorBit ? RowMajor : ColMajor,
-                              Unroll ? InnerSize-1 : Dynamic,
+                              Unroll ? (InnerSize==0 ? 0 : InnerSize-1) : Dynamic,
                               _LhsNested, _RhsNested, PacketScalar, LoadMode>
         ::run(row, col, m_lhs, m_rhs, res);
       return res;
@@ -262,10 +263,7 @@ struct product_coeff_impl<DefaultTraversal, Dynamic, Lhs, Rhs, RetScalar>
   typedef typename Lhs::Index Index;
   static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, RetScalar& res)
   {
-    eigen_assert(lhs.cols()>0 && "you are using a non initialized matrix");
-    res = lhs.coeff(row, 0) * rhs.coeff(0, col);
-      for(Index i = 1; i < lhs.cols(); ++i)
-        res += lhs.coeff(row, i) * rhs.coeff(i, col);
+    res = (lhs.row(row).transpose().cwiseProduct( rhs.col(col) )).sum();
   }
 };
 
diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h
index 3a010ec6a..6d1e6c133 100644
--- a/Eigen/src/Core/util/Macros.h
+++ b/Eigen/src/Core/util/Macros.h
@@ -13,7 +13,7 @@
 
 #define EIGEN_WORLD_VERSION 3
 #define EIGEN_MAJOR_VERSION 2
-#define EIGEN_MINOR_VERSION 2
+#define EIGEN_MINOR_VERSION 4
 
 #define EIGEN_VERSION_AT_LEAST(x,y,z) (EIGEN_WORLD_VERSION>x || (EIGEN_WORLD_VERSION>=x && \
                                       (EIGEN_MAJOR_VERSION>y || (EIGEN_MAJOR_VERSION>=y && \
@@ -96,6 +96,13 @@
 #define EIGEN_DEFAULT_DENSE_INDEX_TYPE std::ptrdiff_t
 #endif
 
+// Cross compiler wrapper around LLVM's __has_builtin
+#ifdef __has_builtin
+#  define EIGEN_HAS_BUILTIN(x) __has_builtin(x)
+#else
+#  define EIGEN_HAS_BUILTIN(x) 0
+#endif
+
 /** Allows to disable some optimizations which might affect the accuracy of the result.
   * Such optimization are enabled by default, and set EIGEN_FAST_MATH to 0 to disable them.
   * They currently include:
@@ -247,7 +254,7 @@ namespace Eigen {
 
 #if !defined(EIGEN_ASM_COMMENT)
   #if (defined __GNUC__) && ( defined(__i386__) || defined(__x86_64__) )
-    #define EIGEN_ASM_COMMENT(X)  asm("#" X)
+    #define EIGEN_ASM_COMMENT(X)  __asm__("#" X)
   #else
     #define EIGEN_ASM_COMMENT(X)
   #endif
@@ -306,7 +313,7 @@ namespace Eigen {
 // just an empty macro !
 #define EIGEN_EMPTY
 
-#if defined(_MSC_VER) && (!defined(__INTEL_COMPILER))
+#if defined(_MSC_VER) && (_MSC_VER < 1900) && (!defined(__INTEL_COMPILER))
 #define EIGEN_INHERIT_ASSIGNMENT_EQUAL_OPERATOR(Derived) \
   using Base::operator =;
 #elif defined(__clang__) // workaround clang bug (see http://forum.kde.org/viewtopic.php?f=74&t=102653)
diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h
index c4011245a..330bcf518 100644
--- a/Eigen/src/Core/util/Memory.h
+++ b/Eigen/src/Core/util/Memory.h
@@ -63,7 +63,7 @@
 // Currently, let's include it only on unix systems:
 #if defined(__unix__) || defined(__unix)
   #include <unistd.h>
-  #if ((defined __QNXNTO__) || (defined _GNU_SOURCE) || ((defined _XOPEN_SOURCE) && (_XOPEN_SOURCE >= 600))) && (defined _POSIX_ADVISORY_INFO) && (_POSIX_ADVISORY_INFO > 0)
+  #if ((defined __QNXNTO__) || (defined _GNU_SOURCE) || (defined __PGI) || ((defined _XOPEN_SOURCE) && (_XOPEN_SOURCE >= 600))) && (defined _POSIX_ADVISORY_INFO) && (_POSIX_ADVISORY_INFO > 0)
     #define EIGEN_HAS_POSIX_MEMALIGN 1
   #endif
 #endif
@@ -417,6 +417,8 @@ template<typename T, bool Align> inline T* conditional_aligned_realloc_new(T* pt
 
 template<typename T, bool Align> inline T* conditional_aligned_new_auto(size_t size)
 {
+  if(size==0)
+    return 0; // short-cut. Also fixes Bug 884
   check_size_for_overflow<T>(size);
   T *result = reinterpret_cast<T*>(conditional_aligned_malloc<Align>(sizeof(T)*size));
   if(NumTraits<T>::RequireInitialization)
@@ -464,9 +466,8 @@ template<typename T, bool Align> inline void conditional_aligned_delete_auto(T *
 template<typename Scalar, typename Index>
 static inline Index first_aligned(const Scalar* array, Index size)
 {
-  enum { PacketSize = packet_traits<Scalar>::size,
-         PacketAlignedMask = PacketSize-1
-  };
+  static const Index PacketSize = packet_traits<Scalar>::size;
+  static const Index PacketAlignedMask = PacketSize-1;
 
   if(PacketSize==1)
   {
@@ -612,7 +613,6 @@ template<typename T> class aligned_stack_memory_handler
       void* operator new(size_t size, const std::nothrow_t&) throw() { \
         try { return Eigen::internal::conditional_aligned_malloc<NeedsToAlign>(size); } \
         catch (...) { return 0; } \
-        return 0; \
       }
   #else
     #define EIGEN_MAKE_ALIGNED_OPERATOR_NEW_NOTHROW(NeedsToAlign) \
diff --git a/Eigen/src/Core/util/StaticAssert.h b/Eigen/src/Core/util/StaticAssert.h
index 8872c5b64..bac5d9fe9 100644
--- a/Eigen/src/Core/util/StaticAssert.h
+++ b/Eigen/src/Core/util/StaticAssert.h
@@ -90,7 +90,9 @@
         YOU_PASSED_A_COLUMN_VECTOR_BUT_A_ROW_VECTOR_WAS_EXPECTED,
         THE_INDEX_TYPE_MUST_BE_A_SIGNED_TYPE,
         THE_STORAGE_ORDER_OF_BOTH_SIDES_MUST_MATCH,
-        OBJECT_ALLOCATED_ON_STACK_IS_TOO_BIG
+        OBJECT_ALLOCATED_ON_STACK_IS_TOO_BIG,
+        IMPLICIT_CONVERSION_TO_SCALAR_IS_FOR_INNER_PRODUCT_ONLY,
+        STORAGE_LAYOUT_DOES_NOT_MATCH
       };
     };
 
diff --git a/Eigen/src/Core/util/XprHelper.h b/Eigen/src/Core/util/XprHelper.h
index 3c4773054..781965d2c 100644
--- a/Eigen/src/Core/util/XprHelper.h
+++ b/Eigen/src/Core/util/XprHelper.h
@@ -341,7 +341,7 @@ template<typename T, int n=1, typename PlainObject = typename eval<T>::type> str
 };
 
 template<typename T>
-T* const_cast_ptr(const T* ptr)
+inline T* const_cast_ptr(const T* ptr)
 {
   return const_cast<T*>(ptr);
 }
diff --git a/Eigen/src/Eigen2Support/LeastSquares.h b/Eigen/src/Eigen2Support/LeastSquares.h
index 0e6fdb488..7992d4944 100644
--- a/Eigen/src/Eigen2Support/LeastSquares.h
+++ b/Eigen/src/Eigen2Support/LeastSquares.h
@@ -147,7 +147,6 @@ void fitHyperplane(int numPoints,
 
   // compute the covariance matrix
   CovMatrixType covMat = CovMatrixType::Zero(size, size);
-  VectorType remean = VectorType::Zero(size);
   for(int i = 0; i < numPoints; ++i)
   {
     VectorType diff = (*(points[i]) - mean).conjugate();
diff --git a/Eigen/src/Eigenvalues/RealQZ.h b/Eigen/src/Eigenvalues/RealQZ.h
index 5706eeebe..4f36091db 100644
--- a/Eigen/src/Eigenvalues/RealQZ.h
+++ b/Eigen/src/Eigenvalues/RealQZ.h
@@ -313,7 +313,7 @@ namespace Eigen {
       using std::abs;
       using std::sqrt;
       const Index dim=m_S.cols();
-      if (abs(m_S.coeff(i+1,i)==Scalar(0)))
+      if (abs(m_S.coeff(i+1,i))==Scalar(0))
         return;
       Index z = findSmallDiagEntry(i,i+1);
       if (z==i-1)
diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
index 3993046a8..be89de4a9 100644
--- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
+++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
@@ -563,7 +563,6 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3
     if(computeEigenvectors)
     {
       Scalar safeNorm2 = Eigen::NumTraits<Scalar>::epsilon();
-      safeNorm2 *= safeNorm2;
       if((eivals(2)-eivals(0))<=Eigen::NumTraits<Scalar>::epsilon())
       {
         eivecs.setIdentity();
@@ -577,7 +576,7 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3
         Scalar d0 = eivals(2) - eivals(1);
         Scalar d1 = eivals(1) - eivals(0);
         int k =  d0 > d1 ? 2 : 0;
-        d0 = d0 > d1 ? d1 : d0;
+        d0 = d0 > d1 ? d0 : d1;
 
         tmp.diagonal().array () -= eivals(k);
         VectorType cross;
@@ -585,19 +584,25 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3
         n = (cross = tmp.row(0).cross(tmp.row(1))).squaredNorm();
 
         if(n>safeNorm2)
+        {
           eivecs.col(k) = cross / sqrt(n);
+        }
         else
         {
           n = (cross = tmp.row(0).cross(tmp.row(2))).squaredNorm();
 
           if(n>safeNorm2)
+          {
             eivecs.col(k) = cross / sqrt(n);
+          }
           else
           {
             n = (cross = tmp.row(1).cross(tmp.row(2))).squaredNorm();
 
             if(n>safeNorm2)
+            {
               eivecs.col(k) = cross / sqrt(n);
+            }
             else
             {
               // the input matrix and/or the eigenvaues probably contains some inf/NaN,
@@ -617,12 +622,16 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3
         tmp.diagonal().array() -= eivals(1);
 
         if(d0<=Eigen::NumTraits<Scalar>::epsilon())
+        {
           eivecs.col(1) = eivecs.col(k).unitOrthogonal();
+        }
         else
         {
-          n = (cross = eivecs.col(k).cross(tmp.row(0).normalized())).squaredNorm();
+          n = (cross = eivecs.col(k).cross(tmp.row(0))).squaredNorm();
           if(n>safeNorm2)
+          {
             eivecs.col(1) = cross / sqrt(n);
+          }
           else
           {
             n = (cross = eivecs.col(k).cross(tmp.row(1))).squaredNorm();
@@ -636,13 +645,14 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3
               else
               {
                 // we should never reach this point,
-                // if so the last two eigenvalues are likely to ve very closed to each other
+                // if so the last two eigenvalues are likely to be very close to each other
                 eivecs.col(1) = eivecs.col(k).unitOrthogonal();
               }
             }
           }
 
           // make sure that eivecs[1] is orthogonal to eivecs[2]
+          // FIXME: this step should not be needed
           Scalar d = eivecs.col(1).dot(eivecs.col(k));
           eivecs.col(1) = (eivecs.col(1) - d * eivecs.col(k)).normalized();
         }
diff --git a/Eigen/src/Geometry/Hyperplane.h b/Eigen/src/Geometry/Hyperplane.h
index aeff43fef..00b7c4300 100644
--- a/Eigen/src/Geometry/Hyperplane.h
+++ b/Eigen/src/Geometry/Hyperplane.h
@@ -100,7 +100,17 @@ public:
   {
     EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(VectorType, 3)
     Hyperplane result(p0.size());
-    result.normal() = (p2 - p0).cross(p1 - p0).normalized();
+    VectorType v0(p2 - p0), v1(p1 - p0);
+    result.normal() = v0.cross(v1);
+    RealScalar norm = result.normal().norm();
+    if(norm <= v0.norm() * v1.norm() * NumTraits<RealScalar>::epsilon())
+    {
+      Matrix<Scalar,2,3> m; m << v0.transpose(), v1.transpose();
+      JacobiSVD<Matrix<Scalar,2,3> > svd(m, ComputeFullV);
+      result.normal() = svd.matrixV().col(2);
+    }
+    else
+      result.normal() /= norm;
     result.offset() = -p0.dot(result.normal());
     return result;
   }
diff --git a/Eigen/src/Geometry/Rotation2D.h b/Eigen/src/Geometry/Rotation2D.h
index 1cac343a5..a2d59fce1 100644
--- a/Eigen/src/Geometry/Rotation2D.h
+++ b/Eigen/src/Geometry/Rotation2D.h
@@ -60,6 +60,9 @@ public:
 
   /** Construct a 2D counter clock wise rotation from the angle \a a in radian. */
   inline Rotation2D(const Scalar& a) : m_angle(a) {}
+  
+  /** Default constructor wihtout initialization. The represented rotation is undefined. */
+  Rotation2D() {}
 
   /** \returns the rotation angle */
   inline Scalar angle() const { return m_angle; }
@@ -81,10 +84,10 @@ public:
   /** Applies the rotation to a 2D vector */
   Vector2 operator* (const Vector2& vec) const
   { return toRotationMatrix() * vec; }
-
+  
   template<typename Derived>
   Rotation2D& fromRotationMatrix(const MatrixBase<Derived>& m);
-  Matrix2 toRotationMatrix(void) const;
+  Matrix2 toRotationMatrix() const;
 
   /** \returns the spherical interpolation between \c *this and \a other using
     * parameter \a t. It is in fact equivalent to a linear interpolation.
diff --git a/Eigen/src/Geometry/Transform.h b/Eigen/src/Geometry/Transform.h
index 56f361566..e786e5356 100644
--- a/Eigen/src/Geometry/Transform.h
+++ b/Eigen/src/Geometry/Transform.h
@@ -62,6 +62,8 @@ struct transform_construct_from_matrix;
 
 template<typename TransformType> struct transform_take_affine_part;
 
+template<int Mode> struct transform_make_affine;
+
 } // end namespace internal
 
 /** \geometry_module \ingroup Geometry_Module
@@ -230,8 +232,7 @@ public:
   inline Transform()
   {
     check_template_params();
-    if (int(Mode)==Affine)
-      makeAffine();
+    internal::transform_make_affine<(int(Mode)==Affine) ? Affine : AffineCompact>::run(m_matrix);
   }
 
   inline Transform(const Transform& other)
@@ -591,11 +592,7 @@ public:
     */
   void makeAffine()
   {
-    if(int(Mode)!=int(AffineCompact))
-    {
-      matrix().template block<1,Dim>(Dim,0).setZero();
-      matrix().coeffRef(Dim,Dim) = Scalar(1);
-    }
+    internal::transform_make_affine<int(Mode)>::run(m_matrix);
   }
 
   /** \internal
@@ -1079,6 +1076,24 @@ Transform<Scalar,Dim,Mode,Options>::fromPositionOrientationScale(const MatrixBas
 
 namespace internal {
 
+template<int Mode>
+struct transform_make_affine
+{
+  template<typename MatrixType>
+  static void run(MatrixType &mat)
+  {
+    static const int Dim = MatrixType::ColsAtCompileTime-1;
+    mat.template block<1,Dim>(Dim,0).setZero();
+    mat.coeffRef(Dim,Dim) = typename MatrixType::Scalar(1);
+  }
+};
+
+template<>
+struct transform_make_affine<AffineCompact>
+{
+  template<typename MatrixType> static void run(MatrixType &) { }
+};
+    
 // selector needed to avoid taking the inverse of a 3x4 matrix
 template<typename TransformType, int Mode=TransformType::Mode>
 struct projective_transform_inverse
diff --git a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h
index 2b9fb7f88..dd135c21f 100644
--- a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h
+++ b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h
@@ -39,7 +39,6 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
   int maxIters = iters;
 
   int n = mat.cols();
-  x = precond.solve(x);
   VectorType r  = rhs - mat * x;
   VectorType r0 = r;
   
@@ -143,7 +142,7 @@ struct traits<BiCGSTAB<_MatrixType,_Preconditioner> >
   * SparseMatrix<double> A(n,n);
   * // fill A and b
   * BiCGSTAB<SparseMatrix<double> > solver;
-  * solver(A);
+  * solver.compute(A);
   * x = solver.solve(b);
   * std::cout << "#iterations:     " << solver.iterations() << std::endl;
   * std::cout << "estimated error: " << solver.error()      << std::endl;
diff --git a/Eigen/src/PardisoSupport/PardisoSupport.h b/Eigen/src/PardisoSupport/PardisoSupport.h
index 1c48f0df7..18cd7d88a 100644
--- a/Eigen/src/PardisoSupport/PardisoSupport.h
+++ b/Eigen/src/PardisoSupport/PardisoSupport.h
@@ -219,7 +219,7 @@ class PardisoImpl
     void pardisoInit(int type)
     {
       m_type = type;
-      bool symmetric = abs(m_type) < 10;
+      bool symmetric = std::abs(m_type) < 10;
       m_iparm[0] = 1;   // No solver default
       m_iparm[1] = 3;   // use Metis for the ordering
       m_iparm[2] = 1;   // Numbers of processors, value of OMP_NUM_THREADS
diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h
index dff9e44eb..c57f2974c 100644
--- a/Eigen/src/SVD/JacobiSVD.h
+++ b/Eigen/src/SVD/JacobiSVD.h
@@ -762,6 +762,7 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
 
     internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows> m_qr_precond_morecols;
     internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> m_qr_precond_morerows;
+    MatrixType m_scaledMatrix;
 };
 
 template<typename MatrixType, int QRPreconditioner>
@@ -808,8 +809,9 @@ void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols, u
                             : 0);
   m_workMatrix.resize(m_diagSize, m_diagSize);
   
-  if(m_cols>m_rows) m_qr_precond_morecols.allocate(*this);
-  if(m_rows>m_cols) m_qr_precond_morerows.allocate(*this);
+  if(m_cols>m_rows)   m_qr_precond_morecols.allocate(*this);
+  if(m_rows>m_cols)   m_qr_precond_morerows.allocate(*this);
+  if(m_cols!=m_cols)  m_scaledMatrix.resize(rows,cols);
 }
 
 template<typename MatrixType, int QRPreconditioner>
@@ -826,21 +828,26 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
   // limit for very small denormal numbers to be considered zero in order to avoid infinite loops (see bug 286)
   const RealScalar considerAsZero = RealScalar(2) * std::numeric_limits<RealScalar>::denorm_min();
 
+  // Scaling factor to reduce over/under-flows
+  RealScalar scale = matrix.cwiseAbs().maxCoeff();
+  if(scale==RealScalar(0)) scale = RealScalar(1);
+  
   /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */
 
-  if(!m_qr_precond_morecols.run(*this, matrix) && !m_qr_precond_morerows.run(*this, matrix))
+  if(m_rows!=m_cols)
   {
-    m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize);
+    m_scaledMatrix = matrix / scale;
+    m_qr_precond_morecols.run(*this, m_scaledMatrix);
+    m_qr_precond_morerows.run(*this, m_scaledMatrix);
+  }
+  else
+  {
+    m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize) / scale;
     if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows);
     if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize);
     if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols);
     if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize);
   }
-  
-  // Scaling factor to reduce over/under-flows
-  RealScalar scale = m_workMatrix.cwiseAbs().maxCoeff();
-  if(scale==RealScalar(0)) scale = RealScalar(1);
-  m_workMatrix /= scale;
 
   /*** step 2. The main Jacobi SVD iteration. ***/
 
@@ -861,7 +868,8 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
         using std::max;
         RealScalar threshold = (max)(considerAsZero, precision * (max)(abs(m_workMatrix.coeff(p,p)),
                                                                        abs(m_workMatrix.coeff(q,q))));
-        if((max)(abs(m_workMatrix.coeff(p,q)),abs(m_workMatrix.coeff(q,p))) > threshold)
+        // We compare both values to threshold instead of calling max to be robust to NaN (See bug 791)
+        if(abs(m_workMatrix.coeff(p,q))>threshold || abs(m_workMatrix.coeff(q,p)) > threshold)
         {
           finished = false;
 
diff --git a/Eigen/src/SparseCore/AmbiVector.h b/Eigen/src/SparseCore/AmbiVector.h
index 17fff96a7..220c6451c 100644
--- a/Eigen/src/SparseCore/AmbiVector.h
+++ b/Eigen/src/SparseCore/AmbiVector.h
@@ -69,7 +69,7 @@ class AmbiVector
       delete[] m_buffer;
       if (size<1000)
       {
-        Index allocSize = (size * sizeof(ListEl))/sizeof(Scalar);
+        Index allocSize = (size * sizeof(ListEl) + sizeof(Scalar) - 1)/sizeof(Scalar);
         m_allocatedElements = (allocSize*sizeof(Scalar))/sizeof(ListEl);
         m_buffer = new Scalar[allocSize];
       }
@@ -88,7 +88,7 @@ class AmbiVector
       Index copyElements = m_allocatedElements;
       m_allocatedElements = (std::min)(Index(m_allocatedElements*1.5),m_size);
       Index allocSize = m_allocatedElements * sizeof(ListEl);
-      allocSize = allocSize/sizeof(Scalar) + (allocSize%sizeof(Scalar)>0?1:0);
+      allocSize = (allocSize + sizeof(Scalar) - 1)/sizeof(Scalar);
       Scalar* newBuffer = new Scalar[allocSize];
       memcpy(newBuffer,  m_buffer,  copyElements * sizeof(ListEl));
       delete[] m_buffer;
diff --git a/Eigen/src/SparseCore/SparseBlock.h b/Eigen/src/SparseCore/SparseBlock.h
index 16a20a574..0ede034ba 100644
--- a/Eigen/src/SparseCore/SparseBlock.h
+++ b/Eigen/src/SparseCore/SparseBlock.h
@@ -68,6 +68,8 @@ public:
     const internal::variable_if_dynamic<Index, OuterSize> m_outerSize;
 
     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(BlockImpl)
+  private:
+    Index nonZeros() const;
 };
 
 
@@ -82,6 +84,7 @@ class BlockImpl<SparseMatrix<_Scalar, _Options, _Index>,BlockRows,BlockCols,true
     typedef SparseMatrix<_Scalar, _Options, _Index> SparseMatrixType;
     typedef typename internal::remove_all<typename SparseMatrixType::Nested>::type _MatrixTypeNested;
     typedef Block<SparseMatrixType, BlockRows, BlockCols, true> BlockType;
+    typedef Block<const SparseMatrixType, BlockRows, BlockCols, true> ConstBlockType;
 public:
     enum { IsRowMajor = internal::traits<BlockType>::IsRowMajor };
     EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType)
@@ -245,6 +248,93 @@ public:
 
 };
 
+
+template<typename _Scalar, int _Options, typename _Index, int BlockRows, int BlockCols>
+class BlockImpl<const SparseMatrix<_Scalar, _Options, _Index>,BlockRows,BlockCols,true,Sparse>
+  : public SparseMatrixBase<Block<const SparseMatrix<_Scalar, _Options, _Index>,BlockRows,BlockCols,true> >
+{
+    typedef SparseMatrix<_Scalar, _Options, _Index> SparseMatrixType;
+    typedef typename internal::remove_all<typename SparseMatrixType::Nested>::type _MatrixTypeNested;
+    typedef Block<const SparseMatrixType, BlockRows, BlockCols, true> BlockType;
+public:
+    enum { IsRowMajor = internal::traits<BlockType>::IsRowMajor };
+    EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType)
+protected:
+    enum { OuterSize = IsRowMajor ? BlockRows : BlockCols };
+public:
+    
+    class InnerIterator: public SparseMatrixType::InnerIterator
+    {
+      public:
+        inline InnerIterator(const BlockType& xpr, Index outer)
+          : SparseMatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer)
+        {}
+        inline Index row() const { return IsRowMajor ? m_outer : this->index(); }
+        inline Index col() const { return IsRowMajor ? this->index() : m_outer; }
+      protected:
+        Index m_outer;
+    };
+    class ReverseInnerIterator: public SparseMatrixType::ReverseInnerIterator
+    {
+      public:
+        inline ReverseInnerIterator(const BlockType& xpr, Index outer)
+          : SparseMatrixType::ReverseInnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer)
+        {}
+        inline Index row() const { return IsRowMajor ? m_outer : this->index(); }
+        inline Index col() const { return IsRowMajor ? this->index() : m_outer; }
+      protected:
+        Index m_outer;
+    };
+
+    inline BlockImpl(const SparseMatrixType& xpr, int i)
+      : m_matrix(xpr), m_outerStart(i), m_outerSize(OuterSize)
+    {}
+
+    inline BlockImpl(const SparseMatrixType& xpr, int startRow, int startCol, int blockRows, int blockCols)
+      : m_matrix(xpr), m_outerStart(IsRowMajor ? startRow : startCol), m_outerSize(IsRowMajor ? blockRows : blockCols)
+    {}
+
+    inline const Scalar* valuePtr() const
+    { return m_matrix.valuePtr() + m_matrix.outerIndexPtr()[m_outerStart]; }
+
+    inline const Index* innerIndexPtr() const
+    { return m_matrix.innerIndexPtr() + m_matrix.outerIndexPtr()[m_outerStart]; }
+
+    inline const Index* outerIndexPtr() const
+    { return m_matrix.outerIndexPtr() + m_outerStart; }
+
+    Index nonZeros() const
+    {
+      if(m_matrix.isCompressed())
+        return  std::size_t(m_matrix.outerIndexPtr()[m_outerStart+m_outerSize.value()])
+              - std::size_t(m_matrix.outerIndexPtr()[m_outerStart]);
+      else if(m_outerSize.value()==0)
+        return 0;
+      else
+        return Map<const Matrix<Index,OuterSize,1> >(m_matrix.innerNonZeroPtr()+m_outerStart, m_outerSize.value()).sum();
+    }
+
+    const Scalar& lastCoeff() const
+    {
+      EIGEN_STATIC_ASSERT_VECTOR_ONLY(BlockImpl);
+      eigen_assert(nonZeros()>0);
+      if(m_matrix.isCompressed())
+        return m_matrix.valuePtr()[m_matrix.outerIndexPtr()[m_outerStart+1]-1];
+      else
+        return m_matrix.valuePtr()[m_matrix.outerIndexPtr()[m_outerStart]+m_matrix.innerNonZeroPtr()[m_outerStart]-1];
+    }
+
+    EIGEN_STRONG_INLINE Index rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); }
+    EIGEN_STRONG_INLINE Index cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); }
+
+  protected:
+
+    typename SparseMatrixType::Nested m_matrix;
+    Index m_outerStart;
+    const internal::variable_if_dynamic<Index, OuterSize> m_outerSize;
+
+};
+
 //----------
 
 /** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this
diff --git a/Eigen/src/SparseCore/SparseDenseProduct.h b/Eigen/src/SparseCore/SparseDenseProduct.h
index 78411db98..a9084192e 100644
--- a/Eigen/src/SparseCore/SparseDenseProduct.h
+++ b/Eigen/src/SparseCore/SparseDenseProduct.h
@@ -306,15 +306,6 @@ class DenseTimeSparseProduct
     DenseTimeSparseProduct& operator=(const DenseTimeSparseProduct&);
 };
 
-// sparse * dense
-template<typename Derived>
-template<typename OtherDerived>
-inline const typename SparseDenseProductReturnType<Derived,OtherDerived>::Type
-SparseMatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const
-{
-  return typename SparseDenseProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
-}
-
 } // end namespace Eigen
 
 #endif // EIGEN_SPARSEDENSEPRODUCT_H
diff --git a/Eigen/src/SparseCore/SparseMatrixBase.h b/Eigen/src/SparseCore/SparseMatrixBase.h
index bbcf7fb1c..485e85bec 100644
--- a/Eigen/src/SparseCore/SparseMatrixBase.h
+++ b/Eigen/src/SparseCore/SparseMatrixBase.h
@@ -358,7 +358,8 @@ template<typename Derived> class SparseMatrixBase : public EigenBase<Derived>
     /** sparse * dense (returns a dense object unless it is an outer product) */
     template<typename OtherDerived>
     const typename SparseDenseProductReturnType<Derived,OtherDerived>::Type
-    operator*(const MatrixBase<OtherDerived> &other) const;
+    operator*(const MatrixBase<OtherDerived> &other) const
+    { return typename SparseDenseProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived()); }
     
      /** \returns an expression of P H P^-1 where H is the matrix represented by \c *this */
     SparseSymmetricPermutationProduct<Derived,Upper|Lower> twistedBy(const PermutationMatrix<Dynamic,Dynamic,Index>& perm) const
diff --git a/Eigen/src/SparseCore/SparsePermutation.h b/Eigen/src/SparseCore/SparsePermutation.h
index b85be93f6..75e210009 100644
--- a/Eigen/src/SparseCore/SparsePermutation.h
+++ b/Eigen/src/SparseCore/SparsePermutation.h
@@ -61,7 +61,7 @@ struct permut_sparsematrix_product_retval
         for(Index j=0; j<m_matrix.outerSize(); ++j)
         {
           Index jp = m_permutation.indices().coeff(j);
-          sizes[((Side==OnTheLeft) ^ Transposed) ? jp : j] = m_matrix.innerVector(((Side==OnTheRight) ^ Transposed) ? jp : j).size();
+          sizes[((Side==OnTheLeft) ^ Transposed) ? jp : j] = m_matrix.innerVector(((Side==OnTheRight) ^ Transposed) ? jp : j).nonZeros();
         }
         tmp.reserve(sizes);
         for(Index j=0; j<m_matrix.outerSize(); ++j)
diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h
index 1d592f2c8..d6f5d59ca 100644
--- a/Eigen/src/SparseLU/SparseLU.h
+++ b/Eigen/src/SparseLU/SparseLU.h
@@ -260,14 +260,13 @@ class SparseLU : public internal::SparseLUImpl<typename _MatrixType::Scalar, typ
       eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
       // Initialize with the determinant of the row matrix
       Scalar det = Scalar(1.);
-      //Note that the diagonal blocks of U are stored in supernodes,
+      // Note that the diagonal blocks of U are stored in supernodes,
       // which are available in the  L part :)
       for (Index j = 0; j < this->cols(); ++j)
       {
         for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
         {
-          if(it.row() < j) continue;
-          if(it.row() == j)
+          if(it.index() == j)
           {
             det *= (std::abs)(it.value());
             break;
diff --git a/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h b/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h
index ad6f2183f..b16afd6a4 100644
--- a/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h
+++ b/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h
@@ -189,8 +189,8 @@ class MappedSuperNodalMatrix<Scalar,Index>::InnerIterator
         m_idval(mat.colIndexPtr()[outer]),
         m_startidval(m_idval),
         m_endidval(mat.colIndexPtr()[outer+1]),
-        m_idrow(mat.rowIndexPtr()[outer]),
-        m_endidrow(mat.rowIndexPtr()[outer+1])
+        m_idrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[outer]]]),
+        m_endidrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[outer]]+1])
     {}
     inline InnerIterator& operator++()
     { 
diff --git a/Eigen/src/SparseQR/SparseQR.h b/Eigen/src/SparseQR/SparseQR.h
index 4c6553bf2..a00bd5db1 100644
--- a/Eigen/src/SparseQR/SparseQR.h
+++ b/Eigen/src/SparseQR/SparseQR.h
@@ -75,7 +75,7 @@ class SparseQR
     typedef Matrix<Scalar, Dynamic, 1> ScalarVector;
     typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
   public:
-    SparseQR () : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false)
+    SparseQR () : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
     { }
     
     /** Construct a QR factorization of the matrix \a mat.
@@ -84,7 +84,7 @@ class SparseQR
       * 
       * \sa compute()
       */
-    SparseQR(const MatrixType& mat) : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false)
+    SparseQR(const MatrixType& mat) : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
     {
       compute(mat);
     }
@@ -262,6 +262,7 @@ class SparseQR
     IndexVector m_etree;            // Column elimination tree
     IndexVector m_firstRowElt;      // First element in each row
     bool m_isQSorted;               // whether Q is sorted or not
+    bool m_isEtreeOk;               // whether the elimination tree match the initial input matrix
     
     template <typename, typename > friend struct SparseQR_QProduct;
     template <typename > friend struct SparseQRMatrixQReturnType;
@@ -281,9 +282,11 @@ template <typename MatrixType, typename OrderingType>
 void SparseQR<MatrixType,OrderingType>::analyzePattern(const MatrixType& mat)
 {
   eigen_assert(mat.isCompressed() && "SparseQR requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to SparseQR");
+  // Copy to a column major matrix if the input is rowmajor
+  typename internal::conditional<MatrixType::IsRowMajor,QRMatrixType,const MatrixType&>::type matCpy(mat);
   // Compute the column fill reducing ordering
   OrderingType ord; 
-  ord(mat, m_perm_c); 
+  ord(matCpy, m_perm_c); 
   Index n = mat.cols();
   Index m = mat.rows();
   Index diagSize = (std::min)(m,n);
@@ -296,7 +299,8 @@ void SparseQR<MatrixType,OrderingType>::analyzePattern(const MatrixType& mat)
   
   // Compute the column elimination tree of the permuted matrix
   m_outputPerm_c = m_perm_c.inverse();
-  internal::coletree(mat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
+  internal::coletree(matCpy, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
+  m_isEtreeOk = true;
   
   m_R.resize(m, n);
   m_Q.resize(m, diagSize);
@@ -330,15 +334,38 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
   Index nzcolR, nzcolQ;                                     // Number of nonzero for the current column of R and Q
   ScalarVector tval(m);                                     // The dense vector used to compute the current column
   RealScalar pivotThreshold = m_threshold;
-    
+  
+  m_R.setZero();
+  m_Q.setZero();
   m_pmat = mat;
+  if(!m_isEtreeOk)
+  {
+    m_outputPerm_c = m_perm_c.inverse();
+    internal::coletree(m_pmat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
+    m_isEtreeOk = true;
+  }
+
   m_pmat.uncompress(); // To have the innerNonZeroPtr allocated
+  
   // Apply the fill-in reducing permutation lazily:
-  for (int i = 0; i < n; i++)
   {
-    Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i;
-    m_pmat.outerIndexPtr()[p] = mat.outerIndexPtr()[i]; 
-    m_pmat.innerNonZeroPtr()[p] = mat.outerIndexPtr()[i+1] - mat.outerIndexPtr()[i]; 
+    // If the input is row major, copy the original column indices,
+    // otherwise directly use the input matrix
+    // 
+    IndexVector originalOuterIndicesCpy;
+    const Index *originalOuterIndices = mat.outerIndexPtr();
+    if(MatrixType::IsRowMajor)
+    {
+      originalOuterIndicesCpy = IndexVector::Map(m_pmat.outerIndexPtr(),n+1);
+      originalOuterIndices = originalOuterIndicesCpy.data();
+    }
+    
+    for (int i = 0; i < n; i++)
+    {
+      Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i;
+      m_pmat.outerIndexPtr()[p] = originalOuterIndices[i]; 
+      m_pmat.innerNonZeroPtr()[p] = originalOuterIndices[i+1] - originalOuterIndices[i]; 
+    }
   }
   
   /* Compute the default threshold as in MatLab, see:
@@ -349,6 +376,8 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
   {
     RealScalar max2Norm = 0.0;
     for (int j = 0; j < n; j++) max2Norm = (max)(max2Norm, m_pmat.col(j).norm());
+    if(max2Norm==RealScalar(0))
+      max2Norm = RealScalar(1);
     pivotThreshold = 20 * (m + n) * max2Norm * NumTraits<RealScalar>::epsilon();
   }
   
@@ -357,7 +386,7 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
   
   Index nonzeroCol = 0; // Record the number of valid pivots
   m_Q.startVec(0);
-  
+
   // Left looking rank-revealing QR factorization: compute a column of R and Q at a time
   for (Index col = 0; col < n; ++col)
   {
@@ -373,7 +402,7 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
     // all the nodes (with indexes lower than rank) reachable through the column elimination tree (etree) rooted at node k.
     // Note: if the diagonal entry does not exist, then its contribution must be explicitly added,
     // thus the trick with found_diag that permits to do one more iteration on the diagonal element if this one has not been found.
-    for (typename MatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp)
+    for (typename QRMatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp)
     {
       Index curIdx = nonzeroCol;
       if(itp) curIdx = itp.row();
@@ -447,7 +476,7 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
       }
     } // End update current column
     
-    Scalar tau;
+    Scalar tau = 0;
     RealScalar beta = 0;
     
     if(nonzeroCol < diagSize)
@@ -461,7 +490,6 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
       for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq)));
       if(sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0))
       {
-        tau = RealScalar(0);
         beta = numext::real(c0);
         tval(Qidx(0)) = 1;
       }
@@ -514,6 +542,7 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
       
       // Recompute the column elimination tree
       internal::coletree(m_pmat, m_etree, m_firstRowElt, m_pivotperm.indices().data());
+      m_isEtreeOk = false;
     }
   }
   
@@ -525,13 +554,13 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
   m_R.finalize();
   m_R.makeCompressed();
   m_isQSorted = false;
-  
+
   m_nonzeropivots = nonzeroCol;
   
   if(nonzeroCol<n)
   {
     // Permute the triangular factor to put the 'dead' columns to the end
-    MatrixType tempR(m_R);
+    QRMatrixType tempR(m_R);
     m_R = tempR * m_pivotperm;
     
     // Update the column permutation
diff --git a/Eigen/src/UmfPackSupport/UmfPackSupport.h b/Eigen/src/UmfPackSupport/UmfPackSupport.h
index 3a48cecf7..29c60c378 100644
--- a/Eigen/src/UmfPackSupport/UmfPackSupport.h
+++ b/Eigen/src/UmfPackSupport/UmfPackSupport.h
@@ -107,6 +107,16 @@ inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *N
   return umfpack_zi_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info);
 }
 
+namespace internal {
+  template<typename T> struct umfpack_helper_is_sparse_plain : false_type {};
+  template<typename Scalar, int Options, typename StorageIndex>
+  struct umfpack_helper_is_sparse_plain<SparseMatrix<Scalar,Options,StorageIndex> >
+    : true_type {};
+  template<typename Scalar, int Options, typename StorageIndex>
+  struct umfpack_helper_is_sparse_plain<MappedSparseMatrix<Scalar,Options,StorageIndex> >
+    : true_type {};
+}
+
 /** \ingroup UmfPackSupport_Module
   * \brief A sparse LU factorization and solver based on UmfPack
   *
@@ -192,10 +202,14 @@ class UmfPackLU : internal::noncopyable
      *  Note that the matrix should be column-major, and in compressed format for best performance.
      *  \sa SparseMatrix::makeCompressed().
      */
-    void compute(const MatrixType& matrix)
+    template<typename InputMatrixType>
+    void compute(const InputMatrixType& matrix)
     {
-      analyzePattern(matrix);
-      factorize(matrix);
+      if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar());
+      if(m_numeric)  umfpack_free_numeric(&m_numeric,Scalar());
+      grapInput(matrix.derived());
+      analyzePattern_impl();
+      factorize_impl();
     }
 
     /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
@@ -230,23 +244,15 @@ class UmfPackLU : internal::noncopyable
       *
       * \sa factorize(), compute()
       */
-    void analyzePattern(const MatrixType& matrix)
+    template<typename InputMatrixType>
+    void analyzePattern(const InputMatrixType& matrix)
     {
-      if(m_symbolic)
-        umfpack_free_symbolic(&m_symbolic,Scalar());
-      if(m_numeric)
-        umfpack_free_numeric(&m_numeric,Scalar());
+      if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar());
+      if(m_numeric)  umfpack_free_numeric(&m_numeric,Scalar());
       
-      grapInput(matrix);
-
-      int errorCode = 0;
-      errorCode = umfpack_symbolic(matrix.rows(), matrix.cols(), m_outerIndexPtr, m_innerIndexPtr, m_valuePtr,
-                                   &m_symbolic, 0, 0);
+      grapInput(matrix.derived());
 
-      m_isInitialized = true;
-      m_info = errorCode ? InvalidInput : Success;
-      m_analysisIsOk = true;
-      m_factorizationIsOk = false;
+      analyzePattern_impl();
     }
 
     /** Performs a numeric decomposition of \a matrix
@@ -255,20 +261,16 @@ class UmfPackLU : internal::noncopyable
       *
       * \sa analyzePattern(), compute()
       */
-    void factorize(const MatrixType& matrix)
+    template<typename InputMatrixType>
+    void factorize(const InputMatrixType& matrix)
     {
       eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
       if(m_numeric)
         umfpack_free_numeric(&m_numeric,Scalar());
 
-      grapInput(matrix);
-
-      int errorCode;
-      errorCode = umfpack_numeric(m_outerIndexPtr, m_innerIndexPtr, m_valuePtr,
-                                  m_symbolic, &m_numeric, 0, 0);
-
-      m_info = errorCode ? NumericalIssue : Success;
-      m_factorizationIsOk = true;
+      grapInput(matrix.derived());
+      
+      factorize_impl();
     }
 
     #ifndef EIGEN_PARSED_BY_DOXYGEN
@@ -283,19 +285,20 @@ class UmfPackLU : internal::noncopyable
 
   protected:
 
-
     void init()
     {
-      m_info = InvalidInput;
-      m_isInitialized = false;
-      m_numeric = 0;
-      m_symbolic = 0;
-      m_outerIndexPtr = 0;
-      m_innerIndexPtr = 0;
-      m_valuePtr      = 0;
+      m_info                  = InvalidInput;
+      m_isInitialized         = false;
+      m_numeric               = 0;
+      m_symbolic              = 0;
+      m_outerIndexPtr         = 0;
+      m_innerIndexPtr         = 0;
+      m_valuePtr              = 0;
+      m_extractedDataAreDirty = true;
     }
     
-    void grapInput(const MatrixType& mat)
+    template<typename InputMatrixType>
+    void grapInput_impl(const InputMatrixType& mat, internal::true_type)
     {
       m_copyMatrix.resize(mat.rows(), mat.cols());
       if( ((MatrixType::Flags&RowMajorBit)==RowMajorBit) || sizeof(typename MatrixType::Index)!=sizeof(int) || !mat.isCompressed() )
@@ -313,6 +316,45 @@ class UmfPackLU : internal::noncopyable
         m_valuePtr      = mat.valuePtr();
       }
     }
+    
+    template<typename InputMatrixType>
+    void grapInput_impl(const InputMatrixType& mat, internal::false_type)
+    {
+      m_copyMatrix = mat;
+      m_outerIndexPtr = m_copyMatrix.outerIndexPtr();
+      m_innerIndexPtr = m_copyMatrix.innerIndexPtr();
+      m_valuePtr      = m_copyMatrix.valuePtr();
+    }
+    
+    template<typename InputMatrixType>
+    void grapInput(const InputMatrixType& mat)
+    {
+      grapInput_impl(mat, internal::umfpack_helper_is_sparse_plain<InputMatrixType>());
+    }
+    
+    void analyzePattern_impl()
+    {
+      int errorCode = 0;
+      errorCode = umfpack_symbolic(m_copyMatrix.rows(), m_copyMatrix.cols(), m_outerIndexPtr, m_innerIndexPtr, m_valuePtr,
+                                   &m_symbolic, 0, 0);
+
+      m_isInitialized = true;
+      m_info = errorCode ? InvalidInput : Success;
+      m_analysisIsOk = true;
+      m_factorizationIsOk = false;
+      m_extractedDataAreDirty = true;
+    }
+    
+    void factorize_impl()
+    {
+      int errorCode;
+      errorCode = umfpack_numeric(m_outerIndexPtr, m_innerIndexPtr, m_valuePtr,
+                                  m_symbolic, &m_numeric, 0, 0);
+
+      m_info = errorCode ? NumericalIssue : Success;
+      m_factorizationIsOk = true;
+      m_extractedDataAreDirty = true;
+    }
 
     // cached data to reduce reallocation, etc.
     mutable LUMatrixType m_l;
diff --git a/cmake/EigenTesting.cmake b/cmake/EigenTesting.cmake
index d9e22ab1a..80b2841df 100644
--- a/cmake/EigenTesting.cmake
+++ b/cmake/EigenTesting.cmake
@@ -452,20 +452,12 @@ macro(ei_set_build_string)
 endmacro(ei_set_build_string)
 
 macro(ei_is_64bit_env VAR)
-
-  file(WRITE "${CMAKE_CURRENT_BINARY_DIR}/is64.cpp"
-      "int main() { return (sizeof(int*) == 8 ? 1 : 0); }
-      ")
-  try_run(run_res compile_res
-         ${CMAKE_CURRENT_BINARY_DIR} "${CMAKE_CURRENT_BINARY_DIR}/is64.cpp"
-          RUN_OUTPUT_VARIABLE run_output)
-
-  if(compile_res AND run_res)
-    set(${VAR} ${run_res})
-  elseif(CMAKE_CL_64)
-    set(${VAR} 1)
-  elseif("$ENV{Platform}" STREQUAL "X64") # nmake 64 bit
+  if(CMAKE_SIZEOF_VOID_P EQUAL 8)
     set(${VAR} 1)
+  elseif(CMAKE_SIZEOF_VOID_P EQUAL 4)
+    set(${VAR} 0)
+  else()
+    message(WARNING "Unsupported pointer size. Please contact the authors.")
   endif()
 endmacro(ei_is_64bit_env)
 
diff --git a/cmake/FindCholmod.cmake b/cmake/FindCholmod.cmake
index 7b3046d45..23239c300 100644
--- a/cmake/FindCholmod.cmake
+++ b/cmake/FindCholmod.cmake
@@ -86,4 +86,4 @@ include(FindPackageHandleStandardArgs)
 find_package_handle_standard_args(CHOLMOD DEFAULT_MSG
                                   CHOLMOD_INCLUDES CHOLMOD_LIBRARIES)
 
-mark_as_advanced(CHOLMOD_INCLUDES CHOLMOD_LIBRARIES AMD_LIBRARY COLAMD_LIBRARY SUITESPARSE_LIBRARY)
+mark_as_advanced(CHOLMOD_INCLUDES CHOLMOD_LIBRARIES AMD_LIBRARY COLAMD_LIBRARY SUITESPARSE_LIBRARY CAMD_LIBRARY CCOLAMD_LIBRARY CHOLMOD_METIS_LIBRARY)
diff --git a/cmake/FindFFTW.cmake b/cmake/FindFFTW.cmake
index a9e9925e7..6c4dc9ab4 100644
--- a/cmake/FindFFTW.cmake
+++ b/cmake/FindFFTW.cmake
@@ -115,5 +115,5 @@ include(FindPackageHandleStandardArgs)
 find_package_handle_standard_args(FFTW DEFAULT_MSG
                                   FFTW_INCLUDES FFTW_LIBRARIES)
 
-mark_as_advanced(FFTW_INCLUDES FFTW_LIBRARIES)
+mark_as_advanced(FFTW_INCLUDES FFTW_LIBRARIES FFTW_LIB FFTWF_LIB FFTWL_LIB)
 
diff --git a/cmake/FindMetis.cmake b/cmake/FindMetis.cmake
index 627c3e9ae..e0040d320 100644
--- a/cmake/FindMetis.cmake
+++ b/cmake/FindMetis.cmake
@@ -10,16 +10,50 @@ find_path(METIS_INCLUDES
   PATHS 
   $ENV{METISDIR} 
   ${INCLUDE_INSTALL_DIR} 
-  PATH_SUFFIXES 
+  PATH_SUFFIXES
+  .
   metis
   include
 )
 
+macro(_metis_check_version)
+  file(READ "${METIS_INCLUDES}/metis.h" _metis_version_header)
+
+  string(REGEX MATCH "define[ \t]+METIS_VER_MAJOR[ \t]+([0-9]+)" _metis_major_version_match "${_metis_version_header}")
+  set(METIS_MAJOR_VERSION "${CMAKE_MATCH_1}")
+  string(REGEX MATCH "define[ \t]+METIS_VER_MINOR[ \t]+([0-9]+)" _metis_minor_version_match "${_metis_version_header}")
+  set(METIS_MINOR_VERSION "${CMAKE_MATCH_1}")
+  string(REGEX MATCH "define[ \t]+METIS_VER_SUBMINOR[ \t]+([0-9]+)" _metis_subminor_version_match "${_metis_version_header}")
+  set(METIS_SUBMINOR_VERSION "${CMAKE_MATCH_1}")
+  if(NOT METIS_MAJOR_VERSION)
+    message(WARNING "Could not determine Metis version. Assuming version 4.0.0")
+    set(METIS_VERSION 4.0.0)
+  else()
+    set(METIS_VERSION ${METIS_MAJOR_VERSION}.${METIS_MINOR_VERSION}.${METIS_SUBMINOR_VERSION})
+  endif()
+  if(${METIS_VERSION} VERSION_LESS ${Metis_FIND_VERSION})
+    set(METIS_VERSION_OK FALSE)
+  else()
+    set(METIS_VERSION_OK TRUE)
+  endif()
+
+  if(NOT METIS_VERSION_OK)
+    message(STATUS "Metis version ${METIS_VERSION} found in ${METIS_INCLUDES}, "
+                   "but at least version ${Metis_FIND_VERSION} is required")
+  endif(NOT METIS_VERSION_OK)
+endmacro(_metis_check_version)
+
+  if(METIS_INCLUDES AND Metis_FIND_VERSION)
+    _metis_check_version()
+  else()
+    set(METIS_VERSION_OK TRUE)
+  endif()
+
 
 find_library(METIS_LIBRARIES metis PATHS $ENV{METISDIR} ${LIB_INSTALL_DIR} PATH_SUFFIXES lib)
 
 include(FindPackageHandleStandardArgs)
 find_package_handle_standard_args(METIS DEFAULT_MSG
-                                  METIS_INCLUDES METIS_LIBRARIES)
+                                  METIS_INCLUDES METIS_LIBRARIES METIS_VERSION_OK)
 
 mark_as_advanced(METIS_INCLUDES METIS_LIBRARIES)
diff --git a/cmake/language_support.cmake b/cmake/language_support.cmake
index d687b71f6..231f7ff70 100644
--- a/cmake/language_support.cmake
+++ b/cmake/language_support.cmake
@@ -33,7 +33,7 @@ function(workaround_9220 language language_works)
   file(WRITE ${CMAKE_BINARY_DIR}/language_tests/${language}/CMakeLists.txt
     ${text})
   execute_process(
-    COMMAND ${CMAKE_COMMAND} .
+    COMMAND ${CMAKE_COMMAND} . -G "${CMAKE_GENERATOR}"
     WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/language_tests/${language}
     RESULT_VARIABLE return_code
     OUTPUT_QUIET
diff --git a/doc/AsciiQuickReference.txt b/doc/AsciiQuickReference.txt
index 1e74e0528..b9f497f87 100644
--- a/doc/AsciiQuickReference.txt
+++ b/doc/AsciiQuickReference.txt
@@ -67,10 +67,10 @@ P.rightCols<cols>()                // P(:, end-cols+1:end)
 P.rightCols(cols)                  // P(:, end-cols+1:end)
 P.topRows<rows>()                  // P(1:rows, :)
 P.topRows(rows)                    // P(1:rows, :)
-P.middleRows<rows>(i)              // P(:, i+1:i+rows)
-P.middleRows(i, rows)              // P(:, i+1:i+rows)
-P.bottomRows<rows>()               // P(:, end-rows+1:end)
-P.bottomRows(rows)                 // P(:, end-rows+1:end)
+P.middleRows<rows>(i)              // P(i+1:i+rows, :)
+P.middleRows(i, rows)              // P(i+1:i+rows, :)
+P.bottomRows<rows>()               // P(end-rows+1:end, :)
+P.bottomRows(rows)                 // P(end-rows+1:end, :)
 P.topLeftCorner(rows, cols)        // P(1:rows, 1:cols)
 P.topRightCorner(rows, cols)       // P(1:rows, end-cols+1:end)
 P.bottomLeftCorner(rows, cols)     // P(end-rows+1:end, 1:cols)
diff --git a/doc/SparseQuickReference.dox b/doc/SparseQuickReference.dox
index 4a33d0cc9..d04ac35c5 100644
--- a/doc/SparseQuickReference.dox
+++ b/doc/SparseQuickReference.dox
@@ -71,11 +71,10 @@ i.e either row major or column major. The default is column major. Most arithmet
 <td> Constant or Random Insertion</td>
 <td>
 \code
-sm1.setZero(); // Set the matrix with zero elements
-sm1.setConstant(val); //Replace all the nonzero values with val
+sm1.setZero();
 \endcode
 </td>
-<td> The matrix sm1 should have been created before ???</td>
+<td>Remove all non-zero coefficients</td>
 </tr>
 </table>
 
diff --git a/doc/examples/CMakeLists.txt b/doc/examples/CMakeLists.txt
index 71b272a15..08cf8efd7 100644
--- a/doc/examples/CMakeLists.txt
+++ b/doc/examples/CMakeLists.txt
@@ -6,12 +6,10 @@ foreach(example_src ${examples_SRCS})
   if(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO)
     target_link_libraries(${example} ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO})
   endif()
-  get_target_property(example_executable
-                      ${example} LOCATION)
   add_custom_command(
     TARGET ${example}
     POST_BUILD
-    COMMAND ${example_executable}
+    COMMAND ${example}
     ARGS >${CMAKE_CURRENT_BINARY_DIR}/${example}.out
   )
   add_dependencies(all_examples ${example})
diff --git a/doc/snippets/CMakeLists.txt b/doc/snippets/CMakeLists.txt
index 92a22ea61..1135900cf 100644
--- a/doc/snippets/CMakeLists.txt
+++ b/doc/snippets/CMakeLists.txt
@@ -14,12 +14,10 @@ foreach(snippet_src ${snippets_SRCS})
   if(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO)
     target_link_libraries(${compile_snippet_target} ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO})
   endif()
-  get_target_property(compile_snippet_executable
-                      ${compile_snippet_target} LOCATION)
   add_custom_command(
     TARGET ${compile_snippet_target}
     POST_BUILD
-    COMMAND ${compile_snippet_executable}
+    COMMAND ${compile_snippet_target}
     ARGS >${CMAKE_CURRENT_BINARY_DIR}/${snippet}.out
   )
   add_dependencies(all_snippets ${compile_snippet_target})
@@ -27,4 +25,4 @@ foreach(snippet_src ${snippets_SRCS})
                               PROPERTIES OBJECT_DEPENDS ${snippet_src})
 endforeach(snippet_src)
 
-ei_add_target_property(compile_tut_arithmetic_transpose_aliasing COMPILE_FLAGS -DEIGEN_NO_DEBUG)
\ No newline at end of file
+ei_add_target_property(compile_tut_arithmetic_transpose_aliasing COMPILE_FLAGS -DEIGEN_NO_DEBUG)
diff --git a/doc/special_examples/CMakeLists.txt b/doc/special_examples/CMakeLists.txt
index 0c9b3c3ba..f8a0148c8 100644
--- a/doc/special_examples/CMakeLists.txt
+++ b/doc/special_examples/CMakeLists.txt
@@ -1,4 +1,3 @@
-
 if(NOT EIGEN_TEST_NOQT)
   find_package(Qt4)
   if(QT4_FOUND)
@@ -6,16 +5,16 @@ if(NOT EIGEN_TEST_NOQT)
   endif()
 endif(NOT EIGEN_TEST_NOQT)
 
-
 if(QT4_FOUND)
   add_executable(Tutorial_sparse_example Tutorial_sparse_example.cpp Tutorial_sparse_example_details.cpp)
   target_link_libraries(Tutorial_sparse_example ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO} ${QT_QTCORE_LIBRARY} ${QT_QTGUI_LIBRARY})
-  
+
   add_custom_command(
       TARGET Tutorial_sparse_example
       POST_BUILD
       COMMAND Tutorial_sparse_example ARGS ${CMAKE_CURRENT_BINARY_DIR}/../html/Tutorial_sparse_example.jpeg
   )
-                     
+
   add_dependencies(all_examples Tutorial_sparse_example)
 endif(QT4_FOUND)
+
diff --git a/failtest/CMakeLists.txt b/failtest/CMakeLists.txt
index 5afa2ac82..5c4860237 100644
--- a/failtest/CMakeLists.txt
+++ b/failtest/CMakeLists.txt
@@ -26,6 +26,12 @@ ei_add_failtest("block_on_const_type_actually_const_1")
 ei_add_failtest("transpose_on_const_type_actually_const")
 ei_add_failtest("diagonal_on_const_type_actually_const")
 
+ei_add_failtest("ref_1")
+ei_add_failtest("ref_2")
+ei_add_failtest("ref_3")
+ei_add_failtest("ref_4")
+ei_add_failtest("ref_5")
+
 if (EIGEN_FAILTEST_FAILURE_COUNT)
   message(FATAL_ERROR
           "${EIGEN_FAILTEST_FAILURE_COUNT} out of ${EIGEN_FAILTEST_COUNT} failtests FAILED. "
diff --git a/failtest/ref_1.cpp b/failtest/ref_1.cpp
new file mode 100644
index 000000000..8b798d53d
--- /dev/null
+++ b/failtest/ref_1.cpp
@@ -0,0 +1,18 @@
+#include "../Eigen/Core"
+
+#ifdef EIGEN_SHOULD_FAIL_TO_BUILD
+#define CV_QUALIFIER const
+#else
+#define CV_QUALIFIER
+#endif
+
+using namespace Eigen;
+
+void call_ref(Ref<VectorXf> a) { }
+
+int main()
+{
+  VectorXf a(10);
+  CV_QUALIFIER VectorXf& ac(a);
+  call_ref(ac);
+}
diff --git a/failtest/ref_2.cpp b/failtest/ref_2.cpp
new file mode 100644
index 000000000..0b779ccf5
--- /dev/null
+++ b/failtest/ref_2.cpp
@@ -0,0 +1,15 @@
+#include "../Eigen/Core"
+
+using namespace Eigen;
+
+void call_ref(Ref<VectorXf> a) { }
+
+int main()
+{
+  MatrixXf A(10,10);
+#ifdef EIGEN_SHOULD_FAIL_TO_BUILD
+  call_ref(A.row(3));
+#else
+  call_ref(A.col(3));
+#endif
+}
diff --git a/failtest/ref_3.cpp b/failtest/ref_3.cpp
new file mode 100644
index 000000000..f46027d48
--- /dev/null
+++ b/failtest/ref_3.cpp
@@ -0,0 +1,15 @@
+#include "../Eigen/Core"
+
+using namespace Eigen;
+
+#ifdef EIGEN_SHOULD_FAIL_TO_BUILD
+void call_ref(Ref<VectorXf> a) { }
+#else
+void call_ref(const Ref<const VectorXf> &a) { }
+#endif
+
+int main()
+{
+  VectorXf a(10);
+  call_ref(a+a);
+}
diff --git a/failtest/ref_4.cpp b/failtest/ref_4.cpp
new file mode 100644
index 000000000..6c11fa4cb
--- /dev/null
+++ b/failtest/ref_4.cpp
@@ -0,0 +1,15 @@
+#include "../Eigen/Core"
+
+using namespace Eigen;
+
+void call_ref(Ref<MatrixXf,0,OuterStride<> > a) {}
+
+int main()
+{
+  MatrixXf A(10,10);
+#ifdef EIGEN_SHOULD_FAIL_TO_BUILD
+  call_ref(A.transpose());
+#else
+  call_ref(A);
+#endif
+}
diff --git a/failtest/ref_5.cpp b/failtest/ref_5.cpp
new file mode 100644
index 000000000..846d52795
--- /dev/null
+++ b/failtest/ref_5.cpp
@@ -0,0 +1,16 @@
+#include "../Eigen/Core"
+
+using namespace Eigen;
+
+void call_ref(Ref<VectorXf> a) { }
+
+int main()
+{
+  VectorXf a(10);
+  DenseBase<VectorXf> &ac(a);
+#ifdef EIGEN_SHOULD_FAIL_TO_BUILD
+  call_ref(ac);
+#else
+  call_ref(ac.derived());
+#endif
+}
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index ccb0fc798..d32451df6 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -66,7 +66,7 @@ endif()
 
 find_package(Pastix)
 find_package(Scotch)
-find_package(Metis)
+find_package(Metis 5.0 REQUIRED)
 if(PASTIX_FOUND AND BLAS_FOUND)
   add_definitions("-DEIGEN_PASTIX_SUPPORT")
   include_directories(${PASTIX_INCLUDES})
@@ -279,6 +279,7 @@ ei_add_property(EIGEN_TESTING_SUMMARY "CXX_FLAGS:         ${CMAKE_CXX_FLAGS}\n")
 ei_add_property(EIGEN_TESTING_SUMMARY "Sparse lib flags:  ${SPARSE_LIBS}\n")
 
 option(EIGEN_TEST_EIGEN2 "Run whole Eigen2 test suite against EIGEN2_SUPPORT" OFF)
+mark_as_advanced(EIGEN_TEST_EIGEN2)
 if(EIGEN_TEST_EIGEN2)
   add_subdirectory(eigen2)
 endif()
diff --git a/test/cholesky.cpp b/test/cholesky.cpp
index 64bcbccc4..56885deb7 100644
--- a/test/cholesky.cpp
+++ b/test/cholesky.cpp
@@ -320,33 +320,35 @@ template<typename MatrixType> void cholesky_definiteness(const MatrixType& m)
 {
   eigen_assert(m.rows() == 2 && m.cols() == 2);
   MatrixType mat;
+  LDLT<MatrixType> ldlt(2);
+  
   {
     mat << 1, 0, 0, -1;
-    LDLT<MatrixType> ldlt(mat);
+    ldlt.compute(mat);
     VERIFY(!ldlt.isNegative());
     VERIFY(!ldlt.isPositive());
   }
   {
     mat << 1, 2, 2, 1;
-    LDLT<MatrixType> ldlt(mat);
+    ldlt.compute(mat);
     VERIFY(!ldlt.isNegative());
     VERIFY(!ldlt.isPositive());
   }
   {
     mat << 0, 0, 0, 0;
-    LDLT<MatrixType> ldlt(mat);
+    ldlt.compute(mat);
     VERIFY(ldlt.isNegative());
     VERIFY(ldlt.isPositive());
   }
   {
     mat << 0, 0, 0, 1;
-    LDLT<MatrixType> ldlt(mat);
+    ldlt.compute(mat);
     VERIFY(!ldlt.isNegative());
     VERIFY(ldlt.isPositive());
   }
   {
     mat << -1, 0, 0, 0;
-    LDLT<MatrixType> ldlt(mat);
+    ldlt.compute(mat);
     VERIFY(ldlt.isNegative());
     VERIFY(!ldlt.isPositive());
   }
diff --git a/test/cwiseop.cpp b/test/cwiseop.cpp
index 247fa2a09..e3361da17 100644
--- a/test/cwiseop.cpp
+++ b/test/cwiseop.cpp
@@ -9,6 +9,8 @@
 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
 #define EIGEN2_SUPPORT
+#define EIGEN_NO_EIGEN2_DEPRECATED_WARNING
+
 #define EIGEN_NO_STATIC_ASSERT
 #include "main.h"
 #include <functional>
diff --git a/test/eigen2/CMakeLists.txt b/test/eigen2/CMakeLists.txt
index 84931e037..9615a6026 100644
--- a/test/eigen2/CMakeLists.txt
+++ b/test/eigen2/CMakeLists.txt
@@ -4,6 +4,7 @@ add_dependencies(eigen2_check eigen2_buildtests)
 add_dependencies(buildtests eigen2_buildtests)
 
 add_definitions("-DEIGEN2_SUPPORT_STAGE10_FULL_EIGEN2_API")
+add_definitions("-DEIGEN_NO_EIGEN2_DEPRECATED_WARNING")
 
 ei_add_test(eigen2_meta)
 ei_add_test(eigen2_sizeof)
diff --git a/test/eigen2/eigen2_adjoint.cpp b/test/eigen2/eigen2_adjoint.cpp
index 8ec9c9202..c0f811459 100644
--- a/test/eigen2/eigen2_adjoint.cpp
+++ b/test/eigen2/eigen2_adjoint.cpp
@@ -29,8 +29,6 @@ template<typename MatrixType> void adjoint(const MatrixType& m)
   MatrixType m1 = MatrixType::Random(rows, cols),
              m2 = MatrixType::Random(rows, cols),
              m3(rows, cols),
-             mzero = MatrixType::Zero(rows, cols),
-             identity = SquareMatrixType::Identity(rows, rows),
              square = SquareMatrixType::Random(rows, rows);
   VectorType v1 = VectorType::Random(rows),
              v2 = VectorType::Random(rows),
diff --git a/test/eigen2/eigen2_basicstuff.cpp b/test/eigen2/eigen2_basicstuff.cpp
index 4fa16d5ae..dd2dec1ef 100644
--- a/test/eigen2/eigen2_basicstuff.cpp
+++ b/test/eigen2/eigen2_basicstuff.cpp
@@ -23,11 +23,8 @@ template<typename MatrixType> void basicStuff(const MatrixType& m)
              m2 = MatrixType::Random(rows, cols),
              m3(rows, cols),
              mzero = MatrixType::Zero(rows, cols),
-             identity = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>
-                              ::Identity(rows, rows),
              square = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>::Random(rows, rows);
   VectorType v1 = VectorType::Random(rows),
-             v2 = VectorType::Random(rows),
              vzero = VectorType::Zero(rows);
 
   Scalar x = ei_random<Scalar>();
diff --git a/test/eigen2/eigen2_cwiseop.cpp b/test/eigen2/eigen2_cwiseop.cpp
index 4391d19b4..22e1cc342 100644
--- a/test/eigen2/eigen2_cwiseop.cpp
+++ b/test/eigen2/eigen2_cwiseop.cpp
@@ -35,11 +35,8 @@ template<typename MatrixType> void cwiseops(const MatrixType& m)
              mzero = MatrixType::Zero(rows, cols),
              mones = MatrixType::Ones(rows, cols),
              identity = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>
-                              ::Identity(rows, rows),
-             square = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>::Random(rows, rows);
-  VectorType v1 = VectorType::Random(rows),
-             v2 = VectorType::Random(rows),
-             vzero = VectorType::Zero(rows),
+                              ::Identity(rows, rows);
+  VectorType vzero = VectorType::Zero(rows),
              vones = VectorType::Ones(rows),
              v3(rows);
 
diff --git a/test/eigen2/eigen2_geometry.cpp b/test/eigen2/eigen2_geometry.cpp
index 70b4ab5f1..514040774 100644
--- a/test/eigen2/eigen2_geometry.cpp
+++ b/test/eigen2/eigen2_geometry.cpp
@@ -392,6 +392,7 @@ template<typename Scalar> void geometry(void)
   #define VERIFY_EULER(I,J,K, X,Y,Z) { \
     Vector3 ea = m.eulerAngles(I,J,K); \
     Matrix3 m1 = Matrix3(AngleAxisx(ea[0], Vector3::Unit##X()) * AngleAxisx(ea[1], Vector3::Unit##Y()) * AngleAxisx(ea[2], Vector3::Unit##Z())); \
+    VERIFY_IS_APPROX(m, m1); \
     VERIFY_IS_APPROX(m,  Matrix3(AngleAxisx(ea[0], Vector3::Unit##X()) * AngleAxisx(ea[1], Vector3::Unit##Y()) * AngleAxisx(ea[2], Vector3::Unit##Z()))); \
   }
   VERIFY_EULER(0,1,2, X,Y,Z);
diff --git a/test/eigen2/eigen2_geometry_with_eigen2_prefix.cpp b/test/eigen2/eigen2_geometry_with_eigen2_prefix.cpp
index f83b57249..12d4a71c3 100644
--- a/test/eigen2/eigen2_geometry_with_eigen2_prefix.cpp
+++ b/test/eigen2/eigen2_geometry_with_eigen2_prefix.cpp
@@ -394,6 +394,7 @@ template<typename Scalar> void geometry(void)
   #define VERIFY_EULER(I,J,K, X,Y,Z) { \
     Vector3 ea = m.eulerAngles(I,J,K); \
     Matrix3 m1 = Matrix3(AngleAxisx(ea[0], Vector3::Unit##X()) * AngleAxisx(ea[1], Vector3::Unit##Y()) * AngleAxisx(ea[2], Vector3::Unit##Z())); \
+    VERIFY_IS_APPROX(m, m1); \
     VERIFY_IS_APPROX(m,  Matrix3(AngleAxisx(ea[0], Vector3::Unit##X()) * AngleAxisx(ea[1], Vector3::Unit##Y()) * AngleAxisx(ea[2], Vector3::Unit##Z()))); \
   }
   VERIFY_EULER(0,1,2, X,Y,Z);
diff --git a/test/eigen2/eigen2_inverse.cpp b/test/eigen2/eigen2_inverse.cpp
index 5de1dfefa..ccd24a194 100644
--- a/test/eigen2/eigen2_inverse.cpp
+++ b/test/eigen2/eigen2_inverse.cpp
@@ -25,7 +25,6 @@ template<typename MatrixType> void inverse(const MatrixType& m)
 
   MatrixType m1 = MatrixType::Random(rows, cols),
              m2(rows, cols),
-             mzero = MatrixType::Zero(rows, cols),
              identity = MatrixType::Identity(rows, rows);
 
   while(ei_abs(m1.determinant()) < RealScalar(0.1) && rows <= 8)
diff --git a/test/eigen2/eigen2_linearstructure.cpp b/test/eigen2/eigen2_linearstructure.cpp
index 22f02c396..488f4c485 100644
--- a/test/eigen2/eigen2_linearstructure.cpp
+++ b/test/eigen2/eigen2_linearstructure.cpp
@@ -25,8 +25,7 @@ template<typename MatrixType> void linearStructure(const MatrixType& m)
   // to test it, hence I consider that we will have tested Random.h
   MatrixType m1 = MatrixType::Random(rows, cols),
              m2 = MatrixType::Random(rows, cols),
-             m3(rows, cols),
-             mzero = MatrixType::Zero(rows, cols);
+             m3(rows, cols);
 
   Scalar s1 = ei_random<Scalar>();
   while (ei_abs(s1)<1e-3) s1 = ei_random<Scalar>();
diff --git a/test/eigen2/eigen2_nomalloc.cpp b/test/eigen2/eigen2_nomalloc.cpp
index e234abe4b..d34a69999 100644
--- a/test/eigen2/eigen2_nomalloc.cpp
+++ b/test/eigen2/eigen2_nomalloc.cpp
@@ -25,22 +25,12 @@ template<typename MatrixType> void nomalloc(const MatrixType& m)
   */
 
   typedef typename MatrixType::Scalar Scalar;
-  typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
 
   int rows = m.rows();
   int cols = m.cols();
 
   MatrixType m1 = MatrixType::Random(rows, cols),
-             m2 = MatrixType::Random(rows, cols),
-             m3(rows, cols),
-             mzero = MatrixType::Zero(rows, cols),
-             identity = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>
-                              ::Identity(rows, rows),
-             square = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>
-                              ::Random(rows, rows);
-  VectorType v1 = VectorType::Random(rows),
-             v2 = VectorType::Random(rows),
-             vzero = VectorType::Zero(rows);
+             m2 = MatrixType::Random(rows, cols);
 
   Scalar s1 = ei_random<Scalar>();
 
diff --git a/test/eigen2/eigen2_submatrices.cpp b/test/eigen2/eigen2_submatrices.cpp
index c5d3f243d..dee970b63 100644
--- a/test/eigen2/eigen2_submatrices.cpp
+++ b/test/eigen2/eigen2_submatrices.cpp
@@ -51,16 +51,10 @@ template<typename MatrixType> void submatrices(const MatrixType& m)
   MatrixType m1 = MatrixType::Random(rows, cols),
              m2 = MatrixType::Random(rows, cols),
              m3(rows, cols),
-             mzero = MatrixType::Zero(rows, cols),
              ones = MatrixType::Ones(rows, cols),
-             identity = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>
-                              ::Identity(rows, rows),
              square = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>
                               ::Random(rows, rows);
-  VectorType v1 = VectorType::Random(rows),
-             v2 = VectorType::Random(rows),
-             v3 = VectorType::Random(rows),
-             vzero = VectorType::Zero(rows);
+  VectorType v1 = VectorType::Random(rows);
 
   Scalar s1 = ei_random<Scalar>();
 
diff --git a/test/eigen2/eigen2_triangular.cpp b/test/eigen2/eigen2_triangular.cpp
index 3748c7d71..6f17b7dff 100644
--- a/test/eigen2/eigen2_triangular.cpp
+++ b/test/eigen2/eigen2_triangular.cpp
@@ -13,7 +13,6 @@ template<typename MatrixType> void triangular(const MatrixType& m)
 {
   typedef typename MatrixType::Scalar Scalar;
   typedef typename NumTraits<Scalar>::Real RealScalar;
-  typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
 
   RealScalar largerEps = 10*test_precision<RealScalar>();
 
@@ -25,16 +24,7 @@ template<typename MatrixType> void triangular(const MatrixType& m)
              m3(rows, cols),
              m4(rows, cols),
              r1(rows, cols),
-             r2(rows, cols),
-             mzero = MatrixType::Zero(rows, cols),
-             mones = MatrixType::Ones(rows, cols),
-             identity = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>
-                              ::Identity(rows, rows),
-             square = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>
-                              ::Random(rows, rows);
-  VectorType v1 = VectorType::Random(rows),
-             v2 = VectorType::Random(rows),
-             vzero = VectorType::Zero(rows);
+             r2(rows, cols);
 
   MatrixType m1up = m1.template part<Eigen::UpperTriangular>();
   MatrixType m2up = m2.template part<Eigen::UpperTriangular>();
diff --git a/test/eigen2/product.h b/test/eigen2/product.h
index 2c9604d9a..ae1b4bae4 100644
--- a/test/eigen2/product.h
+++ b/test/eigen2/product.h
@@ -40,8 +40,7 @@ template<typename MatrixType> void product(const MatrixType& m)
   // to test it, hence I consider that we will have tested Random.h
   MatrixType m1 = MatrixType::Random(rows, cols),
              m2 = MatrixType::Random(rows, cols),
-             m3(rows, cols),
-             mzero = MatrixType::Zero(rows, cols);
+             m3(rows, cols);
   RowSquareMatrixType
              identity = RowSquareMatrixType::Identity(rows, rows),
              square = RowSquareMatrixType::Random(rows, rows),
@@ -49,9 +48,7 @@ template<typename MatrixType> void product(const MatrixType& m)
   ColSquareMatrixType
              square2 = ColSquareMatrixType::Random(cols, cols),
              res2 = ColSquareMatrixType::Random(cols, cols);
-  RowVectorType v1 = RowVectorType::Random(rows),
-             v2 = RowVectorType::Random(rows),
-             vzero = RowVectorType::Zero(rows);
+  RowVectorType v1 = RowVectorType::Random(rows);
   ColVectorType vc2 = ColVectorType::Random(cols), vcres(cols);
   OtherMajorMatrixType tm1 = m1;
 
diff --git a/test/eigen2support.cpp b/test/eigen2support.cpp
index ad1d98091..1fa49a8c8 100644
--- a/test/eigen2support.cpp
+++ b/test/eigen2support.cpp
@@ -8,6 +8,7 @@
 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
 #define EIGEN2_SUPPORT
+#define EIGEN_NO_EIGEN2_DEPRECATED_WARNING
 
 #include "main.h"
 
diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp
index 5c6ecd875..38689cfbf 100644
--- a/test/eigensolver_selfadjoint.cpp
+++ b/test/eigensolver_selfadjoint.cpp
@@ -29,7 +29,21 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
   MatrixType a = MatrixType::Random(rows,cols);
   MatrixType a1 = MatrixType::Random(rows,cols);
   MatrixType symmA =  a.adjoint() * a + a1.adjoint() * a1;
+  MatrixType symmC = symmA;
+  
+  // randomly nullify some rows/columns
+  {
+    Index count = 1;//internal::random<Index>(-cols,cols);
+    for(Index k=0; k<count; ++k)
+    {
+      Index i = internal::random<Index>(0,cols-1);
+      symmA.row(i).setZero();
+      symmA.col(i).setZero();
+    }
+  }
+  
   symmA.template triangularView<StrictlyUpper>().setZero();
+  symmC.template triangularView<StrictlyUpper>().setZero();
 
   MatrixType b = MatrixType::Random(rows,cols);
   MatrixType b1 = MatrixType::Random(rows,cols);
@@ -40,7 +54,7 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
   SelfAdjointEigenSolver<MatrixType> eiDirect;
   eiDirect.computeDirect(symmA);
   // generalized eigen pb
-  GeneralizedSelfAdjointEigenSolver<MatrixType> eiSymmGen(symmA, symmB);
+  GeneralizedSelfAdjointEigenSolver<MatrixType> eiSymmGen(symmC, symmB);
 
   VERIFY_IS_EQUAL(eiSymm.info(), Success);
   VERIFY((symmA.template selfadjointView<Lower>() * eiSymm.eigenvectors()).isApprox(
@@ -57,27 +71,28 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
   VERIFY_IS_APPROX(eiSymm.eigenvalues(), eiSymmNoEivecs.eigenvalues());
   
   // generalized eigen problem Ax = lBx
-  eiSymmGen.compute(symmA, symmB,Ax_lBx);
+  eiSymmGen.compute(symmC, symmB,Ax_lBx);
   VERIFY_IS_EQUAL(eiSymmGen.info(), Success);
-  VERIFY((symmA.template selfadjointView<Lower>() * eiSymmGen.eigenvectors()).isApprox(
+  VERIFY((symmC.template selfadjointView<Lower>() * eiSymmGen.eigenvectors()).isApprox(
           symmB.template selfadjointView<Lower>() * (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
 
   // generalized eigen problem BAx = lx
-  eiSymmGen.compute(symmA, symmB,BAx_lx);
+  eiSymmGen.compute(symmC, symmB,BAx_lx);
   VERIFY_IS_EQUAL(eiSymmGen.info(), Success);
-  VERIFY((symmB.template selfadjointView<Lower>() * (symmA.template selfadjointView<Lower>() * eiSymmGen.eigenvectors())).isApprox(
+  VERIFY((symmB.template selfadjointView<Lower>() * (symmC.template selfadjointView<Lower>() * eiSymmGen.eigenvectors())).isApprox(
          (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
 
   // generalized eigen problem ABx = lx
-  eiSymmGen.compute(symmA, symmB,ABx_lx);
+  eiSymmGen.compute(symmC, symmB,ABx_lx);
   VERIFY_IS_EQUAL(eiSymmGen.info(), Success);
-  VERIFY((symmA.template selfadjointView<Lower>() * (symmB.template selfadjointView<Lower>() * eiSymmGen.eigenvectors())).isApprox(
+  VERIFY((symmC.template selfadjointView<Lower>() * (symmB.template selfadjointView<Lower>() * eiSymmGen.eigenvectors())).isApprox(
          (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
 
 
+  eiSymm.compute(symmC);
   MatrixType sqrtSymmA = eiSymm.operatorSqrt();
-  VERIFY_IS_APPROX(MatrixType(symmA.template selfadjointView<Lower>()), sqrtSymmA*sqrtSymmA);
-  VERIFY_IS_APPROX(sqrtSymmA, symmA.template selfadjointView<Lower>()*eiSymm.operatorInverseSqrt());
+  VERIFY_IS_APPROX(MatrixType(symmC.template selfadjointView<Lower>()), sqrtSymmA*sqrtSymmA);
+  VERIFY_IS_APPROX(sqrtSymmA, symmC.template selfadjointView<Lower>()*eiSymm.operatorInverseSqrt());
 
   MatrixType id = MatrixType::Identity(rows, cols);
   VERIFY_IS_APPROX(id.template selfadjointView<Lower>().operatorNorm(), RealScalar(1));
@@ -95,15 +110,15 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
   VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorInverseSqrt());
 
   // test Tridiagonalization's methods
-  Tridiagonalization<MatrixType> tridiag(symmA);
+  Tridiagonalization<MatrixType> tridiag(symmC);
   // FIXME tridiag.matrixQ().adjoint() does not work
-  VERIFY_IS_APPROX(MatrixType(symmA.template selfadjointView<Lower>()), tridiag.matrixQ() * tridiag.matrixT().eval() * MatrixType(tridiag.matrixQ()).adjoint());
+  VERIFY_IS_APPROX(MatrixType(symmC.template selfadjointView<Lower>()), tridiag.matrixQ() * tridiag.matrixT().eval() * MatrixType(tridiag.matrixQ()).adjoint());
   
   if (rows > 1)
   {
     // Test matrix with NaN
-    symmA(0,0) = std::numeric_limits<typename MatrixType::RealScalar>::quiet_NaN();
-    SelfAdjointEigenSolver<MatrixType> eiSymmNaN(symmA);
+    symmC(0,0) = std::numeric_limits<typename MatrixType::RealScalar>::quiet_NaN();
+    SelfAdjointEigenSolver<MatrixType> eiSymmNaN(symmC);
     VERIFY_IS_EQUAL(eiSymmNaN.info(), NoConvergence);
   }
 }
@@ -113,8 +128,10 @@ void test_eigensolver_selfadjoint()
   int s = 0;
   for(int i = 0; i < g_repeat; i++) {
     // very important to test 3x3 and 2x2 matrices since we provide special paths for them
+    CALL_SUBTEST_1( selfadjointeigensolver(Matrix2f()) );
     CALL_SUBTEST_1( selfadjointeigensolver(Matrix2d()) );
     CALL_SUBTEST_1( selfadjointeigensolver(Matrix3f()) );
+    CALL_SUBTEST_1( selfadjointeigensolver(Matrix3d()) );
     CALL_SUBTEST_2( selfadjointeigensolver(Matrix4d()) );
     s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4);
     CALL_SUBTEST_3( selfadjointeigensolver(MatrixXf(s,s)) );
diff --git a/test/geo_hyperplane.cpp b/test/geo_hyperplane.cpp
index f26fc1329..327537801 100644
--- a/test/geo_hyperplane.cpp
+++ b/test/geo_hyperplane.cpp
@@ -114,6 +114,32 @@ template<typename Scalar> void lines()
   }
 }
 
+template<typename Scalar> void planes()
+{
+  using std::abs;
+  typedef Hyperplane<Scalar, 3> Plane;
+  typedef Matrix<Scalar,3,1> Vector;
+
+  for(int i = 0; i < 10; i++)
+  {
+    Vector v0 = Vector::Random();
+    Vector v1(v0), v2(v0);
+    if(internal::random<double>(0,1)>0.25)
+      v1 += Vector::Random();
+    if(internal::random<double>(0,1)>0.25)
+      v2 += v1 * std::pow(internal::random<Scalar>(0,1),internal::random<int>(1,16));
+    if(internal::random<double>(0,1)>0.25)
+      v2 += Vector::Random() * std::pow(internal::random<Scalar>(0,1),internal::random<int>(1,16));
+
+    Plane p0 = Plane::Through(v0, v1, v2);
+
+    VERIFY_IS_APPROX(p0.normal().norm(), Scalar(1));
+    VERIFY_IS_MUCH_SMALLER_THAN(p0.absDistance(v0), Scalar(1));
+    VERIFY_IS_MUCH_SMALLER_THAN(p0.absDistance(v1), Scalar(1));
+    VERIFY_IS_MUCH_SMALLER_THAN(p0.absDistance(v2), Scalar(1));
+  }
+}
+
 template<typename Scalar> void hyperplane_alignment()
 {
   typedef Hyperplane<Scalar,3,AutoAlign> Plane3a;
@@ -153,5 +179,7 @@ void test_geo_hyperplane()
     CALL_SUBTEST_4( hyperplane(Hyperplane<std::complex<double>,5>()) );
     CALL_SUBTEST_1( lines<float>() );
     CALL_SUBTEST_3( lines<double>() );
+    CALL_SUBTEST_2( planes<float>() );
+    CALL_SUBTEST_5( planes<double>() );
   }
 }
diff --git a/test/geo_transformations.cpp b/test/geo_transformations.cpp
index ee3030b5d..547765714 100644
--- a/test/geo_transformations.cpp
+++ b/test/geo_transformations.cpp
@@ -98,11 +98,19 @@ template<typename Scalar, int Mode, int Options> void transformations()
   Matrix3 matrot1, m;
 
   Scalar a = internal::random<Scalar>(-Scalar(M_PI), Scalar(M_PI));
-  Scalar s0 = internal::random<Scalar>();
+  Scalar s0 = internal::random<Scalar>(),
+         s1 = internal::random<Scalar>();
+  
+  while(v0.norm() < test_precision<Scalar>()) v0 = Vector3::Random();
+  while(v1.norm() < test_precision<Scalar>()) v1 = Vector3::Random();
+    
 
   VERIFY_IS_APPROX(v0, AngleAxisx(a, v0.normalized()) * v0);
   VERIFY_IS_APPROX(-v0, AngleAxisx(Scalar(M_PI), v0.unitOrthogonal()) * v0);
-  VERIFY_IS_APPROX(cos(a)*v0.squaredNorm(), v0.dot(AngleAxisx(a, v0.unitOrthogonal()) * v0));
+  if(abs(cos(a)) > test_precision<Scalar>())
+  {
+    VERIFY_IS_APPROX(cos(a)*v0.squaredNorm(), v0.dot(AngleAxisx(a, v0.unitOrthogonal()) * v0));
+  }
   m = AngleAxisx(a, v0.normalized()).toRotationMatrix().adjoint();
   VERIFY_IS_APPROX(Matrix3::Identity(), m * AngleAxisx(a, v0.normalized()));
   VERIFY_IS_APPROX(Matrix3::Identity(), AngleAxisx(a, v0.normalized()) * m);
@@ -123,11 +131,18 @@ template<typename Scalar, int Mode, int Options> void transformations()
   // angle-axis conversion
   AngleAxisx aa = AngleAxisx(q1);
   VERIFY_IS_APPROX(q1 * v1, Quaternionx(aa) * v1);
-  VERIFY_IS_NOT_APPROX(q1 * v1, Quaternionx(AngleAxisx(aa.angle()*2,aa.axis())) * v1);
+  
+  if(abs(aa.angle()) > NumTraits<Scalar>::dummy_precision())
+  {
+    VERIFY( !(q1 * v1).isApprox(Quaternionx(AngleAxisx(aa.angle()*2,aa.axis())) * v1) );
+  }
 
   aa.fromRotationMatrix(aa.toRotationMatrix());
   VERIFY_IS_APPROX(q1 * v1, Quaternionx(aa) * v1);
-  VERIFY_IS_NOT_APPROX(q1 * v1, Quaternionx(AngleAxisx(aa.angle()*2,aa.axis())) * v1);
+  if(abs(aa.angle()) > NumTraits<Scalar>::dummy_precision())
+  {
+    VERIFY( !(q1 * v1).isApprox(Quaternionx(AngleAxisx(aa.angle()*2,aa.axis())) * v1) );
+  }
 
   // AngleAxis
   VERIFY_IS_APPROX(AngleAxisx(a,v1.normalized()).toRotationMatrix(),
@@ -347,7 +362,9 @@ template<typename Scalar, int Mode, int Options> void transformations()
   // test transform inversion
   t0.setIdentity();
   t0.translate(v0);
-  t0.linear().setRandom();
+  do {
+    t0.linear().setRandom();
+  } while(t0.linear().jacobiSvd().singularValues()(2)<test_precision<Scalar>());
   Matrix4 t044 = Matrix4::Zero();
   t044(3,3) = 1;
   t044.block(0,0,t0.matrix().rows(),4) = t0.matrix();
@@ -397,6 +414,29 @@ template<typename Scalar, int Mode, int Options> void transformations()
   t20 = Translation2(v20) * (Rotation2D<Scalar>(s0) * Eigen::Scaling(s0));
   t21 = Translation2(v20) * Rotation2D<Scalar>(s0) * Eigen::Scaling(s0);
   VERIFY_IS_APPROX(t20,t21);
+  
+  Rotation2D<Scalar> R0(s0), R1(s1);
+  t20 = Translation2(v20) * (R0 * Eigen::Scaling(s0));
+  t21 = Translation2(v20) * R0 * Eigen::Scaling(s0);
+  VERIFY_IS_APPROX(t20,t21);
+  
+  t20 = Translation2(v20) * (R0 * R0.inverse() * Eigen::Scaling(s0));
+  t21 = Translation2(v20) * Eigen::Scaling(s0);
+  VERIFY_IS_APPROX(t20,t21);
+  
+  VERIFY_IS_APPROX(s0, (R0.slerp(0, R1)).angle());
+  VERIFY_IS_APPROX(s1, (R0.slerp(1, R1)).angle());
+  VERIFY_IS_APPROX(s0, (R0.slerp(0.5, R0)).angle());
+  VERIFY_IS_APPROX(Scalar(0), (R0.slerp(0.5, R0.inverse())).angle());
+  
+  // check basic features
+  {
+    Rotation2D<Scalar> r1;           // default ctor
+    r1 = Rotation2D<Scalar>(s0);     // copy assignment
+    VERIFY_IS_APPROX(r1.angle(),s0);
+    Rotation2D<Scalar> r2(r1);       // copy ctor
+    VERIFY_IS_APPROX(r2.angle(),s0);
+  }
 }
 
 template<typename Scalar> void transform_alignment()
diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp
index 7c21f0ab3..12c556e59 100644
--- a/test/jacobisvd.cpp
+++ b/test/jacobisvd.cpp
@@ -321,16 +321,23 @@ void jacobisvd_inf_nan()
   VERIFY(sub(some_inf, some_inf) != sub(some_inf, some_inf));
   svd.compute(MatrixType::Constant(10,10,some_inf), ComputeFullU | ComputeFullV);
 
-  Scalar some_nan = zero<Scalar>() / zero<Scalar>();
-  VERIFY(some_nan != some_nan);
-  svd.compute(MatrixType::Constant(10,10,some_nan), ComputeFullU | ComputeFullV);
+  Scalar nan = std::numeric_limits<Scalar>::quiet_NaN();
+  VERIFY(nan != nan);
+  svd.compute(MatrixType::Constant(10,10,nan), ComputeFullU | ComputeFullV);
 
   MatrixType m = MatrixType::Zero(10,10);
   m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_inf;
   svd.compute(m, ComputeFullU | ComputeFullV);
 
   m = MatrixType::Zero(10,10);
-  m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_nan;
+  m(internal::random<int>(0,9), internal::random<int>(0,9)) = nan;
+  svd.compute(m, ComputeFullU | ComputeFullV);
+  
+  // regression test for bug 791
+  m.resize(3,3);
+  m << 0,    2*NumTraits<Scalar>::epsilon(),  0.5,
+       0,   -0.5,                             0,
+       nan,  0,                               0;
   svd.compute(m, ComputeFullU | ComputeFullV);
 }
 
@@ -434,6 +441,7 @@ void test_jacobisvd()
 
     // Test on inf/nan matrix
     CALL_SUBTEST_7( jacobisvd_inf_nan<MatrixXf>() );
+    CALL_SUBTEST_10( jacobisvd_inf_nan<MatrixXd>() );
   }
 
   CALL_SUBTEST_7(( jacobisvd<MatrixXf>(MatrixXf(internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2), internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2))) ));
diff --git a/test/main.h b/test/main.h
index 14f0d2f78..664204866 100644
--- a/test/main.h
+++ b/test/main.h
@@ -17,13 +17,36 @@
 #include <sstream>
 #include <vector>
 #include <typeinfo>
+
+// The following includes of STL headers have to be done _before_ the
+// definition of macros min() and max().  The reason is that many STL
+// implementations will not work properly as the min and max symbols collide
+// with the STL functions std:min() and std::max().  The STL headers may check
+// for the macro definition of min/max and issue a warning or undefine the
+// macros.
+//
+// Still, Windows defines min() and max() in windef.h as part of the regular
+// Windows system interfaces and many other Windows APIs depend on these
+// macros being available.  To prevent the macro expansion of min/max and to
+// make Eigen compatible with the Windows environment all function calls of
+// std::min() and std::max() have to be written with parenthesis around the
+// function name.
+//
+// All STL headers used by Eigen should be included here.  Because main.h is
+// included before any Eigen header and because the STL headers are guarded
+// against multiple inclusions, no STL header will see our own min/max macro
+// definitions.
 #include <limits>
 #include <algorithm>
-#include <sstream>
 #include <complex>
 #include <deque>
 #include <queue>
+#include <list>
 
+// To test that all calls from Eigen code to std::min() and std::max() are
+// protected by parenthesis against macro expansion, the min()/max() macros
+// are defined here and any not-parenthesized min/max call will cause a
+// compiler error.
 #define min(A,B) please_protect_your_min_with_parentheses
 #define max(A,B) please_protect_your_max_with_parentheses
 
@@ -383,6 +406,26 @@ void randomPermutationVector(PermutationVectorType& v, typename PermutationVecto
   }
 }
 
+template<typename T> bool isNotNaN(const T& x)
+{
+  return x==x;
+}
+
+template<typename T> bool isNaN(const T& x)
+{
+  return x!=x;
+}
+
+template<typename T> bool isInf(const T& x)
+{
+  return x > NumTraits<T>::highest();
+}
+
+template<typename T> bool isMinusInf(const T& x)
+{
+  return x < NumTraits<T>::lowest();
+}
+
 } // end namespace Eigen
 
 template<typename T> struct GetDifferentType;
diff --git a/test/nomalloc.cpp b/test/nomalloc.cpp
index cbd02dd21..8e0402358 100644
--- a/test/nomalloc.cpp
+++ b/test/nomalloc.cpp
@@ -165,6 +165,38 @@ void ctms_decompositions()
   Eigen::JacobiSVD<Matrix> jSVD; jSVD.compute(A, ComputeFullU | ComputeFullV);
 }
 
+void test_zerosized() {
+  // default constructors:
+  Eigen::MatrixXd A;
+  Eigen::VectorXd v;
+  // explicit zero-sized:
+  Eigen::ArrayXXd A0(0,0);
+  Eigen::ArrayXd v0(std::ptrdiff_t(0)); // FIXME ArrayXd(0) is ambiguous
+
+  // assigning empty objects to each other:
+  A=A0;
+  v=v0;
+}
+
+template<typename MatrixType> void test_reference(const MatrixType& m) {
+  typedef typename MatrixType::Scalar Scalar;
+  enum { Flag          =  MatrixType::IsRowMajor ? Eigen::RowMajor : Eigen::ColMajor};
+  enum { TransposeFlag = !MatrixType::IsRowMajor ? Eigen::RowMajor : Eigen::ColMajor};
+  typename MatrixType::Index rows = m.rows(), cols=m.cols();
+  // Dynamic reference:
+  typedef Eigen::Ref<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Flag         > > Ref;
+  typedef Eigen::Ref<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, TransposeFlag> > RefT;
+
+  Ref r1(m);
+  Ref r2(m.block(rows/3, cols/4, rows/2, cols/2));
+  RefT r3(m.transpose());
+  RefT r4(m.topLeftCorner(rows/2, cols/2).transpose());
+
+  VERIFY_RAISES_ASSERT(RefT r5(m));
+  VERIFY_RAISES_ASSERT(Ref r6(m.transpose()));
+  VERIFY_RAISES_ASSERT(Ref r7(Scalar(2) * m));
+}
+
 void test_nomalloc()
 {
   // check that our operator new is indeed called:
@@ -175,5 +207,6 @@ void test_nomalloc()
   
   // Check decomposition modules with dynamic matrices that have a known compile-time max size (ctms)
   CALL_SUBTEST_4(ctms_decompositions<float>());
-
+  CALL_SUBTEST_5(test_zerosized());
+  CALL_SUBTEST_6(test_reference(Matrix<float,32,32>()));
 }
diff --git a/test/nullary.cpp b/test/nullary.cpp
index 5408d88b2..fbc721a1a 100644
--- a/test/nullary.cpp
+++ b/test/nullary.cpp
@@ -80,7 +80,9 @@ void testVectorType(const VectorType& base)
   Matrix<Scalar,1,Dynamic> col_vector(size);
   row_vector.setLinSpaced(size,low,high);
   col_vector.setLinSpaced(size,low,high);
-  VERIFY( row_vector.isApprox(col_vector.transpose(), NumTraits<Scalar>::epsilon()));
+  // when using the extended precision (e.g., FPU) the relative error might exceed 1 bit
+  // when computing the squared sum in isApprox, thus the 2x factor.
+  VERIFY( row_vector.isApprox(col_vector.transpose(), Scalar(2)*NumTraits<Scalar>::epsilon()));
 
   Matrix<Scalar,Dynamic,1> size_changer(size+50);
   size_changer.setLinSpaced(size,low,high);
diff --git a/test/packetmath.cpp b/test/packetmath.cpp
index 2c0519c41..38aa256ce 100644
--- a/test/packetmath.cpp
+++ b/test/packetmath.cpp
@@ -239,6 +239,12 @@ template<typename Scalar> void packetmath_real()
     data2[i] = internal::random<Scalar>(-87,88);
   }
   CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasExp, std::exp, internal::pexp);
+  {
+    data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
+    packet_helper<internal::packet_traits<Scalar>::HasExp,Packet> h;
+    h.store(data2, internal::pexp(h.load(data1))); 
+    VERIFY(isNaN(data2[0]));
+  }
 
   for (int i=0; i<size; ++i)
   {
@@ -247,8 +253,22 @@ template<typename Scalar> void packetmath_real()
   }
   if(internal::random<float>(0,1)<0.1)
     data1[internal::random<int>(0, PacketSize)] = 0;
-  CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasLog, std::log, internal::plog);
   CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasSqrt, std::sqrt, internal::psqrt);
+  CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasLog, std::log, internal::plog);
+  {
+    data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
+    packet_helper<internal::packet_traits<Scalar>::HasLog,Packet> h;
+    h.store(data2, internal::plog(h.load(data1)));
+    VERIFY(isNaN(data2[0]));
+    data1[0] = -1.0f;
+    h.store(data2, internal::plog(h.load(data1)));
+    VERIFY(isNaN(data2[0]));
+#if !EIGEN_FAST_MATH
+    h.store(data2, internal::psqrt(h.load(data1)));
+    VERIFY(isNaN(data2[0]));
+    VERIFY(isNaN(data2[1]));
+#endif
+  }
 }
 
 template<typename Scalar> void packetmath_notcomplex()
diff --git a/test/product.h b/test/product.h
index 856b234ac..0b3abe402 100644
--- a/test/product.h
+++ b/test/product.h
@@ -139,4 +139,12 @@ template<typename MatrixType> void product(const MatrixType& m)
   // inner product
   Scalar x = square2.row(c) * square2.col(c2);
   VERIFY_IS_APPROX(x, square2.row(c).transpose().cwiseProduct(square2.col(c2)).sum());
+  
+  // outer product
+  VERIFY_IS_APPROX(m1.col(c) * m1.row(r), m1.block(0,c,rows,1) * m1.block(r,0,1,cols));
+  VERIFY_IS_APPROX(m1.row(r).transpose() * m1.col(c).transpose(), m1.block(r,0,1,cols).transpose() * m1.block(0,c,rows,1).transpose());
+  VERIFY_IS_APPROX(m1.block(0,c,rows,1) * m1.row(r), m1.block(0,c,rows,1) * m1.block(r,0,1,cols));
+  VERIFY_IS_APPROX(m1.col(c) * m1.block(r,0,1,cols), m1.block(0,c,rows,1) * m1.block(r,0,1,cols));
+  VERIFY_IS_APPROX(m1.leftCols(1) * m1.row(r), m1.block(0,0,rows,1) * m1.block(r,0,1,cols));
+  VERIFY_IS_APPROX(m1.col(c) * m1.topRows(1), m1.block(0,c,rows,1) * m1.block(0,0,1,cols));  
 }
diff --git a/test/ref.cpp b/test/ref.cpp
index 19e81549c..32eb31048 100644
--- a/test/ref.cpp
+++ b/test/ref.cpp
@@ -183,15 +183,15 @@ void call_ref()
 
   VERIFY_EVALUATION_COUNT( call_ref_1(a,a), 0);
   VERIFY_EVALUATION_COUNT( call_ref_1(b,b.transpose()), 0);
-//   call_ref_1(ac);           // does not compile because ac is const
+//   call_ref_1(ac,a<c);           // does not compile because ac is const
   VERIFY_EVALUATION_COUNT( call_ref_1(ab,ab), 0);
   VERIFY_EVALUATION_COUNT( call_ref_1(a.head(4),a.head(4)), 0);
   VERIFY_EVALUATION_COUNT( call_ref_1(abc,abc), 0);
   VERIFY_EVALUATION_COUNT( call_ref_1(A.col(3),A.col(3)), 0);
-//   call_ref_1(A.row(3));    // does not compile because innerstride!=1
+//   call_ref_1(A.row(3),A.row(3));    // does not compile because innerstride!=1
   VERIFY_EVALUATION_COUNT( call_ref_3(A.row(3),A.row(3).transpose()), 0);
   VERIFY_EVALUATION_COUNT( call_ref_4(A.row(3),A.row(3).transpose()), 0);
-//   call_ref_1(a+a);          // does not compile for obvious reason
+//   call_ref_1(a+a, a+a);          // does not compile for obvious reason
 
   MatrixXf tmp = A*A.col(1);
   VERIFY_EVALUATION_COUNT( call_ref_2(A*A.col(1), tmp), 1);     // evaluated into a temp
@@ -212,7 +212,7 @@ void call_ref()
   VERIFY_EVALUATION_COUNT( call_ref_5(a,a), 0);
   VERIFY_EVALUATION_COUNT( call_ref_5(a.head(3),a.head(3)), 0);
   VERIFY_EVALUATION_COUNT( call_ref_5(A,A), 0);
-//   call_ref_5(A.transpose());   // does not compile
+//   call_ref_5(A.transpose(),A.transpose());   // does not compile because storage order does not match
   VERIFY_EVALUATION_COUNT( call_ref_5(A.block(1,1,2,2),A.block(1,1,2,2)), 0);
   VERIFY_EVALUATION_COUNT( call_ref_5(b,b), 0);             // storage order do not match, but this is a degenerate case that should work
   VERIFY_EVALUATION_COUNT( call_ref_5(a.row(3),a.row(3)), 0);
diff --git a/test/sparse_solver.h b/test/sparse_solver.h
index d84aff070..244e81c1b 100644
--- a/test/sparse_solver.h
+++ b/test/sparse_solver.h
@@ -124,7 +124,23 @@ void check_sparse_determinant(Solver& solver, const typename Solver::MatrixType&
   Scalar refDet = dA.determinant();
   VERIFY_IS_APPROX(refDet,solver.determinant());
 }
+template<typename Solver, typename DenseMat>
+void check_sparse_abs_determinant(Solver& solver, const typename Solver::MatrixType& A, const DenseMat& dA)
+{
+  using std::abs;
+  typedef typename Solver::MatrixType Mat;
+  typedef typename Mat::Scalar Scalar;
+  
+  solver.compute(A);
+  if (solver.info() != Success)
+  {
+    std::cerr << "sparse solver testing: factorization failed (check_sparse_abs_determinant)\n";
+    return;
+  }
 
+  Scalar refDet = abs(dA.determinant());
+  VERIFY_IS_APPROX(refDet,solver.absDeterminant());
+}
 
 template<typename Solver, typename DenseMat>
 int generate_sparse_spd_problem(Solver& , typename Solver::MatrixType& A, typename Solver::MatrixType& halfA, DenseMat& dA, int maxSize = 300)
@@ -324,3 +340,20 @@ template<typename Solver> void check_sparse_square_determinant(Solver& solver)
     check_sparse_determinant(solver, A, dA);
   }
 }
+
+template<typename Solver> void check_sparse_square_abs_determinant(Solver& solver)
+{
+  typedef typename Solver::MatrixType Mat;
+  typedef typename Mat::Scalar Scalar;
+  typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
+
+  // generate the problem
+  Mat A;
+  DenseMatrix dA;
+  generate_sparse_square_problem(solver, A, dA, 30);
+  A.makeCompressed();
+  for (int i = 0; i < g_repeat; i++) {
+    check_sparse_abs_determinant(solver, A, dA);
+  }
+}
+
diff --git a/test/sparselu.cpp b/test/sparselu.cpp
index 37980defc..52371cb12 100644
--- a/test/sparselu.cpp
+++ b/test/sparselu.cpp
@@ -44,6 +44,9 @@ template<typename T> void test_sparselu_T()
   check_sparse_square_solving(sparselu_colamd); 
   check_sparse_square_solving(sparselu_amd);
   check_sparse_square_solving(sparselu_natural);
+  
+  check_sparse_square_abs_determinant(sparselu_colamd);
+  check_sparse_square_abs_determinant(sparselu_amd);
 }
 
 void test_sparselu()
diff --git a/test/sparseqr.cpp b/test/sparseqr.cpp
index 1fe4a98ee..451c0e7f8 100644
--- a/test/sparseqr.cpp
+++ b/test/sparseqr.cpp
@@ -10,12 +10,11 @@
 #include <Eigen/SparseQR>
 
 template<typename MatrixType,typename DenseMat>
-int generate_sparse_rectangular_problem(MatrixType& A, DenseMat& dA, int maxRows = 300, int maxCols = 150)
+int generate_sparse_rectangular_problem(MatrixType& A, DenseMat& dA, int maxRows = 300)
 {
-  eigen_assert(maxRows >= maxCols);
   typedef typename MatrixType::Scalar Scalar;
   int rows = internal::random<int>(1,maxRows);
-  int cols = internal::random<int>(1,maxCols);
+  int cols = internal::random<int>(1,rows);
   double density = (std::max)(8./(rows*cols), 0.01);
   
   A.resize(rows,cols);
@@ -54,6 +53,8 @@ template<typename Scalar> void test_sparseqr_scalar()
   
   b = dA * DenseVector::Random(A.cols());
   solver.compute(A);
+  if(internal::random<float>(0,1)>0.5)
+    solver.factorize(A);  // this checks that calling analyzePattern is not needed if the pattern do not change.
   if (solver.info() != Success)
   {
     std::cerr << "sparse QR factorization failed\n";
diff --git a/test/stable_norm.cpp b/test/stable_norm.cpp
index 549f91fbf..231dd9189 100644
--- a/test/stable_norm.cpp
+++ b/test/stable_norm.cpp
@@ -9,11 +9,6 @@
 
 #include "main.h"
 
-template<typename T> bool isNotNaN(const T& x)
-{
-  return x==x;
-}
-
 // workaround aggressive optimization in ICC
 template<typename T> EIGEN_DONT_INLINE  T sub(T a, T b) { return a - b; }
 
diff --git a/unsupported/Eigen/OpenGLSupport b/unsupported/Eigen/OpenGLSupport
index c4090ab11..e2769449c 100644
--- a/unsupported/Eigen/OpenGLSupport
+++ b/unsupported/Eigen/OpenGLSupport
@@ -178,11 +178,11 @@ template<typename Scalar> void glLoadMatrix(const Transform<Scalar,3,Affine>& t)
 template<typename Scalar> void glLoadMatrix(const Transform<Scalar,3,Projective>& t)    { glLoadMatrix(t.matrix()); }
 template<typename Scalar> void glLoadMatrix(const Transform<Scalar,3,AffineCompact>& t) { glLoadMatrix(Transform<Scalar,3,Affine>(t).matrix()); }
 
-static void glRotate(const Rotation2D<float>& rot)
+inline void glRotate(const Rotation2D<float>& rot)
 {
   glRotatef(rot.angle()*180.f/float(M_PI), 0.f, 0.f, 1.f);
 }
-static void glRotate(const Rotation2D<double>& rot)
+inline void glRotate(const Rotation2D<double>& rot)
 {
   glRotated(rot.angle()*180.0/M_PI, 0.0, 0.0, 1.0);
 }
@@ -246,18 +246,18 @@ EIGEN_GL_FUNC1_SPECIALIZATION_MAT(glGet,GLenum,_,double, 4,4,Doublev)
 
 #ifdef GL_VERSION_2_0
 
-static void glUniform2fv_ei  (GLint loc, const float* v)         { glUniform2fv(loc,1,v); }
-static void glUniform2iv_ei  (GLint loc, const int* v)           { glUniform2iv(loc,1,v); }
+inline void glUniform2fv_ei  (GLint loc, const float* v)         { glUniform2fv(loc,1,v); }
+inline void glUniform2iv_ei  (GLint loc, const int* v)           { glUniform2iv(loc,1,v); }
 
-static void glUniform3fv_ei  (GLint loc, const float* v)         { glUniform3fv(loc,1,v); }
-static void glUniform3iv_ei  (GLint loc, const int* v)           { glUniform3iv(loc,1,v); }
+inline void glUniform3fv_ei  (GLint loc, const float* v)         { glUniform3fv(loc,1,v); }
+inline void glUniform3iv_ei  (GLint loc, const int* v)           { glUniform3iv(loc,1,v); }
 
-static void glUniform4fv_ei  (GLint loc, const float* v)         { glUniform4fv(loc,1,v); }
-static void glUniform4iv_ei  (GLint loc, const int* v)           { glUniform4iv(loc,1,v); }
+inline void glUniform4fv_ei  (GLint loc, const float* v)         { glUniform4fv(loc,1,v); }
+inline void glUniform4iv_ei  (GLint loc, const int* v)           { glUniform4iv(loc,1,v); }
 
-static void glUniformMatrix2fv_ei  (GLint loc, const float* v)         { glUniformMatrix2fv(loc,1,false,v); }
-static void glUniformMatrix3fv_ei  (GLint loc, const float* v)         { glUniformMatrix3fv(loc,1,false,v); }
-static void glUniformMatrix4fv_ei  (GLint loc, const float* v)         { glUniformMatrix4fv(loc,1,false,v); }
+inline void glUniformMatrix2fv_ei  (GLint loc, const float* v)         { glUniformMatrix2fv(loc,1,false,v); }
+inline void glUniformMatrix3fv_ei  (GLint loc, const float* v)         { glUniformMatrix3fv(loc,1,false,v); }
+inline void glUniformMatrix4fv_ei  (GLint loc, const float* v)         { glUniformMatrix4fv(loc,1,false,v); }
 
 
 EIGEN_GL_FUNC1_DECLARATION       (glUniform,GLint,const)
@@ -294,9 +294,9 @@ EIGEN_GL_FUNC1_SPECIALIZATION_MAT(glUniform,GLint,const,float,        4,3,Matrix
 
 #ifdef GL_VERSION_3_0
 
-static void glUniform2uiv_ei (GLint loc, const unsigned int* v)  { glUniform2uiv(loc,1,v); }
-static void glUniform3uiv_ei (GLint loc, const unsigned int* v)  { glUniform3uiv(loc,1,v); }
-static void glUniform4uiv_ei (GLint loc, const unsigned int* v)  { glUniform4uiv(loc,1,v); }
+inline void glUniform2uiv_ei (GLint loc, const unsigned int* v)  { glUniform2uiv(loc,1,v); }
+inline void glUniform3uiv_ei (GLint loc, const unsigned int* v)  { glUniform3uiv(loc,1,v); }
+inline void glUniform4uiv_ei (GLint loc, const unsigned int* v)  { glUniform4uiv(loc,1,v); }
 
 EIGEN_GL_FUNC1_SPECIALIZATION_VEC(glUniform,GLint,const,unsigned int, 2,2uiv_ei)
 EIGEN_GL_FUNC1_SPECIALIZATION_VEC(glUniform,GLint,const,unsigned int, 3,3uiv_ei)
@@ -305,9 +305,9 @@ EIGEN_GL_FUNC1_SPECIALIZATION_VEC(glUniform,GLint,const,unsigned int, 4,4uiv_ei)
 #endif
 
 #ifdef GL_ARB_gpu_shader_fp64
-static void glUniform2dv_ei  (GLint loc, const double* v)        { glUniform2dv(loc,1,v); }
-static void glUniform3dv_ei  (GLint loc, const double* v)        { glUniform3dv(loc,1,v); }
-static void glUniform4dv_ei  (GLint loc, const double* v)        { glUniform4dv(loc,1,v); }
+inline void glUniform2dv_ei  (GLint loc, const double* v)        { glUniform2dv(loc,1,v); }
+inline void glUniform3dv_ei  (GLint loc, const double* v)        { glUniform3dv(loc,1,v); }
+inline void glUniform4dv_ei  (GLint loc, const double* v)        { glUniform4dv(loc,1,v); }
 
 EIGEN_GL_FUNC1_SPECIALIZATION_VEC(glUniform,GLint,const,double,       2,2dv_ei)
 EIGEN_GL_FUNC1_SPECIALIZATION_VEC(glUniform,GLint,const,double,       3,3dv_ei)
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
index c32437281..78a307e96 100644
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
@@ -110,7 +110,6 @@ void MatrixPowerAtomic<MatrixType>::compute2x2(MatrixType& res, RealScalar p) co
   using std::abs;
   using std::pow;
   
-  ArrayType logTdiag = m_A.diagonal().array().log();
   res.coeffRef(0,0) = pow(m_A.coeff(0,0), p);
 
   for (Index i=1; i < m_A.cols(); ++i) {
diff --git a/unsupported/doc/examples/CMakeLists.txt b/unsupported/doc/examples/CMakeLists.txt
index 978f9afd8..c47646dfc 100644
--- a/unsupported/doc/examples/CMakeLists.txt
+++ b/unsupported/doc/examples/CMakeLists.txt
@@ -10,12 +10,10 @@ FOREACH(example_src ${examples_SRCS})
   if(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO)
     target_link_libraries(example_${example} ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO})
   endif()
-  GET_TARGET_PROPERTY(example_executable
-                      example_${example} LOCATION)
   ADD_CUSTOM_COMMAND(
     TARGET example_${example}
     POST_BUILD
-    COMMAND ${example_executable}
+    COMMAND example_${example}
     ARGS >${CMAKE_CURRENT_BINARY_DIR}/${example}.out
   )
   ADD_DEPENDENCIES(unsupported_examples example_${example})
diff --git a/unsupported/doc/snippets/CMakeLists.txt b/unsupported/doc/snippets/CMakeLists.txt
index 4a4157933..f0c5cc2a8 100644
--- a/unsupported/doc/snippets/CMakeLists.txt
+++ b/unsupported/doc/snippets/CMakeLists.txt
@@ -14,12 +14,10 @@ FOREACH(snippet_src ${snippets_SRCS})
   if(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO)
     target_link_libraries(${compile_snippet_target} ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO})
   endif()
-  GET_TARGET_PROPERTY(compile_snippet_executable
-                      ${compile_snippet_target} LOCATION)
   ADD_CUSTOM_COMMAND(
     TARGET ${compile_snippet_target}
     POST_BUILD
-    COMMAND ${compile_snippet_executable}
+    COMMAND ${compile_snippet_target}
     ARGS >${CMAKE_CURRENT_BINARY_DIR}/${snippet}.out
   )
   ADD_DEPENDENCIES(unsupported_snippets ${compile_snippet_target})
diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt
index a94a3b5e5..2e4cfdb2e 100644
--- a/unsupported/test/CMakeLists.txt
+++ b/unsupported/test/CMakeLists.txt
@@ -29,11 +29,7 @@ ei_add_test(NonLinearOptimization)
 
 ei_add_test(NumericalDiff)
 ei_add_test(autodiff)
-
-if (NOT CMAKE_CXX_COMPILER MATCHES "clang\\+\\+$")
 ei_add_test(BVH)
-endif()
-
 ei_add_test(matrix_exponential)
 ei_add_test(matrix_function)
 ei_add_test(matrix_power)
@@ -73,8 +69,9 @@ if(NOT EIGEN_TEST_NO_OPENGL)
   find_package(GLUT)
   find_package(GLEW)
   if(OPENGL_FOUND AND GLUT_FOUND AND GLEW_FOUND)
+    include_directories(${OPENGL_INCLUDE_DIR} ${GLUT_INCLUDE_DIR} ${GLEW_INCLUDE_DIRS})
     ei_add_property(EIGEN_TESTED_BACKENDS "OpenGL, ")
-    set(EIGEN_GL_LIB ${GLUT_LIBRARIES} ${GLEW_LIBRARIES})
+    set(EIGEN_GL_LIB ${GLUT_LIBRARIES} ${GLEW_LIBRARIES} ${OPENGL_LIBRARIES})
     ei_add_test(openglsupport  "" "${EIGEN_GL_LIB}" )
   else()
     ei_add_property(EIGEN_MISSING_BACKENDS "OpenGL, ")
diff --git a/unsupported/test/mpreal/mpreal.h b/unsupported/test/mpreal/mpreal.h
index ef0a6a9f0..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-2012 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__

@@ -65,11 +53,21 @@
 #include <stdexcept>

 #include <cfloat>

 #include <cmath>

+#include <cstring>

 #include <limits>

 

 // 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)

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

 #endif

 

+// A Clang feature extension to determine compiler features.

+#ifndef __has_feature

+    #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__) || __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

 #if defined(MPREAL_HAVE_INT64_SUPPORT)

     

     #define MPFR_USE_INTMAX_T                   // Should be defined before mpfr.h

@@ -97,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

@@ -111,7 +135,7 @@
 #endif 

 

 #if defined(MPREAL_HAVE_MSVC_DEBUGVIEW) && defined(_MSC_VER) && defined(_DEBUG)

-#define MPREAL_MSVC_DEBUGVIEW_CODE         DebugView = toString();

+    #define MPREAL_MSVC_DEBUGVIEW_CODE     DebugView = toString();

     #define MPREAL_MSVC_DEBUGVIEW_DATA     std::string DebugView;

 #else

     #define MPREAL_MSVC_DEBUGVIEW_CODE 

@@ -149,7 +173,6 @@ public:
     // Constructors && type conversions

     mpreal();

     mpreal(const mpreal& u);

-    mpreal(const mpfr_t u);    

     mpreal(const mpf_t u);    

     mpreal(const mpz_t u,             mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd());    

     mpreal(const mpq_t u,             mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd());    

@@ -159,6 +182,10 @@ public:
     mpreal(const unsigned int u,      mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd());

     mpreal(const long int u,          mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd());

     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' share the same data & pointers.    

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

 

 #if defined (MPREAL_HAVE_INT64_SUPPORT)

     mpreal(const uint64_t u,          mp_prec_t prec = mpreal::get_default_prec(),  mp_rnd_t mode = mpreal::get_default_rnd());

@@ -170,6 +197,11 @@ public:
 

     ~mpreal();                           

 

+#ifdef MPREAL_HAVE_MOVE_SUPPORT

+    mpreal& operator=(mpreal&& v);

+    mpreal(mpreal&& u);

+#endif

+

     // Operations

     // =

     // +, -, *, /, ++, --, <<, >> 

@@ -292,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

@@ -308,121 +360,125 @@ 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 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);    

     friend const mpreal ldexp(const mpreal& v, mp_exp_t exp);

@@ -433,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);

@@ -524,19 +580,22 @@ 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.

     // Put one of these lines:

     //

-    // mpfr::mpreal=<DebugView>                                ; Show value only

+    // mpfr::mpreal=<DebugView>                              ; Show value only

     // mpfr::mpreal=<DebugView>, <mp[0]._mpfr_prec,u>bits    ; Show value & precision

     // 

     // at the beginning of

     // [Visual Studio Installation Folder]\Common7\Packages\Debugger\autoexp.dat

     MPREAL_MSVC_DEBUGVIEW_DATA

+

+    // "Smart" resources deallocation. Checks if instance initialized before deletion.

+    void clear(::mpfr_ptr);

 };

 

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

@@ -551,64 +610,89 @@ 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;

+}

+

+#ifdef MPREAL_HAVE_MOVE_SUPPORT

+inline mpreal::mpreal(mpreal&& other)

+{

+    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;

+}

+

+inline mpreal& mpreal::operator=(mpreal&& other)

+{

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

 

     MPREAL_MSVC_DEBUGVIEW_CODE;

+    return *this;

 }

+#endif

 

-inline mpreal::mpreal(const mpfr_t u)

+inline mpreal::mpreal(const mpfr_t  u, bool shared)

 {

-    mpfr_init2(mp,mpfr_get_prec(u));

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

+    if(shared)

+    {

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

+    }

+    else

+    {

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

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

+    }

 

     MPREAL_MSVC_DEBUGVIEW_CODE;

 }

 

 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;

@@ -616,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;

 }

@@ -657,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;

 }

@@ -674,23 +758,31 @@ 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)

+{

+#ifdef MPREAL_HAVE_MOVE_SUPPORT

+    if(mpfr_is_initialized(x)) 

+#endif

+    mpfr_clear(x);

+}

+

 inline mpreal::~mpreal() 

 { 

-    mpfr_clear(mp);

+    clear(mpfr_ptr());

 }                           

 

 // internal namespace needed for template magic

@@ -761,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());

@@ -813,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

@@ -860,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){

-			mpfr_clear(mp);

-			mpfr_init2(mp, 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;

     }

@@ -877,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;

@@ -885,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;

@@ -893,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;

@@ -901,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;

@@ -910,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;

@@ -933,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;

@@ -941,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;

@@ -949,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;

@@ -966,15 +1066,15 @@ 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;

     }

 

-    mpfr_clear(t);

+    clear(t);

     return *this;

 }

 

@@ -989,15 +1089,15 @@ 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;

     }

 

-    mpfr_clear(t);

+    clear(t);

     return *this;

 }

 

@@ -1006,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;

 }

@@ -1020,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;

 }

@@ -1042,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

@@ -1053,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;

 }

@@ -1094,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++() 

@@ -1127,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;

 }

@@ -1156,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

@@ -1167,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;

 }

@@ -1196,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)

@@ -1252,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;

 }

@@ -1281,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

@@ -1291,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;

 }

@@ -1357,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

@@ -1367,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()));

@@ -1445,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;

 }

@@ -1543,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;

 }

 

@@ -1551,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; }

@@ -1629,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);

 

@@ -1644,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

 

@@ -1675,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)

     {

@@ -1761,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;

 }

 

@@ -1784,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;

 }

@@ -1839,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 

 

@@ -1851,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;

 }

@@ -1885,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;

 }

 

@@ -1900,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;

     }

 }

 

@@ -1922,37 +2049,37 @@ inline mpreal maxval(mp_prec_t prec)
 

 inline bool isEqualUlps(const mpreal& a, const mpreal& b, int maxUlps)

 {

-  return abs(a - b) <= machine_epsilon((max)(abs(a), abs(b))) * maxUlps;

+    return abs(a - b) <= machine_epsilon((max)(abs(a), abs(b))) * maxUlps;

 }

 

 inline bool isEqualFuzzy(const mpreal& a, const mpreal& b, const mpreal& eps)

 {

-    return abs(a - b) <= (min)(abs(a), abs(b)) * eps;

+    return abs(a - b) <= eps;

 }

 

 inline bool isEqualFuzzy(const mpreal& a, const mpreal& b)

 {

-    return isEqualFuzzy(a, b, machine_epsilon((min)(abs(a), abs(b))));

+    return isEqualFuzzy(a, b, machine_epsilon((max)(1, (min)(abs(a), abs(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;

 }

@@ -2005,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)

 {

@@ -2032,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);

@@ -2048,161 +2178,130 @@ inline const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t r)
 

 inline int cmpabs(const mpreal& a,const mpreal& b)

 {

-    return mpfr_cmpabs(a.mp,b.mp);

+    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.mp,c.mp,v.mp,rnd_mode);

+    return mpfr_sin_cos(s.mpfr_ptr(), c.mpfr_ptr(), v.mpfr_srcptr(), 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)

-{

-    mpreal a;

-    mp_prec_t yp, xp;

-

-    yp = y.get_prec(); 

-    xp = x.get_prec(); 

-

-    a.set_prec(yp>xp?yp:xp);

-

-    mpfr_atan2(a.mp, y.mp, x.mp, 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;

-    mp_prec_t yp, xp;

-

-    yp = y.get_prec(); 

-    xp = x.get_prec(); 

-

-    a.set_prec(yp>xp?yp:xp);

-

-    mpfr_hypot(a.mp, x.mp, y.mp, rnd_mode);

-

+    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;

-    mp_prec_t yp, xp;

-

-    yp = y.get_prec(); 

-    xp = x.get_prec(); 

-

-    a.set_prec(yp>xp?yp:xp);

-

-    mpfr_remainder(a.mp, x.mp, y.mp, rnd_mode);

-

+    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;

-    mp_prec_t yp, xp;

-

-    yp = y.get_prec(); 

-    xp = x.get_prec(); 

-

-    a.set_prec(yp>xp?yp:xp);

-

-    mpfr_remquo(a.mp,q, x.mp, y.mp, rnd_mode);

-

+    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.mp,v,rnd_mode);

+    mpfr_fac_ui(x.mpfr_ptr(),v,rnd_mode);

     return x;

 }

 

 

-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;

@@ -2217,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;

@@ -2232,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;

@@ -2247,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;

@@ -2264,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;

     

@@ -2305,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;

@@ -2320,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);

@@ -2331,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);

@@ -2402,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

@@ -2415,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);

@@ -2457,14 +2556,22 @@ 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 = mpreal::get_default_rnd())

+{

+    mpreal x;

+    mpfr_grandom(x.mp, NULL, state, rnd_mode);

+    return x;

+}

+

 #endif 

 

 #if (MPFR_VERSION <= MPFR_VERSION_NUM(2,4,2))

@@ -2480,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))

@@ -2504,6 +2611,25 @@ inline const mpreal random(unsigned int seed)
 

 }

 

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

+inline const mpreal grandom(unsigned int seed = 0)

+{

+    static gmp_randstate_t state;

+    static bool isFirstTime = true;

+

+    if(isFirstTime)

+    {

+        gmp_randinit_default(state);

+        gmp_randseed_ui(state,0);

+        isFirstTime = false;

+    }

+

+    if(seed != 0) gmp_randseed_ui(state,seed);

+

+    return mpfr::grandom(state);

+}

+#endif

+

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

 // Set/Get global properties

 inline void mpreal::set_default_prec(mp_prec_t prec)

@@ -2523,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);

@@ -2549,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);

@@ -2571,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);

@@ -2586,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)

@@ -2621,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)

@@ -2824,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) 

     { 

@@ -2852,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 MPFR_RNDN: return round_to_nearest;

-                case MPFR_RNDZ: return round_toward_zero; 

-                case MPFR_RNDU: return round_toward_infinity; 

-                case MPFR_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);  }

@@ -2873,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);  }

 

@@ -2881,7 +2993,7 @@ namespace std
         {

             mp_rnd_t r = mpfr::mpreal::get_default_rnd();

 

-            if(r == MPFR_RNDN) return mpfr::mpreal(0.5, precision); 

+            if(r == GMP_RNDN)  return mpfr::mpreal(0.5, precision); 

             else               return mpfr::mpreal(1.0, precision);    

         }

 

@@ -2896,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())

         {

@@ -2915,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/polynomialsolver.cpp b/unsupported/test/polynomialsolver.cpp
index 13f92169e..de79f1538 100644
--- a/unsupported/test/polynomialsolver.cpp
+++ b/unsupported/test/polynomialsolver.cpp
@@ -104,9 +104,7 @@ void evalSolverSugarFunction( const POLYNOMIAL& pols, const ROOTS& roots, const
     // 1) the roots found are correct
     // 2) the roots have distinct moduli
 
-    typedef typename POLYNOMIAL::Scalar                 Scalar;
     typedef typename REAL_ROOTS::Scalar                 Real;
-    typedef PolynomialSolver<Scalar, Deg >              PolynomialSolverType;
 
     //Test realRoots
     std::vector< Real > calc_realRoots;