diff options
| -rw-r--r-- | CMakeLists.txt | 7 | ||||
| -rw-r--r-- | Eigen/src/Cholesky/LDLT.h | 8 | ||||
| -rw-r--r-- | Eigen/src/Core/PermutationMatrix.h | 60 | ||||
| -rw-r--r-- | Eigen/src/Core/Transpositions.h | 77 | ||||
| -rwxr-xr-x | Eigen/src/Core/arch/AltiVec/PacketMath.h | 26 | ||||
| -rw-r--r-- | Eigen/src/Core/util/Meta.h | 23 | ||||
| -rw-r--r-- | Eigen/src/LU/PartialPivLU.h | 6 | ||||
| -rw-r--r-- | Eigen/src/SVD/JacobiSVD.h | 12 | ||||
| -rw-r--r-- | Eigen/src/SparseCore/SparseAssign.h | 16 | ||||
| -rw-r--r-- | Eigen/src/SparseCore/SparseCwiseBinaryOp.h | 11 | ||||
| -rw-r--r-- | Eigen/src/SparseCore/SparseDenseProduct.h | 18 | ||||
| -rw-r--r-- | Eigen/src/SparseCore/SparseDiagonalProduct.h | 14 | ||||
| -rw-r--r-- | Eigen/src/SparseCore/SparseMatrix.h | 6 | ||||
| -rw-r--r-- | Eigen/src/SparseCore/SparseMatrixBase.h | 2 | ||||
| -rw-r--r-- | Eigen/src/SparseCore/SparseSelfAdjointView.h | 2 | ||||
| -rw-r--r-- | test/jacobisvd.cpp | 88 | ||||
| -rw-r--r-- | test/meta.cpp | 24 | ||||
| -rw-r--r-- | test/product_large.cpp | 9 | ||||
| -rw-r--r-- | test/product_trsolve.cpp | 18 | ||||
| -rw-r--r-- | test/sparse_product.cpp | 43 | ||||
| -rw-r--r-- | unsupported/Eigen/MPRealSupport | 48 | ||||
| -rw-r--r-- | unsupported/test/mpreal/mpreal.h | 1056 | ||||
| -rw-r--r-- | unsupported/test/mpreal_support.cpp | 3 |
23 files changed, 899 insertions, 678 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index a719e47fd..470095680 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -108,7 +108,8 @@ endif() set(EIGEN_TEST_MAX_SIZE "320" CACHE STRING "Maximal matrix/vector size, default is 320") macro(ei_add_cxx_compiler_flag FLAG) - string(REGEX REPLACE "-" "" SFLAG ${FLAG}) + string(REGEX REPLACE "-" "" SFLAG1 ${FLAG}) + string(REGEX REPLACE "\\+" "p" SFLAG ${SFLAG1}) check_cxx_compiler_flag(${FLAG} COMPILER_SUPPORT_${SFLAG}) if(COMPILER_SUPPORT_${SFLAG}) set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${FLAG}") @@ -142,6 +143,9 @@ if(NOT MSVC) ei_add_cxx_compiler_flag("-Wpointer-arith") ei_add_cxx_compiler_flag("-Wwrite-strings") ei_add_cxx_compiler_flag("-Wformat-security") + ei_add_cxx_compiler_flag("-Wshorten-64-to-32") + ei_add_cxx_compiler_flag("-Wenum-conversion") + ei_add_cxx_compiler_flag("-Wc++11-extensions") ei_add_cxx_compiler_flag("-Wno-psabi") ei_add_cxx_compiler_flag("-Wno-variadic-macros") @@ -153,6 +157,7 @@ if(NOT MSVC) ei_add_cxx_compiler_flag("-wd981") # disable ICC's "operands are evaluated in unspecified order" remark ei_add_cxx_compiler_flag("-wd2304") # disbale ICC's "warning #2304: non-explicit constructor with single argument may cause implicit type conversion" produced by -Wnon-virtual-dtor + # The -ansi flag must be added last, otherwise it is also used as a linker flag by check_cxx_compiler_flag making it fails # Moreover we should not set both -strict-ansi and -ansi check_cxx_compiler_flag("-strict-ansi" COMPILER_SUPPORT_STRICTANSI) diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index 7415b826a..f621323a3 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -264,6 +264,7 @@ template<> struct ldlt_inplace<Lower> typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; typedef typename MatrixType::Index Index; + typedef typename TranspositionType::StorageIndexType IndexType; eigen_assert(mat.rows()==mat.cols()); const Index size = mat.rows(); @@ -283,7 +284,7 @@ template<> struct ldlt_inplace<Lower> mat.diagonal().tail(size-k).cwiseAbs().maxCoeff(&index_of_biggest_in_corner); index_of_biggest_in_corner += k; - transpositions.coeffRef(k) = index_of_biggest_in_corner; + transpositions.coeffRef(k) = IndexType(index_of_biggest_in_corner); if(k != index_of_biggest_in_corner) { // apply the transposition while taking care to consider only @@ -292,7 +293,7 @@ template<> struct ldlt_inplace<Lower> mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k)); mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s)); std::swap(mat.coeffRef(k,k),mat.coeffRef(index_of_biggest_in_corner,index_of_biggest_in_corner)); - for(int i=k+1;i<index_of_biggest_in_corner;++i) + for(Index i=k+1;i<index_of_biggest_in_corner;++i) { Scalar tmp = mat.coeffRef(i,k); mat.coeffRef(i,k) = numext::conj(mat.coeffRef(index_of_biggest_in_corner,i)); @@ -460,6 +461,7 @@ template<typename MatrixType, int _UpLo> template<typename Derived> LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::rankUpdate(const MatrixBase<Derived>& w, const typename NumTraits<typename MatrixType::Scalar>::Real& sigma) { + typedef typename TranspositionType::StorageIndexType IndexType; const Index size = w.rows(); if (m_isInitialized) { @@ -471,7 +473,7 @@ LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::rankUpdate(const MatrixBase<Deri m_matrix.setZero(); m_transpositions.resize(size); for (Index i = 0; i < size; i++) - m_transpositions.coeffRef(i) = i; + m_transpositions.coeffRef(i) = IndexType(i); m_temporary.resize(size); m_sign = sigma>=0 ? internal::PositiveSemiDef : internal::NegativeSemiDef; m_isInitialized = true; diff --git a/Eigen/src/Core/PermutationMatrix.h b/Eigen/src/Core/PermutationMatrix.h index 658473397..5d648c68c 100644 --- a/Eigen/src/Core/PermutationMatrix.h +++ b/Eigen/src/Core/PermutationMatrix.h @@ -68,11 +68,11 @@ class PermutationBase : public EigenBase<Derived> MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime, MaxColsAtCompileTime = Traits::MaxColsAtCompileTime }; - typedef typename Traits::Scalar Scalar; + typedef typename Traits::StorageIndexType StorageIndexType; typedef typename Traits::Index Index; - typedef Matrix<Scalar,RowsAtCompileTime,ColsAtCompileTime,0,MaxRowsAtCompileTime,MaxColsAtCompileTime> + typedef Matrix<StorageIndexType,RowsAtCompileTime,ColsAtCompileTime,0,MaxRowsAtCompileTime,MaxColsAtCompileTime> DenseMatrixType; - typedef PermutationMatrix<IndicesType::SizeAtCompileTime,IndicesType::MaxSizeAtCompileTime,Index> + typedef PermutationMatrix<IndicesType::SizeAtCompileTime,IndicesType::MaxSizeAtCompileTime,StorageIndexType> PlainPermutationType; using Base::derived; #endif @@ -149,7 +149,7 @@ class PermutationBase : public EigenBase<Derived> /** Sets *this to be the identity permutation matrix */ void setIdentity() { - for(Index i = 0; i < size(); ++i) + for(StorageIndexType i = 0; i < size(); ++i) indices().coeffRef(i) = i; } @@ -175,8 +175,8 @@ class PermutationBase : public EigenBase<Derived> eigen_assert(i>=0 && j>=0 && i<size() && j<size()); for(Index k = 0; k < size(); ++k) { - if(indices().coeff(k) == i) indices().coeffRef(k) = j; - else if(indices().coeff(k) == j) indices().coeffRef(k) = i; + if(indices().coeff(k) == i) indices().coeffRef(k) = StorageIndexType(j); + else if(indices().coeff(k) == j) indices().coeffRef(k) = StorageIndexType(i); } return derived(); } @@ -264,7 +264,7 @@ class PermutationBase : public EigenBase<Derived> * * \param SizeAtCompileTime the number of rows/cols, or Dynamic * \param MaxSizeAtCompileTime the maximum number of rows/cols, or Dynamic. This optional parameter defaults to SizeAtCompileTime. Most of the time, you should not have to specify it. - * \param IndexType the interger type of the indices + * \param StorageIndexType the interger type of the indices * * This class represents a permutation matrix, internally stored as a vector of integers. * @@ -272,17 +272,18 @@ class PermutationBase : public EigenBase<Derived> */ namespace internal { -template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType> -struct traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType> > - : traits<Matrix<IndexType,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> > +template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndexType> +struct traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndexType> > + : traits<Matrix<_StorageIndexType,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> > { - typedef IndexType Index; - typedef Matrix<IndexType, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType; + typedef Matrix<_StorageIndexType, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType; + typedef typename IndicesType::Index Index; + typedef _StorageIndexType StorageIndexType; }; } -template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType> -class PermutationMatrix : public PermutationBase<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType> > +template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndexType> +class PermutationMatrix : public PermutationBase<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndexType> > { typedef PermutationBase<PermutationMatrix> Base; typedef internal::traits<PermutationMatrix> Traits; @@ -294,6 +295,8 @@ class PermutationMatrix : public PermutationBase<PermutationMatrix<SizeAtCompile #ifndef EIGEN_PARSED_BY_DOXYGEN typedef typename Traits::IndicesType IndicesType; + typedef typename Traits::StorageIndexType StorageIndexType; + typedef typename Traits::Index Index; #endif inline PermutationMatrix() @@ -301,7 +304,7 @@ class PermutationMatrix : public PermutationBase<PermutationMatrix<SizeAtCompile /** Constructs an uninitialized permutation matrix of given size. */ - inline PermutationMatrix(int size) : m_indices(size) + inline PermutationMatrix(Index size) : m_indices(size) {} /** Copy constructor. */ @@ -390,18 +393,19 @@ class PermutationMatrix : public PermutationBase<PermutationMatrix<SizeAtCompile namespace internal { -template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess> -struct traits<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess> > - : traits<Matrix<IndexType,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> > +template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndexType, int _PacketAccess> +struct traits<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndexType>,_PacketAccess> > + : traits<Matrix<_StorageIndexType,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> > { - typedef IndexType Index; - typedef Map<const Matrix<IndexType, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1>, _PacketAccess> IndicesType; + typedef Map<const Matrix<_StorageIndexType, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1>, _PacketAccess> IndicesType; + typedef typename IndicesType::Index Index; + typedef _StorageIndexType StorageIndexType; }; } -template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess> -class Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess> - : public PermutationBase<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess> > +template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndexType, int _PacketAccess> +class Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndexType>,_PacketAccess> + : public PermutationBase<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndexType>,_PacketAccess> > { typedef PermutationBase<Map> Base; typedef internal::traits<Map> Traits; @@ -409,14 +413,15 @@ class Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>, #ifndef EIGEN_PARSED_BY_DOXYGEN typedef typename Traits::IndicesType IndicesType; - typedef typename IndicesType::Scalar Index; + typedef typename IndicesType::Scalar StorageIndexType; + typedef typename IndicesType::Index Index; #endif - inline Map(const Index* indicesPtr) + inline Map(const StorageIndexType* indicesPtr) : m_indices(indicesPtr) {} - inline Map(const Index* indicesPtr, Index size) + inline Map(const StorageIndexType* indicesPtr, Index size) : m_indices(indicesPtr,size) {} @@ -472,7 +477,8 @@ struct traits<PermutationWrapper<_IndicesType> > { typedef PermutationStorage StorageKind; typedef typename _IndicesType::Scalar Scalar; - typedef typename _IndicesType::Scalar Index; + typedef typename _IndicesType::Scalar StorageIndexType; + typedef typename _IndicesType::Index Index; typedef _IndicesType IndicesType; enum { RowsAtCompileTime = _IndicesType::SizeAtCompileTime, diff --git a/Eigen/src/Core/Transpositions.h b/Eigen/src/Core/Transpositions.h index e4ba0756f..92261118f 100644 --- a/Eigen/src/Core/Transpositions.h +++ b/Eigen/src/Core/Transpositions.h @@ -53,7 +53,8 @@ class TranspositionsBase public: typedef typename Traits::IndicesType IndicesType; - typedef typename IndicesType::Scalar Index; + typedef typename IndicesType::Scalar StorageIndexType; + typedef typename IndicesType::Index Index; Derived& derived() { return *static_cast<Derived*>(this); } const Derived& derived() const { return *static_cast<const Derived*>(this); } @@ -81,17 +82,17 @@ class TranspositionsBase inline Index size() const { return indices().size(); } /** Direct access to the underlying index vector */ - inline const Index& coeff(Index i) const { return indices().coeff(i); } + inline const StorageIndexType& coeff(Index i) const { return indices().coeff(i); } /** Direct access to the underlying index vector */ - inline Index& coeffRef(Index i) { return indices().coeffRef(i); } + inline StorageIndexType& coeffRef(Index i) { return indices().coeffRef(i); } /** Direct access to the underlying index vector */ - inline const Index& operator()(Index i) const { return indices()(i); } + inline const StorageIndexType& operator()(Index i) const { return indices()(i); } /** Direct access to the underlying index vector */ - inline Index& operator()(Index i) { return indices()(i); } + inline StorageIndexType& operator()(Index i) { return indices()(i); } /** Direct access to the underlying index vector */ - inline const Index& operator[](Index i) const { return indices()(i); } + inline const StorageIndexType& operator[](Index i) const { return indices()(i); } /** Direct access to the underlying index vector */ - inline Index& operator[](Index i) { return indices()(i); } + inline StorageIndexType& operator[](Index i) { return indices()(i); } /** const version of indices(). */ const IndicesType& indices() const { return derived().indices(); } @@ -99,7 +100,7 @@ class TranspositionsBase IndicesType& indices() { return derived().indices(); } /** Resizes to given size. */ - inline void resize(int newSize) + inline void resize(Index newSize) { indices().resize(newSize); } @@ -107,7 +108,7 @@ class TranspositionsBase /** Sets \c *this to represents an identity transformation */ void setIdentity() { - for(int i = 0; i < indices().size(); ++i) + for(StorageIndexType i = 0; i < indices().size(); ++i) coeffRef(i) = i; } @@ -144,23 +145,26 @@ class TranspositionsBase }; namespace internal { -template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType> -struct traits<Transpositions<SizeAtCompileTime,MaxSizeAtCompileTime,IndexType> > +template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndexType> +struct traits<Transpositions<SizeAtCompileTime,MaxSizeAtCompileTime,_StorageIndexType> > { - typedef IndexType Index; - typedef Matrix<Index, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType; + typedef Matrix<_StorageIndexType, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType; + typedef typename IndicesType::Index Index; + typedef _StorageIndexType StorageIndexType; }; } -template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType> -class Transpositions : public TranspositionsBase<Transpositions<SizeAtCompileTime,MaxSizeAtCompileTime,IndexType> > +template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndexType> +class Transpositions : public TranspositionsBase<Transpositions<SizeAtCompileTime,MaxSizeAtCompileTime,_StorageIndexType> > { typedef internal::traits<Transpositions> Traits; public: typedef TranspositionsBase<Transpositions> Base; typedef typename Traits::IndicesType IndicesType; - typedef typename IndicesType::Scalar Index; + typedef typename IndicesType::Scalar StorageIndexType; + typedef typename IndicesType::Index Index; + inline Transpositions() {} @@ -215,30 +219,32 @@ class Transpositions : public TranspositionsBase<Transpositions<SizeAtCompileTim namespace internal { -template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess> -struct traits<Map<Transpositions<SizeAtCompileTime,MaxSizeAtCompileTime,IndexType>,_PacketAccess> > +template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndexType, int _PacketAccess> +struct traits<Map<Transpositions<SizeAtCompileTime,MaxSizeAtCompileTime,_StorageIndexType>,_PacketAccess> > { - typedef IndexType Index; - typedef Map<const Matrix<Index,SizeAtCompileTime,1,0,MaxSizeAtCompileTime,1>, _PacketAccess> IndicesType; + typedef Map<const Matrix<_StorageIndexType,SizeAtCompileTime,1,0,MaxSizeAtCompileTime,1>, _PacketAccess> IndicesType; + typedef typename IndicesType::Index Index; + typedef _StorageIndexType StorageIndexType; }; } -template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int PacketAccess> -class Map<Transpositions<SizeAtCompileTime,MaxSizeAtCompileTime,IndexType>,PacketAccess> - : public TranspositionsBase<Map<Transpositions<SizeAtCompileTime,MaxSizeAtCompileTime,IndexType>,PacketAccess> > +template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndexType, int PacketAccess> +class Map<Transpositions<SizeAtCompileTime,MaxSizeAtCompileTime,_StorageIndexType>,PacketAccess> + : public TranspositionsBase<Map<Transpositions<SizeAtCompileTime,MaxSizeAtCompileTime,_StorageIndexType>,PacketAccess> > { typedef internal::traits<Map> Traits; public: typedef TranspositionsBase<Map> Base; typedef typename Traits::IndicesType IndicesType; - typedef typename IndicesType::Scalar Index; + typedef typename IndicesType::Scalar StorageIndexType; + typedef typename IndicesType::Index Index; - inline Map(const Index* indicesPtr) + inline Map(const StorageIndexType* indicesPtr) : m_indices(indicesPtr) {} - inline Map(const Index* indicesPtr, Index size) + inline Map(const StorageIndexType* indicesPtr, Index size) : m_indices(indicesPtr,size) {} @@ -275,7 +281,8 @@ namespace internal { template<typename _IndicesType> struct traits<TranspositionsWrapper<_IndicesType> > { - typedef typename _IndicesType::Scalar Index; + typedef typename _IndicesType::Scalar StorageIndexType; + typedef typename _IndicesType::Index Index; typedef _IndicesType IndicesType; }; } @@ -289,7 +296,8 @@ class TranspositionsWrapper typedef TranspositionsBase<TranspositionsWrapper> Base; typedef typename Traits::IndicesType IndicesType; - typedef typename IndicesType::Scalar Index; + typedef typename IndicesType::Scalar StorageIndexType; + typedef typename IndicesType::Index Index; inline TranspositionsWrapper(IndicesType& a_indices) : m_indices(a_indices) @@ -363,24 +371,25 @@ struct transposition_matrix_product_retval { typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned; typedef typename TranspositionType::Index Index; + typedef typename TranspositionType::StorageIndexType StorageIndexType; transposition_matrix_product_retval(const TranspositionType& tr, const MatrixType& matrix) : m_transpositions(tr), m_matrix(matrix) {} - inline int rows() const { return m_matrix.rows(); } - inline int cols() const { return m_matrix.cols(); } + inline Index rows() const { return m_matrix.rows(); } + inline Index cols() const { return m_matrix.cols(); } template<typename Dest> inline void evalTo(Dest& dst) const { - const int size = m_transpositions.size(); - Index j = 0; + const Index size = m_transpositions.size(); + StorageIndexType j = 0; if(!(is_same<MatrixTypeNestedCleaned,Dest>::value && extract_data(dst) == extract_data(m_matrix))) dst = m_matrix; - for(int k=(Transposed?size-1:0) ; Transposed?k>=0:k<size ; Transposed?--k:++k) - if((j=m_transpositions.coeff(k))!=k) + for(Index k=(Transposed?size-1:0) ; Transposed?k>=0:k<size ; Transposed?--k:++k) + if(Index(j=m_transpositions.coeff(k))!=k) { if(Side==OnTheLeft) dst.row(k).swap(dst.row(j)); diff --git a/Eigen/src/Core/arch/AltiVec/PacketMath.h b/Eigen/src/Core/arch/AltiVec/PacketMath.h index aadfc17be..b43e8ace3 100755 --- a/Eigen/src/Core/arch/AltiVec/PacketMath.h +++ b/Eigen/src/Core/arch/AltiVec/PacketMath.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008 Konstantinos Margaritis <markos@codex.gr> +// Copyright (C) 2008-2014 Konstantinos Margaritis <markos@freevec.org> // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed @@ -62,17 +62,17 @@ typedef __vector unsigned char Packet16uc; // Define global static constants: static Packet4f p4f_COUNTDOWN = { 0.0, 1.0, 2.0, 3.0 }; static Packet4i p4i_COUNTDOWN = { 0, 1, 2, 3 }; -static Packet16uc p16uc_REVERSE = {12,13,14,15, 8,9,10,11, 4,5,6,7, 0,1,2,3}; -static Packet16uc p16uc_FORWARD = vec_lvsl(0, (float*)0); -static Packet16uc p16uc_DUPLICATE = {0,1,2,3, 0,1,2,3, 4,5,6,7, 4,5,6,7}; - -static _EIGEN_DECLARE_CONST_FAST_Packet4f(ZERO, 0); -static _EIGEN_DECLARE_CONST_FAST_Packet4i(ZERO, 0); -static _EIGEN_DECLARE_CONST_FAST_Packet4i(ONE,1); -static _EIGEN_DECLARE_CONST_FAST_Packet4i(MINUS16,-16); -static _EIGEN_DECLARE_CONST_FAST_Packet4i(MINUS1,-1); -static Packet4f p4f_ONE = vec_ctf(p4i_ONE, 0); -static Packet4f p4f_ZERO_ = (Packet4f) vec_sl((Packet4ui)p4i_MINUS1, (Packet4ui)p4i_MINUS1); +static Packet16uc p16uc_REVERSE = { 12,13,14,15, 8,9,10,11, 4,5,6,7, 0,1,2,3}; +static Packet16uc p16uc_FORWARD = vec_lvsl(0, (float*)0); //{ 0,1,2,3, 4,5,6,7, 8,9,10,11, 12,13,14,15} +static Packet16uc p16uc_DUPLICATE = { 0,1,2,3, 0,1,2,3, 4,5,6,7, 4,5,6,7}; + +static _EIGEN_DECLARE_CONST_FAST_Packet4f(ZERO, 0); //{ 0.0, 0.0, 0.0, 0.0} +static _EIGEN_DECLARE_CONST_FAST_Packet4i(ZERO, 0); //{ 0, 0, 0, 0,} +static _EIGEN_DECLARE_CONST_FAST_Packet4i(ONE,1); //{ 1, 1, 1, 1} +static _EIGEN_DECLARE_CONST_FAST_Packet4i(MINUS16,-16); //{ -16, -16, -16, -16} +static _EIGEN_DECLARE_CONST_FAST_Packet4i(MINUS1,-1); //{ -1, -1, -1, -1} +static Packet4f p4f_ONE = vec_ctf(p4i_ONE, 0); //{ 1.0, 1.0, 1.0, 1.0} +static Packet4f p4f_ZERO_ = (Packet4f) vec_sl((Packet4ui)p4i_MINUS1, (Packet4ui)p4i_MINUS1); //{ 0x80000000, 0x80000000, 0x80000000, 0x80000000} template<> struct packet_traits<float> : default_packet_traits { @@ -82,8 +82,10 @@ template<> struct packet_traits<float> : default_packet_traits Vectorizable = 1, AlignedOnScalar = 1, size=4, + HasHalfPacket=0, // FIXME check the Has* + HasDiv = 1, HasSin = 0, HasCos = 0, HasLog = 0, diff --git a/Eigen/src/Core/util/Meta.h b/Eigen/src/Core/util/Meta.h index c578b6c86..f3bafd5af 100644 --- a/Eigen/src/Core/util/Meta.h +++ b/Eigen/src/Core/util/Meta.h @@ -82,22 +82,31 @@ template<typename T> struct add_const_on_value_type<T const* const> { typedef T template<typename From, typename To> -struct is_convertible +struct is_convertible_impl { private: + struct any_conversion + { + template <typename T> any_conversion(const volatile T&); + template <typename T> any_conversion(T&); + }; struct yes {int a[1];}; struct no {int a[2];}; - - template<typename T> - static yes test (const T&) {} - - template<typename> static no test (...) {} + + static yes test(const To&, int); + static no test(any_conversion, ...); public: static From ms_from; - enum { value = sizeof(test<To>(ms_from))==sizeof(yes) }; + enum { value = sizeof(test(ms_from, 0))==sizeof(yes) }; }; +template<typename From, typename To> +struct is_convertible +{ + enum { value = is_convertible_impl<typename remove_all<From>::type, + typename remove_all<To >::type>::value }; +}; /** \internal Allows to enable/disable an overload * according to a compile time condition. diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index 02db81e9a..57076f3a4 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -427,13 +427,13 @@ struct partial_lu_impl /** \internal performs the LU decomposition with partial pivoting in-place. */ template<typename MatrixType, typename TranspositionType> -void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, typename TranspositionType::Index& nb_transpositions) +void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, typename TranspositionType::StorageIndexType& nb_transpositions) { eigen_assert(lu.cols() == row_transpositions.size()); eigen_assert((&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1); partial_lu_impl - <typename MatrixType::Scalar, MatrixType::Flags&RowMajorBit?RowMajor:ColMajor, typename TranspositionType::Index> + <typename MatrixType::Scalar, MatrixType::Flags&RowMajorBit?RowMajor:ColMajor, typename TranspositionType::StorageIndexType> ::blocked_lu(lu.rows(), lu.cols(), &lu.coeffRef(0,0), lu.outerStride(), &row_transpositions.coeffRef(0), nb_transpositions); } @@ -452,7 +452,7 @@ PartialPivLU<MatrixType>& PartialPivLU<MatrixType>::compute(const MatrixType& ma m_rowsTranspositions.resize(size); - typename TranspositionType::Index nb_transpositions; + typename TranspositionType::StorageIndexType nb_transpositions; internal::partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions); m_det_p = (nb_transpositions%2) ? -1 : 1; diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index a0c0c9172..d1b63b607 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -375,17 +375,19 @@ struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true> Scalar z; JacobiRotation<Scalar> rot; RealScalar n = sqrt(numext::abs2(work_matrix.coeff(p,p)) + numext::abs2(work_matrix.coeff(q,p))); + if(n==0) { z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); work_matrix.row(p) *= z; if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z); if(work_matrix.coeff(q,q)!=Scalar(0)) + { z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); - else - z = Scalar(0); - work_matrix.row(q) *= z; - if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z); + work_matrix.row(q) *= z; + if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z); + } + // otherwise the second row is already zero, so we have nothing to do. } else { @@ -852,7 +854,7 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize); } - // Scaling factor to reducover/under-flows + // Scaling factor to reduce over/under-flows RealScalar scale = m_workMatrix.cwiseAbs().maxCoeff(); if(scale==RealScalar(0)) scale = RealScalar(1); m_workMatrix /= scale; diff --git a/Eigen/src/SparseCore/SparseAssign.h b/Eigen/src/SparseCore/SparseAssign.h index c99f4b74e..066bc2de1 100644 --- a/Eigen/src/SparseCore/SparseAssign.h +++ b/Eigen/src/SparseCore/SparseAssign.h @@ -51,20 +51,20 @@ template<typename OtherDerived> inline Derived& SparseMatrixBase<Derived>::assign(const OtherDerived& other) { const bool transpose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit); - const Index outerSize = (int(OtherDerived::Flags) & RowMajorBit) ? other.rows() : other.cols(); + const Index outerSize = (int(OtherDerived::Flags) & RowMajorBit) ? Index(other.rows()) : Index(other.cols()); if ((!transpose) && other.isRValue()) { // eval without temporary - derived().resize(other.rows(), other.cols()); + derived().resize(Index(other.rows()), Index(other.cols())); derived().setZero(); derived().reserve((std::max)(this->rows(),this->cols())*2); for (Index j=0; j<outerSize; ++j) { derived().startVec(j); - for (typename OtherDerived::InnerIterator it(other, j); it; ++it) + for (typename OtherDerived::InnerIterator it(other, typename OtherDerived::Index(j)); it; ++it) { Scalar v = it.value(); - derived().insertBackByOuterInner(j,it.index()) = v; + derived().insertBackByOuterInner(j,Index(it.index())) = v; } } derived().finalize(); @@ -87,19 +87,19 @@ inline void SparseMatrixBase<Derived>::assignGeneric(const OtherDerived& other) enum { Flip = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit) }; - const Index outerSize = other.outerSize(); + const Index outerSize = Index(other.outerSize()); //typedef typename internal::conditional<transpose, LinkedVectorMatrix<Scalar,Flags&RowMajorBit>, Derived>::type TempType; // thanks to shallow copies, we always eval to a tempary - Derived temp(other.rows(), other.cols()); + Derived temp(Index(other.rows()), Index(other.cols())); temp.reserve((std::max)(this->rows(),this->cols())*2); for (Index j=0; j<outerSize; ++j) { temp.startVec(j); - for (typename OtherDerived::InnerIterator it(other.derived(), j); it; ++it) + for (typename OtherDerived::InnerIterator it(other.derived(), typename OtherDerived::Index(j)); it; ++it) { Scalar v = it.value(); - temp.insertBackByOuterInner(Flip?it.index():j,Flip?j:it.index()) = v; + temp.insertBackByOuterInner(Flip?Index(it.index()):j,Flip?j:Index(it.index())) = v; } } temp.finalize(); diff --git a/Eigen/src/SparseCore/SparseCwiseBinaryOp.h b/Eigen/src/SparseCore/SparseCwiseBinaryOp.h index 68ea6cb11..3a7e72cd2 100644 --- a/Eigen/src/SparseCore/SparseCwiseBinaryOp.h +++ b/Eigen/src/SparseCore/SparseCwiseBinaryOp.h @@ -65,7 +65,6 @@ class CwiseBinaryOpImpl<BinaryOp,Lhs,Rhs,Sparse>::InnerIterator : public internal::sparse_cwise_binary_op_inner_iterator_selector<BinaryOp,Lhs,Rhs,typename CwiseBinaryOpImpl<BinaryOp,Lhs,Rhs,Sparse>::InnerIterator> { public: - typedef typename Lhs::Index Index; typedef internal::sparse_cwise_binary_op_inner_iterator_selector< BinaryOp,Lhs,Rhs, InnerIterator> Base; @@ -91,11 +90,11 @@ class sparse_cwise_binary_op_inner_iterator_selector<BinaryOp, Lhs, Rhs, Derived { typedef CwiseBinaryOp<BinaryOp, Lhs, Rhs> CwiseBinaryXpr; typedef typename traits<CwiseBinaryXpr>::Scalar Scalar; + typedef typename traits<CwiseBinaryXpr>::Index Index; typedef typename traits<CwiseBinaryXpr>::_LhsNested _LhsNested; typedef typename traits<CwiseBinaryXpr>::_RhsNested _RhsNested; typedef typename _LhsNested::InnerIterator LhsIterator; typedef typename _RhsNested::InnerIterator RhsIterator; - typedef typename Lhs::Index Index; public: @@ -157,11 +156,11 @@ class sparse_cwise_binary_op_inner_iterator_selector<scalar_product_op<T>, Lhs, typedef scalar_product_op<T> BinaryFunc; typedef CwiseBinaryOp<BinaryFunc, Lhs, Rhs> CwiseBinaryXpr; typedef typename CwiseBinaryXpr::Scalar Scalar; + typedef typename CwiseBinaryXpr::Index Index; typedef typename traits<CwiseBinaryXpr>::_LhsNested _LhsNested; typedef typename _LhsNested::InnerIterator LhsIterator; typedef typename traits<CwiseBinaryXpr>::_RhsNested _RhsNested; typedef typename _RhsNested::InnerIterator RhsIterator; - typedef typename Lhs::Index Index; public: EIGEN_STRONG_INLINE sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, Index outer) @@ -211,15 +210,15 @@ class sparse_cwise_binary_op_inner_iterator_selector<scalar_product_op<T>, Lhs, typedef scalar_product_op<T> BinaryFunc; typedef CwiseBinaryOp<BinaryFunc, Lhs, Rhs> CwiseBinaryXpr; typedef typename CwiseBinaryXpr::Scalar Scalar; + typedef typename CwiseBinaryXpr::Index Index; typedef typename traits<CwiseBinaryXpr>::_LhsNested _LhsNested; typedef typename traits<CwiseBinaryXpr>::RhsNested RhsNested; typedef typename _LhsNested::InnerIterator LhsIterator; - typedef typename Lhs::Index Index; enum { IsRowMajor = (int(Lhs::Flags)&RowMajorBit)==RowMajorBit }; public: EIGEN_STRONG_INLINE sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, Index outer) - : m_rhs(xpr.rhs()), m_lhsIter(xpr.lhs(),outer), m_functor(xpr.functor()), m_outer(outer) + : m_rhs(xpr.rhs()), m_lhsIter(xpr.lhs(),typename _LhsNested::Index(outer)), m_functor(xpr.functor()), m_outer(outer) {} EIGEN_STRONG_INLINE Derived& operator++() @@ -252,9 +251,9 @@ class sparse_cwise_binary_op_inner_iterator_selector<scalar_product_op<T>, Lhs, typedef scalar_product_op<T> BinaryFunc; typedef CwiseBinaryOp<BinaryFunc, Lhs, Rhs> CwiseBinaryXpr; typedef typename CwiseBinaryXpr::Scalar Scalar; + typedef typename CwiseBinaryXpr::Index Index; typedef typename traits<CwiseBinaryXpr>::_RhsNested _RhsNested; typedef typename _RhsNested::InnerIterator RhsIterator; - typedef typename Lhs::Index Index; enum { IsRowMajor = (int(Rhs::Flags)&RowMajorBit)==RowMajorBit }; public: diff --git a/Eigen/src/SparseCore/SparseDenseProduct.h b/Eigen/src/SparseCore/SparseDenseProduct.h index 36d9f637c..116edd62e 100644 --- a/Eigen/src/SparseCore/SparseDenseProduct.h +++ b/Eigen/src/SparseCore/SparseDenseProduct.h @@ -248,8 +248,8 @@ class SparseDenseOuterProduct EIGEN_STATIC_ASSERT(Tr,YOU_MADE_A_PROGRAMMING_MISTAKE); } - EIGEN_STRONG_INLINE Index rows() const { return Tr ? m_rhs.rows() : m_lhs.rows(); } - EIGEN_STRONG_INLINE Index cols() const { return Tr ? m_lhs.cols() : m_rhs.cols(); } + EIGEN_STRONG_INLINE Index rows() const { return Tr ? Index(m_rhs.rows()) : m_lhs.rows(); } + EIGEN_STRONG_INLINE Index cols() const { return Tr ? m_lhs.cols() : Index(m_rhs.cols()); } EIGEN_STRONG_INLINE const _LhsNested& lhs() const { return m_lhs; } EIGEN_STRONG_INLINE const _RhsNested& rhs() const { return m_rhs; } @@ -266,31 +266,33 @@ class SparseDenseOuterProduct<Lhs,Rhs,Transpose>::InnerIterator : public _LhsNes typedef typename SparseDenseOuterProduct::Index Index; public: EIGEN_STRONG_INLINE InnerIterator(const SparseDenseOuterProduct& prod, Index outer) - : Base(prod.lhs(), 0), m_outer(outer), m_factor(get(prod.rhs(), outer, typename internal::traits<Rhs>::StorageKind() )) - { } + : Base(prod.lhs(), 0), m_outer(outer), m_empty(false), m_factor(get(prod.rhs(), outer, typename internal::traits<Rhs>::StorageKind() )) + {} inline Index outer() const { return m_outer; } inline Index row() const { return Transpose ? m_outer : Base::index(); } inline Index col() const { return Transpose ? Base::index() : m_outer; } inline Scalar value() const { return Base::value() * m_factor; } + inline operator bool() const { return Base::operator bool() && !m_empty; } protected: - static Scalar get(const _RhsNested &rhs, Index outer, Dense = Dense()) + Scalar get(const _RhsNested &rhs, Index outer, Dense = Dense()) const { return rhs.coeff(outer); } - static Scalar get(const _RhsNested &rhs, Index outer, Sparse = Sparse()) + Scalar get(const _RhsNested &rhs, Index outer, Sparse = Sparse()) { typename Traits::_RhsNested::InnerIterator it(rhs, outer); - if (it && it.index()==0) + if (it && it.index()==0 && it.value()!=Scalar(0)) return it.value(); - + m_empty = true; return Scalar(0); } Index m_outer; + bool m_empty; Scalar m_factor; }; diff --git a/Eigen/src/SparseCore/SparseDiagonalProduct.h b/Eigen/src/SparseCore/SparseDiagonalProduct.h index 0998feba3..4c51881e0 100644 --- a/Eigen/src/SparseCore/SparseDiagonalProduct.h +++ b/Eigen/src/SparseCore/SparseDiagonalProduct.h @@ -34,8 +34,10 @@ struct traits<SparseDiagonalProduct<Lhs, Rhs> > typedef typename remove_all<Lhs>::type _Lhs; typedef typename remove_all<Rhs>::type _Rhs; typedef typename _Lhs::Scalar Scalar; - typedef typename promote_index_type<typename traits<Lhs>::Index, - typename traits<Rhs>::Index>::type Index; + // propagate the index type of the sparse matrix + typedef typename conditional< is_diagonal<_Lhs>::ret, + typename traits<Rhs>::Index, + typename traits<Lhs>::Index>::type Index; typedef Sparse StorageKind; typedef MatrixXpr XprKind; enum { @@ -92,8 +94,8 @@ class SparseDiagonalProduct eigen_assert(lhs.cols() == rhs.rows() && "invalid sparse matrix * diagonal matrix product"); } - EIGEN_STRONG_INLINE Index rows() const { return m_lhs.rows(); } - EIGEN_STRONG_INLINE Index cols() const { return m_rhs.cols(); } + EIGEN_STRONG_INLINE Index rows() const { return Index(m_lhs.rows()); } + EIGEN_STRONG_INLINE Index cols() const { return Index(m_rhs.cols()); } EIGEN_STRONG_INLINE const _LhsNested& lhs() const { return m_lhs; } EIGEN_STRONG_INLINE const _RhsNested& rhs() const { return m_rhs; } @@ -116,7 +118,7 @@ class sparse_diagonal_product_inner_iterator_selector : public CwiseUnaryOp<scalar_multiple_op<typename Lhs::Scalar>,const Rhs>::InnerIterator { typedef typename CwiseUnaryOp<scalar_multiple_op<typename Lhs::Scalar>,const Rhs>::InnerIterator Base; - typedef typename Lhs::Index Index; + typedef typename Rhs::Index Index; public: inline sparse_diagonal_product_inner_iterator_selector( const SparseDiagonalProductType& expr, Index outer) @@ -136,7 +138,7 @@ class sparse_diagonal_product_inner_iterator_selector scalar_product_op<typename Lhs::Scalar>, const typename Rhs::ConstInnerVectorReturnType, const typename Lhs::DiagonalVectorType>::InnerIterator Base; - typedef typename Lhs::Index Index; + typedef typename Rhs::Index Index; Index m_outer; public: inline sparse_diagonal_product_inner_iterator_selector( diff --git a/Eigen/src/SparseCore/SparseMatrix.h b/Eigen/src/SparseCore/SparseMatrix.h index db39d901b..9d0d6c3ab 100644 --- a/Eigen/src/SparseCore/SparseMatrix.h +++ b/Eigen/src/SparseCore/SparseMatrix.h @@ -808,7 +808,9 @@ protected: template<typename Other> void initAssignment(const Other& other) { - resize(other.rows(), other.cols()); + eigen_assert( other.rows() == typename Other::Index(Index(other.rows())) + && other.cols() == typename Other::Index(Index(other.cols())) ); + resize(Index(other.rows()), Index(other.cols())); if(m_innerNonZeros) { std::free(m_innerNonZeros); @@ -948,7 +950,7 @@ void set_from_triplets(const InputIterator& begin, const InputIterator& end, Spa enum { IsRowMajor = SparseMatrixType::IsRowMajor }; typedef typename SparseMatrixType::Scalar Scalar; typedef typename SparseMatrixType::Index Index; - SparseMatrix<Scalar,IsRowMajor?ColMajor:RowMajor> trMat(mat.rows(),mat.cols()); + SparseMatrix<Scalar,IsRowMajor?ColMajor:RowMajor,Index> trMat(mat.rows(),mat.cols()); if(begin!=end) { diff --git a/Eigen/src/SparseCore/SparseMatrixBase.h b/Eigen/src/SparseCore/SparseMatrixBase.h index 22c4040e9..a6452fa4e 100644 --- a/Eigen/src/SparseCore/SparseMatrixBase.h +++ b/Eigen/src/SparseCore/SparseMatrixBase.h @@ -366,7 +366,7 @@ template<typename Derived> class SparseMatrixBase : public EigenBase<Derived> { dst.setZero(); for (Index j=0; j<outerSize(); ++j) - for (typename Derived::InnerIterator i(derived(),j); i; ++i) + for (typename Derived::InnerIterator i(derived(),typename Derived::Index(j)); i; ++i) dst.coeffRef(i.row(),i.col()) = i.value(); } #endif diff --git a/Eigen/src/SparseCore/SparseSelfAdjointView.h b/Eigen/src/SparseCore/SparseSelfAdjointView.h index 0eda96bc4..56c922929 100644 --- a/Eigen/src/SparseCore/SparseSelfAdjointView.h +++ b/Eigen/src/SparseCore/SparseSelfAdjointView.h @@ -246,7 +246,7 @@ class SparseSelfAdjointTimeDenseProduct || ( (UpLo&Lower) && LhsIsRowMajor), ProcessSecondHalf = !ProcessFirstHalf }; - for (Index j=0; j<m_lhs.outerSize(); ++j) + for (typename _Lhs::Index j=0; j<m_lhs.outerSize(); ++j) { LhsInnerIterator i(m_lhs,j); if (ProcessSecondHalf) diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp index d441a6eca..36721b496 100644 --- a/test/jacobisvd.cpp +++ b/test/jacobisvd.cpp @@ -67,6 +67,7 @@ template<typename MatrixType, int QRPreconditioner> void jacobisvd_solve(const MatrixType& m, unsigned int computationOptions) { typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; typedef typename MatrixType::Index Index; Index rows = m.rows(); Index cols = m.cols(); @@ -81,9 +82,37 @@ void jacobisvd_solve(const MatrixType& m, unsigned int computationOptions) RhsType rhs = RhsType::Random(rows, internal::random<Index>(1, cols)); JacobiSVD<MatrixType, QRPreconditioner> svd(m, computationOptions); + + if(internal::is_same<RealScalar,double>::value) svd.setThreshold(1e-8); + else if(internal::is_same<RealScalar,float>::value) svd.setThreshold(1e-4); + SolutionType x = svd.solve(rhs); + + RealScalar residual = (m*x-rhs).norm(); + // Check that there is no significantly better solution in the neighborhood of x + if(!test_isMuchSmallerThan(residual,rhs.norm())) + { + // If the residual is very small, then we have an exact solution, so we are already good. + for(int k=0;k<x.rows();++k) + { + SolutionType y(x); + y.row(k).array() += 2*NumTraits<RealScalar>::epsilon(); + RealScalar residual_y = (m*y-rhs).norm(); + VERIFY( test_isApprox(residual_y,residual) || residual < residual_y ); + + y.row(k) = x.row(k).array() - 2*NumTraits<RealScalar>::epsilon(); + residual_y = (m*y-rhs).norm(); + VERIFY( test_isApprox(residual_y,residual) || residual < residual_y ); + } + } + // evaluate normal equation which works also for least-squares solutions - VERIFY_IS_APPROX(m.adjoint()*m*x,m.adjoint()*rhs); + if(internal::is_same<RealScalar,double>::value) + { + // This test is not stable with single precision. + // This is probably because squaring m signicantly affects the precision. + VERIFY_IS_APPROX(m.adjoint()*m*x,m.adjoint()*rhs); + } // check minimal norm solutions { @@ -139,10 +168,9 @@ void jacobisvd_test_all_computation_options(const MatrixType& m) if (QRPreconditioner == NoQRPreconditioner && m.rows() != m.cols()) return; JacobiSVD<MatrixType, QRPreconditioner> fullSvd(m, ComputeFullU|ComputeFullV); - - jacobisvd_check_full(m, fullSvd); - jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeFullU | ComputeFullV); - + CALL_SUBTEST(( jacobisvd_check_full(m, fullSvd) )); + CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeFullU | ComputeFullV) )); + #if defined __INTEL_COMPILER // remark #111: statement is unreachable #pragma warning disable 111 @@ -150,20 +178,20 @@ void jacobisvd_test_all_computation_options(const MatrixType& m) if(QRPreconditioner == FullPivHouseholderQRPreconditioner) return; - jacobisvd_compare_to_full(m, ComputeFullU, fullSvd); - jacobisvd_compare_to_full(m, ComputeFullV, fullSvd); - jacobisvd_compare_to_full(m, 0, fullSvd); + CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeFullU, fullSvd) )); + CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeFullV, fullSvd) )); + CALL_SUBTEST(( jacobisvd_compare_to_full(m, 0, fullSvd) )); if (MatrixType::ColsAtCompileTime == Dynamic) { // thin U/V are only available with dynamic number of columns - jacobisvd_compare_to_full(m, ComputeFullU|ComputeThinV, fullSvd); - jacobisvd_compare_to_full(m, ComputeThinV, fullSvd); - jacobisvd_compare_to_full(m, ComputeThinU|ComputeFullV, fullSvd); - jacobisvd_compare_to_full(m, ComputeThinU , fullSvd); - jacobisvd_compare_to_full(m, ComputeThinU|ComputeThinV, fullSvd); - jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeFullU | ComputeThinV); - jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeThinU | ComputeFullV); - jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeThinU | ComputeThinV); + CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeFullU|ComputeThinV, fullSvd) )); + CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeThinV, fullSvd) )); + CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeThinU|ComputeFullV, fullSvd) )); + CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeThinU , fullSvd) )); + CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeThinU|ComputeThinV, fullSvd) )); + CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeFullU | ComputeThinV) )); + CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeThinU | ComputeFullV) )); + CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeThinU | ComputeThinV) )); // test reconstruction typedef typename MatrixType::Index Index; @@ -176,12 +204,29 @@ void jacobisvd_test_all_computation_options(const MatrixType& m) template<typename MatrixType> void jacobisvd(const MatrixType& a = MatrixType(), bool pickrandom = true) { - MatrixType m = pickrandom ? MatrixType::Random(a.rows(), a.cols()) : a; + MatrixType m = a; + if(pickrandom) + { + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef typename MatrixType::Index Index; + Index diagSize = (std::min)(a.rows(), a.cols()); + RealScalar s = std::numeric_limits<RealScalar>::max_exponent10/4; + s = internal::random<RealScalar>(1,s); + Matrix<RealScalar,Dynamic,1> d = Matrix<RealScalar,Dynamic,1>::Random(diagSize); + for(Index k=0; k<diagSize; ++k) + d(k) = d(k)*std::pow(RealScalar(10),internal::random<RealScalar>(-s,s)); + m = Matrix<Scalar,Dynamic,Dynamic>::Random(a.rows(),diagSize) * d.asDiagonal() * Matrix<Scalar,Dynamic,Dynamic>::Random(diagSize,a.cols()); + // cancel some coeffs + Index n = internal::random<Index>(0,m.size()-1); + for(Index i=0; i<n; ++i) + m(internal::random<Index>(0,m.rows()-1), internal::random<Index>(0,m.cols()-1)) = Scalar(0); + } - jacobisvd_test_all_computation_options<MatrixType, FullPivHouseholderQRPreconditioner>(m); - jacobisvd_test_all_computation_options<MatrixType, ColPivHouseholderQRPreconditioner>(m); - jacobisvd_test_all_computation_options<MatrixType, HouseholderQRPreconditioner>(m); - jacobisvd_test_all_computation_options<MatrixType, NoQRPreconditioner>(m); + CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, FullPivHouseholderQRPreconditioner>(m) )); + CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, ColPivHouseholderQRPreconditioner>(m) )); + CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, HouseholderQRPreconditioner>(m) )); + CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, NoQRPreconditioner>(m) )); } template<typename MatrixType> void jacobisvd_verify_assert(const MatrixType& m) @@ -384,6 +429,7 @@ void test_jacobisvd() TEST_SET_BUT_UNUSED_VARIABLE(r) TEST_SET_BUT_UNUSED_VARIABLE(c) + CALL_SUBTEST_10(( jacobisvd<MatrixXd>(MatrixXd(r,c)) )); CALL_SUBTEST_7(( jacobisvd<MatrixXf>(MatrixXf(r,c)) )); CALL_SUBTEST_8(( jacobisvd<MatrixXcd>(MatrixXcd(r,c)) )); (void) r; diff --git a/test/meta.cpp b/test/meta.cpp index 3302c5887..b8dea68e8 100644 --- a/test/meta.cpp +++ b/test/meta.cpp @@ -9,6 +9,12 @@ #include "main.h" +template<typename From, typename To> +bool check_is_convertible(const From&, const To&) +{ + return internal::is_convertible<From,To>::value; +} + void test_meta() { VERIFY((internal::conditional<(3<4),internal::true_type, internal::false_type>::type::value)); @@ -52,6 +58,24 @@ void test_meta() VERIFY(( internal::is_same<const float,internal::remove_pointer<const float*>::type >::value)); VERIFY(( internal::is_same<float,internal::remove_pointer<float* const >::type >::value)); + VERIFY(( internal::is_convertible<float,double>::value )); + VERIFY(( internal::is_convertible<int,double>::value )); + VERIFY(( internal::is_convertible<double,int>::value )); + VERIFY((!internal::is_convertible<std::complex<double>,double>::value )); + VERIFY(( internal::is_convertible<Array33f,Matrix3f>::value )); +// VERIFY((!internal::is_convertible<Matrix3f,Matrix3d>::value )); //does not work because the conversion is prevented by a static assertion + VERIFY((!internal::is_convertible<Array33f,int>::value )); + VERIFY((!internal::is_convertible<MatrixXf,float>::value )); + { + float f; + MatrixXf A, B; + VectorXf a, b; + VERIFY(( check_is_convertible(a.dot(b), f) )); + VERIFY(( check_is_convertible(a.transpose()*b, f) )); + VERIFY((!check_is_convertible(A*B, f) )); + VERIFY(( check_is_convertible(A*B, A) )); + } + VERIFY(internal::meta_sqrt<1>::ret == 1); #define VERIFY_META_SQRT(X) VERIFY(internal::meta_sqrt<X>::ret == int(std::sqrt(double(X)))) VERIFY_META_SQRT(2); diff --git a/test/product_large.cpp b/test/product_large.cpp index 03d7bd8ed..11531aa1d 100644 --- a/test/product_large.cpp +++ b/test/product_large.cpp @@ -61,4 +61,13 @@ void test_product_large() VERIFY_IS_APPROX(r2, (mat1.row(2)*mat2).eval()); } #endif + + // Regression test for bug 714: +#ifdef EIGEN_HAS_OPENMP + std::cout << "Testing omp_set_dynamic(1)\n"; + omp_set_dynamic(1); + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_6( product(Matrix<float,Dynamic,Dynamic>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) ); + } +#endif } diff --git a/test/product_trsolve.cpp b/test/product_trsolve.cpp index 69892b3a8..4b97fa9d6 100644 --- a/test/product_trsolve.cpp +++ b/test/product_trsolve.cpp @@ -84,10 +84,18 @@ void test_product_trsolve() CALL_SUBTEST_4((trsolve<std::complex<double>,Dynamic,Dynamic>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2),internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2)))); // vectors - CALL_SUBTEST_1((trsolve<float,Dynamic,1>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE)))); - CALL_SUBTEST_5((trsolve<std::complex<double>,Dynamic,1>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE)))); - CALL_SUBTEST_6((trsolve<float,1,1>())); - CALL_SUBTEST_7((trsolve<float,1,2>())); - CALL_SUBTEST_8((trsolve<std::complex<float>,4,1>())); + CALL_SUBTEST_5((trsolve<float,Dynamic,1>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE)))); + CALL_SUBTEST_6((trsolve<double,Dynamic,1>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE)))); + CALL_SUBTEST_7((trsolve<std::complex<float>,Dynamic,1>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE)))); + CALL_SUBTEST_8((trsolve<std::complex<double>,Dynamic,1>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE)))); + + // meta-unrollers + CALL_SUBTEST_9((trsolve<float,4,1>())); + CALL_SUBTEST_10((trsolve<double,4,1>())); + CALL_SUBTEST_11((trsolve<std::complex<float>,4,1>())); + CALL_SUBTEST_12((trsolve<float,1,1>())); + CALL_SUBTEST_13((trsolve<float,1,2>())); + CALL_SUBTEST_14((trsolve<float,3,1>())); + } } diff --git a/test/sparse_product.cpp b/test/sparse_product.cpp index 6d32132bc..27bc548f8 100644 --- a/test/sparse_product.cpp +++ b/test/sparse_product.cpp @@ -108,28 +108,39 @@ template<typename SparseMatrixType> void sparse_product() Index r = internal::random<Index>(0,rows-1); Index c1 = internal::random<Index>(0,cols-1); Index r1 = internal::random<Index>(0,depth-1); + DenseMatrix dm5 = DenseMatrix::Random(depth, cols); - VERIFY_IS_APPROX( m4=m2.col(c)*refMat3.col(c1).transpose(), refMat4=refMat2.col(c)*refMat3.col(c1).transpose()); - VERIFY_IS_APPROX( m4=m2.middleCols(c,1)*refMat3.col(c1).transpose(), refMat4=refMat2.col(c)*refMat3.col(c1).transpose()); - VERIFY_IS_APPROX(dm4=m2.col(c)*refMat3.col(c1).transpose(), refMat4=refMat2.col(c)*refMat3.col(c1).transpose()); + VERIFY_IS_APPROX( m4=m2.col(c)*dm5.col(c1).transpose(), refMat4=refMat2.col(c)*dm5.col(c1).transpose()); + VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array()!=0).count()); + VERIFY_IS_APPROX( m4=m2.middleCols(c,1)*dm5.col(c1).transpose(), refMat4=refMat2.col(c)*dm5.col(c1).transpose()); + VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array()!=0).count()); + VERIFY_IS_APPROX(dm4=m2.col(c)*dm5.col(c1).transpose(), refMat4=refMat2.col(c)*dm5.col(c1).transpose()); - VERIFY_IS_APPROX(m4=refMat3.col(c1)*m2.col(c).transpose(), refMat4=refMat3.col(c1)*refMat2.col(c).transpose()); - VERIFY_IS_APPROX(m4=refMat3.col(c1)*m2.middleCols(c,1).transpose(), refMat4=refMat3.col(c1)*refMat2.col(c).transpose()); - VERIFY_IS_APPROX(dm4=refMat3.col(c1)*m2.col(c).transpose(), refMat4=refMat3.col(c1)*refMat2.col(c).transpose()); + VERIFY_IS_APPROX(m4=dm5.col(c1)*m2.col(c).transpose(), refMat4=dm5.col(c1)*refMat2.col(c).transpose()); + VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array()!=0).count()); + VERIFY_IS_APPROX(m4=dm5.col(c1)*m2.middleCols(c,1).transpose(), refMat4=dm5.col(c1)*refMat2.col(c).transpose()); + VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array()!=0).count()); + VERIFY_IS_APPROX(dm4=dm5.col(c1)*m2.col(c).transpose(), refMat4=dm5.col(c1)*refMat2.col(c).transpose()); - VERIFY_IS_APPROX( m4=refMat3.row(r1).transpose()*m2.col(c).transpose(), refMat4=refMat3.row(r1).transpose()*refMat2.col(c).transpose()); - VERIFY_IS_APPROX(dm4=refMat3.row(r1).transpose()*m2.col(c).transpose(), refMat4=refMat3.row(r1).transpose()*refMat2.col(c).transpose()); + VERIFY_IS_APPROX( m4=dm5.row(r1).transpose()*m2.col(c).transpose(), refMat4=dm5.row(r1).transpose()*refMat2.col(c).transpose()); + VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array()!=0).count()); + VERIFY_IS_APPROX(dm4=dm5.row(r1).transpose()*m2.col(c).transpose(), refMat4=dm5.row(r1).transpose()*refMat2.col(c).transpose()); - VERIFY_IS_APPROX( m4=m2.row(r).transpose()*refMat3.col(c1).transpose(), refMat4=refMat2.row(r).transpose()*refMat3.col(c1).transpose()); - VERIFY_IS_APPROX( m4=m2.middleRows(r,1).transpose()*refMat3.col(c1).transpose(), refMat4=refMat2.row(r).transpose()*refMat3.col(c1).transpose()); - VERIFY_IS_APPROX(dm4=m2.row(r).transpose()*refMat3.col(c1).transpose(), refMat4=refMat2.row(r).transpose()*refMat3.col(c1).transpose()); + VERIFY_IS_APPROX( m4=m2.row(r).transpose()*dm5.col(c1).transpose(), refMat4=refMat2.row(r).transpose()*dm5.col(c1).transpose()); + VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array()!=0).count()); + VERIFY_IS_APPROX( m4=m2.middleRows(r,1).transpose()*dm5.col(c1).transpose(), refMat4=refMat2.row(r).transpose()*dm5.col(c1).transpose()); + VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array()!=0).count()); + VERIFY_IS_APPROX(dm4=m2.row(r).transpose()*dm5.col(c1).transpose(), refMat4=refMat2.row(r).transpose()*dm5.col(c1).transpose()); - VERIFY_IS_APPROX( m4=refMat3.col(c1)*m2.row(r), refMat4=refMat3.col(c1)*refMat2.row(r)); - VERIFY_IS_APPROX( m4=refMat3.col(c1)*m2.middleRows(r,1), refMat4=refMat3.col(c1)*refMat2.row(r)); - VERIFY_IS_APPROX(dm4=refMat3.col(c1)*m2.row(r), refMat4=refMat3.col(c1)*refMat2.row(r)); + VERIFY_IS_APPROX( m4=dm5.col(c1)*m2.row(r), refMat4=dm5.col(c1)*refMat2.row(r)); + VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array()!=0).count()); + VERIFY_IS_APPROX( m4=dm5.col(c1)*m2.middleRows(r,1), refMat4=dm5.col(c1)*refMat2.row(r)); + VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array()!=0).count()); + VERIFY_IS_APPROX(dm4=dm5.col(c1)*m2.row(r), refMat4=dm5.col(c1)*refMat2.row(r)); - VERIFY_IS_APPROX( m4=refMat3.row(r1).transpose()*m2.row(r), refMat4=refMat3.row(r1).transpose()*refMat2.row(r)); - VERIFY_IS_APPROX(dm4=refMat3.row(r1).transpose()*m2.row(r), refMat4=refMat3.row(r1).transpose()*refMat2.row(r)); + VERIFY_IS_APPROX( m4=dm5.row(r1).transpose()*m2.row(r), refMat4=dm5.row(r1).transpose()*refMat2.row(r)); + VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array()!=0).count()); + VERIFY_IS_APPROX(dm4=dm5.row(r1).transpose()*m2.row(r), refMat4=dm5.row(r1).transpose()*refMat2.row(r)); } VERIFY_IS_APPROX(m6=m6*m6, refMat6=refMat6*refMat6); diff --git a/unsupported/Eigen/MPRealSupport b/unsupported/Eigen/MPRealSupport index d4b03647d..35d77e5bd 100644 --- a/unsupported/Eigen/MPRealSupport +++ b/unsupported/Eigen/MPRealSupport @@ -27,6 +27,8 @@ namespace Eigen { * via the <a href="http://www.holoborodko.com/pavel/mpfr">MPFR C++</a> * library which itself is built upon <a href="http://www.mpfr.org/">MPFR</a>/<a href="http://gmplib.org/">GMP</a>. * + * \warning MPFR C++ is licensed under the GPL. + * * You can find a copy of MPFR C++ that is known to be compatible in the unsupported/test/mpreal folder. * * Here is an example: @@ -139,64 +141,50 @@ int main() public: typedef mpfr::mpreal ResScalar; enum { - nr = 2, // must be 2 for proper packing... + nr = 1, mr = 1, - WorkSpaceFactor = nr, LhsProgress = 1, RhsProgress = 1 }; }; - template<typename Index, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs> - struct gebp_kernel<mpfr::mpreal,mpfr::mpreal,Index,mr,nr,ConjugateLhs,ConjugateRhs> + template<typename Index, bool ConjugateLhs, bool ConjugateRhs> + struct gebp_kernel<mpfr::mpreal,mpfr::mpreal,Index,1,1,ConjugateLhs,ConjugateRhs> { typedef mpfr::mpreal mpreal; EIGEN_DONT_INLINE void operator()(mpreal* res, Index resStride, const mpreal* blockA, const mpreal* blockB, Index rows, Index depth, Index cols, mpreal alpha, - Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0, mpreal* /*unpackedB*/ = 0) + Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0) { - mpreal acc1, acc2, tmp; + mpreal acc1(0,mpfr_get_prec(blockA[0].mpfr_srcptr())), + tmp (0,mpfr_get_prec(blockA[0].mpfr_srcptr())); if(strideA==-1) strideA = depth; if(strideB==-1) strideB = depth; - for(Index j=0; j<cols; j+=nr) + for(Index i=0; i<rows; ++i) { - Index actual_nr = (std::min<Index>)(nr,cols-j); - mpreal *C1 = res + j*resStride; - mpreal *C2 = res + (j+1)*resStride; - for(Index i=0; i<rows; i++) + for(Index j=0; j<cols; ++j) { - mpreal *B = const_cast<mpreal*>(blockB) + j*strideB + offsetB*actual_nr; - mpreal *A = const_cast<mpreal*>(blockA) + i*strideA + offsetA; + mpreal *C1 = res + j*resStride; + + const mpreal *A = blockA + i*strideA + offsetA; + const mpreal *B = blockB + j*strideB + offsetB; + acc1 = 0; - acc2 = 0; for(Index k=0; k<depth; k++) { - mpfr_mul(tmp.mpfr_ptr(), A[k].mpfr_ptr(), B[0].mpfr_ptr(), mpreal::get_default_rnd()); + mpfr_mul(tmp.mpfr_ptr(), A[k].mpfr_srcptr(), B[k].mpfr_srcptr(), mpreal::get_default_rnd()); mpfr_add(acc1.mpfr_ptr(), acc1.mpfr_ptr(), tmp.mpfr_ptr(), mpreal::get_default_rnd()); - - if(actual_nr==2) { - mpfr_mul(tmp.mpfr_ptr(), A[k].mpfr_ptr(), B[1].mpfr_ptr(), mpreal::get_default_rnd()); - mpfr_add(acc2.mpfr_ptr(), acc2.mpfr_ptr(), tmp.mpfr_ptr(), mpreal::get_default_rnd()); - } - - B+=actual_nr; } - mpfr_mul(acc1.mpfr_ptr(), acc1.mpfr_ptr(), alpha.mpfr_ptr(), mpreal::get_default_rnd()); - mpfr_add(C1[i].mpfr_ptr(), C1[i].mpfr_ptr(), acc1.mpfr_ptr(), mpreal::get_default_rnd()); - - if(actual_nr==2) { - mpfr_mul(acc2.mpfr_ptr(), acc2.mpfr_ptr(), alpha.mpfr_ptr(), mpreal::get_default_rnd()); - mpfr_add(C2[i].mpfr_ptr(), C2[i].mpfr_ptr(), acc2.mpfr_ptr(), mpreal::get_default_rnd()); - } + mpfr_mul(acc1.mpfr_ptr(), acc1.mpfr_srcptr(), alpha.mpfr_srcptr(), mpreal::get_default_rnd()); + mpfr_add(C1[i].mpfr_ptr(), C1[i].mpfr_srcptr(), acc1.mpfr_srcptr(), mpreal::get_default_rnd()); } } } }; - } // end namespace internal } diff --git a/unsupported/test/mpreal/mpreal.h b/unsupported/test/mpreal/mpreal.h index 70d37ab71..4c29d2ac2 100644 --- a/unsupported/test/mpreal/mpreal.h +++ b/unsupported/test/mpreal/mpreal.h @@ -1,59 +1,47 @@ /*
- Multi-precision floating point number class for C++.
+ MPFR C++: Multi-precision floating point number class for C++.
Based on MPFR library: http://mpfr.org
Project homepage: http://www.holoborodko.com/pavel/mpfr
Contact e-mail: pavel@holoborodko.com
- Copyright (c) 2008-2013 Pavel Holoborodko
+ Copyright (c) 2008-2014 Pavel Holoborodko
Contributors:
Dmitriy Gubanov, Konstantin Holoborodko, Brian Gladman,
Helmut Jarausch, Fokko Beekhof, Ulrich Mutze, Heinz van Saanen,
Pere Constans, Peter van Hoof, Gael Guennebaud, Tsai Chia Cheng,
- Alexei Zubanov, Jauhien Piatlicki, Victor Berger, John Westwood.
+ Alexei Zubanov, Jauhien Piatlicki, Victor Berger, John Westwood,
+ Petr Aleksandrov, Orion Poplawski, Charles Karney.
- ****************************************************************************
- This library is free software; you can redistribute it and/or
- modify it under the terms of the GNU Lesser General Public
- License as published by the Free Software Foundation; either
- version 2.1 of the License, or (at your option) any later version.
-
- This library is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- Lesser General Public License for more details.
+ Licensing:
+ (A) MPFR C++ is under GNU General Public License ("GPL").
+
+ (B) Non-free licenses may also be purchased from the author, for users who
+ do not want their programs protected by the GPL.
- You should have received a copy of the GNU Lesser General Public
- License along with this library; if not, write to the Free Software
- Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ The non-free licenses are for users that wish to use MPFR C++ in
+ their products but are unwilling to release their software
+ under the GPL (which would require them to release source code
+ and allow free redistribution).
- ****************************************************************************
- Redistribution and use in source and binary forms, with or without
- modification, are permitted provided that the following conditions
- are met:
-
- 1. Redistributions of source code must retain the above copyright
- notice, this list of conditions and the following disclaimer.
-
- 2. Redistributions in binary form must reproduce the above copyright
- notice, this list of conditions and the following disclaimer in the
- documentation and/or other materials provided with the distribution.
+ Such users can purchase an unlimited-use license from the author.
+ Contact us for more details.
- 3. The name of the author may not be used to endorse or promote products
- derived from this software without specific prior written permission.
-
- THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
- ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
- FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- SUCH DAMAGE.
+ GNU General Public License ("GPL") copyright permissions statement:
+ **************************************************************************
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef __MPREAL_H__
@@ -71,6 +59,15 @@ // Options
#define MPREAL_HAVE_INT64_SUPPORT // Enable int64_t support if possible. Available only for MSVC 2010 & GCC.
#define MPREAL_HAVE_MSVC_DEBUGVIEW // Enable Debugger Visualizer for "Debug" builds in MSVC.
+#define MPREAL_HAVE_DYNAMIC_STD_NUMERIC_LIMITS // Enable extended std::numeric_limits<mpfr::mpreal> specialization.
+ // Meaning that "digits", "round_style" and similar members are defined as functions, not constants.
+ // See std::numeric_limits<mpfr::mpreal> at the end of the file for more information.
+
+// Library version
+#define MPREAL_VERSION_MAJOR 3
+#define MPREAL_VERSION_MINOR 5
+#define MPREAL_VERSION_PATCHLEVEL 9
+#define MPREAL_VERSION_STRING "3.5.9"
// Detect compiler using signatures from http://predef.sourceforge.net/
#if defined(__GNUC__) && defined(__INTEL_COMPILER)
@@ -83,17 +80,29 @@ #define IsInf(x) std::isinf(x) // GNU C/C++ (and/or other compilers), just hope for C99 conformance
#endif
-// Detect support for r-value references (move semantic). Borrowed from Eigen.
// A Clang feature extension to determine compiler features.
-// We use it to determine 'cxx_rvalue_references'
#ifndef __has_feature
-# define __has_feature(x) 0
+ #define __has_feature(x) 0
#endif
+// Detect support for r-value references (move semantic). Borrowed from Eigen.
#if (__has_feature(cxx_rvalue_references) || \
- defined(__GXX_EXPERIMENTAL_CXX0X__) || \
- (defined(_MSC_VER) && _MSC_VER >= 1600))
-#define MPREAL_HAVE_MOVE_SUPPORT
+ defined(__GXX_EXPERIMENTAL_CXX0X__) || __cplusplus >= 201103L || \
+ (defined(_MSC_VER) && _MSC_VER >= 1600))
+
+ #define MPREAL_HAVE_MOVE_SUPPORT
+
+ // Use fields in mpfr_t structure to check if it was initialized / set dummy initialization
+ #define mpfr_is_initialized(x) (0 != (x)->_mpfr_d)
+ #define mpfr_set_uninitialized(x) ((x)->_mpfr_d = 0 )
+#endif
+
+// Detect support for explicit converters.
+#if (__has_feature(cxx_explicit_conversions) || \
+ defined(__GXX_EXPERIMENTAL_CXX0X__) || __cplusplus >= 201103L || \
+ (defined(_MSC_VER) && _MSC_VER >= 1800))
+
+ #define MPREAL_HAVE_EXPLICIT_CONVERTERS
#endif
// Detect available 64-bit capabilities
@@ -112,7 +121,7 @@ #endif
#elif defined (__GNUC__) && defined(__linux__)
- #if defined(__amd64__) || defined(__amd64) || defined(__x86_64__) || defined(__x86_64) || defined(__ia64) || defined(__itanium__) || defined(_M_IA64)
+ #if defined(__amd64__) || defined(__amd64) || defined(__x86_64__) || defined(__x86_64) || defined(__ia64) || defined(__itanium__) || defined(_M_IA64) || defined (__PPC64__)
#undef MPREAL_HAVE_INT64_SUPPORT // Remove all shaman dances for x64 builds since
#undef MPFR_USE_INTMAX_T // GCC already supports x64 as of "long int" is 64-bit integer, nothing left to do
#else
@@ -175,7 +184,7 @@ public: mpreal(const int u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd());
// Construct mpreal from mpfr_t structure.
- // shared = true allows to avoid deep copy, so that mpreal and 'u' shared same data & pointers.
+ // shared = true allows to avoid deep copy, so that mpreal and 'u' share the same data & pointers.
mpreal(const mpfr_t u, bool shared = false);
#if defined (MPREAL_HAVE_INT64_SUPPORT)
@@ -315,14 +324,34 @@ public: friend bool operator == (const mpreal& a, const double b);
// Type Conversion operators
+ bool toBool (mp_rnd_t mode = GMP_RNDZ) const;
long toLong (mp_rnd_t mode = GMP_RNDZ) const;
unsigned long toULong (mp_rnd_t mode = GMP_RNDZ) const;
+ float toFloat (mp_rnd_t mode = GMP_RNDN) const;
double toDouble (mp_rnd_t mode = GMP_RNDN) const;
long double toLDouble (mp_rnd_t mode = GMP_RNDN) const;
+#if defined (MPREAL_HAVE_EXPLICIT_CONVERTERS)
+ explicit operator bool () const { return toBool(); }
+ explicit operator int () const { return toLong(); }
+ explicit operator long () const { return toLong(); }
+ explicit operator long long () const { return toLong(); }
+ explicit operator unsigned () const { return toULong(); }
+ explicit operator unsigned long () const { return toULong(); }
+ explicit operator unsigned long long () const { return toULong(); }
+ explicit operator float () const { return toFloat(); }
+ explicit operator double () const { return toDouble(); }
+ explicit operator long double () const { return toLDouble(); }
+#endif
+
#if defined (MPREAL_HAVE_INT64_SUPPORT)
int64_t toInt64 (mp_rnd_t mode = GMP_RNDZ) const;
uint64_t toUInt64 (mp_rnd_t mode = GMP_RNDZ) const;
+
+ #if defined (MPREAL_HAVE_EXPLICIT_CONVERTERS)
+ explicit operator int64_t () const { return toInt64(); }
+ explicit operator uint64_t () const { return toUInt64(); }
+ #endif
#endif
// Get raw pointers so that mpreal can be directly used in raw mpfr_* functions
@@ -331,122 +360,124 @@ public: ::mpfr_srcptr mpfr_srcptr() const;
// Convert mpreal to string with n significant digits in base b
- // n = 0 -> convert with the maximum available digits
- std::string toString(int n = 0, int b = 10, mp_rnd_t mode = mpreal::get_default_rnd()) const;
+ // n = -1 -> convert with the maximum available digits
+ std::string toString(int n = -1, int b = 10, mp_rnd_t mode = mpreal::get_default_rnd()) const;
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
- std::string toString(const std::string& format) const;
+ std::string toString(const std::string& format) const;
#endif
+ std::ostream& output(std::ostream& os) const;
+
// Math Functions
- friend const mpreal sqr (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal sqrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal sqrt(const unsigned long int v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal cbrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal root(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal pow (const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal pow (const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal pow (const mpreal& a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal pow (const mpreal& a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal pow (const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal pow (const unsigned long int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal fabs(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
-
- friend const mpreal abs(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+ friend const mpreal sqr (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal sqrt(const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal sqrt(const unsigned long int v, mp_rnd_t rnd_mode);
+ friend const mpreal cbrt(const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal root(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode);
+ friend const mpreal pow (const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode);
+ friend const mpreal pow (const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode);
+ friend const mpreal pow (const mpreal& a, const unsigned long int b, mp_rnd_t rnd_mode);
+ friend const mpreal pow (const mpreal& a, const long int b, mp_rnd_t rnd_mode);
+ friend const mpreal pow (const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode);
+ friend const mpreal pow (const unsigned long int a, const unsigned long int b, mp_rnd_t rnd_mode);
+ friend const mpreal fabs(const mpreal& v, mp_rnd_t rnd_mode);
+
+ friend const mpreal abs(const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode);
+ friend inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode);
+ friend inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode);
+ friend inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode);
+ friend inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode);
friend int cmpabs(const mpreal& a,const mpreal& b);
- friend const mpreal log (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal log2 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal log10(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal exp (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal exp2 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal exp10(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal log1p(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal expm1(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
-
- friend const mpreal cos(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal sin(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal tan(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal sec(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal csc(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal cot(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
-
- friend const mpreal acos (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal asin (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal atan (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal acot (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal asec (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal acsc (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
-
- friend const mpreal cosh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal sinh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal tanh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal sech (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal csch (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal coth (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal acosh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal asinh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal atanh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal acoth (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal asech (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal acsch (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
-
- friend const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
-
- friend const mpreal fac_ui (unsigned long int v, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal eint (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
-
- friend const mpreal gamma (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal lngamma (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal lgamma (const mpreal& v, int *signp = 0, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal zeta (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal erf (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal erfc (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal besselj0 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal besselj1 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal besseljn (long n, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal bessely0 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal bessely1 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal besselyn (long n, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+ friend const mpreal log (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal log2 (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal log10(const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal exp (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal exp2 (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal exp10(const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal log1p(const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal expm1(const mpreal& v, mp_rnd_t rnd_mode);
+
+ friend const mpreal cos(const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal sin(const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal tan(const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal sec(const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal csc(const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal cot(const mpreal& v, mp_rnd_t rnd_mode);
+ friend int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode);
+
+ friend const mpreal acos (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal asin (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal atan (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode);
+ friend const mpreal acot (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal asec (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal acsc (const mpreal& v, mp_rnd_t rnd_mode);
+
+ friend const mpreal cosh (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal sinh (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal tanh (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal sech (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal csch (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal coth (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal acosh (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal asinh (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal atanh (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal acoth (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal asech (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal acsch (const mpreal& v, mp_rnd_t rnd_mode);
+
+ friend const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode);
+
+ friend const mpreal fac_ui (unsigned long int v, mp_prec_t prec, mp_rnd_t rnd_mode);
+ friend const mpreal eint (const mpreal& v, mp_rnd_t rnd_mode);
+
+ friend const mpreal gamma (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal lngamma (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal lgamma (const mpreal& v, int *signp, mp_rnd_t rnd_mode);
+ friend const mpreal zeta (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal erf (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal erfc (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal besselj0 (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal besselj1 (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal besseljn (long n, const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal bessely0 (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal bessely1 (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal besselyn (long n, const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode);
+ friend const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode);
+ friend const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode);
+ friend const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode);
friend int sgn(const mpreal& v); // returns -1 or +1
// MPFR 2.4.0 Specifics
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
- friend int sinh_cosh (mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal li2 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal rec_sqrt (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+ friend int sinh_cosh (mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal li2 (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode);
+ friend const mpreal rec_sqrt (const mpreal& v, mp_rnd_t rnd_mode);
// MATLAB's semantic equivalents
- friend const mpreal rem (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); // Remainder after division
- friend const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); // Modulus after division
+ friend const mpreal rem (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); // Remainder after division
+ friend const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); // Modulus after division
#endif
// MPFR 3.0.0 Specifics
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
- friend const mpreal digamma (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal ai (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal urandom (gmp_randstate_t& state, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); // use gmp_randinit_default() to init state, gmp_randclear() to clear
- friend const mpreal grandom (gmp_randstate_t& state, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); // use gmp_randinit_default() to init state, gmp_randclear() to clear
- friend const mpreal grandom (unsigned int seed = 0);
+ friend const mpreal digamma (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal ai (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal urandom (gmp_randstate_t& state, mp_rnd_t rnd_mode); // use gmp_randinit_default() to init state, gmp_randclear() to clear
+ friend const mpreal grandom (gmp_randstate_t& state, mp_rnd_t rnd_mode); // use gmp_randinit_default() to init state, gmp_randclear() to clear
+ friend const mpreal grandom (unsigned int seed);
#endif
// Uniformly distributed random number generation in [0,1] using
// Mersenne-Twister algorithm by default.
// Use parameter to setup seed, e.g.: random((unsigned)time(NULL))
// Check urandom() for more precise control.
- friend const mpreal random(unsigned int seed = 0);
+ friend const mpreal random(unsigned int seed);
// Exponent and mantissa manipulation
friend const mpreal frexp(const mpreal& v, mp_exp_t* exp);
@@ -458,31 +489,31 @@ public: // Constants
// don't forget to call mpfr_free_cache() for every thread where you are using const-functions
- friend const mpreal const_log2 (mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal const_pi (mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal const_euler (mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal const_catalan (mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+ friend const mpreal const_log2 (mp_prec_t prec, mp_rnd_t rnd_mode);
+ friend const mpreal const_pi (mp_prec_t prec, mp_rnd_t rnd_mode);
+ friend const mpreal const_euler (mp_prec_t prec, mp_rnd_t rnd_mode);
+ friend const mpreal const_catalan (mp_prec_t prec, mp_rnd_t rnd_mode);
// returns +inf iff sign>=0 otherwise -inf
- friend const mpreal const_infinity(int sign = 1, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+ friend const mpreal const_infinity(int sign, mp_prec_t prec);
// Output/ Input
friend std::ostream& operator<<(std::ostream& os, const mpreal& v);
friend std::istream& operator>>(std::istream& is, mpreal& v);
// Integer Related Functions
- friend const mpreal rint (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+ friend const mpreal rint (const mpreal& v, mp_rnd_t rnd_mode);
friend const mpreal ceil (const mpreal& v);
friend const mpreal floor(const mpreal& v);
friend const mpreal round(const mpreal& v);
friend const mpreal trunc(const mpreal& v);
- friend const mpreal rint_ceil (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal rint_floor (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal rint_round (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal rint_trunc (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal frac (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal remainder ( const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+ friend const mpreal rint_ceil (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal rint_floor (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal rint_round (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal rint_trunc (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal frac (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal remainder ( const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode);
+ friend const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode);
// Miscellaneous Functions
friend const mpreal nexttoward (const mpreal& x, const mpreal& y);
@@ -549,8 +580,8 @@ public: // Efficient swapping of two mpreal values - needed for std algorithms
friend void swap(mpreal& x, mpreal& y);
- friend const mpreal fmax(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal fmin(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+ friend const mpreal fmax(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode);
+ friend const mpreal fmin(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode);
private:
// Human friendly Debug Preview in Visual Studio.
@@ -579,16 +610,16 @@ public: // Default constructor: creates mp number and initializes it to 0.
inline mpreal::mpreal()
{
- mpfr_init2(mp,mpreal::get_default_prec());
- mpfr_set_ui(mp,0,mpreal::get_default_rnd());
+ mpfr_init2 (mpfr_ptr(), mpreal::get_default_prec());
+ mpfr_set_ui(mpfr_ptr(), 0, mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mpreal::mpreal(const mpreal& u)
{
- mpfr_init2(mp,mpfr_get_prec(u.mp));
- mpfr_set(mp,u.mp,mpreal::get_default_rnd());
+ mpfr_init2(mpfr_ptr(),mpfr_get_prec(u.mpfr_srcptr()));
+ mpfr_set (mpfr_ptr(),u.mpfr_srcptr(),mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
}
@@ -596,8 +627,7 @@ inline mpreal::mpreal(const mpreal& u) #ifdef MPREAL_HAVE_MOVE_SUPPORT
inline mpreal::mpreal(mpreal&& other)
{
- mpfr_ptr()->_mpfr_d = 0; // make sure "other" holds no pinter to actual data
-
+ mpfr_set_uninitialized(mpfr_ptr()); // make sure "other" holds no pinter to actual data
mpfr_swap(mpfr_ptr(), other.mpfr_ptr());
MPREAL_MSVC_DEBUGVIEW_CODE;
@@ -616,12 +646,12 @@ inline mpreal::mpreal(const mpfr_t u, bool shared) {
if(shared)
{
- std::memcpy(mp, u, sizeof(mpfr_t));
+ std::memcpy(mpfr_ptr(), u, sizeof(mpfr_t));
}
else
{
- mpfr_init2(mp, mpfr_get_prec(u));
- mpfr_set (mp, u, mpreal::get_default_rnd());
+ mpfr_init2(mpfr_ptr(), mpfr_get_prec(u));
+ mpfr_set (mpfr_ptr(), u, mpreal::get_default_rnd());
}
MPREAL_MSVC_DEBUGVIEW_CODE;
@@ -629,40 +659,40 @@ inline mpreal::mpreal(const mpfr_t u, bool shared) inline mpreal::mpreal(const mpf_t u)
{
- mpfr_init2(mp,(mp_prec_t) mpf_get_prec(u)); // (gmp: mp_bitcnt_t) unsigned long -> long (mpfr: mp_prec_t)
- mpfr_set_f(mp,u,mpreal::get_default_rnd());
+ mpfr_init2(mpfr_ptr(),(mp_prec_t) mpf_get_prec(u)); // (gmp: mp_bitcnt_t) unsigned long -> long (mpfr: mp_prec_t)
+ mpfr_set_f(mpfr_ptr(),u,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mpreal::mpreal(const mpz_t u, mp_prec_t prec, mp_rnd_t mode)
{
- mpfr_init2(mp,prec);
- mpfr_set_z(mp,u,mode);
+ mpfr_init2(mpfr_ptr(), prec);
+ mpfr_set_z(mpfr_ptr(), u, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mpreal::mpreal(const mpq_t u, mp_prec_t prec, mp_rnd_t mode)
{
- mpfr_init2(mp,prec);
- mpfr_set_q(mp,u,mode);
+ mpfr_init2(mpfr_ptr(), prec);
+ mpfr_set_q(mpfr_ptr(), u, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode)
{
- mpfr_init2(mp, prec);
+ mpfr_init2(mpfr_ptr(), prec);
#if (MPREAL_DOUBLE_BITS_OVERFLOW > -1)
- if(fits_in_bits(u, MPREAL_DOUBLE_BITS_OVERFLOW))
- {
- mpfr_set_d(mp, u, mode);
- }else
- throw conversion_overflow();
+ if(fits_in_bits(u, MPREAL_DOUBLE_BITS_OVERFLOW))
+ {
+ mpfr_set_d(mpfr_ptr(), u, mode);
+ }else
+ throw conversion_overflow();
#else
- mpfr_set_d(mp, u, mode);
+ mpfr_set_d(mpfr_ptr(), u, mode);
#endif
MPREAL_MSVC_DEBUGVIEW_CODE;
@@ -670,40 +700,40 @@ inline mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode) inline mpreal::mpreal(const long double u, mp_prec_t prec, mp_rnd_t mode)
{
- mpfr_init2(mp,prec);
- mpfr_set_ld(mp,u,mode);
+ mpfr_init2 (mpfr_ptr(), prec);
+ mpfr_set_ld(mpfr_ptr(), u, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mpreal::mpreal(const unsigned long int u, mp_prec_t prec, mp_rnd_t mode)
{
- mpfr_init2(mp,prec);
- mpfr_set_ui(mp,u,mode);
+ mpfr_init2 (mpfr_ptr(), prec);
+ mpfr_set_ui(mpfr_ptr(), u, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mpreal::mpreal(const unsigned int u, mp_prec_t prec, mp_rnd_t mode)
{
- mpfr_init2(mp,prec);
- mpfr_set_ui(mp,u,mode);
+ mpfr_init2 (mpfr_ptr(), prec);
+ mpfr_set_ui(mpfr_ptr(), u, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mpreal::mpreal(const long int u, mp_prec_t prec, mp_rnd_t mode)
{
- mpfr_init2(mp,prec);
- mpfr_set_si(mp,u,mode);
+ mpfr_init2 (mpfr_ptr(), prec);
+ mpfr_set_si(mpfr_ptr(), u, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mpreal::mpreal(const int u, mp_prec_t prec, mp_rnd_t mode)
{
- mpfr_init2(mp,prec);
- mpfr_set_si(mp,u,mode);
+ mpfr_init2 (mpfr_ptr(), prec);
+ mpfr_set_si(mpfr_ptr(), u, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
@@ -711,16 +741,16 @@ inline mpreal::mpreal(const int u, mp_prec_t prec, mp_rnd_t mode) #if defined (MPREAL_HAVE_INT64_SUPPORT)
inline mpreal::mpreal(const uint64_t u, mp_prec_t prec, mp_rnd_t mode)
{
- mpfr_init2(mp,prec);
- mpfr_set_uj(mp, u, mode);
+ mpfr_init2 (mpfr_ptr(), prec);
+ mpfr_set_uj(mpfr_ptr(), u, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mpreal::mpreal(const int64_t u, mp_prec_t prec, mp_rnd_t mode)
{
- mpfr_init2(mp,prec);
- mpfr_set_sj(mp, u, mode);
+ mpfr_init2 (mpfr_ptr(), prec);
+ mpfr_set_sj(mpfr_ptr(), u, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
@@ -728,24 +758,26 @@ inline mpreal::mpreal(const int64_t u, mp_prec_t prec, mp_rnd_t mode) inline mpreal::mpreal(const char* s, mp_prec_t prec, int base, mp_rnd_t mode)
{
- mpfr_init2(mp, prec);
- mpfr_set_str(mp, s, base, mode);
+ mpfr_init2 (mpfr_ptr(), prec);
+ mpfr_set_str(mpfr_ptr(), s, base, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mpreal::mpreal(const std::string& s, mp_prec_t prec, int base, mp_rnd_t mode)
{
- mpfr_init2(mp, prec);
- mpfr_set_str(mp, s.c_str(), base, mode);
+ mpfr_init2 (mpfr_ptr(), prec);
+ mpfr_set_str(mpfr_ptr(), s.c_str(), base, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline void mpreal::clear(::mpfr_ptr x)
{
- // Use undocumented field in mpfr_t structure to check if it was initialized
- if(0 != x->_mpfr_d) mpfr_clear(x);
+#ifdef MPREAL_HAVE_MOVE_SUPPORT
+ if(mpfr_is_initialized(x))
+#endif
+ mpfr_clear(x);
}
inline mpreal::~mpreal()
@@ -821,6 +853,9 @@ const mpreal sqrt(const int v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); const mpreal sqrt(const long double v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
const mpreal sqrt(const double v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+// abs
+inline const mpreal abs(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd());
+
//////////////////////////////////////////////////////////////////////////
// pow
const mpreal pow(const mpreal& a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
@@ -873,6 +908,11 @@ const mpreal pow(const double a, const unsigned int b, mp_rnd_t rnd_mode = mprea const mpreal pow(const double a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
const mpreal pow(const double a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+
//////////////////////////////////////////////////////////////////////////
// Estimate machine epsilon for the given precision
// Returns smallest eps such that 1.0 + eps != 1.0
@@ -920,15 +960,15 @@ inline mpreal& mpreal::operator=(const mpreal& v) {
if (this != &v)
{
- mp_prec_t tp = mpfr_get_prec(mp);
- mp_prec_t vp = mpfr_get_prec(v.mp);
+ mp_prec_t tp = mpfr_get_prec( mpfr_srcptr());
+ mp_prec_t vp = mpfr_get_prec(v.mpfr_srcptr());
- if(tp != vp){
- clear(mpfr_ptr());
- mpfr_init2(mpfr_ptr(), vp);
- }
+ if(tp != vp){
+ clear(mpfr_ptr());
+ mpfr_init2(mpfr_ptr(), vp);
+ }
- mpfr_set(mp, v.mp, mpreal::get_default_rnd());
+ mpfr_set(mpfr_ptr(), v.mpfr_srcptr(), mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
}
@@ -937,7 +977,7 @@ inline mpreal& mpreal::operator=(const mpreal& v) inline mpreal& mpreal::operator=(const mpf_t v)
{
- mpfr_set_f(mp, v, mpreal::get_default_rnd());
+ mpfr_set_f(mpfr_ptr(), v, mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
@@ -945,7 +985,7 @@ inline mpreal& mpreal::operator=(const mpf_t v) inline mpreal& mpreal::operator=(const mpz_t v)
{
- mpfr_set_z(mp, v, mpreal::get_default_rnd());
+ mpfr_set_z(mpfr_ptr(), v, mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
@@ -953,7 +993,7 @@ inline mpreal& mpreal::operator=(const mpz_t v) inline mpreal& mpreal::operator=(const mpq_t v)
{
- mpfr_set_q(mp, v, mpreal::get_default_rnd());
+ mpfr_set_q(mpfr_ptr(), v, mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
@@ -961,7 +1001,7 @@ inline mpreal& mpreal::operator=(const mpq_t v) inline mpreal& mpreal::operator=(const long double v)
{
- mpfr_set_ld(mp, v, mpreal::get_default_rnd());
+ mpfr_set_ld(mpfr_ptr(), v, mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
@@ -970,22 +1010,22 @@ inline mpreal& mpreal::operator=(const long double v) inline mpreal& mpreal::operator=(const double v)
{
#if (MPREAL_DOUBLE_BITS_OVERFLOW > -1)
- if(fits_in_bits(v, MPREAL_DOUBLE_BITS_OVERFLOW))
- {
- mpfr_set_d(mp,v,mpreal::get_default_rnd());
- }else
- throw conversion_overflow();
+ if(fits_in_bits(v, MPREAL_DOUBLE_BITS_OVERFLOW))
+ {
+ mpfr_set_d(mpfr_ptr(),v,mpreal::get_default_rnd());
+ }else
+ throw conversion_overflow();
#else
- mpfr_set_d(mp,v,mpreal::get_default_rnd());
+ mpfr_set_d(mpfr_ptr(),v,mpreal::get_default_rnd());
#endif
- MPREAL_MSVC_DEBUGVIEW_CODE;
+ MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator=(const unsigned long int v)
{
- mpfr_set_ui(mp, v, mpreal::get_default_rnd());
+ mpfr_set_ui(mpfr_ptr(), v, mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
@@ -993,7 +1033,7 @@ inline mpreal& mpreal::operator=(const unsigned long int v) inline mpreal& mpreal::operator=(const unsigned int v)
{
- mpfr_set_ui(mp, v, mpreal::get_default_rnd());
+ mpfr_set_ui(mpfr_ptr(), v, mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
@@ -1001,7 +1041,7 @@ inline mpreal& mpreal::operator=(const unsigned int v) inline mpreal& mpreal::operator=(const long int v)
{
- mpfr_set_si(mp, v, mpreal::get_default_rnd());
+ mpfr_set_si(mpfr_ptr(), v, mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
@@ -1009,7 +1049,7 @@ inline mpreal& mpreal::operator=(const long int v) inline mpreal& mpreal::operator=(const int v)
{
- mpfr_set_si(mp, v, mpreal::get_default_rnd());
+ mpfr_set_si(mpfr_ptr(), v, mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
@@ -1026,11 +1066,11 @@ inline mpreal& mpreal::operator=(const char* s) mpfr_t t;
- mpfr_init2(t, mpfr_get_prec(mp));
+ mpfr_init2(t, mpfr_get_prec(mpfr_srcptr()));
if(0 == mpfr_set_str(t, s, 10, mpreal::get_default_rnd()))
{
- mpfr_set(mp, t, mpreal::get_default_rnd());
+ mpfr_set(mpfr_ptr(), t, mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
}
@@ -1049,11 +1089,11 @@ inline mpreal& mpreal::operator=(const std::string& s) mpfr_t t;
- mpfr_init2(t, mpfr_get_prec(mp));
+ mpfr_init2(t, mpfr_get_prec(mpfr_srcptr()));
if(0 == mpfr_set_str(t, s.c_str(), 10, mpreal::get_default_rnd()))
{
- mpfr_set(mp, t, mpreal::get_default_rnd());
+ mpfr_set(mpfr_ptr(), t, mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
}
@@ -1066,7 +1106,7 @@ inline mpreal& mpreal::operator=(const std::string& s) // + Addition
inline mpreal& mpreal::operator+=(const mpreal& v)
{
- mpfr_add(mp,mp,v.mp,mpreal::get_default_rnd());
+ mpfr_add(mpfr_ptr(), mpfr_srcptr(), v.mpfr_srcptr(), mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
@@ -1080,14 +1120,14 @@ inline mpreal& mpreal::operator+=(const mpf_t u) inline mpreal& mpreal::operator+=(const mpz_t u)
{
- mpfr_add_z(mp,mp,u,mpreal::get_default_rnd());
+ mpfr_add_z(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator+=(const mpq_t u)
{
- mpfr_add_q(mp,mp,u,mpreal::get_default_rnd());
+ mpfr_add_q(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
@@ -1102,7 +1142,7 @@ inline mpreal& mpreal::operator+= (const long double u) inline mpreal& mpreal::operator+= (const double u)
{
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
- mpfr_add_d(mp,mp,u,mpreal::get_default_rnd());
+ mpfr_add_d(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
#else
*this += mpreal(u);
#endif
@@ -1113,28 +1153,28 @@ inline mpreal& mpreal::operator+= (const double u) inline mpreal& mpreal::operator+=(const unsigned long int u)
{
- mpfr_add_ui(mp,mp,u,mpreal::get_default_rnd());
+ mpfr_add_ui(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator+=(const unsigned int u)
{
- mpfr_add_ui(mp,mp,u,mpreal::get_default_rnd());
+ mpfr_add_ui(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator+=(const long int u)
{
- mpfr_add_si(mp,mp,u,mpreal::get_default_rnd());
+ mpfr_add_si(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator+=(const int u)
{
- mpfr_add_si(mp,mp,u,mpreal::get_default_rnd());
+ mpfr_add_si(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
@@ -1154,9 +1194,9 @@ inline const mpreal mpreal::operator+()const { return mpreal(*this); } inline const mpreal operator+(const mpreal& a, const mpreal& b)
{
- mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));
- mpfr_add(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
- return c;
+ mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));
+ mpfr_add(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
+ return c;
}
inline mpreal& mpreal::operator++()
@@ -1187,21 +1227,21 @@ inline const mpreal mpreal::operator-- (int) // - Subtraction
inline mpreal& mpreal::operator-=(const mpreal& v)
{
- mpfr_sub(mp,mp,v.mp,mpreal::get_default_rnd());
+ mpfr_sub(mpfr_ptr(),mpfr_srcptr(),v.mpfr_srcptr(),mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator-=(const mpz_t v)
{
- mpfr_sub_z(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_sub_z(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator-=(const mpq_t v)
{
- mpfr_sub_q(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_sub_q(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
@@ -1216,7 +1256,7 @@ inline mpreal& mpreal::operator-=(const long double v) inline mpreal& mpreal::operator-=(const double v)
{
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
- mpfr_sub_d(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_sub_d(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
#else
*this -= mpreal(v);
#endif
@@ -1227,28 +1267,28 @@ inline mpreal& mpreal::operator-=(const double v) inline mpreal& mpreal::operator-=(const unsigned long int v)
{
- mpfr_sub_ui(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_sub_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator-=(const unsigned int v)
{
- mpfr_sub_ui(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_sub_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator-=(const long int v)
{
- mpfr_sub_si(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_sub_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator-=(const int v)
{
- mpfr_sub_si(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_sub_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
@@ -1256,15 +1296,15 @@ inline mpreal& mpreal::operator-=(const int v) inline const mpreal mpreal::operator-()const
{
mpreal u(*this);
- mpfr_neg(u.mp,u.mp,mpreal::get_default_rnd());
+ mpfr_neg(u.mpfr_ptr(),u.mpfr_srcptr(),mpreal::get_default_rnd());
return u;
}
inline const mpreal operator-(const mpreal& a, const mpreal& b)
{
- mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));
- mpfr_sub(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
- return c;
+ mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));
+ mpfr_sub(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
+ return c;
}
inline const mpreal operator-(const double b, const mpreal& a)
@@ -1312,21 +1352,21 @@ inline const mpreal operator-(const int b, const mpreal& a) // * Multiplication
inline mpreal& mpreal::operator*= (const mpreal& v)
{
- mpfr_mul(mp,mp,v.mp,mpreal::get_default_rnd());
+ mpfr_mul(mpfr_ptr(),mpfr_srcptr(),v.mpfr_srcptr(),mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator*=(const mpz_t v)
{
- mpfr_mul_z(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_mul_z(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator*=(const mpq_t v)
{
- mpfr_mul_q(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_mul_q(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
@@ -1341,7 +1381,7 @@ inline mpreal& mpreal::operator*=(const long double v) inline mpreal& mpreal::operator*=(const double v)
{
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
- mpfr_mul_d(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_mul_d(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
#else
*this *= mpreal(v);
#endif
@@ -1351,58 +1391,58 @@ inline mpreal& mpreal::operator*=(const double v) inline mpreal& mpreal::operator*=(const unsigned long int v)
{
- mpfr_mul_ui(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_mul_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator*=(const unsigned int v)
{
- mpfr_mul_ui(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_mul_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator*=(const long int v)
{
- mpfr_mul_si(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_mul_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator*=(const int v)
{
- mpfr_mul_si(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_mul_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline const mpreal operator*(const mpreal& a, const mpreal& b)
{
- mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));
- mpfr_mul(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
- return c;
+ mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));
+ mpfr_mul(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
+ return c;
}
//////////////////////////////////////////////////////////////////////////
// / Division
inline mpreal& mpreal::operator/=(const mpreal& v)
{
- mpfr_div(mp,mp,v.mp,mpreal::get_default_rnd());
+ mpfr_div(mpfr_ptr(),mpfr_srcptr(),v.mpfr_srcptr(),mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator/=(const mpz_t v)
{
- mpfr_div_z(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_div_z(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator/=(const mpq_t v)
{
- mpfr_div_q(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_div_q(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
@@ -1417,7 +1457,7 @@ inline mpreal& mpreal::operator/=(const long double v) inline mpreal& mpreal::operator/=(const double v)
{
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
- mpfr_div_d(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_div_d(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
#else
*this /= mpreal(v);
#endif
@@ -1427,72 +1467,72 @@ inline mpreal& mpreal::operator/=(const double v) inline mpreal& mpreal::operator/=(const unsigned long int v)
{
- mpfr_div_ui(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_div_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator/=(const unsigned int v)
{
- mpfr_div_ui(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_div_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator/=(const long int v)
{
- mpfr_div_si(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_div_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator/=(const int v)
{
- mpfr_div_si(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_div_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline const mpreal operator/(const mpreal& a, const mpreal& b)
{
- mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));
- mpfr_div(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
- return c;
+ mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_srcptr()), mpfr_get_prec(b.mpfr_srcptr())));
+ mpfr_div(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
+ return c;
}
inline const mpreal operator/(const unsigned long int b, const mpreal& a)
{
- mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));
+ mpreal x(0, mpfr_get_prec(a.mpfr_srcptr()));
mpfr_ui_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());
return x;
}
inline const mpreal operator/(const unsigned int b, const mpreal& a)
{
- mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));
+ mpreal x(0, mpfr_get_prec(a.mpfr_srcptr()));
mpfr_ui_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());
return x;
}
inline const mpreal operator/(const long int b, const mpreal& a)
{
- mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));
- mpfr_si_div(x.mpfr_ptr(), b, a.mpfr_srcptr(),mpreal::get_default_rnd());
+ mpreal x(0, mpfr_get_prec(a.mpfr_srcptr()));
+ mpfr_si_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());
return x;
}
inline const mpreal operator/(const int b, const mpreal& a)
{
- mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));
- mpfr_si_div(x.mpfr_ptr(), b, a.mpfr_srcptr(),mpreal::get_default_rnd());
+ mpreal x(0, mpfr_get_prec(a.mpfr_srcptr()));
+ mpfr_si_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());
return x;
}
inline const mpreal operator/(const double b, const mpreal& a)
{
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
- mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));
- mpfr_d_div(x.mpfr_ptr(), b, a.mpfr_srcptr(),mpreal::get_default_rnd());
+ mpreal x(0, mpfr_get_prec(a.mpfr_srcptr()));
+ mpfr_d_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());
return x;
#else
mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));
@@ -1505,56 +1545,56 @@ inline const mpreal operator/(const double b, const mpreal& a) // Shifts operators - Multiplication/Division by power of 2
inline mpreal& mpreal::operator<<=(const unsigned long int u)
{
- mpfr_mul_2ui(mp,mp,u,mpreal::get_default_rnd());
+ mpfr_mul_2ui(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator<<=(const unsigned int u)
{
- mpfr_mul_2ui(mp,mp,static_cast<unsigned long int>(u),mpreal::get_default_rnd());
+ mpfr_mul_2ui(mpfr_ptr(),mpfr_srcptr(),static_cast<unsigned long int>(u),mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator<<=(const long int u)
{
- mpfr_mul_2si(mp,mp,u,mpreal::get_default_rnd());
+ mpfr_mul_2si(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator<<=(const int u)
{
- mpfr_mul_2si(mp,mp,static_cast<long int>(u),mpreal::get_default_rnd());
+ mpfr_mul_2si(mpfr_ptr(),mpfr_srcptr(),static_cast<long int>(u),mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator>>=(const unsigned long int u)
{
- mpfr_div_2ui(mp,mp,u,mpreal::get_default_rnd());
+ mpfr_div_2ui(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator>>=(const unsigned int u)
{
- mpfr_div_2ui(mp,mp,static_cast<unsigned long int>(u),mpreal::get_default_rnd());
+ mpfr_div_2ui(mpfr_ptr(),mpfr_srcptr(),static_cast<unsigned long int>(u),mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator>>=(const long int u)
{
- mpfr_div_2si(mp,mp,u,mpreal::get_default_rnd());
+ mpfr_div_2si(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator>>=(const int u)
{
- mpfr_div_2si(mp,mp,static_cast<long int>(u),mpreal::get_default_rnd());
+ mpfr_div_2si(mpfr_ptr(),mpfr_srcptr(),static_cast<long int>(u),mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
@@ -1603,7 +1643,7 @@ inline const mpreal operator>>(const mpreal& v, const int k) inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode)
{
mpreal x(v);
- mpfr_mul_2ui(x.mp,v.mp,k,rnd_mode);
+ mpfr_mul_2ui(x.mpfr_ptr(),v.mpfr_srcptr(),k,rnd_mode);
return x;
}
@@ -1611,61 +1651,63 @@ inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_m inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode)
{
mpreal x(v);
- mpfr_mul_2si(x.mp,v.mp,k,rnd_mode);
+ mpfr_mul_2si(x.mpfr_ptr(),v.mpfr_srcptr(),k,rnd_mode);
return x;
}
inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode)
{
mpreal x(v);
- mpfr_div_2ui(x.mp,v.mp,k,rnd_mode);
+ mpfr_div_2ui(x.mpfr_ptr(),v.mpfr_srcptr(),k,rnd_mode);
return x;
}
inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode)
{
mpreal x(v);
- mpfr_div_2si(x.mp,v.mp,k,rnd_mode);
+ mpfr_div_2si(x.mpfr_ptr(),v.mpfr_srcptr(),k,rnd_mode);
return x;
}
//////////////////////////////////////////////////////////////////////////
//Boolean operators
-inline bool operator > (const mpreal& a, const mpreal& b){ return (mpfr_greater_p(a.mp,b.mp) !=0); }
-inline bool operator >= (const mpreal& a, const mpreal& b){ return (mpfr_greaterequal_p(a.mp,b.mp) !=0); }
-inline bool operator < (const mpreal& a, const mpreal& b){ return (mpfr_less_p(a.mp,b.mp) !=0); }
-inline bool operator <= (const mpreal& a, const mpreal& b){ return (mpfr_lessequal_p(a.mp,b.mp) !=0); }
-inline bool operator == (const mpreal& a, const mpreal& b){ return (mpfr_equal_p(a.mp,b.mp) !=0); }
-inline bool operator != (const mpreal& a, const mpreal& b){ return (mpfr_lessgreater_p(a.mp,b.mp) !=0); }
-
-inline bool operator == (const mpreal& a, const unsigned long int b ){ return (mpfr_cmp_ui(a.mp,b) == 0); }
-inline bool operator == (const mpreal& a, const unsigned int b ){ return (mpfr_cmp_ui(a.mp,b) == 0); }
-inline bool operator == (const mpreal& a, const long int b ){ return (mpfr_cmp_si(a.mp,b) == 0); }
-inline bool operator == (const mpreal& a, const int b ){ return (mpfr_cmp_si(a.mp,b) == 0); }
-inline bool operator == (const mpreal& a, const long double b ){ return (mpfr_cmp_ld(a.mp,b) == 0); }
-inline bool operator == (const mpreal& a, const double b ){ return (mpfr_cmp_d(a.mp,b) == 0); }
-
-
-inline bool isnan (const mpreal& v){ return (mpfr_nan_p(v.mp) != 0); }
-inline bool isinf (const mpreal& v){ return (mpfr_inf_p(v.mp) != 0); }
-inline bool isfinite (const mpreal& v){ return (mpfr_number_p(v.mp) != 0); }
-inline bool iszero (const mpreal& v){ return (mpfr_zero_p(v.mp) != 0); }
-inline bool isint (const mpreal& v){ return (mpfr_integer_p(v.mp) != 0); }
+inline bool operator > (const mpreal& a, const mpreal& b){ return (mpfr_greater_p (a.mpfr_srcptr(),b.mpfr_srcptr()) !=0 ); }
+inline bool operator >= (const mpreal& a, const mpreal& b){ return (mpfr_greaterequal_p (a.mpfr_srcptr(),b.mpfr_srcptr()) !=0 ); }
+inline bool operator < (const mpreal& a, const mpreal& b){ return (mpfr_less_p (a.mpfr_srcptr(),b.mpfr_srcptr()) !=0 ); }
+inline bool operator <= (const mpreal& a, const mpreal& b){ return (mpfr_lessequal_p (a.mpfr_srcptr(),b.mpfr_srcptr()) !=0 ); }
+inline bool operator == (const mpreal& a, const mpreal& b){ return (mpfr_equal_p (a.mpfr_srcptr(),b.mpfr_srcptr()) !=0 ); }
+inline bool operator != (const mpreal& a, const mpreal& b){ return (mpfr_lessgreater_p (a.mpfr_srcptr(),b.mpfr_srcptr()) !=0 ); }
+
+inline bool operator == (const mpreal& a, const unsigned long int b ){ return (mpfr_cmp_ui(a.mpfr_srcptr(),b) == 0 ); }
+inline bool operator == (const mpreal& a, const unsigned int b ){ return (mpfr_cmp_ui(a.mpfr_srcptr(),b) == 0 ); }
+inline bool operator == (const mpreal& a, const long int b ){ return (mpfr_cmp_si(a.mpfr_srcptr(),b) == 0 ); }
+inline bool operator == (const mpreal& a, const int b ){ return (mpfr_cmp_si(a.mpfr_srcptr(),b) == 0 ); }
+inline bool operator == (const mpreal& a, const long double b ){ return (mpfr_cmp_ld(a.mpfr_srcptr(),b) == 0 ); }
+inline bool operator == (const mpreal& a, const double b ){ return (mpfr_cmp_d (a.mpfr_srcptr(),b) == 0 ); }
+
+
+inline bool isnan (const mpreal& op){ return (mpfr_nan_p (op.mpfr_srcptr()) != 0 ); }
+inline bool isinf (const mpreal& op){ return (mpfr_inf_p (op.mpfr_srcptr()) != 0 ); }
+inline bool isfinite (const mpreal& op){ return (mpfr_number_p (op.mpfr_srcptr()) != 0 ); }
+inline bool iszero (const mpreal& op){ return (mpfr_zero_p (op.mpfr_srcptr()) != 0 ); }
+inline bool isint (const mpreal& op){ return (mpfr_integer_p(op.mpfr_srcptr()) != 0 ); }
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
-inline bool isregular(const mpreal& v){ return (mpfr_regular_p(v.mp));}
+inline bool isregular(const mpreal& op){ return (mpfr_regular_p(op.mpfr_srcptr()));}
#endif
//////////////////////////////////////////////////////////////////////////
// Type Converters
-inline long mpreal::toLong (mp_rnd_t mode) const { return mpfr_get_si(mp, mode); }
-inline unsigned long mpreal::toULong (mp_rnd_t mode) const { return mpfr_get_ui(mp, mode); }
-inline double mpreal::toDouble (mp_rnd_t mode) const { return mpfr_get_d (mp, mode); }
-inline long double mpreal::toLDouble(mp_rnd_t mode) const { return mpfr_get_ld(mp, mode); }
+inline bool mpreal::toBool (mp_rnd_t mode) const { return mpfr_zero_p (mpfr_srcptr()) == 0; }
+inline long mpreal::toLong (mp_rnd_t mode) const { return mpfr_get_si (mpfr_srcptr(), mode); }
+inline unsigned long mpreal::toULong (mp_rnd_t mode) const { return mpfr_get_ui (mpfr_srcptr(), mode); }
+inline float mpreal::toFloat (mp_rnd_t mode) const { return mpfr_get_flt(mpfr_srcptr(), mode); }
+inline double mpreal::toDouble (mp_rnd_t mode) const { return mpfr_get_d (mpfr_srcptr(), mode); }
+inline long double mpreal::toLDouble(mp_rnd_t mode) const { return mpfr_get_ld (mpfr_srcptr(), mode); }
#if defined (MPREAL_HAVE_INT64_SUPPORT)
-inline int64_t mpreal::toInt64 (mp_rnd_t mode) const{ return mpfr_get_sj(mp, mode); }
-inline uint64_t mpreal::toUInt64(mp_rnd_t mode) const{ return mpfr_get_uj(mp, mode); }
+inline int64_t mpreal::toInt64 (mp_rnd_t mode) const{ return mpfr_get_sj(mpfr_srcptr(), mode); }
+inline uint64_t mpreal::toUInt64(mp_rnd_t mode) const{ return mpfr_get_uj(mpfr_srcptr(), mode); }
#endif
inline ::mpfr_ptr mpreal::mpfr_ptr() { return mp; }
@@ -1689,7 +1731,7 @@ inline std::string mpreal::toString(const std::string& format) const if( !format.empty() )
{
- if(!(mpfr_asprintf(&s,format.c_str(),mp) < 0))
+ if(!(mpfr_asprintf(&s, format.c_str(), mpfr_srcptr()) < 0))
{
out = std::string(s);
@@ -1704,20 +1746,19 @@ inline std::string mpreal::toString(const std::string& format) const inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const
{
+ // TODO: Add extended format specification (f, e, rounding mode) as it done in output operator
(void)b;
(void)mode;
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
- // Use MPFR native function for output
- char format[128];
- int digits;
+ std::ostringstream format;
- digits = n > 0 ? n : bits2digits(mpfr_get_prec(mp));
-
- sprintf(format,"%%.%dRNg",digits); // Default format
+ int digits = (n >= 0) ? n : bits2digits(mpfr_get_prec(mpfr_srcptr()));
+
+ format << "%." << digits << "RNg";
- return toString(std::string(format));
+ return toString(format.str());
#else
@@ -1735,8 +1776,8 @@ inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const if(mpfr_zero_p(mp)) return "0";
if(mpfr_nan_p(mp)) return "NaN";
- s = mpfr_get_str(NULL,&exp,b,0,mp,mode);
- ns = mpfr_get_str(NULL,&exp,b,n,mp,mode);
+ s = mpfr_get_str(NULL, &exp, b, 0, mp, mode);
+ ns = mpfr_get_str(NULL, &exp, b, (std::max)(0,n), mp, mode);
if(s!=NULL && ns!=NULL)
{
@@ -1821,17 +1862,43 @@ inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const //////////////////////////////////////////////////////////////////////////
// I/O
+inline std::ostream& mpreal::output(std::ostream& os) const
+{
+ std::ostringstream format;
+ const std::ios::fmtflags flags = os.flags();
+
+ format << ((flags & std::ios::showpos) ? "%+" : "%");
+ if (os.precision() >= 0)
+ format << '.' << os.precision() << "R*"
+ << ((flags & std::ios::floatfield) == std::ios::fixed ? 'f' :
+ (flags & std::ios::floatfield) == std::ios::scientific ? 'e' :
+ 'g');
+ else
+ format << "R*e";
+
+ char *s = NULL;
+ if(!(mpfr_asprintf(&s, format.str().c_str(),
+ mpfr::mpreal::get_default_rnd(),
+ mpfr_srcptr())
+ < 0))
+ {
+ os << std::string(s);
+ mpfr_free_str(s);
+ }
+ return os;
+}
+
inline std::ostream& operator<<(std::ostream& os, const mpreal& v)
{
- return os<<v.toString(static_cast<int>(os.precision()));
+ return v.output(os);
}
inline std::istream& operator>>(std::istream &is, mpreal& v)
{
- // ToDo, use cout::hexfloat and other flags to setup base
+ // TODO: use cout::hexfloat and other flags to setup base
std::string tmp;
is >> tmp;
- mpfr_set_str(v.mp, tmp.c_str(), 10, mpreal::get_default_rnd());
+ mpfr_set_str(v.mpfr_ptr(), tmp.c_str(), 10, mpreal::get_default_rnd());
return is;
}
@@ -1844,53 +1911,53 @@ inline mp_prec_t digits2bits(int d) {
const double LOG2_10 = 3.3219280948873624;
- return (mp_prec_t) std::ceil( d * LOG2_10 );
+ return mp_prec_t(std::ceil( d * LOG2_10 ));
}
inline int bits2digits(mp_prec_t b)
{
const double LOG10_2 = 0.30102999566398119;
- return (int) std::floor( b * LOG10_2 );
+ return int(std::floor( b * LOG10_2 ));
}
//////////////////////////////////////////////////////////////////////////
// Set/Get number properties
-inline int sgn(const mpreal& v)
+inline int sgn(const mpreal& op)
{
- int r = mpfr_signbit(v.mp);
- return (r>0?-1:1);
+ int r = mpfr_signbit(op.mpfr_srcptr());
+ return (r > 0? -1 : 1);
}
inline mpreal& mpreal::setSign(int sign, mp_rnd_t RoundingMode)
{
- mpfr_setsign(mp,mp,(sign < 0 ? 1 : 0),RoundingMode);
+ mpfr_setsign(mpfr_ptr(), mpfr_srcptr(), (sign < 0 ? 1 : 0), RoundingMode);
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline int mpreal::getPrecision() const
{
- return mpfr_get_prec(mp);
+ return int(mpfr_get_prec(mpfr_srcptr()));
}
inline mpreal& mpreal::setPrecision(int Precision, mp_rnd_t RoundingMode)
{
- mpfr_prec_round(mp, Precision, RoundingMode);
+ mpfr_prec_round(mpfr_ptr(), Precision, RoundingMode);
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::setInf(int sign)
{
- mpfr_set_inf(mp,sign);
+ mpfr_set_inf(mpfr_ptr(), sign);
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::setNan()
{
- mpfr_set_nan(mp);
+ mpfr_set_nan(mpfr_ptr());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
@@ -1899,9 +1966,9 @@ inline mpreal& mpreal::setZero(int sign) {
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
- mpfr_set_zero(mp, sign);
+ mpfr_set_zero(mpfr_ptr(), sign);
#else
- mpfr_set_si(mp, 0, (mpfr_get_default_rounding_mode)());
+ mpfr_set_si(mpfr_ptr(), 0, (mpfr_get_default_rounding_mode)());
setSign(sign);
#endif
@@ -1911,23 +1978,23 @@ inline mpreal& mpreal::setZero(int sign) inline mp_prec_t mpreal::get_prec() const
{
- return mpfr_get_prec(mp);
+ return mpfr_get_prec(mpfr_srcptr());
}
inline void mpreal::set_prec(mp_prec_t prec, mp_rnd_t rnd_mode)
{
- mpfr_prec_round(mp,prec,rnd_mode);
+ mpfr_prec_round(mpfr_ptr(),prec,rnd_mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mp_exp_t mpreal::get_exp ()
{
- return mpfr_get_exp(mp);
+ return mpfr_get_exp(mpfr_srcptr());
}
inline int mpreal::set_exp (mp_exp_t e)
{
- int x = mpfr_set_exp(mp, e);
+ int x = mpfr_set_exp(mpfr_ptr(), e);
MPREAL_MSVC_DEBUGVIEW_CODE;
return x;
}
@@ -1945,7 +2012,7 @@ inline const mpreal ldexp(const mpreal& v, mp_exp_t exp) mpreal x(v);
// rounding is not important since we just increasing the exponent
- mpfr_mul_2si(x.mp,x.mp,exp,mpreal::get_default_rnd());
+ mpfr_mul_2si(x.mpfr_ptr(), x.mpfr_srcptr(), exp, mpreal::get_default_rnd());
return x;
}
@@ -1960,9 +2027,9 @@ inline mpreal machine_epsilon(const mpreal& x) /* the smallest eps such that x + eps != x */
if( x < 0)
{
- return nextabove(-x)+x;
+ return nextabove(-x) + x;
}else{
- return nextabove(x)-x;
+ return nextabove( x) - x;
}
}
@@ -1997,22 +2064,22 @@ inline bool isEqualFuzzy(const mpreal& a, const mpreal& b) inline const mpreal modf(const mpreal& v, mpreal& n)
{
- mpreal frac(v);
+ mpreal f(v);
// rounding is not important since we are using the same number
- mpfr_frac(frac.mp,frac.mp,mpreal::get_default_rnd());
- mpfr_trunc(n.mp,v.mp);
- return frac;
+ mpfr_frac (f.mpfr_ptr(),f.mpfr_srcptr(),mpreal::get_default_rnd());
+ mpfr_trunc(n.mpfr_ptr(),v.mpfr_srcptr());
+ return f;
}
inline int mpreal::check_range (int t, mp_rnd_t rnd_mode)
{
- return mpfr_check_range(mp,t,rnd_mode);
+ return mpfr_check_range(mpfr_ptr(),t,rnd_mode);
}
inline int mpreal::subnormalize (int t,mp_rnd_t rnd_mode)
{
- int r = mpfr_subnormalize(mp,t,rnd_mode);
+ int r = mpfr_subnormalize(mpfr_ptr(),t,rnd_mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
return r;
}
@@ -2065,8 +2132,11 @@ inline mp_exp_t mpreal::get_emax_max (void) mpfr_##f(y.mpfr_ptr(), x.mpfr_srcptr(), r); \
return y;
-inline const mpreal sqr (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(sqr ); }
-inline const mpreal sqrt (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(sqrt); }
+inline const mpreal sqr (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd())
+{ MPREAL_UNARY_MATH_FUNCTION_BODY(sqr ); }
+
+inline const mpreal sqrt (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd())
+{ MPREAL_UNARY_MATH_FUNCTION_BODY(sqrt); }
inline const mpreal sqrt(const unsigned long int x, mp_rnd_t r)
{
@@ -2092,14 +2162,14 @@ inline const mpreal sqrt(const int v, mp_rnd_t rnd_mode) else return mpreal().setNan(); // NaN
}
-inline const mpreal root(const mpreal& x, unsigned long int k, mp_rnd_t r)
+inline const mpreal root(const mpreal& x, unsigned long int k, mp_rnd_t r = mpreal::get_default_rnd())
{
mpreal y(0, mpfr_get_prec(x.mpfr_srcptr()));
mpfr_root(y.mpfr_ptr(), x.mpfr_srcptr(), k, r);
return y;
}
-inline const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t r)
+inline const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t r = mpreal::get_default_rnd())
{
mpreal y(0, mpfr_get_prec(a.mpfr_srcptr()));
mpfr_dim(y.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), r);
@@ -2111,7 +2181,7 @@ inline int cmpabs(const mpreal& a,const mpreal& b) return mpfr_cmpabs(a.mpfr_ptr(), b.mpfr_srcptr());
}
-inline int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode)
+inline int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
return mpfr_sin_cos(s.mpfr_ptr(), c.mpfr_ptr(), v.mpfr_srcptr(), rnd_mode);
}
@@ -2119,84 +2189,85 @@ inline int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode) inline const mpreal sqrt (const long double v, mp_rnd_t rnd_mode) { return sqrt(mpreal(v),rnd_mode); }
inline const mpreal sqrt (const double v, mp_rnd_t rnd_mode) { return sqrt(mpreal(v),rnd_mode); }
-inline const mpreal cbrt (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(cbrt ); }
-inline const mpreal fabs (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(abs ); }
-inline const mpreal abs (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(abs ); }
-inline const mpreal log (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(log ); }
-inline const mpreal log2 (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(log2 ); }
-inline const mpreal log10 (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(log10); }
-inline const mpreal exp (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(exp ); }
-inline const mpreal exp2 (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(exp2 ); }
-inline const mpreal exp10 (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(exp10); }
-inline const mpreal cos (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(cos ); }
-inline const mpreal sin (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(sin ); }
-inline const mpreal tan (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(tan ); }
-inline const mpreal sec (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(sec ); }
-inline const mpreal csc (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(csc ); }
-inline const mpreal cot (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(cot ); }
-inline const mpreal acos (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(acos ); }
-inline const mpreal asin (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(asin ); }
-inline const mpreal atan (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(atan ); }
-
-inline const mpreal acot (const mpreal& v, mp_rnd_t r) { return atan (1/v, r); }
-inline const mpreal asec (const mpreal& v, mp_rnd_t r) { return acos (1/v, r); }
-inline const mpreal acsc (const mpreal& v, mp_rnd_t r) { return asin (1/v, r); }
-inline const mpreal acoth (const mpreal& v, mp_rnd_t r) { return atanh(1/v, r); }
-inline const mpreal asech (const mpreal& v, mp_rnd_t r) { return acosh(1/v, r); }
-inline const mpreal acsch (const mpreal& v, mp_rnd_t r) { return asinh(1/v, r); }
-
-inline const mpreal cosh (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(cosh ); }
-inline const mpreal sinh (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(sinh ); }
-inline const mpreal tanh (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(tanh ); }
-inline const mpreal sech (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(sech ); }
-inline const mpreal csch (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(csch ); }
-inline const mpreal coth (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(coth ); }
-inline const mpreal acosh (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(acosh); }
-inline const mpreal asinh (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(asinh); }
-inline const mpreal atanh (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(atanh); }
-
-inline const mpreal log1p (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(log1p ); }
-inline const mpreal expm1 (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(expm1 ); }
-inline const mpreal eint (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(eint ); }
-inline const mpreal gamma (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(gamma ); }
-inline const mpreal lngamma (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(lngamma); }
-inline const mpreal zeta (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(zeta ); }
-inline const mpreal erf (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(erf ); }
-inline const mpreal erfc (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(erfc ); }
-inline const mpreal besselj0(const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(j0 ); }
-inline const mpreal besselj1(const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(j1 ); }
-inline const mpreal bessely0(const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(y0 ); }
-inline const mpreal bessely1(const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(y1 ); }
-
-inline const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode)
+inline const mpreal cbrt (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(cbrt ); }
+inline const mpreal fabs (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(abs ); }
+inline const mpreal abs (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(abs ); }
+inline const mpreal log (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(log ); }
+inline const mpreal log2 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(log2 ); }
+inline const mpreal log10 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(log10); }
+inline const mpreal exp (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(exp ); }
+inline const mpreal exp2 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(exp2 ); }
+inline const mpreal exp10 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(exp10); }
+inline const mpreal cos (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(cos ); }
+inline const mpreal sin (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(sin ); }
+inline const mpreal tan (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(tan ); }
+inline const mpreal sec (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(sec ); }
+inline const mpreal csc (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(csc ); }
+inline const mpreal cot (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(cot ); }
+inline const mpreal acos (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(acos ); }
+inline const mpreal asin (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(asin ); }
+inline const mpreal atan (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(atan ); }
+
+inline const mpreal acot (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return atan (1/v, r); }
+inline const mpreal asec (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return acos (1/v, r); }
+inline const mpreal acsc (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return asin (1/v, r); }
+inline const mpreal acoth (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return atanh(1/v, r); }
+inline const mpreal asech (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return acosh(1/v, r); }
+inline const mpreal acsch (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return asinh(1/v, r); }
+
+inline const mpreal cosh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(cosh ); }
+inline const mpreal sinh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(sinh ); }
+inline const mpreal tanh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(tanh ); }
+inline const mpreal sech (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(sech ); }
+inline const mpreal csch (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(csch ); }
+inline const mpreal coth (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(coth ); }
+inline const mpreal acosh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(acosh); }
+inline const mpreal asinh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(asinh); }
+inline const mpreal atanh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(atanh); }
+
+inline const mpreal log1p (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(log1p ); }
+inline const mpreal expm1 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(expm1 ); }
+inline const mpreal eint (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(eint ); }
+inline const mpreal gamma (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(gamma ); }
+inline const mpreal lngamma (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(lngamma); }
+inline const mpreal zeta (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(zeta ); }
+inline const mpreal erf (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(erf ); }
+inline const mpreal erfc (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(erfc ); }
+inline const mpreal besselj0(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(j0 ); }
+inline const mpreal besselj1(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(j1 ); }
+inline const mpreal bessely0(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(y0 ); }
+inline const mpreal bessely1(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(y1 ); }
+
+inline const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal a(0,(std::max)(y.getPrecision(), x.getPrecision()));
mpfr_atan2(a.mpfr_ptr(), y.mpfr_srcptr(), x.mpfr_srcptr(), rnd_mode);
return a;
}
-inline const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
+inline const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal a(0,(std::max)(y.getPrecision(), x.getPrecision()));
mpfr_hypot(a.mpfr_ptr(), x.mpfr_srcptr(), y.mpfr_srcptr(), rnd_mode);
return a;
}
-inline const mpreal remainder (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
+inline const mpreal remainder (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal a(0,(std::max)(y.getPrecision(), x.getPrecision()));
mpfr_remainder(a.mpfr_ptr(), x.mpfr_srcptr(), y.mpfr_srcptr(), rnd_mode);
return a;
}
-inline const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
+inline const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal a(0,(std::max)(y.getPrecision(), x.getPrecision()));
mpfr_remquo(a.mpfr_ptr(),q, x.mpfr_srcptr(), y.mpfr_srcptr(), rnd_mode);
return a;
}
-inline const mpreal fac_ui (unsigned long int v, mp_prec_t prec, mp_rnd_t rnd_mode)
+inline const mpreal fac_ui (unsigned long int v, mp_prec_t prec = mpreal::get_default_prec(),
+ mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal x(0, prec);
mpfr_fac_ui(x.mpfr_ptr(),v,rnd_mode);
@@ -2204,33 +2275,33 @@ inline const mpreal fac_ui (unsigned long int v, mp_prec_t prec, mp_rnd_t rnd_mo }
-inline const mpreal lgamma (const mpreal& v, int *signp, mp_rnd_t rnd_mode)
+inline const mpreal lgamma (const mpreal& v, int *signp = 0, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal x(v);
int tsignp;
- if(signp) mpfr_lgamma(x.mp, signp,v.mp,rnd_mode);
- else mpfr_lgamma(x.mp,&tsignp,v.mp,rnd_mode);
+ if(signp) mpfr_lgamma(x.mpfr_ptr(), signp,v.mpfr_srcptr(),rnd_mode);
+ else mpfr_lgamma(x.mpfr_ptr(),&tsignp,v.mpfr_srcptr(),rnd_mode);
return x;
}
-inline const mpreal besseljn (long n, const mpreal& x, mp_rnd_t r)
+inline const mpreal besseljn (long n, const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd())
{
- mpreal y(0, mpfr_get_prec(x.mpfr_srcptr()));
+ mpreal y(0, x.getPrecision());
mpfr_jn(y.mpfr_ptr(), n, x.mpfr_srcptr(), r);
return y;
}
-inline const mpreal besselyn (long n, const mpreal& x, mp_rnd_t r)
+inline const mpreal besselyn (long n, const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd())
{
- mpreal y(0, mpfr_get_prec(x.mpfr_srcptr()));
+ mpreal y(0, x.getPrecision());
mpfr_yn(y.mpfr_ptr(), n, x.mpfr_srcptr(), r);
return y;
}
-inline const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode)
+inline const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal a;
mp_prec_t p1, p2, p3;
@@ -2245,7 +2316,7 @@ inline const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, m return a;
}
-inline const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode)
+inline const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal a;
mp_prec_t p1, p2, p3;
@@ -2260,7 +2331,7 @@ inline const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, m return a;
}
-inline const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode)
+inline const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal a;
mp_prec_t p1, p2;
@@ -2275,7 +2346,7 @@ inline const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode) return a;
}
-inline const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode)
+inline const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal x;
mpfr_ptr* t;
@@ -2292,23 +2363,23 @@ inline const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_m // MPFR 2.4.0 Specifics
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
-inline int sinh_cosh(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode)
+inline int sinh_cosh(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
return mpfr_sinh_cosh(s.mp,c.mp,v.mp,rnd_mode);
}
-inline const mpreal li2 (const mpreal& x, mp_rnd_t r)
+inline const mpreal li2 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd())
{
MPREAL_UNARY_MATH_FUNCTION_BODY(li2);
}
-inline const mpreal rem (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
+inline const mpreal rem (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
/* R = rem(X,Y) if Y != 0, returns X - n * Y where n = trunc(X/Y). */
return fmod(x, y, rnd_mode);
}
-inline const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
+inline const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
(void)rnd_mode;
@@ -2333,7 +2404,7 @@ inline const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode) return m;
}
-inline const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
+inline const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal a;
mp_prec_t yp, xp;
@@ -2348,7 +2419,7 @@ inline const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode) return a;
}
-inline const mpreal rec_sqrt(const mpreal& v, mp_rnd_t rnd_mode)
+inline const mpreal rec_sqrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal x(v);
mpfr_rec_sqrt(x.mp,v.mp,rnd_mode);
@@ -2359,41 +2430,41 @@ inline const mpreal rec_sqrt(const mpreal& v, mp_rnd_t rnd_mode) //////////////////////////////////////////////////////////////////////////
// MPFR 3.0.0 Specifics
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
-inline const mpreal digamma (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(digamma); }
-inline const mpreal ai (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(ai); }
+inline const mpreal digamma (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(digamma); }
+inline const mpreal ai (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(ai); }
#endif // MPFR 3.0.0 Specifics
//////////////////////////////////////////////////////////////////////////
// Constants
-inline const mpreal const_log2 (mp_prec_t p, mp_rnd_t r)
+inline const mpreal const_log2 (mp_prec_t p = mpreal::get_default_prec(), mp_rnd_t r = mpreal::get_default_rnd())
{
mpreal x(0, p);
mpfr_const_log2(x.mpfr_ptr(), r);
return x;
}
-inline const mpreal const_pi (mp_prec_t p, mp_rnd_t r)
+inline const mpreal const_pi (mp_prec_t p = mpreal::get_default_prec(), mp_rnd_t r = mpreal::get_default_rnd())
{
mpreal x(0, p);
mpfr_const_pi(x.mpfr_ptr(), r);
return x;
}
-inline const mpreal const_euler (mp_prec_t p, mp_rnd_t r)
+inline const mpreal const_euler (mp_prec_t p = mpreal::get_default_prec(), mp_rnd_t r = mpreal::get_default_rnd())
{
mpreal x(0, p);
mpfr_const_euler(x.mpfr_ptr(), r);
return x;
}
-inline const mpreal const_catalan (mp_prec_t p, mp_rnd_t r)
+inline const mpreal const_catalan (mp_prec_t p = mpreal::get_default_prec(), mp_rnd_t r = mpreal::get_default_rnd())
{
mpreal x(0, p);
mpfr_const_catalan(x.mpfr_ptr(), r);
return x;
}
-inline const mpreal const_infinity (int sign, mp_prec_t p, mp_rnd_t /*r*/)
+inline const mpreal const_infinity (int sign = 1, mp_prec_t p = mpreal::get_default_prec())
{
mpreal x(0, p);
mpfr_set_inf(x.mpfr_ptr(), sign);
@@ -2430,12 +2501,12 @@ inline const mpreal trunc(const mpreal& v) return x;
}
-inline const mpreal rint (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint ); }
-inline const mpreal rint_ceil (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint_ceil ); }
-inline const mpreal rint_floor (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint_floor); }
-inline const mpreal rint_round (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint_round); }
-inline const mpreal rint_trunc (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint_trunc); }
-inline const mpreal frac (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(frac ); }
+inline const mpreal rint (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint ); }
+inline const mpreal rint_ceil (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint_ceil ); }
+inline const mpreal rint_floor (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint_floor); }
+inline const mpreal rint_round (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint_round); }
+inline const mpreal rint_trunc (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint_trunc); }
+inline const mpreal frac (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(frac ); }
//////////////////////////////////////////////////////////////////////////
// Miscellaneous Functions
@@ -2443,14 +2514,14 @@ inline void swap (mpreal& a, mpreal& b) { mpfr_swap(a.mp,b inline const mpreal (max)(const mpreal& x, const mpreal& y){ return (x>y?x:y); }
inline const mpreal (min)(const mpreal& x, const mpreal& y){ return (x<y?x:y); }
-inline const mpreal fmax(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
+inline const mpreal fmax(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal a;
mpfr_max(a.mp,x.mp,y.mp,rnd_mode);
return a;
}
-inline const mpreal fmin(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
+inline const mpreal fmin(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal a;
mpfr_min(a.mp,x.mp,y.mp,rnd_mode);
@@ -2485,16 +2556,16 @@ inline const mpreal urandomb (gmp_randstate_t& state) return x;
}
-#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
+#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,1,0))
// use gmp_randinit_default() to init state, gmp_randclear() to clear
-inline const mpreal urandom (gmp_randstate_t& state, mp_rnd_t rnd_mode)
+inline const mpreal urandom (gmp_randstate_t& state, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal x;
mpfr_urandom(x.mp,state,rnd_mode);
return x;
}
-inline const mpreal grandom (gmp_randstate_t& state, mp_rnd_t rnd_mode)
+inline const mpreal grandom (gmp_randstate_t& state, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal x;
mpfr_grandom(x.mp, NULL, state, rnd_mode);
@@ -2516,7 +2587,7 @@ inline const mpreal random2 (mp_size_t size, mp_exp_t exp) // a = random(seed); <- initialization & first random number generation
// a = random(); <- next random numbers generation
// seed != 0
-inline const mpreal random(unsigned int seed)
+inline const mpreal random(unsigned int seed = 0)
{
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
@@ -2541,7 +2612,7 @@ inline const mpreal random(unsigned int seed) }
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
-inline const mpreal grandom(unsigned int seed)
+inline const mpreal grandom(unsigned int seed = 0)
{
static gmp_randstate_t state;
static bool isFirstTime = true;
@@ -2578,21 +2649,21 @@ inline bool mpreal::fits_in_bits(double x, int n) return IsInf(x) || (std::modf ( std::ldexp ( std::frexp ( x, &i ), n ), &t ) == 0.0);
}
-inline const mpreal pow(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode)
+inline const mpreal pow(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal x(a);
mpfr_pow(x.mp,x.mp,b.mp,rnd_mode);
return x;
}
-inline const mpreal pow(const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode)
+inline const mpreal pow(const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal x(a);
mpfr_pow_z(x.mp,x.mp,b,rnd_mode);
return x;
}
-inline const mpreal pow(const mpreal& a, const unsigned long int b, mp_rnd_t rnd_mode)
+inline const mpreal pow(const mpreal& a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal x(a);
mpfr_pow_ui(x.mp,x.mp,b,rnd_mode);
@@ -2604,7 +2675,7 @@ inline const mpreal pow(const mpreal& a, const unsigned int b, mp_rnd_t rnd_mode return pow(a,static_cast<unsigned long int>(b),rnd_mode);
}
-inline const mpreal pow(const mpreal& a, const long int b, mp_rnd_t rnd_mode)
+inline const mpreal pow(const mpreal& a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal x(a);
mpfr_pow_si(x.mp,x.mp,b,rnd_mode);
@@ -2626,7 +2697,7 @@ inline const mpreal pow(const mpreal& a, const double b, mp_rnd_t rnd_mode) return pow(a,mpreal(b),rnd_mode);
}
-inline const mpreal pow(const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode)
+inline const mpreal pow(const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal x(a);
mpfr_ui_pow(x.mp,a,b.mp,rnd_mode);
@@ -2641,13 +2712,13 @@ inline const mpreal pow(const unsigned int a, const mpreal& b, mp_rnd_t rnd_mode inline const mpreal pow(const long int a, const mpreal& b, mp_rnd_t rnd_mode)
{
if (a>=0) return pow(static_cast<unsigned long int>(a),b,rnd_mode);
- else return pow(mpreal(a),b,rnd_mode);
+ else return pow(mpreal(a),b,rnd_mode);
}
inline const mpreal pow(const int a, const mpreal& b, mp_rnd_t rnd_mode)
{
if (a>=0) return pow(static_cast<unsigned long int>(a),b,rnd_mode);
- else return pow(mpreal(a),b,rnd_mode);
+ else return pow(mpreal(a),b,rnd_mode);
}
inline const mpreal pow(const long double a, const mpreal& b, mp_rnd_t rnd_mode)
@@ -2676,13 +2747,13 @@ inline const mpreal pow(const unsigned long int a, const unsigned int b, mp_rnd_ inline const mpreal pow(const unsigned long int a, const long int b, mp_rnd_t rnd_mode)
{
if(b>0) return pow(a,static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui
- else return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow
+ else return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow
}
inline const mpreal pow(const unsigned long int a, const int b, mp_rnd_t rnd_mode)
{
if(b>0) return pow(a,static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui
- else return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow
+ else return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow
}
inline const mpreal pow(const unsigned long int a, const long double b, mp_rnd_t rnd_mode)
@@ -2879,7 +2950,7 @@ inline const mpreal pow(const double a, const int b, mp_rnd_t rnd_mode) // Non-throwing swap C++ idiom: http://en.wikibooks.org/wiki/More_C%2B%2B_Idioms/Non-throwing_swap
namespace std
{
- // only allowed to extend namespace std with specializations
+ // we are allowed to extend namespace std with specializations only
template <>
inline void swap(mpfr::mpreal& x, mpfr::mpreal& y)
{
@@ -2907,20 +2978,6 @@ namespace std static const bool tinyness_before = true;
static const float_denorm_style has_denorm = denorm_absent;
-
- inline static float_round_style round_style()
- {
- mp_rnd_t r = mpfr::mpreal::get_default_rnd();
-
- switch (r)
- {
- case GMP_RNDN: return round_to_nearest;
- case GMP_RNDZ: return round_toward_zero;
- case GMP_RNDU: return round_toward_infinity;
- case GMP_RNDD: return round_toward_neg_infinity;
- default: return round_indeterminate;
- }
- }
inline static mpfr::mpreal (min) (mp_prec_t precision = mpfr::mpreal::get_default_prec()) { return mpfr::minval(precision); }
inline static mpfr::mpreal (max) (mp_prec_t precision = mpfr::mpreal::get_default_prec()) { return mpfr::maxval(precision); }
@@ -2928,7 +2985,7 @@ namespace std // Returns smallest eps such that 1 + eps != 1 (classic machine epsilon)
inline static mpfr::mpreal epsilon(mp_prec_t precision = mpfr::mpreal::get_default_prec()) { return mpfr::machine_epsilon(precision); }
-
+
// Returns smallest eps such that x + eps != x (relative machine epsilon)
inline static mpfr::mpreal epsilon(const mpfr::mpreal& x) { return mpfr::machine_epsilon(x); }
@@ -2951,10 +3008,30 @@ namespace std MPREAL_PERMISSIVE_EXPR static const int min_exponent10 = (int) (MPFR_EMIN_DEFAULT * 0.3010299956639811);
MPREAL_PERMISSIVE_EXPR static const int max_exponent10 = (int) (MPFR_EMAX_DEFAULT * 0.3010299956639811);
- // Should be constant according to standard, but 'digits' depends on precision in MPFR
+#ifdef MPREAL_HAVE_DYNAMIC_STD_NUMERIC_LIMITS
+
+ // Following members should be constant according to standard, but they can be variable in MPFR
+ // So we define them as functions here.
+ //
+ // This is preferable way for std::numeric_limits<mpfr::mpreal> specialization.
+ // But it is incompatible with standard std::numeric_limits and might not work with other libraries, e.g. boost.
+ // See below for compatible implementation.
+ inline static float_round_style round_style()
+ {
+ mp_rnd_t r = mpfr::mpreal::get_default_rnd();
+
+ switch (r)
+ {
+ case GMP_RNDN: return round_to_nearest;
+ case GMP_RNDZ: return round_toward_zero;
+ case GMP_RNDU: return round_toward_infinity;
+ case GMP_RNDD: return round_toward_neg_infinity;
+ default: return round_indeterminate;
+ }
+ }
- inline static int digits() { return mpfr::mpreal::get_default_prec(); }
- inline static int digits(const mpfr::mpreal& x) { return x.getPrecision(); }
+ inline static int digits() { return int(mpfr::mpreal::get_default_prec()); }
+ inline static int digits(const mpfr::mpreal& x) { return x.getPrecision(); }
inline static int digits10(mp_prec_t precision = mpfr::mpreal::get_default_prec())
{
@@ -2970,8 +3047,27 @@ 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
};
}
-#endif /* __MPREAL_H__ */
+#endif /* __MPREAL_H__ */
\ No newline at end of file diff --git a/unsupported/test/mpreal_support.cpp b/unsupported/test/mpreal_support.cpp index bc00382be..1aa9e786a 100644 --- a/unsupported/test/mpreal_support.cpp +++ b/unsupported/test/mpreal_support.cpp @@ -32,12 +32,11 @@ void test_mpreal_support() VERIFY_IS_APPROX(A.array().abs2().sqrt(), A.array().abs()); VERIFY_IS_APPROX(A.array().sin(), sin(A.array())); VERIFY_IS_APPROX(A.array().cos(), cos(A.array())); - // Cholesky X = S.selfadjointView<Lower>().llt().solve(B); VERIFY_IS_APPROX((S.selfadjointView<Lower>()*X).eval(),B); - + // partial LU X = A.lu().solve(B); VERIFY_IS_APPROX((A*X).eval(),B); |