diff options
92 files changed, 1533 insertions, 839 deletions
diff --git a/Android.mk b/Android.mk index 4078ec1b8..272cd956e 100644 --- a/Android.mk +++ b/Android.mk @@ -11,6 +11,3 @@ # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. - -# This empty file is here solely for the purpose of optimizing the Android build -# Please keep it there, and empty, thanks. diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index c52b7d1a6..02ab93880 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -442,6 +442,7 @@ LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const MatrixType& a) m_transpositions.resize(size); m_isInitialized = false; m_temporary.resize(size); + m_sign = internal::ZeroSign; internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, m_sign); @@ -502,7 +503,6 @@ struct solve_retval<LDLT<_MatrixType,_UpLo>, Rhs> using std::abs; using std::max; typedef typename LDLTType::MatrixType MatrixType; - typedef typename LDLTType::Scalar Scalar; typedef typename LDLTType::RealScalar RealScalar; const typename Diagonal<const MatrixType>::RealReturnType vectorD(dec().vectorD()); // In some previous versions, tolerance was set to the max of 1/highest and the maximal diagonal entry * epsilon diff --git a/Eigen/src/Core/ArrayWrapper.h b/Eigen/src/Core/ArrayWrapper.h index a791bc358..b4641e2a0 100644 --- a/Eigen/src/Core/ArrayWrapper.h +++ b/Eigen/src/Core/ArrayWrapper.h @@ -29,6 +29,11 @@ struct traits<ArrayWrapper<ExpressionType> > : public traits<typename remove_all<typename ExpressionType::Nested>::type > { typedef ArrayXpr XprKind; + // Let's remove NestByRefBit + enum { + Flags0 = traits<typename remove_all<typename ExpressionType::Nested>::type >::Flags, + Flags = Flags0 & ~NestByRefBit + }; }; } @@ -149,6 +154,11 @@ struct traits<MatrixWrapper<ExpressionType> > : public traits<typename remove_all<typename ExpressionType::Nested>::type > { typedef MatrixXpr XprKind; + // Let's remove NestByRefBit + enum { + Flags0 = traits<typename remove_all<typename ExpressionType::Nested>::type >::Flags, + Flags = Flags0 & ~NestByRefBit + }; }; } diff --git a/Eigen/src/Core/DenseBase.h b/Eigen/src/Core/DenseBase.h index c5800f6c8..04862f374 100644 --- a/Eigen/src/Core/DenseBase.h +++ b/Eigen/src/Core/DenseBase.h @@ -462,8 +462,10 @@ template<typename Derived> class DenseBase template<int p> RealScalar lpNorm() const; template<int RowFactor, int ColFactor> - const Replicate<Derived,RowFactor,ColFactor> replicate() const; - const Replicate<Derived,Dynamic,Dynamic> replicate(Index rowFacor,Index colFactor) const; + inline const Replicate<Derived,RowFactor,ColFactor> replicate() const; + + typedef Replicate<Derived,Dynamic,Dynamic> ReplicateReturnType; + inline const ReplicateReturnType replicate(Index rowFacor,Index colFactor) const; typedef Reverse<Derived, BothDirections> ReverseReturnType; typedef const Reverse<const Derived, BothDirections> ConstReverseReturnType; diff --git a/Eigen/src/Core/Diagonal.h b/Eigen/src/Core/Diagonal.h index aab8007b3..68cf6d4b0 100644 --- a/Eigen/src/Core/Diagonal.h +++ b/Eigen/src/Core/Diagonal.h @@ -190,18 +190,18 @@ MatrixBase<Derived>::diagonal() const * * \sa MatrixBase::diagonal(), class Diagonal */ template<typename Derived> -inline typename MatrixBase<Derived>::template DiagonalIndexReturnType<DynamicIndex>::Type +inline typename MatrixBase<Derived>::DiagonalDynamicIndexReturnType MatrixBase<Derived>::diagonal(Index index) { - return typename DiagonalIndexReturnType<DynamicIndex>::Type(derived(), index); + return DiagonalDynamicIndexReturnType(derived(), index); } /** This is the const version of diagonal(Index). */ template<typename Derived> -inline typename MatrixBase<Derived>::template ConstDiagonalIndexReturnType<DynamicIndex>::Type +inline typename MatrixBase<Derived>::ConstDiagonalDynamicIndexReturnType MatrixBase<Derived>::diagonal(Index index) const { - return typename ConstDiagonalIndexReturnType<DynamicIndex>::Type(derived(), index); + return ConstDiagonalDynamicIndexReturnType(derived(), index); } /** \returns an expression of the \a DiagIndex-th sub or super diagonal of the matrix \c *this diff --git a/Eigen/src/Core/GeneralProduct.h b/Eigen/src/Core/GeneralProduct.h index 2a59d9464..9e805a80f 100644 --- a/Eigen/src/Core/GeneralProduct.h +++ b/Eigen/src/Core/GeneralProduct.h @@ -232,7 +232,7 @@ EIGEN_DONT_INLINE void outer_product_selector_run(const ProductType& prod, Dest& // FIXME not very good if rhs is real and lhs complex while alpha is real too const Index cols = dest.cols(); for (Index j=0; j<cols; ++j) - func(dest.col(j), prod.rhs().coeff(j) * prod.lhs()); + func(dest.col(j), prod.rhs().coeff(0,j) * prod.lhs()); } // Row major @@ -243,7 +243,7 @@ EIGEN_DONT_INLINE void outer_product_selector_run(const ProductType& prod, Dest& // FIXME not very good if lhs is real and rhs complex while alpha is real too const Index rows = dest.rows(); for (Index i=0; i<rows; ++i) - func(dest.row(i), prod.lhs().coeff(i) * prod.rhs()); + func(dest.row(i), prod.lhs().coeff(i,0) * prod.rhs()); } template<typename Lhs, typename Rhs> diff --git a/Eigen/src/Core/MapBase.h b/Eigen/src/Core/MapBase.h index ab50c9b81..cebed2bb6 100644 --- a/Eigen/src/Core/MapBase.h +++ b/Eigen/src/Core/MapBase.h @@ -168,6 +168,7 @@ template<typename Derived> class MapBase<Derived, ReadOnlyAccessors> template<typename Derived> class MapBase<Derived, WriteAccessors> : public MapBase<Derived, ReadOnlyAccessors> { + typedef MapBase<Derived, ReadOnlyAccessors> ReadOnlyMapBase; public: typedef MapBase<Derived, ReadOnlyAccessors> Base; @@ -230,11 +231,13 @@ template<typename Derived> class MapBase<Derived, WriteAccessors> Derived& operator=(const MapBase& other) { - Base::Base::operator=(other); + ReadOnlyMapBase::Base::operator=(other); return derived(); } - using Base::Base::operator=; + // In theory we could simply refer to Base:Base::operator=, but MSVC does not like Base::Base, + // see bugs 821 and 920. + using ReadOnlyMapBase::Base::operator=; }; #undef EIGEN_STATIC_ASSERT_INDEX_BASED_ACCESS diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 344b38f2f..cc3279746 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -215,7 +215,7 @@ template<typename Derived> class MatrixBase typedef Diagonal<Derived> DiagonalReturnType; DiagonalReturnType diagonal(); - typedef typename internal::add_const<Diagonal<const Derived> >::type ConstDiagonalReturnType; + typedef typename internal::add_const<Diagonal<const Derived> >::type ConstDiagonalReturnType; ConstDiagonalReturnType diagonal() const; template<int Index> struct DiagonalIndexReturnType { typedef Diagonal<Derived,Index> Type; }; @@ -223,16 +223,12 @@ template<typename Derived> class MatrixBase template<int Index> typename DiagonalIndexReturnType<Index>::Type diagonal(); template<int Index> typename ConstDiagonalIndexReturnType<Index>::Type diagonal() const; + + typedef Diagonal<Derived,DynamicIndex> DiagonalDynamicIndexReturnType; + typedef typename internal::add_const<Diagonal<const Derived,DynamicIndex> >::type ConstDiagonalDynamicIndexReturnType; - // Note: The "MatrixBase::" prefixes are added to help MSVC9 to match these declarations with the later implementations. - // On the other hand they confuse MSVC8... - #if (defined _MSC_VER) && (_MSC_VER >= 1500) // 2008 or later - typename MatrixBase::template DiagonalIndexReturnType<DynamicIndex>::Type diagonal(Index index); - typename MatrixBase::template ConstDiagonalIndexReturnType<DynamicIndex>::Type diagonal(Index index) const; - #else - typename DiagonalIndexReturnType<DynamicIndex>::Type diagonal(Index index); - typename ConstDiagonalIndexReturnType<DynamicIndex>::Type diagonal(Index index) const; - #endif + DiagonalDynamicIndexReturnType diagonal(Index index); + ConstDiagonalDynamicIndexReturnType diagonal(Index index) const; #ifdef EIGEN2_SUPPORT template<unsigned int Mode> typename internal::eigen2_part_return_type<Derived, Mode>::type part(); diff --git a/Eigen/src/Core/PermutationMatrix.h b/Eigen/src/Core/PermutationMatrix.h index 1297b8413..f26f3e5cc 100644 --- a/Eigen/src/Core/PermutationMatrix.h +++ b/Eigen/src/Core/PermutationMatrix.h @@ -555,7 +555,10 @@ struct permut_matrix_product_retval const Index n = Side==OnTheLeft ? rows() : cols(); // FIXME we need an is_same for expression that is not sensitive to constness. For instance // is_same_xpr<Block<const Matrix>, Block<Matrix> >::value should be true. - if(is_same<MatrixTypeNestedCleaned,Dest>::value && extract_data(dst) == extract_data(m_matrix)) + if( is_same<MatrixTypeNestedCleaned,Dest>::value + && blas_traits<MatrixTypeNestedCleaned>::HasUsableDirectAccess + && blas_traits<Dest>::HasUsableDirectAccess + && extract_data(dst) == extract_data(m_matrix)) { // apply the permutation inplace Matrix<bool,PermutationType::RowsAtCompileTime,1,0,PermutationType::MaxRowsAtCompileTime> mask(m_permutation.size()); diff --git a/Eigen/src/Core/ProductBase.h b/Eigen/src/Core/ProductBase.h index a494b5f87..cf74470a9 100644 --- a/Eigen/src/Core/ProductBase.h +++ b/Eigen/src/Core/ProductBase.h @@ -85,7 +85,14 @@ class ProductBase : public MatrixBase<Derived> public: +#ifndef EIGEN_NO_MALLOC + typedef typename Base::PlainObject BasePlainObject; + typedef Matrix<Scalar,RowsAtCompileTime==1?1:Dynamic,ColsAtCompileTime==1?1:Dynamic,BasePlainObject::Options> DynPlainObject; + typedef typename internal::conditional<(BasePlainObject::SizeAtCompileTime==Dynamic) || (BasePlainObject::SizeAtCompileTime*int(sizeof(Scalar)) < int(EIGEN_STACK_ALLOCATION_LIMIT)), + BasePlainObject, DynPlainObject>::type PlainObject; +#else typedef typename Base::PlainObject PlainObject; +#endif ProductBase(const Lhs& a_lhs, const Rhs& a_rhs) : m_lhs(a_lhs), m_rhs(a_rhs) @@ -180,7 +187,12 @@ namespace internal { template<typename Lhs, typename Rhs, int Mode, int N, typename PlainObject> struct nested<GeneralProduct<Lhs,Rhs,Mode>, N, PlainObject> { - typedef PlainObject const& type; + typedef typename GeneralProduct<Lhs,Rhs,Mode>::PlainObject const& type; +}; +template<typename Lhs, typename Rhs, int Mode, int N, typename PlainObject> +struct nested<const GeneralProduct<Lhs,Rhs,Mode>, N, PlainObject> +{ + typedef typename GeneralProduct<Lhs,Rhs,Mode>::PlainObject const& type; }; } diff --git a/Eigen/src/Core/Ref.h b/Eigen/src/Core/Ref.h index cd6d949c4..df87b1af4 100644 --- a/Eigen/src/Core/Ref.h +++ b/Eigen/src/Core/Ref.h @@ -188,6 +188,8 @@ template<typename PlainObjectType, int Options, typename StrideType> class Ref : public RefBase<Ref<PlainObjectType, Options, StrideType> > { typedef internal::traits<Ref> Traits; + template<typename Derived> + inline Ref(const PlainObjectBase<Derived>& expr); public: typedef RefBase<Ref> Base; @@ -196,20 +198,21 @@ template<typename PlainObjectType, int Options, typename StrideType> class Ref #ifndef EIGEN_PARSED_BY_DOXYGEN template<typename Derived> - inline Ref(PlainObjectBase<Derived>& expr, - typename internal::enable_if<bool(Traits::template match<Derived>::MatchAtCompileTime),Derived>::type* = 0) + inline Ref(PlainObjectBase<Derived>& expr) { - Base::construct(expr); + EIGEN_STATIC_ASSERT(static_cast<bool>(Traits::template match<Derived>::MatchAtCompileTime), STORAGE_LAYOUT_DOES_NOT_MATCH); + Base::construct(expr.derived()); } template<typename Derived> - inline Ref(const DenseBase<Derived>& expr, - typename internal::enable_if<bool(internal::is_lvalue<Derived>::value&&bool(Traits::template match<Derived>::MatchAtCompileTime)),Derived>::type* = 0, - int = Derived::ThisConstantIsPrivateInPlainObjectBase) + inline Ref(const DenseBase<Derived>& expr) #else template<typename Derived> inline Ref(DenseBase<Derived>& expr) #endif { + EIGEN_STATIC_ASSERT(static_cast<bool>(internal::is_lvalue<Derived>::value), THIS_EXPRESSION_IS_NOT_A_LVALUE__IT_IS_READ_ONLY); + EIGEN_STATIC_ASSERT(static_cast<bool>(Traits::template match<Derived>::MatchAtCompileTime), STORAGE_LAYOUT_DOES_NOT_MATCH); + enum { THIS_EXPRESSION_IS_NOT_A_LVALUE__IT_IS_READ_ONLY = Derived::ThisConstantIsPrivateInPlainObjectBase}; Base::construct(expr.const_cast_derived()); } diff --git a/Eigen/src/Core/Replicate.h b/Eigen/src/Core/Replicate.h index dde86a834..ac4537c14 100644 --- a/Eigen/src/Core/Replicate.h +++ b/Eigen/src/Core/Replicate.h @@ -135,7 +135,7 @@ template<typename MatrixType,int RowFactor,int ColFactor> class Replicate */ template<typename Derived> template<int RowFactor, int ColFactor> -inline const Replicate<Derived,RowFactor,ColFactor> +const Replicate<Derived,RowFactor,ColFactor> DenseBase<Derived>::replicate() const { return Replicate<Derived,RowFactor,ColFactor>(derived()); @@ -150,7 +150,7 @@ DenseBase<Derived>::replicate() const * \sa VectorwiseOp::replicate(), DenseBase::replicate<int,int>(), class Replicate */ template<typename Derived> -inline const Replicate<Derived,Dynamic,Dynamic> +const typename DenseBase<Derived>::ReplicateReturnType DenseBase<Derived>::replicate(Index rowFactor,Index colFactor) const { return Replicate<Derived,Dynamic,Dynamic>(derived(),rowFactor,colFactor); diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index 845ae1aec..4d65392c6 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -380,19 +380,19 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView EIGEN_STRONG_INLINE TriangularView& operator=(const ProductBase<ProductDerived, Lhs,Rhs>& other) { setZero(); - return assignProduct(other,1); + return assignProduct(other.derived(),1); } template<typename ProductDerived, typename Lhs, typename Rhs> EIGEN_STRONG_INLINE TriangularView& operator+=(const ProductBase<ProductDerived, Lhs,Rhs>& other) { - return assignProduct(other,1); + return assignProduct(other.derived(),1); } template<typename ProductDerived, typename Lhs, typename Rhs> EIGEN_STRONG_INLINE TriangularView& operator-=(const ProductBase<ProductDerived, Lhs,Rhs>& other) { - return assignProduct(other,-1); + return assignProduct(other.derived(),-1); } @@ -400,25 +400,34 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView EIGEN_STRONG_INLINE TriangularView& operator=(const ScaledProduct<ProductDerived>& other) { setZero(); - return assignProduct(other,other.alpha()); + return assignProduct(other.derived(),other.alpha()); } template<typename ProductDerived> EIGEN_STRONG_INLINE TriangularView& operator+=(const ScaledProduct<ProductDerived>& other) { - return assignProduct(other,other.alpha()); + return assignProduct(other.derived(),other.alpha()); } template<typename ProductDerived> EIGEN_STRONG_INLINE TriangularView& operator-=(const ScaledProduct<ProductDerived>& other) { - return assignProduct(other,-other.alpha()); + return assignProduct(other.derived(),-other.alpha()); } protected: template<typename ProductDerived, typename Lhs, typename Rhs> EIGEN_STRONG_INLINE TriangularView& assignProduct(const ProductBase<ProductDerived, Lhs,Rhs>& prod, const Scalar& alpha); + + template<int Mode, bool LhsIsTriangular, + typename Lhs, bool LhsIsVector, + typename Rhs, bool RhsIsVector> + EIGEN_STRONG_INLINE TriangularView& assignProduct(const TriangularProduct<Mode, LhsIsTriangular, Lhs, LhsIsVector, Rhs, RhsIsVector>& prod, const Scalar& alpha) + { + lazyAssign(alpha*prod.eval()); + return *this; + } MatrixTypeNested m_matrix; }; diff --git a/Eigen/src/Core/arch/NEON/Complex.h b/Eigen/src/Core/arch/NEON/Complex.h index f183d31de..8d9255eef 100644 --- a/Eigen/src/Core/arch/NEON/Complex.h +++ b/Eigen/src/Core/arch/NEON/Complex.h @@ -110,7 +110,7 @@ template<> EIGEN_STRONG_INLINE Packet2cf ploaddup<Packet2cf>(const std::complex< template<> EIGEN_STRONG_INLINE void pstore <std::complex<float> >(std::complex<float> * to, const Packet2cf& from) { EIGEN_DEBUG_ALIGNED_STORE pstore((float*)to, from.v); } template<> EIGEN_STRONG_INLINE void pstoreu<std::complex<float> >(std::complex<float> * to, const Packet2cf& from) { EIGEN_DEBUG_UNALIGNED_STORE pstoreu((float*)to, from.v); } -template<> EIGEN_STRONG_INLINE void prefetch<std::complex<float> >(const std::complex<float> * addr) { __pld((float *)addr); } +template<> EIGEN_STRONG_INLINE void prefetch<std::complex<float> >(const std::complex<float> * addr) { EIGEN_ARM_PREFETCH((float *)addr); } template<> EIGEN_STRONG_INLINE std::complex<float> pfirst<Packet2cf>(const Packet2cf& a) { diff --git a/Eigen/src/Core/arch/NEON/PacketMath.h b/Eigen/src/Core/arch/NEON/PacketMath.h index 163bac215..94dfab330 100644 --- a/Eigen/src/Core/arch/NEON/PacketMath.h +++ b/Eigen/src/Core/arch/NEON/PacketMath.h @@ -48,9 +48,18 @@ typedef uint32x4_t Packet4ui; #define EIGEN_INIT_NEON_PACKET2(X, Y) {X, Y} #define EIGEN_INIT_NEON_PACKET4(X, Y, Z, W) {X, Y, Z, W} #endif - -#ifndef __pld -#define __pld(x) asm volatile ( " pld [%[addr]]\n" :: [addr] "r" (x) : "cc" ); + +// arm64 does have the pld instruction. If available, let's trust the __builtin_prefetch built-in function +// which available on LLVM and GCC (at least) +#if EIGEN_HAS_BUILTIN(__builtin_prefetch) || defined(__GNUC__) + #define EIGEN_ARM_PREFETCH(ADDR) __builtin_prefetch(ADDR); +#elif defined __pld + #define EIGEN_ARM_PREFETCH(ADDR) __pld(ADDR) +#elif !defined(__aarch64__) + #define EIGEN_ARM_PREFETCH(ADDR) __asm__ __volatile__ ( " pld [%[addr]]\n" :: [addr] "r" (ADDR) : "cc" ); +#else + // by default no explicit prefetching + #define EIGEN_ARM_PREFETCH(ADDR) #endif template<> struct packet_traits<float> : default_packet_traits @@ -209,8 +218,8 @@ template<> EIGEN_STRONG_INLINE void pstore<int>(int* to, const Packet4i& f template<> EIGEN_STRONG_INLINE void pstoreu<float>(float* to, const Packet4f& from) { EIGEN_DEBUG_UNALIGNED_STORE vst1q_f32(to, from); } template<> EIGEN_STRONG_INLINE void pstoreu<int>(int* to, const Packet4i& from) { EIGEN_DEBUG_UNALIGNED_STORE vst1q_s32(to, from); } -template<> EIGEN_STRONG_INLINE void prefetch<float>(const float* addr) { __pld(addr); } -template<> EIGEN_STRONG_INLINE void prefetch<int>(const int* addr) { __pld(addr); } +template<> EIGEN_STRONG_INLINE void prefetch<float>(const float* addr) { EIGEN_ARM_PREFETCH(addr); } +template<> EIGEN_STRONG_INLINE void prefetch<int>(const int* addr) { EIGEN_ARM_PREFETCH(addr); } // FIXME only store the 2 first elements ? template<> EIGEN_STRONG_INLINE float pfirst<Packet4f>(const Packet4f& a) { float EIGEN_ALIGN16 x[4]; vst1q_f32(x, a); return x[0]; } diff --git a/Eigen/src/Core/arch/SSE/MathFunctions.h b/Eigen/src/Core/arch/SSE/MathFunctions.h index 99cbd0d95..d16f30bb0 100644 --- a/Eigen/src/Core/arch/SSE/MathFunctions.h +++ b/Eigen/src/Core/arch/SSE/MathFunctions.h @@ -52,7 +52,7 @@ Packet4f plog<Packet4f>(const Packet4f& _x) Packet4i emm0; - Packet4f invalid_mask = _mm_cmplt_ps(x, _mm_setzero_ps()); + Packet4f invalid_mask = _mm_cmpnge_ps(x, _mm_setzero_ps()); // not greater equal is true if x is NaN Packet4f iszero_mask = _mm_cmpeq_ps(x, _mm_setzero_ps()); x = pmax(x, p4f_min_norm_pos); /* cut off denormalized stuff */ @@ -166,7 +166,7 @@ Packet4f pexp<Packet4f>(const Packet4f& _x) emm0 = _mm_cvttps_epi32(fx); emm0 = _mm_add_epi32(emm0, p4i_0x7f); emm0 = _mm_slli_epi32(emm0, 23); - return pmul(y, _mm_castsi128_ps(emm0)); + return pmax(pmul(y, Packet4f(_mm_castsi128_ps(emm0))), _x); } template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet2d pexp<Packet2d>(const Packet2d& _x) @@ -239,7 +239,7 @@ Packet2d pexp<Packet2d>(const Packet2d& _x) emm0 = _mm_add_epi32(emm0, p4i_1023_0); emm0 = _mm_slli_epi32(emm0, 20); emm0 = _mm_shuffle_epi32(emm0, _MM_SHUFFLE(1,2,0,3)); - return pmul(x, _mm_castsi128_pd(emm0)); + return pmax(pmul(x, Packet2d(_mm_castsi128_pd(emm0))), _x); } /* evaluation of 4 sines at onces, using SSE2 intrinsics. diff --git a/Eigen/src/Core/products/CoeffBasedProduct.h b/Eigen/src/Core/products/CoeffBasedProduct.h index c06a0df1c..421f925e1 100644 --- a/Eigen/src/Core/products/CoeffBasedProduct.h +++ b/Eigen/src/Core/products/CoeffBasedProduct.h @@ -90,6 +90,7 @@ struct traits<CoeffBasedProduct<LhsNested,RhsNested,NestingFlags> > | (SameType && (CanVectorizeLhs || CanVectorizeRhs) ? PacketAccessBit : 0), CoeffReadCost = InnerSize == Dynamic ? Dynamic + : InnerSize == 0 ? 0 : InnerSize * (NumTraits<Scalar>::MulCost + LhsCoeffReadCost + RhsCoeffReadCost) + (InnerSize - 1) * NumTraits<Scalar>::AddCost, @@ -133,7 +134,7 @@ class CoeffBasedProduct }; typedef internal::product_coeff_impl<CanVectorizeInner ? InnerVectorizedTraversal : DefaultTraversal, - Unroll ? InnerSize-1 : Dynamic, + Unroll ? (InnerSize==0 ? 0 : InnerSize-1) : Dynamic, _LhsNested, _RhsNested, Scalar> ScalarCoeffImpl; typedef CoeffBasedProduct<LhsNested,RhsNested,NestByRefBit> LazyCoeffBasedProductType; @@ -184,7 +185,7 @@ class CoeffBasedProduct { PacketScalar res; internal::product_packet_impl<Flags&RowMajorBit ? RowMajor : ColMajor, - Unroll ? InnerSize-1 : Dynamic, + Unroll ? (InnerSize==0 ? 0 : InnerSize-1) : Dynamic, _LhsNested, _RhsNested, PacketScalar, LoadMode> ::run(row, col, m_lhs, m_rhs, res); return res; @@ -262,10 +263,7 @@ struct product_coeff_impl<DefaultTraversal, Dynamic, Lhs, Rhs, RetScalar> typedef typename Lhs::Index Index; static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, RetScalar& res) { - eigen_assert(lhs.cols()>0 && "you are using a non initialized matrix"); - res = lhs.coeff(row, 0) * rhs.coeff(0, col); - for(Index i = 1; i < lhs.cols(); ++i) - res += lhs.coeff(row, i) * rhs.coeff(i, col); + res = (lhs.row(row).transpose().cwiseProduct( rhs.col(col) )).sum(); } }; diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h index 3a010ec6a..6d1e6c133 100644 --- a/Eigen/src/Core/util/Macros.h +++ b/Eigen/src/Core/util/Macros.h @@ -13,7 +13,7 @@ #define EIGEN_WORLD_VERSION 3 #define EIGEN_MAJOR_VERSION 2 -#define EIGEN_MINOR_VERSION 2 +#define EIGEN_MINOR_VERSION 4 #define EIGEN_VERSION_AT_LEAST(x,y,z) (EIGEN_WORLD_VERSION>x || (EIGEN_WORLD_VERSION>=x && \ (EIGEN_MAJOR_VERSION>y || (EIGEN_MAJOR_VERSION>=y && \ @@ -96,6 +96,13 @@ #define EIGEN_DEFAULT_DENSE_INDEX_TYPE std::ptrdiff_t #endif +// Cross compiler wrapper around LLVM's __has_builtin +#ifdef __has_builtin +# define EIGEN_HAS_BUILTIN(x) __has_builtin(x) +#else +# define EIGEN_HAS_BUILTIN(x) 0 +#endif + /** Allows to disable some optimizations which might affect the accuracy of the result. * Such optimization are enabled by default, and set EIGEN_FAST_MATH to 0 to disable them. * They currently include: @@ -247,7 +254,7 @@ namespace Eigen { #if !defined(EIGEN_ASM_COMMENT) #if (defined __GNUC__) && ( defined(__i386__) || defined(__x86_64__) ) - #define EIGEN_ASM_COMMENT(X) asm("#" X) + #define EIGEN_ASM_COMMENT(X) __asm__("#" X) #else #define EIGEN_ASM_COMMENT(X) #endif @@ -306,7 +313,7 @@ namespace Eigen { // just an empty macro ! #define EIGEN_EMPTY -#if defined(_MSC_VER) && (!defined(__INTEL_COMPILER)) +#if defined(_MSC_VER) && (_MSC_VER < 1900) && (!defined(__INTEL_COMPILER)) #define EIGEN_INHERIT_ASSIGNMENT_EQUAL_OPERATOR(Derived) \ using Base::operator =; #elif defined(__clang__) // workaround clang bug (see http://forum.kde.org/viewtopic.php?f=74&t=102653) diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h index c4011245a..330bcf518 100644 --- a/Eigen/src/Core/util/Memory.h +++ b/Eigen/src/Core/util/Memory.h @@ -63,7 +63,7 @@ // Currently, let's include it only on unix systems: #if defined(__unix__) || defined(__unix) #include <unistd.h> - #if ((defined __QNXNTO__) || (defined _GNU_SOURCE) || ((defined _XOPEN_SOURCE) && (_XOPEN_SOURCE >= 600))) && (defined _POSIX_ADVISORY_INFO) && (_POSIX_ADVISORY_INFO > 0) + #if ((defined __QNXNTO__) || (defined _GNU_SOURCE) || (defined __PGI) || ((defined _XOPEN_SOURCE) && (_XOPEN_SOURCE >= 600))) && (defined _POSIX_ADVISORY_INFO) && (_POSIX_ADVISORY_INFO > 0) #define EIGEN_HAS_POSIX_MEMALIGN 1 #endif #endif @@ -417,6 +417,8 @@ template<typename T, bool Align> inline T* conditional_aligned_realloc_new(T* pt template<typename T, bool Align> inline T* conditional_aligned_new_auto(size_t size) { + if(size==0) + return 0; // short-cut. Also fixes Bug 884 check_size_for_overflow<T>(size); T *result = reinterpret_cast<T*>(conditional_aligned_malloc<Align>(sizeof(T)*size)); if(NumTraits<T>::RequireInitialization) @@ -464,9 +466,8 @@ template<typename T, bool Align> inline void conditional_aligned_delete_auto(T * template<typename Scalar, typename Index> static inline Index first_aligned(const Scalar* array, Index size) { - enum { PacketSize = packet_traits<Scalar>::size, - PacketAlignedMask = PacketSize-1 - }; + static const Index PacketSize = packet_traits<Scalar>::size; + static const Index PacketAlignedMask = PacketSize-1; if(PacketSize==1) { @@ -612,7 +613,6 @@ template<typename T> class aligned_stack_memory_handler void* operator new(size_t size, const std::nothrow_t&) throw() { \ try { return Eigen::internal::conditional_aligned_malloc<NeedsToAlign>(size); } \ catch (...) { return 0; } \ - return 0; \ } #else #define EIGEN_MAKE_ALIGNED_OPERATOR_NEW_NOTHROW(NeedsToAlign) \ diff --git a/Eigen/src/Core/util/StaticAssert.h b/Eigen/src/Core/util/StaticAssert.h index 8872c5b64..bac5d9fe9 100644 --- a/Eigen/src/Core/util/StaticAssert.h +++ b/Eigen/src/Core/util/StaticAssert.h @@ -90,7 +90,9 @@ YOU_PASSED_A_COLUMN_VECTOR_BUT_A_ROW_VECTOR_WAS_EXPECTED, THE_INDEX_TYPE_MUST_BE_A_SIGNED_TYPE, THE_STORAGE_ORDER_OF_BOTH_SIDES_MUST_MATCH, - OBJECT_ALLOCATED_ON_STACK_IS_TOO_BIG + OBJECT_ALLOCATED_ON_STACK_IS_TOO_BIG, + IMPLICIT_CONVERSION_TO_SCALAR_IS_FOR_INNER_PRODUCT_ONLY, + STORAGE_LAYOUT_DOES_NOT_MATCH }; }; diff --git a/Eigen/src/Core/util/XprHelper.h b/Eigen/src/Core/util/XprHelper.h index 3c4773054..781965d2c 100644 --- a/Eigen/src/Core/util/XprHelper.h +++ b/Eigen/src/Core/util/XprHelper.h @@ -341,7 +341,7 @@ template<typename T, int n=1, typename PlainObject = typename eval<T>::type> str }; template<typename T> -T* const_cast_ptr(const T* ptr) +inline T* const_cast_ptr(const T* ptr) { return const_cast<T*>(ptr); } diff --git a/Eigen/src/Eigen2Support/LeastSquares.h b/Eigen/src/Eigen2Support/LeastSquares.h index 0e6fdb488..7992d4944 100644 --- a/Eigen/src/Eigen2Support/LeastSquares.h +++ b/Eigen/src/Eigen2Support/LeastSquares.h @@ -147,7 +147,6 @@ void fitHyperplane(int numPoints, // compute the covariance matrix CovMatrixType covMat = CovMatrixType::Zero(size, size); - VectorType remean = VectorType::Zero(size); for(int i = 0; i < numPoints; ++i) { VectorType diff = (*(points[i]) - mean).conjugate(); diff --git a/Eigen/src/Eigenvalues/RealQZ.h b/Eigen/src/Eigenvalues/RealQZ.h index 5706eeebe..4f36091db 100644 --- a/Eigen/src/Eigenvalues/RealQZ.h +++ b/Eigen/src/Eigenvalues/RealQZ.h @@ -313,7 +313,7 @@ namespace Eigen { using std::abs; using std::sqrt; const Index dim=m_S.cols(); - if (abs(m_S.coeff(i+1,i)==Scalar(0))) + if (abs(m_S.coeff(i+1,i))==Scalar(0)) return; Index z = findSmallDiagEntry(i,i+1); if (z==i-1) diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 3993046a8..be89de4a9 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -563,7 +563,6 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3 if(computeEigenvectors) { Scalar safeNorm2 = Eigen::NumTraits<Scalar>::epsilon(); - safeNorm2 *= safeNorm2; if((eivals(2)-eivals(0))<=Eigen::NumTraits<Scalar>::epsilon()) { eivecs.setIdentity(); @@ -577,7 +576,7 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3 Scalar d0 = eivals(2) - eivals(1); Scalar d1 = eivals(1) - eivals(0); int k = d0 > d1 ? 2 : 0; - d0 = d0 > d1 ? d1 : d0; + d0 = d0 > d1 ? d0 : d1; tmp.diagonal().array () -= eivals(k); VectorType cross; @@ -585,19 +584,25 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3 n = (cross = tmp.row(0).cross(tmp.row(1))).squaredNorm(); if(n>safeNorm2) + { eivecs.col(k) = cross / sqrt(n); + } else { n = (cross = tmp.row(0).cross(tmp.row(2))).squaredNorm(); if(n>safeNorm2) + { eivecs.col(k) = cross / sqrt(n); + } else { n = (cross = tmp.row(1).cross(tmp.row(2))).squaredNorm(); if(n>safeNorm2) + { eivecs.col(k) = cross / sqrt(n); + } else { // the input matrix and/or the eigenvaues probably contains some inf/NaN, @@ -617,12 +622,16 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3 tmp.diagonal().array() -= eivals(1); if(d0<=Eigen::NumTraits<Scalar>::epsilon()) + { eivecs.col(1) = eivecs.col(k).unitOrthogonal(); + } else { - n = (cross = eivecs.col(k).cross(tmp.row(0).normalized())).squaredNorm(); + n = (cross = eivecs.col(k).cross(tmp.row(0))).squaredNorm(); if(n>safeNorm2) + { eivecs.col(1) = cross / sqrt(n); + } else { n = (cross = eivecs.col(k).cross(tmp.row(1))).squaredNorm(); @@ -636,13 +645,14 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3 else { // we should never reach this point, - // if so the last two eigenvalues are likely to ve very closed to each other + // if so the last two eigenvalues are likely to be very close to each other eivecs.col(1) = eivecs.col(k).unitOrthogonal(); } } } // make sure that eivecs[1] is orthogonal to eivecs[2] + // FIXME: this step should not be needed Scalar d = eivecs.col(1).dot(eivecs.col(k)); eivecs.col(1) = (eivecs.col(1) - d * eivecs.col(k)).normalized(); } diff --git a/Eigen/src/Geometry/Hyperplane.h b/Eigen/src/Geometry/Hyperplane.h index aeff43fef..00b7c4300 100644 --- a/Eigen/src/Geometry/Hyperplane.h +++ b/Eigen/src/Geometry/Hyperplane.h @@ -100,7 +100,17 @@ public: { EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(VectorType, 3) Hyperplane result(p0.size()); - result.normal() = (p2 - p0).cross(p1 - p0).normalized(); + VectorType v0(p2 - p0), v1(p1 - p0); + result.normal() = v0.cross(v1); + RealScalar norm = result.normal().norm(); + if(norm <= v0.norm() * v1.norm() * NumTraits<RealScalar>::epsilon()) + { + Matrix<Scalar,2,3> m; m << v0.transpose(), v1.transpose(); + JacobiSVD<Matrix<Scalar,2,3> > svd(m, ComputeFullV); + result.normal() = svd.matrixV().col(2); + } + else + result.normal() /= norm; result.offset() = -p0.dot(result.normal()); return result; } diff --git a/Eigen/src/Geometry/Rotation2D.h b/Eigen/src/Geometry/Rotation2D.h index 1cac343a5..a2d59fce1 100644 --- a/Eigen/src/Geometry/Rotation2D.h +++ b/Eigen/src/Geometry/Rotation2D.h @@ -60,6 +60,9 @@ public: /** Construct a 2D counter clock wise rotation from the angle \a a in radian. */ inline Rotation2D(const Scalar& a) : m_angle(a) {} + + /** Default constructor wihtout initialization. The represented rotation is undefined. */ + Rotation2D() {} /** \returns the rotation angle */ inline Scalar angle() const { return m_angle; } @@ -81,10 +84,10 @@ public: /** Applies the rotation to a 2D vector */ Vector2 operator* (const Vector2& vec) const { return toRotationMatrix() * vec; } - + template<typename Derived> Rotation2D& fromRotationMatrix(const MatrixBase<Derived>& m); - Matrix2 toRotationMatrix(void) const; + Matrix2 toRotationMatrix() const; /** \returns the spherical interpolation between \c *this and \a other using * parameter \a t. It is in fact equivalent to a linear interpolation. diff --git a/Eigen/src/Geometry/Transform.h b/Eigen/src/Geometry/Transform.h index 56f361566..e786e5356 100644 --- a/Eigen/src/Geometry/Transform.h +++ b/Eigen/src/Geometry/Transform.h @@ -62,6 +62,8 @@ struct transform_construct_from_matrix; template<typename TransformType> struct transform_take_affine_part; +template<int Mode> struct transform_make_affine; + } // end namespace internal /** \geometry_module \ingroup Geometry_Module @@ -230,8 +232,7 @@ public: inline Transform() { check_template_params(); - if (int(Mode)==Affine) - makeAffine(); + internal::transform_make_affine<(int(Mode)==Affine) ? Affine : AffineCompact>::run(m_matrix); } inline Transform(const Transform& other) @@ -591,11 +592,7 @@ public: */ void makeAffine() { - if(int(Mode)!=int(AffineCompact)) - { - matrix().template block<1,Dim>(Dim,0).setZero(); - matrix().coeffRef(Dim,Dim) = Scalar(1); - } + internal::transform_make_affine<int(Mode)>::run(m_matrix); } /** \internal @@ -1079,6 +1076,24 @@ Transform<Scalar,Dim,Mode,Options>::fromPositionOrientationScale(const MatrixBas namespace internal { +template<int Mode> +struct transform_make_affine +{ + template<typename MatrixType> + static void run(MatrixType &mat) + { + static const int Dim = MatrixType::ColsAtCompileTime-1; + mat.template block<1,Dim>(Dim,0).setZero(); + mat.coeffRef(Dim,Dim) = typename MatrixType::Scalar(1); + } +}; + +template<> +struct transform_make_affine<AffineCompact> +{ + template<typename MatrixType> static void run(MatrixType &) { } +}; + // selector needed to avoid taking the inverse of a 3x4 matrix template<typename TransformType, int Mode=TransformType::Mode> struct projective_transform_inverse diff --git a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h index 2b9fb7f88..dd135c21f 100644 --- a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h +++ b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h @@ -39,7 +39,6 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x, int maxIters = iters; int n = mat.cols(); - x = precond.solve(x); VectorType r = rhs - mat * x; VectorType r0 = r; @@ -143,7 +142,7 @@ struct traits<BiCGSTAB<_MatrixType,_Preconditioner> > * SparseMatrix<double> A(n,n); * // fill A and b * BiCGSTAB<SparseMatrix<double> > solver; - * solver(A); + * solver.compute(A); * x = solver.solve(b); * std::cout << "#iterations: " << solver.iterations() << std::endl; * std::cout << "estimated error: " << solver.error() << std::endl; diff --git a/Eigen/src/PardisoSupport/PardisoSupport.h b/Eigen/src/PardisoSupport/PardisoSupport.h index 1c48f0df7..18cd7d88a 100644 --- a/Eigen/src/PardisoSupport/PardisoSupport.h +++ b/Eigen/src/PardisoSupport/PardisoSupport.h @@ -219,7 +219,7 @@ class PardisoImpl void pardisoInit(int type) { m_type = type; - bool symmetric = abs(m_type) < 10; + bool symmetric = std::abs(m_type) < 10; m_iparm[0] = 1; // No solver default m_iparm[1] = 3; // use Metis for the ordering m_iparm[2] = 1; // Numbers of processors, value of OMP_NUM_THREADS diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index dff9e44eb..c57f2974c 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -762,6 +762,7 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows> m_qr_precond_morecols; internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> m_qr_precond_morerows; + MatrixType m_scaledMatrix; }; template<typename MatrixType, int QRPreconditioner> @@ -808,8 +809,9 @@ void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols, u : 0); m_workMatrix.resize(m_diagSize, m_diagSize); - if(m_cols>m_rows) m_qr_precond_morecols.allocate(*this); - if(m_rows>m_cols) m_qr_precond_morerows.allocate(*this); + if(m_cols>m_rows) m_qr_precond_morecols.allocate(*this); + if(m_rows>m_cols) m_qr_precond_morerows.allocate(*this); + if(m_cols!=m_cols) m_scaledMatrix.resize(rows,cols); } template<typename MatrixType, int QRPreconditioner> @@ -826,21 +828,26 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig // limit for very small denormal numbers to be considered zero in order to avoid infinite loops (see bug 286) const RealScalar considerAsZero = RealScalar(2) * std::numeric_limits<RealScalar>::denorm_min(); + // Scaling factor to reduce over/under-flows + RealScalar scale = matrix.cwiseAbs().maxCoeff(); + if(scale==RealScalar(0)) scale = RealScalar(1); + /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */ - if(!m_qr_precond_morecols.run(*this, matrix) && !m_qr_precond_morerows.run(*this, matrix)) + if(m_rows!=m_cols) { - m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize); + m_scaledMatrix = matrix / scale; + m_qr_precond_morecols.run(*this, m_scaledMatrix); + m_qr_precond_morerows.run(*this, m_scaledMatrix); + } + else + { + m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize) / scale; if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows); if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize); if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols); if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize); } - - // Scaling factor to reduce over/under-flows - RealScalar scale = m_workMatrix.cwiseAbs().maxCoeff(); - if(scale==RealScalar(0)) scale = RealScalar(1); - m_workMatrix /= scale; /*** step 2. The main Jacobi SVD iteration. ***/ @@ -861,7 +868,8 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig using std::max; RealScalar threshold = (max)(considerAsZero, precision * (max)(abs(m_workMatrix.coeff(p,p)), abs(m_workMatrix.coeff(q,q)))); - if((max)(abs(m_workMatrix.coeff(p,q)),abs(m_workMatrix.coeff(q,p))) > threshold) + // We compare both values to threshold instead of calling max to be robust to NaN (See bug 791) + if(abs(m_workMatrix.coeff(p,q))>threshold || abs(m_workMatrix.coeff(q,p)) > threshold) { finished = false; diff --git a/Eigen/src/SparseCore/AmbiVector.h b/Eigen/src/SparseCore/AmbiVector.h index 17fff96a7..220c6451c 100644 --- a/Eigen/src/SparseCore/AmbiVector.h +++ b/Eigen/src/SparseCore/AmbiVector.h @@ -69,7 +69,7 @@ class AmbiVector delete[] m_buffer; if (size<1000) { - Index allocSize = (size * sizeof(ListEl))/sizeof(Scalar); + Index allocSize = (size * sizeof(ListEl) + sizeof(Scalar) - 1)/sizeof(Scalar); m_allocatedElements = (allocSize*sizeof(Scalar))/sizeof(ListEl); m_buffer = new Scalar[allocSize]; } @@ -88,7 +88,7 @@ class AmbiVector Index copyElements = m_allocatedElements; m_allocatedElements = (std::min)(Index(m_allocatedElements*1.5),m_size); Index allocSize = m_allocatedElements * sizeof(ListEl); - allocSize = allocSize/sizeof(Scalar) + (allocSize%sizeof(Scalar)>0?1:0); + allocSize = (allocSize + sizeof(Scalar) - 1)/sizeof(Scalar); Scalar* newBuffer = new Scalar[allocSize]; memcpy(newBuffer, m_buffer, copyElements * sizeof(ListEl)); delete[] m_buffer; diff --git a/Eigen/src/SparseCore/SparseBlock.h b/Eigen/src/SparseCore/SparseBlock.h index 16a20a574..0ede034ba 100644 --- a/Eigen/src/SparseCore/SparseBlock.h +++ b/Eigen/src/SparseCore/SparseBlock.h @@ -68,6 +68,8 @@ public: const internal::variable_if_dynamic<Index, OuterSize> m_outerSize; EIGEN_INHERIT_ASSIGNMENT_OPERATORS(BlockImpl) + private: + Index nonZeros() const; }; @@ -82,6 +84,7 @@ class BlockImpl<SparseMatrix<_Scalar, _Options, _Index>,BlockRows,BlockCols,true typedef SparseMatrix<_Scalar, _Options, _Index> SparseMatrixType; typedef typename internal::remove_all<typename SparseMatrixType::Nested>::type _MatrixTypeNested; typedef Block<SparseMatrixType, BlockRows, BlockCols, true> BlockType; + typedef Block<const SparseMatrixType, BlockRows, BlockCols, true> ConstBlockType; public: enum { IsRowMajor = internal::traits<BlockType>::IsRowMajor }; EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType) @@ -245,6 +248,93 @@ public: }; + +template<typename _Scalar, int _Options, typename _Index, int BlockRows, int BlockCols> +class BlockImpl<const SparseMatrix<_Scalar, _Options, _Index>,BlockRows,BlockCols,true,Sparse> + : public SparseMatrixBase<Block<const SparseMatrix<_Scalar, _Options, _Index>,BlockRows,BlockCols,true> > +{ + typedef SparseMatrix<_Scalar, _Options, _Index> SparseMatrixType; + typedef typename internal::remove_all<typename SparseMatrixType::Nested>::type _MatrixTypeNested; + typedef Block<const SparseMatrixType, BlockRows, BlockCols, true> BlockType; +public: + enum { IsRowMajor = internal::traits<BlockType>::IsRowMajor }; + EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType) +protected: + enum { OuterSize = IsRowMajor ? BlockRows : BlockCols }; +public: + + class InnerIterator: public SparseMatrixType::InnerIterator + { + public: + inline InnerIterator(const BlockType& xpr, Index outer) + : SparseMatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) + {} + inline Index row() const { return IsRowMajor ? m_outer : this->index(); } + inline Index col() const { return IsRowMajor ? this->index() : m_outer; } + protected: + Index m_outer; + }; + class ReverseInnerIterator: public SparseMatrixType::ReverseInnerIterator + { + public: + inline ReverseInnerIterator(const BlockType& xpr, Index outer) + : SparseMatrixType::ReverseInnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) + {} + inline Index row() const { return IsRowMajor ? m_outer : this->index(); } + inline Index col() const { return IsRowMajor ? this->index() : m_outer; } + protected: + Index m_outer; + }; + + inline BlockImpl(const SparseMatrixType& xpr, int i) + : m_matrix(xpr), m_outerStart(i), m_outerSize(OuterSize) + {} + + inline BlockImpl(const SparseMatrixType& xpr, int startRow, int startCol, int blockRows, int blockCols) + : m_matrix(xpr), m_outerStart(IsRowMajor ? startRow : startCol), m_outerSize(IsRowMajor ? blockRows : blockCols) + {} + + inline const Scalar* valuePtr() const + { return m_matrix.valuePtr() + m_matrix.outerIndexPtr()[m_outerStart]; } + + inline const Index* innerIndexPtr() const + { return m_matrix.innerIndexPtr() + m_matrix.outerIndexPtr()[m_outerStart]; } + + inline const Index* outerIndexPtr() const + { return m_matrix.outerIndexPtr() + m_outerStart; } + + Index nonZeros() const + { + if(m_matrix.isCompressed()) + return std::size_t(m_matrix.outerIndexPtr()[m_outerStart+m_outerSize.value()]) + - std::size_t(m_matrix.outerIndexPtr()[m_outerStart]); + else if(m_outerSize.value()==0) + return 0; + else + return Map<const Matrix<Index,OuterSize,1> >(m_matrix.innerNonZeroPtr()+m_outerStart, m_outerSize.value()).sum(); + } + + const Scalar& lastCoeff() const + { + EIGEN_STATIC_ASSERT_VECTOR_ONLY(BlockImpl); + eigen_assert(nonZeros()>0); + if(m_matrix.isCompressed()) + return m_matrix.valuePtr()[m_matrix.outerIndexPtr()[m_outerStart+1]-1]; + else + return m_matrix.valuePtr()[m_matrix.outerIndexPtr()[m_outerStart]+m_matrix.innerNonZeroPtr()[m_outerStart]-1]; + } + + EIGEN_STRONG_INLINE Index rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); } + EIGEN_STRONG_INLINE Index cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); } + + protected: + + typename SparseMatrixType::Nested m_matrix; + Index m_outerStart; + const internal::variable_if_dynamic<Index, OuterSize> m_outerSize; + +}; + //---------- /** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this diff --git a/Eigen/src/SparseCore/SparseDenseProduct.h b/Eigen/src/SparseCore/SparseDenseProduct.h index 78411db98..a9084192e 100644 --- a/Eigen/src/SparseCore/SparseDenseProduct.h +++ b/Eigen/src/SparseCore/SparseDenseProduct.h @@ -306,15 +306,6 @@ class DenseTimeSparseProduct DenseTimeSparseProduct& operator=(const DenseTimeSparseProduct&); }; -// sparse * dense -template<typename Derived> -template<typename OtherDerived> -inline const typename SparseDenseProductReturnType<Derived,OtherDerived>::Type -SparseMatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const -{ - return typename SparseDenseProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived()); -} - } // end namespace Eigen #endif // EIGEN_SPARSEDENSEPRODUCT_H diff --git a/Eigen/src/SparseCore/SparseMatrixBase.h b/Eigen/src/SparseCore/SparseMatrixBase.h index bbcf7fb1c..485e85bec 100644 --- a/Eigen/src/SparseCore/SparseMatrixBase.h +++ b/Eigen/src/SparseCore/SparseMatrixBase.h @@ -358,7 +358,8 @@ template<typename Derived> class SparseMatrixBase : public EigenBase<Derived> /** sparse * dense (returns a dense object unless it is an outer product) */ template<typename OtherDerived> const typename SparseDenseProductReturnType<Derived,OtherDerived>::Type - operator*(const MatrixBase<OtherDerived> &other) const; + operator*(const MatrixBase<OtherDerived> &other) const + { return typename SparseDenseProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived()); } /** \returns an expression of P H P^-1 where H is the matrix represented by \c *this */ SparseSymmetricPermutationProduct<Derived,Upper|Lower> twistedBy(const PermutationMatrix<Dynamic,Dynamic,Index>& perm) const diff --git a/Eigen/src/SparseCore/SparsePermutation.h b/Eigen/src/SparseCore/SparsePermutation.h index b85be93f6..75e210009 100644 --- a/Eigen/src/SparseCore/SparsePermutation.h +++ b/Eigen/src/SparseCore/SparsePermutation.h @@ -61,7 +61,7 @@ struct permut_sparsematrix_product_retval for(Index j=0; j<m_matrix.outerSize(); ++j) { Index jp = m_permutation.indices().coeff(j); - sizes[((Side==OnTheLeft) ^ Transposed) ? jp : j] = m_matrix.innerVector(((Side==OnTheRight) ^ Transposed) ? jp : j).size(); + sizes[((Side==OnTheLeft) ^ Transposed) ? jp : j] = m_matrix.innerVector(((Side==OnTheRight) ^ Transposed) ? jp : j).nonZeros(); } tmp.reserve(sizes); for(Index j=0; j<m_matrix.outerSize(); ++j) diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h index 1d592f2c8..d6f5d59ca 100644 --- a/Eigen/src/SparseLU/SparseLU.h +++ b/Eigen/src/SparseLU/SparseLU.h @@ -260,14 +260,13 @@ class SparseLU : public internal::SparseLUImpl<typename _MatrixType::Scalar, typ eigen_assert(m_factorizationIsOk && "The matrix should be factorized first."); // Initialize with the determinant of the row matrix Scalar det = Scalar(1.); - //Note that the diagonal blocks of U are stored in supernodes, + // Note that the diagonal blocks of U are stored in supernodes, // which are available in the L part :) for (Index j = 0; j < this->cols(); ++j) { for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it) { - if(it.row() < j) continue; - if(it.row() == j) + if(it.index() == j) { det *= (std::abs)(it.value()); break; diff --git a/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h b/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h index ad6f2183f..b16afd6a4 100644 --- a/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h +++ b/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h @@ -189,8 +189,8 @@ class MappedSuperNodalMatrix<Scalar,Index>::InnerIterator m_idval(mat.colIndexPtr()[outer]), m_startidval(m_idval), m_endidval(mat.colIndexPtr()[outer+1]), - m_idrow(mat.rowIndexPtr()[outer]), - m_endidrow(mat.rowIndexPtr()[outer+1]) + m_idrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[outer]]]), + m_endidrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[outer]]+1]) {} inline InnerIterator& operator++() { diff --git a/Eigen/src/SparseQR/SparseQR.h b/Eigen/src/SparseQR/SparseQR.h index 4c6553bf2..a00bd5db1 100644 --- a/Eigen/src/SparseQR/SparseQR.h +++ b/Eigen/src/SparseQR/SparseQR.h @@ -75,7 +75,7 @@ class SparseQR typedef Matrix<Scalar, Dynamic, 1> ScalarVector; typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType; public: - SparseQR () : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false) + SparseQR () : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false) { } /** Construct a QR factorization of the matrix \a mat. @@ -84,7 +84,7 @@ class SparseQR * * \sa compute() */ - SparseQR(const MatrixType& mat) : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false) + SparseQR(const MatrixType& mat) : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false) { compute(mat); } @@ -262,6 +262,7 @@ class SparseQR IndexVector m_etree; // Column elimination tree IndexVector m_firstRowElt; // First element in each row bool m_isQSorted; // whether Q is sorted or not + bool m_isEtreeOk; // whether the elimination tree match the initial input matrix template <typename, typename > friend struct SparseQR_QProduct; template <typename > friend struct SparseQRMatrixQReturnType; @@ -281,9 +282,11 @@ template <typename MatrixType, typename OrderingType> void SparseQR<MatrixType,OrderingType>::analyzePattern(const MatrixType& mat) { eigen_assert(mat.isCompressed() && "SparseQR requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to SparseQR"); + // Copy to a column major matrix if the input is rowmajor + typename internal::conditional<MatrixType::IsRowMajor,QRMatrixType,const MatrixType&>::type matCpy(mat); // Compute the column fill reducing ordering OrderingType ord; - ord(mat, m_perm_c); + ord(matCpy, m_perm_c); Index n = mat.cols(); Index m = mat.rows(); Index diagSize = (std::min)(m,n); @@ -296,7 +299,8 @@ void SparseQR<MatrixType,OrderingType>::analyzePattern(const MatrixType& mat) // Compute the column elimination tree of the permuted matrix m_outputPerm_c = m_perm_c.inverse(); - internal::coletree(mat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data()); + internal::coletree(matCpy, m_etree, m_firstRowElt, m_outputPerm_c.indices().data()); + m_isEtreeOk = true; m_R.resize(m, n); m_Q.resize(m, diagSize); @@ -330,15 +334,38 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat) Index nzcolR, nzcolQ; // Number of nonzero for the current column of R and Q ScalarVector tval(m); // The dense vector used to compute the current column RealScalar pivotThreshold = m_threshold; - + + m_R.setZero(); + m_Q.setZero(); m_pmat = mat; + if(!m_isEtreeOk) + { + m_outputPerm_c = m_perm_c.inverse(); + internal::coletree(m_pmat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data()); + m_isEtreeOk = true; + } + m_pmat.uncompress(); // To have the innerNonZeroPtr allocated + // Apply the fill-in reducing permutation lazily: - for (int i = 0; i < n; i++) { - Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i; - m_pmat.outerIndexPtr()[p] = mat.outerIndexPtr()[i]; - m_pmat.innerNonZeroPtr()[p] = mat.outerIndexPtr()[i+1] - mat.outerIndexPtr()[i]; + // If the input is row major, copy the original column indices, + // otherwise directly use the input matrix + // + IndexVector originalOuterIndicesCpy; + const Index *originalOuterIndices = mat.outerIndexPtr(); + if(MatrixType::IsRowMajor) + { + originalOuterIndicesCpy = IndexVector::Map(m_pmat.outerIndexPtr(),n+1); + originalOuterIndices = originalOuterIndicesCpy.data(); + } + + for (int i = 0; i < n; i++) + { + Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i; + m_pmat.outerIndexPtr()[p] = originalOuterIndices[i]; + m_pmat.innerNonZeroPtr()[p] = originalOuterIndices[i+1] - originalOuterIndices[i]; + } } /* Compute the default threshold as in MatLab, see: @@ -349,6 +376,8 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat) { RealScalar max2Norm = 0.0; for (int j = 0; j < n; j++) max2Norm = (max)(max2Norm, m_pmat.col(j).norm()); + if(max2Norm==RealScalar(0)) + max2Norm = RealScalar(1); pivotThreshold = 20 * (m + n) * max2Norm * NumTraits<RealScalar>::epsilon(); } @@ -357,7 +386,7 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat) Index nonzeroCol = 0; // Record the number of valid pivots m_Q.startVec(0); - + // Left looking rank-revealing QR factorization: compute a column of R and Q at a time for (Index col = 0; col < n; ++col) { @@ -373,7 +402,7 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat) // all the nodes (with indexes lower than rank) reachable through the column elimination tree (etree) rooted at node k. // Note: if the diagonal entry does not exist, then its contribution must be explicitly added, // thus the trick with found_diag that permits to do one more iteration on the diagonal element if this one has not been found. - for (typename MatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp) + for (typename QRMatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp) { Index curIdx = nonzeroCol; if(itp) curIdx = itp.row(); @@ -447,7 +476,7 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat) } } // End update current column - Scalar tau; + Scalar tau = 0; RealScalar beta = 0; if(nonzeroCol < diagSize) @@ -461,7 +490,6 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat) for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq))); if(sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0)) { - tau = RealScalar(0); beta = numext::real(c0); tval(Qidx(0)) = 1; } @@ -514,6 +542,7 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat) // Recompute the column elimination tree internal::coletree(m_pmat, m_etree, m_firstRowElt, m_pivotperm.indices().data()); + m_isEtreeOk = false; } } @@ -525,13 +554,13 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat) m_R.finalize(); m_R.makeCompressed(); m_isQSorted = false; - + m_nonzeropivots = nonzeroCol; if(nonzeroCol<n) { // Permute the triangular factor to put the 'dead' columns to the end - MatrixType tempR(m_R); + QRMatrixType tempR(m_R); m_R = tempR * m_pivotperm; // Update the column permutation diff --git a/Eigen/src/UmfPackSupport/UmfPackSupport.h b/Eigen/src/UmfPackSupport/UmfPackSupport.h index 3a48cecf7..29c60c378 100644 --- a/Eigen/src/UmfPackSupport/UmfPackSupport.h +++ b/Eigen/src/UmfPackSupport/UmfPackSupport.h @@ -107,6 +107,16 @@ inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *N return umfpack_zi_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info); } +namespace internal { + template<typename T> struct umfpack_helper_is_sparse_plain : false_type {}; + template<typename Scalar, int Options, typename StorageIndex> + struct umfpack_helper_is_sparse_plain<SparseMatrix<Scalar,Options,StorageIndex> > + : true_type {}; + template<typename Scalar, int Options, typename StorageIndex> + struct umfpack_helper_is_sparse_plain<MappedSparseMatrix<Scalar,Options,StorageIndex> > + : true_type {}; +} + /** \ingroup UmfPackSupport_Module * \brief A sparse LU factorization and solver based on UmfPack * @@ -192,10 +202,14 @@ class UmfPackLU : internal::noncopyable * Note that the matrix should be column-major, and in compressed format for best performance. * \sa SparseMatrix::makeCompressed(). */ - void compute(const MatrixType& matrix) + template<typename InputMatrixType> + void compute(const InputMatrixType& matrix) { - analyzePattern(matrix); - factorize(matrix); + if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar()); + if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar()); + grapInput(matrix.derived()); + analyzePattern_impl(); + factorize_impl(); } /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. @@ -230,23 +244,15 @@ class UmfPackLU : internal::noncopyable * * \sa factorize(), compute() */ - void analyzePattern(const MatrixType& matrix) + template<typename InputMatrixType> + void analyzePattern(const InputMatrixType& matrix) { - if(m_symbolic) - umfpack_free_symbolic(&m_symbolic,Scalar()); - if(m_numeric) - umfpack_free_numeric(&m_numeric,Scalar()); + if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar()); + if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar()); - grapInput(matrix); - - int errorCode = 0; - errorCode = umfpack_symbolic(matrix.rows(), matrix.cols(), m_outerIndexPtr, m_innerIndexPtr, m_valuePtr, - &m_symbolic, 0, 0); + grapInput(matrix.derived()); - m_isInitialized = true; - m_info = errorCode ? InvalidInput : Success; - m_analysisIsOk = true; - m_factorizationIsOk = false; + analyzePattern_impl(); } /** Performs a numeric decomposition of \a matrix @@ -255,20 +261,16 @@ class UmfPackLU : internal::noncopyable * * \sa analyzePattern(), compute() */ - void factorize(const MatrixType& matrix) + template<typename InputMatrixType> + void factorize(const InputMatrixType& matrix) { eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()"); if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar()); - grapInput(matrix); - - int errorCode; - errorCode = umfpack_numeric(m_outerIndexPtr, m_innerIndexPtr, m_valuePtr, - m_symbolic, &m_numeric, 0, 0); - - m_info = errorCode ? NumericalIssue : Success; - m_factorizationIsOk = true; + grapInput(matrix.derived()); + + factorize_impl(); } #ifndef EIGEN_PARSED_BY_DOXYGEN @@ -283,19 +285,20 @@ class UmfPackLU : internal::noncopyable protected: - void init() { - m_info = InvalidInput; - m_isInitialized = false; - m_numeric = 0; - m_symbolic = 0; - m_outerIndexPtr = 0; - m_innerIndexPtr = 0; - m_valuePtr = 0; + m_info = InvalidInput; + m_isInitialized = false; + m_numeric = 0; + m_symbolic = 0; + m_outerIndexPtr = 0; + m_innerIndexPtr = 0; + m_valuePtr = 0; + m_extractedDataAreDirty = true; } - void grapInput(const MatrixType& mat) + template<typename InputMatrixType> + void grapInput_impl(const InputMatrixType& mat, internal::true_type) { m_copyMatrix.resize(mat.rows(), mat.cols()); if( ((MatrixType::Flags&RowMajorBit)==RowMajorBit) || sizeof(typename MatrixType::Index)!=sizeof(int) || !mat.isCompressed() ) @@ -313,6 +316,45 @@ class UmfPackLU : internal::noncopyable m_valuePtr = mat.valuePtr(); } } + + template<typename InputMatrixType> + void grapInput_impl(const InputMatrixType& mat, internal::false_type) + { + m_copyMatrix = mat; + m_outerIndexPtr = m_copyMatrix.outerIndexPtr(); + m_innerIndexPtr = m_copyMatrix.innerIndexPtr(); + m_valuePtr = m_copyMatrix.valuePtr(); + } + + template<typename InputMatrixType> + void grapInput(const InputMatrixType& mat) + { + grapInput_impl(mat, internal::umfpack_helper_is_sparse_plain<InputMatrixType>()); + } + + void analyzePattern_impl() + { + int errorCode = 0; + errorCode = umfpack_symbolic(m_copyMatrix.rows(), m_copyMatrix.cols(), m_outerIndexPtr, m_innerIndexPtr, m_valuePtr, + &m_symbolic, 0, 0); + + m_isInitialized = true; + m_info = errorCode ? InvalidInput : Success; + m_analysisIsOk = true; + m_factorizationIsOk = false; + m_extractedDataAreDirty = true; + } + + void factorize_impl() + { + int errorCode; + errorCode = umfpack_numeric(m_outerIndexPtr, m_innerIndexPtr, m_valuePtr, + m_symbolic, &m_numeric, 0, 0); + + m_info = errorCode ? NumericalIssue : Success; + m_factorizationIsOk = true; + m_extractedDataAreDirty = true; + } // cached data to reduce reallocation, etc. mutable LUMatrixType m_l; diff --git a/cmake/EigenTesting.cmake b/cmake/EigenTesting.cmake index d9e22ab1a..80b2841df 100644 --- a/cmake/EigenTesting.cmake +++ b/cmake/EigenTesting.cmake @@ -452,20 +452,12 @@ macro(ei_set_build_string) endmacro(ei_set_build_string) macro(ei_is_64bit_env VAR) - - file(WRITE "${CMAKE_CURRENT_BINARY_DIR}/is64.cpp" - "int main() { return (sizeof(int*) == 8 ? 1 : 0); } - ") - try_run(run_res compile_res - ${CMAKE_CURRENT_BINARY_DIR} "${CMAKE_CURRENT_BINARY_DIR}/is64.cpp" - RUN_OUTPUT_VARIABLE run_output) - - if(compile_res AND run_res) - set(${VAR} ${run_res}) - elseif(CMAKE_CL_64) - set(${VAR} 1) - elseif("$ENV{Platform}" STREQUAL "X64") # nmake 64 bit + if(CMAKE_SIZEOF_VOID_P EQUAL 8) set(${VAR} 1) + elseif(CMAKE_SIZEOF_VOID_P EQUAL 4) + set(${VAR} 0) + else() + message(WARNING "Unsupported pointer size. Please contact the authors.") endif() endmacro(ei_is_64bit_env) diff --git a/cmake/FindCholmod.cmake b/cmake/FindCholmod.cmake index 7b3046d45..23239c300 100644 --- a/cmake/FindCholmod.cmake +++ b/cmake/FindCholmod.cmake @@ -86,4 +86,4 @@ include(FindPackageHandleStandardArgs) find_package_handle_standard_args(CHOLMOD DEFAULT_MSG CHOLMOD_INCLUDES CHOLMOD_LIBRARIES) -mark_as_advanced(CHOLMOD_INCLUDES CHOLMOD_LIBRARIES AMD_LIBRARY COLAMD_LIBRARY SUITESPARSE_LIBRARY) +mark_as_advanced(CHOLMOD_INCLUDES CHOLMOD_LIBRARIES AMD_LIBRARY COLAMD_LIBRARY SUITESPARSE_LIBRARY CAMD_LIBRARY CCOLAMD_LIBRARY CHOLMOD_METIS_LIBRARY) diff --git a/cmake/FindFFTW.cmake b/cmake/FindFFTW.cmake index a9e9925e7..6c4dc9ab4 100644 --- a/cmake/FindFFTW.cmake +++ b/cmake/FindFFTW.cmake @@ -115,5 +115,5 @@ include(FindPackageHandleStandardArgs) find_package_handle_standard_args(FFTW DEFAULT_MSG FFTW_INCLUDES FFTW_LIBRARIES) -mark_as_advanced(FFTW_INCLUDES FFTW_LIBRARIES) +mark_as_advanced(FFTW_INCLUDES FFTW_LIBRARIES FFTW_LIB FFTWF_LIB FFTWL_LIB) diff --git a/cmake/FindMetis.cmake b/cmake/FindMetis.cmake index 627c3e9ae..e0040d320 100644 --- a/cmake/FindMetis.cmake +++ b/cmake/FindMetis.cmake @@ -10,16 +10,50 @@ find_path(METIS_INCLUDES PATHS $ENV{METISDIR} ${INCLUDE_INSTALL_DIR} - PATH_SUFFIXES + PATH_SUFFIXES + . metis include ) +macro(_metis_check_version) + file(READ "${METIS_INCLUDES}/metis.h" _metis_version_header) + + string(REGEX MATCH "define[ \t]+METIS_VER_MAJOR[ \t]+([0-9]+)" _metis_major_version_match "${_metis_version_header}") + set(METIS_MAJOR_VERSION "${CMAKE_MATCH_1}") + string(REGEX MATCH "define[ \t]+METIS_VER_MINOR[ \t]+([0-9]+)" _metis_minor_version_match "${_metis_version_header}") + set(METIS_MINOR_VERSION "${CMAKE_MATCH_1}") + string(REGEX MATCH "define[ \t]+METIS_VER_SUBMINOR[ \t]+([0-9]+)" _metis_subminor_version_match "${_metis_version_header}") + set(METIS_SUBMINOR_VERSION "${CMAKE_MATCH_1}") + if(NOT METIS_MAJOR_VERSION) + message(WARNING "Could not determine Metis version. Assuming version 4.0.0") + set(METIS_VERSION 4.0.0) + else() + set(METIS_VERSION ${METIS_MAJOR_VERSION}.${METIS_MINOR_VERSION}.${METIS_SUBMINOR_VERSION}) + endif() + if(${METIS_VERSION} VERSION_LESS ${Metis_FIND_VERSION}) + set(METIS_VERSION_OK FALSE) + else() + set(METIS_VERSION_OK TRUE) + endif() + + if(NOT METIS_VERSION_OK) + message(STATUS "Metis version ${METIS_VERSION} found in ${METIS_INCLUDES}, " + "but at least version ${Metis_FIND_VERSION} is required") + endif(NOT METIS_VERSION_OK) +endmacro(_metis_check_version) + + if(METIS_INCLUDES AND Metis_FIND_VERSION) + _metis_check_version() + else() + set(METIS_VERSION_OK TRUE) + endif() + find_library(METIS_LIBRARIES metis PATHS $ENV{METISDIR} ${LIB_INSTALL_DIR} PATH_SUFFIXES lib) include(FindPackageHandleStandardArgs) find_package_handle_standard_args(METIS DEFAULT_MSG - METIS_INCLUDES METIS_LIBRARIES) + METIS_INCLUDES METIS_LIBRARIES METIS_VERSION_OK) mark_as_advanced(METIS_INCLUDES METIS_LIBRARIES) diff --git a/cmake/language_support.cmake b/cmake/language_support.cmake index d687b71f6..231f7ff70 100644 --- a/cmake/language_support.cmake +++ b/cmake/language_support.cmake @@ -33,7 +33,7 @@ function(workaround_9220 language language_works) file(WRITE ${CMAKE_BINARY_DIR}/language_tests/${language}/CMakeLists.txt ${text}) execute_process( - COMMAND ${CMAKE_COMMAND} . + COMMAND ${CMAKE_COMMAND} . -G "${CMAKE_GENERATOR}" WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/language_tests/${language} RESULT_VARIABLE return_code OUTPUT_QUIET diff --git a/doc/AsciiQuickReference.txt b/doc/AsciiQuickReference.txt index 1e74e0528..b9f497f87 100644 --- a/doc/AsciiQuickReference.txt +++ b/doc/AsciiQuickReference.txt @@ -67,10 +67,10 @@ P.rightCols<cols>() // P(:, end-cols+1:end) P.rightCols(cols) // P(:, end-cols+1:end) P.topRows<rows>() // P(1:rows, :) P.topRows(rows) // P(1:rows, :) -P.middleRows<rows>(i) // P(:, i+1:i+rows) -P.middleRows(i, rows) // P(:, i+1:i+rows) -P.bottomRows<rows>() // P(:, end-rows+1:end) -P.bottomRows(rows) // P(:, end-rows+1:end) +P.middleRows<rows>(i) // P(i+1:i+rows, :) +P.middleRows(i, rows) // P(i+1:i+rows, :) +P.bottomRows<rows>() // P(end-rows+1:end, :) +P.bottomRows(rows) // P(end-rows+1:end, :) P.topLeftCorner(rows, cols) // P(1:rows, 1:cols) P.topRightCorner(rows, cols) // P(1:rows, end-cols+1:end) P.bottomLeftCorner(rows, cols) // P(end-rows+1:end, 1:cols) diff --git a/doc/SparseQuickReference.dox b/doc/SparseQuickReference.dox index 4a33d0cc9..d04ac35c5 100644 --- a/doc/SparseQuickReference.dox +++ b/doc/SparseQuickReference.dox @@ -71,11 +71,10 @@ i.e either row major or column major. The default is column major. Most arithmet <td> Constant or Random Insertion</td> <td> \code -sm1.setZero(); // Set the matrix with zero elements -sm1.setConstant(val); //Replace all the nonzero values with val +sm1.setZero(); \endcode </td> -<td> The matrix sm1 should have been created before ???</td> +<td>Remove all non-zero coefficients</td> </tr> </table> diff --git a/doc/examples/CMakeLists.txt b/doc/examples/CMakeLists.txt index 71b272a15..08cf8efd7 100644 --- a/doc/examples/CMakeLists.txt +++ b/doc/examples/CMakeLists.txt @@ -6,12 +6,10 @@ foreach(example_src ${examples_SRCS}) if(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO) target_link_libraries(${example} ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO}) endif() - get_target_property(example_executable - ${example} LOCATION) add_custom_command( TARGET ${example} POST_BUILD - COMMAND ${example_executable} + COMMAND ${example} ARGS >${CMAKE_CURRENT_BINARY_DIR}/${example}.out ) add_dependencies(all_examples ${example}) diff --git a/doc/snippets/CMakeLists.txt b/doc/snippets/CMakeLists.txt index 92a22ea61..1135900cf 100644 --- a/doc/snippets/CMakeLists.txt +++ b/doc/snippets/CMakeLists.txt @@ -14,12 +14,10 @@ foreach(snippet_src ${snippets_SRCS}) if(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO) target_link_libraries(${compile_snippet_target} ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO}) endif() - get_target_property(compile_snippet_executable - ${compile_snippet_target} LOCATION) add_custom_command( TARGET ${compile_snippet_target} POST_BUILD - COMMAND ${compile_snippet_executable} + COMMAND ${compile_snippet_target} ARGS >${CMAKE_CURRENT_BINARY_DIR}/${snippet}.out ) add_dependencies(all_snippets ${compile_snippet_target}) @@ -27,4 +25,4 @@ foreach(snippet_src ${snippets_SRCS}) PROPERTIES OBJECT_DEPENDS ${snippet_src}) endforeach(snippet_src) -ei_add_target_property(compile_tut_arithmetic_transpose_aliasing COMPILE_FLAGS -DEIGEN_NO_DEBUG)
\ No newline at end of file +ei_add_target_property(compile_tut_arithmetic_transpose_aliasing COMPILE_FLAGS -DEIGEN_NO_DEBUG) diff --git a/doc/special_examples/CMakeLists.txt b/doc/special_examples/CMakeLists.txt index 0c9b3c3ba..f8a0148c8 100644 --- a/doc/special_examples/CMakeLists.txt +++ b/doc/special_examples/CMakeLists.txt @@ -1,4 +1,3 @@ - if(NOT EIGEN_TEST_NOQT) find_package(Qt4) if(QT4_FOUND) @@ -6,16 +5,16 @@ if(NOT EIGEN_TEST_NOQT) endif() endif(NOT EIGEN_TEST_NOQT) - if(QT4_FOUND) add_executable(Tutorial_sparse_example Tutorial_sparse_example.cpp Tutorial_sparse_example_details.cpp) target_link_libraries(Tutorial_sparse_example ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO} ${QT_QTCORE_LIBRARY} ${QT_QTGUI_LIBRARY}) - + add_custom_command( TARGET Tutorial_sparse_example POST_BUILD COMMAND Tutorial_sparse_example ARGS ${CMAKE_CURRENT_BINARY_DIR}/../html/Tutorial_sparse_example.jpeg ) - + add_dependencies(all_examples Tutorial_sparse_example) endif(QT4_FOUND) + diff --git a/failtest/CMakeLists.txt b/failtest/CMakeLists.txt index 5afa2ac82..5c4860237 100644 --- a/failtest/CMakeLists.txt +++ b/failtest/CMakeLists.txt @@ -26,6 +26,12 @@ ei_add_failtest("block_on_const_type_actually_const_1") ei_add_failtest("transpose_on_const_type_actually_const") ei_add_failtest("diagonal_on_const_type_actually_const") +ei_add_failtest("ref_1") +ei_add_failtest("ref_2") +ei_add_failtest("ref_3") +ei_add_failtest("ref_4") +ei_add_failtest("ref_5") + if (EIGEN_FAILTEST_FAILURE_COUNT) message(FATAL_ERROR "${EIGEN_FAILTEST_FAILURE_COUNT} out of ${EIGEN_FAILTEST_COUNT} failtests FAILED. " diff --git a/failtest/ref_1.cpp b/failtest/ref_1.cpp new file mode 100644 index 000000000..8b798d53d --- /dev/null +++ b/failtest/ref_1.cpp @@ -0,0 +1,18 @@ +#include "../Eigen/Core" + +#ifdef EIGEN_SHOULD_FAIL_TO_BUILD +#define CV_QUALIFIER const +#else +#define CV_QUALIFIER +#endif + +using namespace Eigen; + +void call_ref(Ref<VectorXf> a) { } + +int main() +{ + VectorXf a(10); + CV_QUALIFIER VectorXf& ac(a); + call_ref(ac); +} diff --git a/failtest/ref_2.cpp b/failtest/ref_2.cpp new file mode 100644 index 000000000..0b779ccf5 --- /dev/null +++ b/failtest/ref_2.cpp @@ -0,0 +1,15 @@ +#include "../Eigen/Core" + +using namespace Eigen; + +void call_ref(Ref<VectorXf> a) { } + +int main() +{ + MatrixXf A(10,10); +#ifdef EIGEN_SHOULD_FAIL_TO_BUILD + call_ref(A.row(3)); +#else + call_ref(A.col(3)); +#endif +} diff --git a/failtest/ref_3.cpp b/failtest/ref_3.cpp new file mode 100644 index 000000000..f46027d48 --- /dev/null +++ b/failtest/ref_3.cpp @@ -0,0 +1,15 @@ +#include "../Eigen/Core" + +using namespace Eigen; + +#ifdef EIGEN_SHOULD_FAIL_TO_BUILD +void call_ref(Ref<VectorXf> a) { } +#else +void call_ref(const Ref<const VectorXf> &a) { } +#endif + +int main() +{ + VectorXf a(10); + call_ref(a+a); +} diff --git a/failtest/ref_4.cpp b/failtest/ref_4.cpp new file mode 100644 index 000000000..6c11fa4cb --- /dev/null +++ b/failtest/ref_4.cpp @@ -0,0 +1,15 @@ +#include "../Eigen/Core" + +using namespace Eigen; + +void call_ref(Ref<MatrixXf,0,OuterStride<> > a) {} + +int main() +{ + MatrixXf A(10,10); +#ifdef EIGEN_SHOULD_FAIL_TO_BUILD + call_ref(A.transpose()); +#else + call_ref(A); +#endif +} diff --git a/failtest/ref_5.cpp b/failtest/ref_5.cpp new file mode 100644 index 000000000..846d52795 --- /dev/null +++ b/failtest/ref_5.cpp @@ -0,0 +1,16 @@ +#include "../Eigen/Core" + +using namespace Eigen; + +void call_ref(Ref<VectorXf> a) { } + +int main() +{ + VectorXf a(10); + DenseBase<VectorXf> &ac(a); +#ifdef EIGEN_SHOULD_FAIL_TO_BUILD + call_ref(ac); +#else + call_ref(ac.derived()); +#endif +} diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index ccb0fc798..d32451df6 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -66,7 +66,7 @@ endif() find_package(Pastix) find_package(Scotch) -find_package(Metis) +find_package(Metis 5.0 REQUIRED) if(PASTIX_FOUND AND BLAS_FOUND) add_definitions("-DEIGEN_PASTIX_SUPPORT") include_directories(${PASTIX_INCLUDES}) @@ -279,6 +279,7 @@ ei_add_property(EIGEN_TESTING_SUMMARY "CXX_FLAGS: ${CMAKE_CXX_FLAGS}\n") ei_add_property(EIGEN_TESTING_SUMMARY "Sparse lib flags: ${SPARSE_LIBS}\n") option(EIGEN_TEST_EIGEN2 "Run whole Eigen2 test suite against EIGEN2_SUPPORT" OFF) +mark_as_advanced(EIGEN_TEST_EIGEN2) if(EIGEN_TEST_EIGEN2) add_subdirectory(eigen2) endif() diff --git a/test/cholesky.cpp b/test/cholesky.cpp index 64bcbccc4..56885deb7 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -320,33 +320,35 @@ template<typename MatrixType> void cholesky_definiteness(const MatrixType& m) { eigen_assert(m.rows() == 2 && m.cols() == 2); MatrixType mat; + LDLT<MatrixType> ldlt(2); + { mat << 1, 0, 0, -1; - LDLT<MatrixType> ldlt(mat); + ldlt.compute(mat); VERIFY(!ldlt.isNegative()); VERIFY(!ldlt.isPositive()); } { mat << 1, 2, 2, 1; - LDLT<MatrixType> ldlt(mat); + ldlt.compute(mat); VERIFY(!ldlt.isNegative()); VERIFY(!ldlt.isPositive()); } { mat << 0, 0, 0, 0; - LDLT<MatrixType> ldlt(mat); + ldlt.compute(mat); VERIFY(ldlt.isNegative()); VERIFY(ldlt.isPositive()); } { mat << 0, 0, 0, 1; - LDLT<MatrixType> ldlt(mat); + ldlt.compute(mat); VERIFY(!ldlt.isNegative()); VERIFY(ldlt.isPositive()); } { mat << -1, 0, 0, 0; - LDLT<MatrixType> ldlt(mat); + ldlt.compute(mat); VERIFY(ldlt.isNegative()); VERIFY(!ldlt.isPositive()); } diff --git a/test/cwiseop.cpp b/test/cwiseop.cpp index 247fa2a09..e3361da17 100644 --- a/test/cwiseop.cpp +++ b/test/cwiseop.cpp @@ -9,6 +9,8 @@ // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. #define EIGEN2_SUPPORT +#define EIGEN_NO_EIGEN2_DEPRECATED_WARNING + #define EIGEN_NO_STATIC_ASSERT #include "main.h" #include <functional> diff --git a/test/eigen2/CMakeLists.txt b/test/eigen2/CMakeLists.txt index 84931e037..9615a6026 100644 --- a/test/eigen2/CMakeLists.txt +++ b/test/eigen2/CMakeLists.txt @@ -4,6 +4,7 @@ add_dependencies(eigen2_check eigen2_buildtests) add_dependencies(buildtests eigen2_buildtests) add_definitions("-DEIGEN2_SUPPORT_STAGE10_FULL_EIGEN2_API") +add_definitions("-DEIGEN_NO_EIGEN2_DEPRECATED_WARNING") ei_add_test(eigen2_meta) ei_add_test(eigen2_sizeof) diff --git a/test/eigen2/eigen2_adjoint.cpp b/test/eigen2/eigen2_adjoint.cpp index 8ec9c9202..c0f811459 100644 --- a/test/eigen2/eigen2_adjoint.cpp +++ b/test/eigen2/eigen2_adjoint.cpp @@ -29,8 +29,6 @@ template<typename MatrixType> void adjoint(const MatrixType& m) MatrixType m1 = MatrixType::Random(rows, cols), m2 = MatrixType::Random(rows, cols), m3(rows, cols), - mzero = MatrixType::Zero(rows, cols), - identity = SquareMatrixType::Identity(rows, rows), square = SquareMatrixType::Random(rows, rows); VectorType v1 = VectorType::Random(rows), v2 = VectorType::Random(rows), diff --git a/test/eigen2/eigen2_basicstuff.cpp b/test/eigen2/eigen2_basicstuff.cpp index 4fa16d5ae..dd2dec1ef 100644 --- a/test/eigen2/eigen2_basicstuff.cpp +++ b/test/eigen2/eigen2_basicstuff.cpp @@ -23,11 +23,8 @@ template<typename MatrixType> void basicStuff(const MatrixType& m) m2 = MatrixType::Random(rows, cols), m3(rows, cols), mzero = MatrixType::Zero(rows, cols), - identity = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> - ::Identity(rows, rows), square = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>::Random(rows, rows); VectorType v1 = VectorType::Random(rows), - v2 = VectorType::Random(rows), vzero = VectorType::Zero(rows); Scalar x = ei_random<Scalar>(); diff --git a/test/eigen2/eigen2_cwiseop.cpp b/test/eigen2/eigen2_cwiseop.cpp index 4391d19b4..22e1cc342 100644 --- a/test/eigen2/eigen2_cwiseop.cpp +++ b/test/eigen2/eigen2_cwiseop.cpp @@ -35,11 +35,8 @@ template<typename MatrixType> void cwiseops(const MatrixType& m) mzero = MatrixType::Zero(rows, cols), mones = MatrixType::Ones(rows, cols), identity = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> - ::Identity(rows, rows), - square = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>::Random(rows, rows); - VectorType v1 = VectorType::Random(rows), - v2 = VectorType::Random(rows), - vzero = VectorType::Zero(rows), + ::Identity(rows, rows); + VectorType vzero = VectorType::Zero(rows), vones = VectorType::Ones(rows), v3(rows); diff --git a/test/eigen2/eigen2_geometry.cpp b/test/eigen2/eigen2_geometry.cpp index 70b4ab5f1..514040774 100644 --- a/test/eigen2/eigen2_geometry.cpp +++ b/test/eigen2/eigen2_geometry.cpp @@ -392,6 +392,7 @@ template<typename Scalar> void geometry(void) #define VERIFY_EULER(I,J,K, X,Y,Z) { \ Vector3 ea = m.eulerAngles(I,J,K); \ Matrix3 m1 = Matrix3(AngleAxisx(ea[0], Vector3::Unit##X()) * AngleAxisx(ea[1], Vector3::Unit##Y()) * AngleAxisx(ea[2], Vector3::Unit##Z())); \ + VERIFY_IS_APPROX(m, m1); \ VERIFY_IS_APPROX(m, Matrix3(AngleAxisx(ea[0], Vector3::Unit##X()) * AngleAxisx(ea[1], Vector3::Unit##Y()) * AngleAxisx(ea[2], Vector3::Unit##Z()))); \ } VERIFY_EULER(0,1,2, X,Y,Z); diff --git a/test/eigen2/eigen2_geometry_with_eigen2_prefix.cpp b/test/eigen2/eigen2_geometry_with_eigen2_prefix.cpp index f83b57249..12d4a71c3 100644 --- a/test/eigen2/eigen2_geometry_with_eigen2_prefix.cpp +++ b/test/eigen2/eigen2_geometry_with_eigen2_prefix.cpp @@ -394,6 +394,7 @@ template<typename Scalar> void geometry(void) #define VERIFY_EULER(I,J,K, X,Y,Z) { \ Vector3 ea = m.eulerAngles(I,J,K); \ Matrix3 m1 = Matrix3(AngleAxisx(ea[0], Vector3::Unit##X()) * AngleAxisx(ea[1], Vector3::Unit##Y()) * AngleAxisx(ea[2], Vector3::Unit##Z())); \ + VERIFY_IS_APPROX(m, m1); \ VERIFY_IS_APPROX(m, Matrix3(AngleAxisx(ea[0], Vector3::Unit##X()) * AngleAxisx(ea[1], Vector3::Unit##Y()) * AngleAxisx(ea[2], Vector3::Unit##Z()))); \ } VERIFY_EULER(0,1,2, X,Y,Z); diff --git a/test/eigen2/eigen2_inverse.cpp b/test/eigen2/eigen2_inverse.cpp index 5de1dfefa..ccd24a194 100644 --- a/test/eigen2/eigen2_inverse.cpp +++ b/test/eigen2/eigen2_inverse.cpp @@ -25,7 +25,6 @@ template<typename MatrixType> void inverse(const MatrixType& m) MatrixType m1 = MatrixType::Random(rows, cols), m2(rows, cols), - mzero = MatrixType::Zero(rows, cols), identity = MatrixType::Identity(rows, rows); while(ei_abs(m1.determinant()) < RealScalar(0.1) && rows <= 8) diff --git a/test/eigen2/eigen2_linearstructure.cpp b/test/eigen2/eigen2_linearstructure.cpp index 22f02c396..488f4c485 100644 --- a/test/eigen2/eigen2_linearstructure.cpp +++ b/test/eigen2/eigen2_linearstructure.cpp @@ -25,8 +25,7 @@ template<typename MatrixType> void linearStructure(const MatrixType& m) // to test it, hence I consider that we will have tested Random.h MatrixType m1 = MatrixType::Random(rows, cols), m2 = MatrixType::Random(rows, cols), - m3(rows, cols), - mzero = MatrixType::Zero(rows, cols); + m3(rows, cols); Scalar s1 = ei_random<Scalar>(); while (ei_abs(s1)<1e-3) s1 = ei_random<Scalar>(); diff --git a/test/eigen2/eigen2_nomalloc.cpp b/test/eigen2/eigen2_nomalloc.cpp index e234abe4b..d34a69999 100644 --- a/test/eigen2/eigen2_nomalloc.cpp +++ b/test/eigen2/eigen2_nomalloc.cpp @@ -25,22 +25,12 @@ template<typename MatrixType> void nomalloc(const MatrixType& m) */ typedef typename MatrixType::Scalar Scalar; - typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; int rows = m.rows(); int cols = m.cols(); MatrixType m1 = MatrixType::Random(rows, cols), - m2 = MatrixType::Random(rows, cols), - m3(rows, cols), - mzero = MatrixType::Zero(rows, cols), - identity = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> - ::Identity(rows, rows), - square = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> - ::Random(rows, rows); - VectorType v1 = VectorType::Random(rows), - v2 = VectorType::Random(rows), - vzero = VectorType::Zero(rows); + m2 = MatrixType::Random(rows, cols); Scalar s1 = ei_random<Scalar>(); diff --git a/test/eigen2/eigen2_submatrices.cpp b/test/eigen2/eigen2_submatrices.cpp index c5d3f243d..dee970b63 100644 --- a/test/eigen2/eigen2_submatrices.cpp +++ b/test/eigen2/eigen2_submatrices.cpp @@ -51,16 +51,10 @@ template<typename MatrixType> void submatrices(const MatrixType& m) MatrixType m1 = MatrixType::Random(rows, cols), m2 = MatrixType::Random(rows, cols), m3(rows, cols), - mzero = MatrixType::Zero(rows, cols), ones = MatrixType::Ones(rows, cols), - identity = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> - ::Identity(rows, rows), square = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> ::Random(rows, rows); - VectorType v1 = VectorType::Random(rows), - v2 = VectorType::Random(rows), - v3 = VectorType::Random(rows), - vzero = VectorType::Zero(rows); + VectorType v1 = VectorType::Random(rows); Scalar s1 = ei_random<Scalar>(); diff --git a/test/eigen2/eigen2_triangular.cpp b/test/eigen2/eigen2_triangular.cpp index 3748c7d71..6f17b7dff 100644 --- a/test/eigen2/eigen2_triangular.cpp +++ b/test/eigen2/eigen2_triangular.cpp @@ -13,7 +13,6 @@ template<typename MatrixType> void triangular(const MatrixType& m) { typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; - typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; RealScalar largerEps = 10*test_precision<RealScalar>(); @@ -25,16 +24,7 @@ template<typename MatrixType> void triangular(const MatrixType& m) m3(rows, cols), m4(rows, cols), r1(rows, cols), - r2(rows, cols), - mzero = MatrixType::Zero(rows, cols), - mones = MatrixType::Ones(rows, cols), - identity = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> - ::Identity(rows, rows), - square = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> - ::Random(rows, rows); - VectorType v1 = VectorType::Random(rows), - v2 = VectorType::Random(rows), - vzero = VectorType::Zero(rows); + r2(rows, cols); MatrixType m1up = m1.template part<Eigen::UpperTriangular>(); MatrixType m2up = m2.template part<Eigen::UpperTriangular>(); diff --git a/test/eigen2/product.h b/test/eigen2/product.h index 2c9604d9a..ae1b4bae4 100644 --- a/test/eigen2/product.h +++ b/test/eigen2/product.h @@ -40,8 +40,7 @@ template<typename MatrixType> void product(const MatrixType& m) // to test it, hence I consider that we will have tested Random.h MatrixType m1 = MatrixType::Random(rows, cols), m2 = MatrixType::Random(rows, cols), - m3(rows, cols), - mzero = MatrixType::Zero(rows, cols); + m3(rows, cols); RowSquareMatrixType identity = RowSquareMatrixType::Identity(rows, rows), square = RowSquareMatrixType::Random(rows, rows), @@ -49,9 +48,7 @@ template<typename MatrixType> void product(const MatrixType& m) ColSquareMatrixType square2 = ColSquareMatrixType::Random(cols, cols), res2 = ColSquareMatrixType::Random(cols, cols); - RowVectorType v1 = RowVectorType::Random(rows), - v2 = RowVectorType::Random(rows), - vzero = RowVectorType::Zero(rows); + RowVectorType v1 = RowVectorType::Random(rows); ColVectorType vc2 = ColVectorType::Random(cols), vcres(cols); OtherMajorMatrixType tm1 = m1; diff --git a/test/eigen2support.cpp b/test/eigen2support.cpp index ad1d98091..1fa49a8c8 100644 --- a/test/eigen2support.cpp +++ b/test/eigen2support.cpp @@ -8,6 +8,7 @@ // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. #define EIGEN2_SUPPORT +#define EIGEN_NO_EIGEN2_DEPRECATED_WARNING #include "main.h" diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp index 5c6ecd875..38689cfbf 100644 --- a/test/eigensolver_selfadjoint.cpp +++ b/test/eigensolver_selfadjoint.cpp @@ -29,7 +29,21 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m) MatrixType a = MatrixType::Random(rows,cols); MatrixType a1 = MatrixType::Random(rows,cols); MatrixType symmA = a.adjoint() * a + a1.adjoint() * a1; + MatrixType symmC = symmA; + + // randomly nullify some rows/columns + { + Index count = 1;//internal::random<Index>(-cols,cols); + for(Index k=0; k<count; ++k) + { + Index i = internal::random<Index>(0,cols-1); + symmA.row(i).setZero(); + symmA.col(i).setZero(); + } + } + symmA.template triangularView<StrictlyUpper>().setZero(); + symmC.template triangularView<StrictlyUpper>().setZero(); MatrixType b = MatrixType::Random(rows,cols); MatrixType b1 = MatrixType::Random(rows,cols); @@ -40,7 +54,7 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m) SelfAdjointEigenSolver<MatrixType> eiDirect; eiDirect.computeDirect(symmA); // generalized eigen pb - GeneralizedSelfAdjointEigenSolver<MatrixType> eiSymmGen(symmA, symmB); + GeneralizedSelfAdjointEigenSolver<MatrixType> eiSymmGen(symmC, symmB); VERIFY_IS_EQUAL(eiSymm.info(), Success); VERIFY((symmA.template selfadjointView<Lower>() * eiSymm.eigenvectors()).isApprox( @@ -57,27 +71,28 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m) VERIFY_IS_APPROX(eiSymm.eigenvalues(), eiSymmNoEivecs.eigenvalues()); // generalized eigen problem Ax = lBx - eiSymmGen.compute(symmA, symmB,Ax_lBx); + eiSymmGen.compute(symmC, symmB,Ax_lBx); VERIFY_IS_EQUAL(eiSymmGen.info(), Success); - VERIFY((symmA.template selfadjointView<Lower>() * eiSymmGen.eigenvectors()).isApprox( + VERIFY((symmC.template selfadjointView<Lower>() * eiSymmGen.eigenvectors()).isApprox( symmB.template selfadjointView<Lower>() * (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps)); // generalized eigen problem BAx = lx - eiSymmGen.compute(symmA, symmB,BAx_lx); + eiSymmGen.compute(symmC, symmB,BAx_lx); VERIFY_IS_EQUAL(eiSymmGen.info(), Success); - VERIFY((symmB.template selfadjointView<Lower>() * (symmA.template selfadjointView<Lower>() * eiSymmGen.eigenvectors())).isApprox( + VERIFY((symmB.template selfadjointView<Lower>() * (symmC.template selfadjointView<Lower>() * eiSymmGen.eigenvectors())).isApprox( (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps)); // generalized eigen problem ABx = lx - eiSymmGen.compute(symmA, symmB,ABx_lx); + eiSymmGen.compute(symmC, symmB,ABx_lx); VERIFY_IS_EQUAL(eiSymmGen.info(), Success); - VERIFY((symmA.template selfadjointView<Lower>() * (symmB.template selfadjointView<Lower>() * eiSymmGen.eigenvectors())).isApprox( + VERIFY((symmC.template selfadjointView<Lower>() * (symmB.template selfadjointView<Lower>() * eiSymmGen.eigenvectors())).isApprox( (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps)); + eiSymm.compute(symmC); MatrixType sqrtSymmA = eiSymm.operatorSqrt(); - VERIFY_IS_APPROX(MatrixType(symmA.template selfadjointView<Lower>()), sqrtSymmA*sqrtSymmA); - VERIFY_IS_APPROX(sqrtSymmA, symmA.template selfadjointView<Lower>()*eiSymm.operatorInverseSqrt()); + VERIFY_IS_APPROX(MatrixType(symmC.template selfadjointView<Lower>()), sqrtSymmA*sqrtSymmA); + VERIFY_IS_APPROX(sqrtSymmA, symmC.template selfadjointView<Lower>()*eiSymm.operatorInverseSqrt()); MatrixType id = MatrixType::Identity(rows, cols); VERIFY_IS_APPROX(id.template selfadjointView<Lower>().operatorNorm(), RealScalar(1)); @@ -95,15 +110,15 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m) VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorInverseSqrt()); // test Tridiagonalization's methods - Tridiagonalization<MatrixType> tridiag(symmA); + Tridiagonalization<MatrixType> tridiag(symmC); // FIXME tridiag.matrixQ().adjoint() does not work - VERIFY_IS_APPROX(MatrixType(symmA.template selfadjointView<Lower>()), tridiag.matrixQ() * tridiag.matrixT().eval() * MatrixType(tridiag.matrixQ()).adjoint()); + VERIFY_IS_APPROX(MatrixType(symmC.template selfadjointView<Lower>()), tridiag.matrixQ() * tridiag.matrixT().eval() * MatrixType(tridiag.matrixQ()).adjoint()); if (rows > 1) { // Test matrix with NaN - symmA(0,0) = std::numeric_limits<typename MatrixType::RealScalar>::quiet_NaN(); - SelfAdjointEigenSolver<MatrixType> eiSymmNaN(symmA); + symmC(0,0) = std::numeric_limits<typename MatrixType::RealScalar>::quiet_NaN(); + SelfAdjointEigenSolver<MatrixType> eiSymmNaN(symmC); VERIFY_IS_EQUAL(eiSymmNaN.info(), NoConvergence); } } @@ -113,8 +128,10 @@ void test_eigensolver_selfadjoint() int s = 0; for(int i = 0; i < g_repeat; i++) { // very important to test 3x3 and 2x2 matrices since we provide special paths for them + CALL_SUBTEST_1( selfadjointeigensolver(Matrix2f()) ); CALL_SUBTEST_1( selfadjointeigensolver(Matrix2d()) ); CALL_SUBTEST_1( selfadjointeigensolver(Matrix3f()) ); + CALL_SUBTEST_1( selfadjointeigensolver(Matrix3d()) ); CALL_SUBTEST_2( selfadjointeigensolver(Matrix4d()) ); s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4); CALL_SUBTEST_3( selfadjointeigensolver(MatrixXf(s,s)) ); diff --git a/test/geo_hyperplane.cpp b/test/geo_hyperplane.cpp index f26fc1329..327537801 100644 --- a/test/geo_hyperplane.cpp +++ b/test/geo_hyperplane.cpp @@ -114,6 +114,32 @@ template<typename Scalar> void lines() } } +template<typename Scalar> void planes() +{ + using std::abs; + typedef Hyperplane<Scalar, 3> Plane; + typedef Matrix<Scalar,3,1> Vector; + + for(int i = 0; i < 10; i++) + { + Vector v0 = Vector::Random(); + Vector v1(v0), v2(v0); + if(internal::random<double>(0,1)>0.25) + v1 += Vector::Random(); + if(internal::random<double>(0,1)>0.25) + v2 += v1 * std::pow(internal::random<Scalar>(0,1),internal::random<int>(1,16)); + if(internal::random<double>(0,1)>0.25) + v2 += Vector::Random() * std::pow(internal::random<Scalar>(0,1),internal::random<int>(1,16)); + + Plane p0 = Plane::Through(v0, v1, v2); + + VERIFY_IS_APPROX(p0.normal().norm(), Scalar(1)); + VERIFY_IS_MUCH_SMALLER_THAN(p0.absDistance(v0), Scalar(1)); + VERIFY_IS_MUCH_SMALLER_THAN(p0.absDistance(v1), Scalar(1)); + VERIFY_IS_MUCH_SMALLER_THAN(p0.absDistance(v2), Scalar(1)); + } +} + template<typename Scalar> void hyperplane_alignment() { typedef Hyperplane<Scalar,3,AutoAlign> Plane3a; @@ -153,5 +179,7 @@ void test_geo_hyperplane() CALL_SUBTEST_4( hyperplane(Hyperplane<std::complex<double>,5>()) ); CALL_SUBTEST_1( lines<float>() ); CALL_SUBTEST_3( lines<double>() ); + CALL_SUBTEST_2( planes<float>() ); + CALL_SUBTEST_5( planes<double>() ); } } diff --git a/test/geo_transformations.cpp b/test/geo_transformations.cpp index ee3030b5d..547765714 100644 --- a/test/geo_transformations.cpp +++ b/test/geo_transformations.cpp @@ -98,11 +98,19 @@ template<typename Scalar, int Mode, int Options> void transformations() Matrix3 matrot1, m; Scalar a = internal::random<Scalar>(-Scalar(M_PI), Scalar(M_PI)); - Scalar s0 = internal::random<Scalar>(); + Scalar s0 = internal::random<Scalar>(), + s1 = internal::random<Scalar>(); + + while(v0.norm() < test_precision<Scalar>()) v0 = Vector3::Random(); + while(v1.norm() < test_precision<Scalar>()) v1 = Vector3::Random(); + VERIFY_IS_APPROX(v0, AngleAxisx(a, v0.normalized()) * v0); VERIFY_IS_APPROX(-v0, AngleAxisx(Scalar(M_PI), v0.unitOrthogonal()) * v0); - VERIFY_IS_APPROX(cos(a)*v0.squaredNorm(), v0.dot(AngleAxisx(a, v0.unitOrthogonal()) * v0)); + if(abs(cos(a)) > test_precision<Scalar>()) + { + VERIFY_IS_APPROX(cos(a)*v0.squaredNorm(), v0.dot(AngleAxisx(a, v0.unitOrthogonal()) * v0)); + } m = AngleAxisx(a, v0.normalized()).toRotationMatrix().adjoint(); VERIFY_IS_APPROX(Matrix3::Identity(), m * AngleAxisx(a, v0.normalized())); VERIFY_IS_APPROX(Matrix3::Identity(), AngleAxisx(a, v0.normalized()) * m); @@ -123,11 +131,18 @@ template<typename Scalar, int Mode, int Options> void transformations() // angle-axis conversion AngleAxisx aa = AngleAxisx(q1); VERIFY_IS_APPROX(q1 * v1, Quaternionx(aa) * v1); - VERIFY_IS_NOT_APPROX(q1 * v1, Quaternionx(AngleAxisx(aa.angle()*2,aa.axis())) * v1); + + if(abs(aa.angle()) > NumTraits<Scalar>::dummy_precision()) + { + VERIFY( !(q1 * v1).isApprox(Quaternionx(AngleAxisx(aa.angle()*2,aa.axis())) * v1) ); + } aa.fromRotationMatrix(aa.toRotationMatrix()); VERIFY_IS_APPROX(q1 * v1, Quaternionx(aa) * v1); - VERIFY_IS_NOT_APPROX(q1 * v1, Quaternionx(AngleAxisx(aa.angle()*2,aa.axis())) * v1); + if(abs(aa.angle()) > NumTraits<Scalar>::dummy_precision()) + { + VERIFY( !(q1 * v1).isApprox(Quaternionx(AngleAxisx(aa.angle()*2,aa.axis())) * v1) ); + } // AngleAxis VERIFY_IS_APPROX(AngleAxisx(a,v1.normalized()).toRotationMatrix(), @@ -347,7 +362,9 @@ template<typename Scalar, int Mode, int Options> void transformations() // test transform inversion t0.setIdentity(); t0.translate(v0); - t0.linear().setRandom(); + do { + t0.linear().setRandom(); + } while(t0.linear().jacobiSvd().singularValues()(2)<test_precision<Scalar>()); Matrix4 t044 = Matrix4::Zero(); t044(3,3) = 1; t044.block(0,0,t0.matrix().rows(),4) = t0.matrix(); @@ -397,6 +414,29 @@ template<typename Scalar, int Mode, int Options> void transformations() t20 = Translation2(v20) * (Rotation2D<Scalar>(s0) * Eigen::Scaling(s0)); t21 = Translation2(v20) * Rotation2D<Scalar>(s0) * Eigen::Scaling(s0); VERIFY_IS_APPROX(t20,t21); + + Rotation2D<Scalar> R0(s0), R1(s1); + t20 = Translation2(v20) * (R0 * Eigen::Scaling(s0)); + t21 = Translation2(v20) * R0 * Eigen::Scaling(s0); + VERIFY_IS_APPROX(t20,t21); + + t20 = Translation2(v20) * (R0 * R0.inverse() * Eigen::Scaling(s0)); + t21 = Translation2(v20) * Eigen::Scaling(s0); + VERIFY_IS_APPROX(t20,t21); + + VERIFY_IS_APPROX(s0, (R0.slerp(0, R1)).angle()); + VERIFY_IS_APPROX(s1, (R0.slerp(1, R1)).angle()); + VERIFY_IS_APPROX(s0, (R0.slerp(0.5, R0)).angle()); + VERIFY_IS_APPROX(Scalar(0), (R0.slerp(0.5, R0.inverse())).angle()); + + // check basic features + { + Rotation2D<Scalar> r1; // default ctor + r1 = Rotation2D<Scalar>(s0); // copy assignment + VERIFY_IS_APPROX(r1.angle(),s0); + Rotation2D<Scalar> r2(r1); // copy ctor + VERIFY_IS_APPROX(r2.angle(),s0); + } } template<typename Scalar> void transform_alignment() diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp index 7c21f0ab3..12c556e59 100644 --- a/test/jacobisvd.cpp +++ b/test/jacobisvd.cpp @@ -321,16 +321,23 @@ void jacobisvd_inf_nan() VERIFY(sub(some_inf, some_inf) != sub(some_inf, some_inf)); svd.compute(MatrixType::Constant(10,10,some_inf), ComputeFullU | ComputeFullV); - Scalar some_nan = zero<Scalar>() / zero<Scalar>(); - VERIFY(some_nan != some_nan); - svd.compute(MatrixType::Constant(10,10,some_nan), ComputeFullU | ComputeFullV); + Scalar nan = std::numeric_limits<Scalar>::quiet_NaN(); + VERIFY(nan != nan); + svd.compute(MatrixType::Constant(10,10,nan), ComputeFullU | ComputeFullV); MatrixType m = MatrixType::Zero(10,10); m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_inf; svd.compute(m, ComputeFullU | ComputeFullV); m = MatrixType::Zero(10,10); - m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_nan; + m(internal::random<int>(0,9), internal::random<int>(0,9)) = nan; + svd.compute(m, ComputeFullU | ComputeFullV); + + // regression test for bug 791 + m.resize(3,3); + m << 0, 2*NumTraits<Scalar>::epsilon(), 0.5, + 0, -0.5, 0, + nan, 0, 0; svd.compute(m, ComputeFullU | ComputeFullV); } @@ -434,6 +441,7 @@ void test_jacobisvd() // Test on inf/nan matrix CALL_SUBTEST_7( jacobisvd_inf_nan<MatrixXf>() ); + CALL_SUBTEST_10( jacobisvd_inf_nan<MatrixXd>() ); } CALL_SUBTEST_7(( jacobisvd<MatrixXf>(MatrixXf(internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2), internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2))) )); diff --git a/test/main.h b/test/main.h index 14f0d2f78..664204866 100644 --- a/test/main.h +++ b/test/main.h @@ -17,13 +17,36 @@ #include <sstream> #include <vector> #include <typeinfo> + +// The following includes of STL headers have to be done _before_ the +// definition of macros min() and max(). The reason is that many STL +// implementations will not work properly as the min and max symbols collide +// with the STL functions std:min() and std::max(). The STL headers may check +// for the macro definition of min/max and issue a warning or undefine the +// macros. +// +// Still, Windows defines min() and max() in windef.h as part of the regular +// Windows system interfaces and many other Windows APIs depend on these +// macros being available. To prevent the macro expansion of min/max and to +// make Eigen compatible with the Windows environment all function calls of +// std::min() and std::max() have to be written with parenthesis around the +// function name. +// +// All STL headers used by Eigen should be included here. Because main.h is +// included before any Eigen header and because the STL headers are guarded +// against multiple inclusions, no STL header will see our own min/max macro +// definitions. #include <limits> #include <algorithm> -#include <sstream> #include <complex> #include <deque> #include <queue> +#include <list> +// To test that all calls from Eigen code to std::min() and std::max() are +// protected by parenthesis against macro expansion, the min()/max() macros +// are defined here and any not-parenthesized min/max call will cause a +// compiler error. #define min(A,B) please_protect_your_min_with_parentheses #define max(A,B) please_protect_your_max_with_parentheses @@ -383,6 +406,26 @@ void randomPermutationVector(PermutationVectorType& v, typename PermutationVecto } } +template<typename T> bool isNotNaN(const T& x) +{ + return x==x; +} + +template<typename T> bool isNaN(const T& x) +{ + return x!=x; +} + +template<typename T> bool isInf(const T& x) +{ + return x > NumTraits<T>::highest(); +} + +template<typename T> bool isMinusInf(const T& x) +{ + return x < NumTraits<T>::lowest(); +} + } // end namespace Eigen template<typename T> struct GetDifferentType; diff --git a/test/nomalloc.cpp b/test/nomalloc.cpp index cbd02dd21..8e0402358 100644 --- a/test/nomalloc.cpp +++ b/test/nomalloc.cpp @@ -165,6 +165,38 @@ void ctms_decompositions() Eigen::JacobiSVD<Matrix> jSVD; jSVD.compute(A, ComputeFullU | ComputeFullV); } +void test_zerosized() { + // default constructors: + Eigen::MatrixXd A; + Eigen::VectorXd v; + // explicit zero-sized: + Eigen::ArrayXXd A0(0,0); + Eigen::ArrayXd v0(std::ptrdiff_t(0)); // FIXME ArrayXd(0) is ambiguous + + // assigning empty objects to each other: + A=A0; + v=v0; +} + +template<typename MatrixType> void test_reference(const MatrixType& m) { + typedef typename MatrixType::Scalar Scalar; + enum { Flag = MatrixType::IsRowMajor ? Eigen::RowMajor : Eigen::ColMajor}; + enum { TransposeFlag = !MatrixType::IsRowMajor ? Eigen::RowMajor : Eigen::ColMajor}; + typename MatrixType::Index rows = m.rows(), cols=m.cols(); + // Dynamic reference: + typedef Eigen::Ref<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Flag > > Ref; + typedef Eigen::Ref<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, TransposeFlag> > RefT; + + Ref r1(m); + Ref r2(m.block(rows/3, cols/4, rows/2, cols/2)); + RefT r3(m.transpose()); + RefT r4(m.topLeftCorner(rows/2, cols/2).transpose()); + + VERIFY_RAISES_ASSERT(RefT r5(m)); + VERIFY_RAISES_ASSERT(Ref r6(m.transpose())); + VERIFY_RAISES_ASSERT(Ref r7(Scalar(2) * m)); +} + void test_nomalloc() { // check that our operator new is indeed called: @@ -175,5 +207,6 @@ void test_nomalloc() // Check decomposition modules with dynamic matrices that have a known compile-time max size (ctms) CALL_SUBTEST_4(ctms_decompositions<float>()); - + CALL_SUBTEST_5(test_zerosized()); + CALL_SUBTEST_6(test_reference(Matrix<float,32,32>())); } diff --git a/test/nullary.cpp b/test/nullary.cpp index 5408d88b2..fbc721a1a 100644 --- a/test/nullary.cpp +++ b/test/nullary.cpp @@ -80,7 +80,9 @@ void testVectorType(const VectorType& base) Matrix<Scalar,1,Dynamic> col_vector(size); row_vector.setLinSpaced(size,low,high); col_vector.setLinSpaced(size,low,high); - VERIFY( row_vector.isApprox(col_vector.transpose(), NumTraits<Scalar>::epsilon())); + // when using the extended precision (e.g., FPU) the relative error might exceed 1 bit + // when computing the squared sum in isApprox, thus the 2x factor. + VERIFY( row_vector.isApprox(col_vector.transpose(), Scalar(2)*NumTraits<Scalar>::epsilon())); Matrix<Scalar,Dynamic,1> size_changer(size+50); size_changer.setLinSpaced(size,low,high); diff --git a/test/packetmath.cpp b/test/packetmath.cpp index 2c0519c41..38aa256ce 100644 --- a/test/packetmath.cpp +++ b/test/packetmath.cpp @@ -239,6 +239,12 @@ template<typename Scalar> void packetmath_real() data2[i] = internal::random<Scalar>(-87,88); } CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasExp, std::exp, internal::pexp); + { + data1[0] = std::numeric_limits<Scalar>::quiet_NaN(); + packet_helper<internal::packet_traits<Scalar>::HasExp,Packet> h; + h.store(data2, internal::pexp(h.load(data1))); + VERIFY(isNaN(data2[0])); + } for (int i=0; i<size; ++i) { @@ -247,8 +253,22 @@ template<typename Scalar> void packetmath_real() } if(internal::random<float>(0,1)<0.1) data1[internal::random<int>(0, PacketSize)] = 0; - CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasLog, std::log, internal::plog); CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasSqrt, std::sqrt, internal::psqrt); + CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasLog, std::log, internal::plog); + { + data1[0] = std::numeric_limits<Scalar>::quiet_NaN(); + packet_helper<internal::packet_traits<Scalar>::HasLog,Packet> h; + h.store(data2, internal::plog(h.load(data1))); + VERIFY(isNaN(data2[0])); + data1[0] = -1.0f; + h.store(data2, internal::plog(h.load(data1))); + VERIFY(isNaN(data2[0])); +#if !EIGEN_FAST_MATH + h.store(data2, internal::psqrt(h.load(data1))); + VERIFY(isNaN(data2[0])); + VERIFY(isNaN(data2[1])); +#endif + } } template<typename Scalar> void packetmath_notcomplex() diff --git a/test/product.h b/test/product.h index 856b234ac..0b3abe402 100644 --- a/test/product.h +++ b/test/product.h @@ -139,4 +139,12 @@ template<typename MatrixType> void product(const MatrixType& m) // inner product Scalar x = square2.row(c) * square2.col(c2); VERIFY_IS_APPROX(x, square2.row(c).transpose().cwiseProduct(square2.col(c2)).sum()); + + // outer product + VERIFY_IS_APPROX(m1.col(c) * m1.row(r), m1.block(0,c,rows,1) * m1.block(r,0,1,cols)); + VERIFY_IS_APPROX(m1.row(r).transpose() * m1.col(c).transpose(), m1.block(r,0,1,cols).transpose() * m1.block(0,c,rows,1).transpose()); + VERIFY_IS_APPROX(m1.block(0,c,rows,1) * m1.row(r), m1.block(0,c,rows,1) * m1.block(r,0,1,cols)); + VERIFY_IS_APPROX(m1.col(c) * m1.block(r,0,1,cols), m1.block(0,c,rows,1) * m1.block(r,0,1,cols)); + VERIFY_IS_APPROX(m1.leftCols(1) * m1.row(r), m1.block(0,0,rows,1) * m1.block(r,0,1,cols)); + VERIFY_IS_APPROX(m1.col(c) * m1.topRows(1), m1.block(0,c,rows,1) * m1.block(0,0,1,cols)); } diff --git a/test/ref.cpp b/test/ref.cpp index 19e81549c..32eb31048 100644 --- a/test/ref.cpp +++ b/test/ref.cpp @@ -183,15 +183,15 @@ void call_ref() VERIFY_EVALUATION_COUNT( call_ref_1(a,a), 0); VERIFY_EVALUATION_COUNT( call_ref_1(b,b.transpose()), 0); -// call_ref_1(ac); // does not compile because ac is const +// call_ref_1(ac,a<c); // does not compile because ac is const VERIFY_EVALUATION_COUNT( call_ref_1(ab,ab), 0); VERIFY_EVALUATION_COUNT( call_ref_1(a.head(4),a.head(4)), 0); VERIFY_EVALUATION_COUNT( call_ref_1(abc,abc), 0); VERIFY_EVALUATION_COUNT( call_ref_1(A.col(3),A.col(3)), 0); -// call_ref_1(A.row(3)); // does not compile because innerstride!=1 +// call_ref_1(A.row(3),A.row(3)); // does not compile because innerstride!=1 VERIFY_EVALUATION_COUNT( call_ref_3(A.row(3),A.row(3).transpose()), 0); VERIFY_EVALUATION_COUNT( call_ref_4(A.row(3),A.row(3).transpose()), 0); -// call_ref_1(a+a); // does not compile for obvious reason +// call_ref_1(a+a, a+a); // does not compile for obvious reason MatrixXf tmp = A*A.col(1); VERIFY_EVALUATION_COUNT( call_ref_2(A*A.col(1), tmp), 1); // evaluated into a temp @@ -212,7 +212,7 @@ void call_ref() VERIFY_EVALUATION_COUNT( call_ref_5(a,a), 0); VERIFY_EVALUATION_COUNT( call_ref_5(a.head(3),a.head(3)), 0); VERIFY_EVALUATION_COUNT( call_ref_5(A,A), 0); -// call_ref_5(A.transpose()); // does not compile +// call_ref_5(A.transpose(),A.transpose()); // does not compile because storage order does not match VERIFY_EVALUATION_COUNT( call_ref_5(A.block(1,1,2,2),A.block(1,1,2,2)), 0); VERIFY_EVALUATION_COUNT( call_ref_5(b,b), 0); // storage order do not match, but this is a degenerate case that should work VERIFY_EVALUATION_COUNT( call_ref_5(a.row(3),a.row(3)), 0); diff --git a/test/sparse_solver.h b/test/sparse_solver.h index d84aff070..244e81c1b 100644 --- a/test/sparse_solver.h +++ b/test/sparse_solver.h @@ -124,7 +124,23 @@ void check_sparse_determinant(Solver& solver, const typename Solver::MatrixType& Scalar refDet = dA.determinant(); VERIFY_IS_APPROX(refDet,solver.determinant()); } +template<typename Solver, typename DenseMat> +void check_sparse_abs_determinant(Solver& solver, const typename Solver::MatrixType& A, const DenseMat& dA) +{ + using std::abs; + typedef typename Solver::MatrixType Mat; + typedef typename Mat::Scalar Scalar; + + solver.compute(A); + if (solver.info() != Success) + { + std::cerr << "sparse solver testing: factorization failed (check_sparse_abs_determinant)\n"; + return; + } + Scalar refDet = abs(dA.determinant()); + VERIFY_IS_APPROX(refDet,solver.absDeterminant()); +} template<typename Solver, typename DenseMat> int generate_sparse_spd_problem(Solver& , typename Solver::MatrixType& A, typename Solver::MatrixType& halfA, DenseMat& dA, int maxSize = 300) @@ -324,3 +340,20 @@ template<typename Solver> void check_sparse_square_determinant(Solver& solver) check_sparse_determinant(solver, A, dA); } } + +template<typename Solver> void check_sparse_square_abs_determinant(Solver& solver) +{ + typedef typename Solver::MatrixType Mat; + typedef typename Mat::Scalar Scalar; + typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; + + // generate the problem + Mat A; + DenseMatrix dA; + generate_sparse_square_problem(solver, A, dA, 30); + A.makeCompressed(); + for (int i = 0; i < g_repeat; i++) { + check_sparse_abs_determinant(solver, A, dA); + } +} + diff --git a/test/sparselu.cpp b/test/sparselu.cpp index 37980defc..52371cb12 100644 --- a/test/sparselu.cpp +++ b/test/sparselu.cpp @@ -44,6 +44,9 @@ template<typename T> void test_sparselu_T() check_sparse_square_solving(sparselu_colamd); check_sparse_square_solving(sparselu_amd); check_sparse_square_solving(sparselu_natural); + + check_sparse_square_abs_determinant(sparselu_colamd); + check_sparse_square_abs_determinant(sparselu_amd); } void test_sparselu() diff --git a/test/sparseqr.cpp b/test/sparseqr.cpp index 1fe4a98ee..451c0e7f8 100644 --- a/test/sparseqr.cpp +++ b/test/sparseqr.cpp @@ -10,12 +10,11 @@ #include <Eigen/SparseQR> template<typename MatrixType,typename DenseMat> -int generate_sparse_rectangular_problem(MatrixType& A, DenseMat& dA, int maxRows = 300, int maxCols = 150) +int generate_sparse_rectangular_problem(MatrixType& A, DenseMat& dA, int maxRows = 300) { - eigen_assert(maxRows >= maxCols); typedef typename MatrixType::Scalar Scalar; int rows = internal::random<int>(1,maxRows); - int cols = internal::random<int>(1,maxCols); + int cols = internal::random<int>(1,rows); double density = (std::max)(8./(rows*cols), 0.01); A.resize(rows,cols); @@ -54,6 +53,8 @@ template<typename Scalar> void test_sparseqr_scalar() b = dA * DenseVector::Random(A.cols()); solver.compute(A); + if(internal::random<float>(0,1)>0.5) + solver.factorize(A); // this checks that calling analyzePattern is not needed if the pattern do not change. if (solver.info() != Success) { std::cerr << "sparse QR factorization failed\n"; diff --git a/test/stable_norm.cpp b/test/stable_norm.cpp index 549f91fbf..231dd9189 100644 --- a/test/stable_norm.cpp +++ b/test/stable_norm.cpp @@ -9,11 +9,6 @@ #include "main.h" -template<typename T> bool isNotNaN(const T& x) -{ - return x==x; -} - // workaround aggressive optimization in ICC template<typename T> EIGEN_DONT_INLINE T sub(T a, T b) { return a - b; } diff --git a/unsupported/Eigen/OpenGLSupport b/unsupported/Eigen/OpenGLSupport index c4090ab11..e2769449c 100644 --- a/unsupported/Eigen/OpenGLSupport +++ b/unsupported/Eigen/OpenGLSupport @@ -178,11 +178,11 @@ template<typename Scalar> void glLoadMatrix(const Transform<Scalar,3,Affine>& t) template<typename Scalar> void glLoadMatrix(const Transform<Scalar,3,Projective>& t) { glLoadMatrix(t.matrix()); } template<typename Scalar> void glLoadMatrix(const Transform<Scalar,3,AffineCompact>& t) { glLoadMatrix(Transform<Scalar,3,Affine>(t).matrix()); } -static void glRotate(const Rotation2D<float>& rot) +inline void glRotate(const Rotation2D<float>& rot) { glRotatef(rot.angle()*180.f/float(M_PI), 0.f, 0.f, 1.f); } -static void glRotate(const Rotation2D<double>& rot) +inline void glRotate(const Rotation2D<double>& rot) { glRotated(rot.angle()*180.0/M_PI, 0.0, 0.0, 1.0); } @@ -246,18 +246,18 @@ EIGEN_GL_FUNC1_SPECIALIZATION_MAT(glGet,GLenum,_,double, 4,4,Doublev) #ifdef GL_VERSION_2_0 -static void glUniform2fv_ei (GLint loc, const float* v) { glUniform2fv(loc,1,v); } -static void glUniform2iv_ei (GLint loc, const int* v) { glUniform2iv(loc,1,v); } +inline void glUniform2fv_ei (GLint loc, const float* v) { glUniform2fv(loc,1,v); } +inline void glUniform2iv_ei (GLint loc, const int* v) { glUniform2iv(loc,1,v); } -static void glUniform3fv_ei (GLint loc, const float* v) { glUniform3fv(loc,1,v); } -static void glUniform3iv_ei (GLint loc, const int* v) { glUniform3iv(loc,1,v); } +inline void glUniform3fv_ei (GLint loc, const float* v) { glUniform3fv(loc,1,v); } +inline void glUniform3iv_ei (GLint loc, const int* v) { glUniform3iv(loc,1,v); } -static void glUniform4fv_ei (GLint loc, const float* v) { glUniform4fv(loc,1,v); } -static void glUniform4iv_ei (GLint loc, const int* v) { glUniform4iv(loc,1,v); } +inline void glUniform4fv_ei (GLint loc, const float* v) { glUniform4fv(loc,1,v); } +inline void glUniform4iv_ei (GLint loc, const int* v) { glUniform4iv(loc,1,v); } -static void glUniformMatrix2fv_ei (GLint loc, const float* v) { glUniformMatrix2fv(loc,1,false,v); } -static void glUniformMatrix3fv_ei (GLint loc, const float* v) { glUniformMatrix3fv(loc,1,false,v); } -static void glUniformMatrix4fv_ei (GLint loc, const float* v) { glUniformMatrix4fv(loc,1,false,v); } +inline void glUniformMatrix2fv_ei (GLint loc, const float* v) { glUniformMatrix2fv(loc,1,false,v); } +inline void glUniformMatrix3fv_ei (GLint loc, const float* v) { glUniformMatrix3fv(loc,1,false,v); } +inline void glUniformMatrix4fv_ei (GLint loc, const float* v) { glUniformMatrix4fv(loc,1,false,v); } EIGEN_GL_FUNC1_DECLARATION (glUniform,GLint,const) @@ -294,9 +294,9 @@ EIGEN_GL_FUNC1_SPECIALIZATION_MAT(glUniform,GLint,const,float, 4,3,Matrix #ifdef GL_VERSION_3_0 -static void glUniform2uiv_ei (GLint loc, const unsigned int* v) { glUniform2uiv(loc,1,v); } -static void glUniform3uiv_ei (GLint loc, const unsigned int* v) { glUniform3uiv(loc,1,v); } -static void glUniform4uiv_ei (GLint loc, const unsigned int* v) { glUniform4uiv(loc,1,v); } +inline void glUniform2uiv_ei (GLint loc, const unsigned int* v) { glUniform2uiv(loc,1,v); } +inline void glUniform3uiv_ei (GLint loc, const unsigned int* v) { glUniform3uiv(loc,1,v); } +inline void glUniform4uiv_ei (GLint loc, const unsigned int* v) { glUniform4uiv(loc,1,v); } EIGEN_GL_FUNC1_SPECIALIZATION_VEC(glUniform,GLint,const,unsigned int, 2,2uiv_ei) EIGEN_GL_FUNC1_SPECIALIZATION_VEC(glUniform,GLint,const,unsigned int, 3,3uiv_ei) @@ -305,9 +305,9 @@ EIGEN_GL_FUNC1_SPECIALIZATION_VEC(glUniform,GLint,const,unsigned int, 4,4uiv_ei) #endif #ifdef GL_ARB_gpu_shader_fp64 -static void glUniform2dv_ei (GLint loc, const double* v) { glUniform2dv(loc,1,v); } -static void glUniform3dv_ei (GLint loc, const double* v) { glUniform3dv(loc,1,v); } -static void glUniform4dv_ei (GLint loc, const double* v) { glUniform4dv(loc,1,v); } +inline void glUniform2dv_ei (GLint loc, const double* v) { glUniform2dv(loc,1,v); } +inline void glUniform3dv_ei (GLint loc, const double* v) { glUniform3dv(loc,1,v); } +inline void glUniform4dv_ei (GLint loc, const double* v) { glUniform4dv(loc,1,v); } EIGEN_GL_FUNC1_SPECIALIZATION_VEC(glUniform,GLint,const,double, 2,2dv_ei) EIGEN_GL_FUNC1_SPECIALIZATION_VEC(glUniform,GLint,const,double, 3,3dv_ei) diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index c32437281..78a307e96 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -110,7 +110,6 @@ void MatrixPowerAtomic<MatrixType>::compute2x2(MatrixType& res, RealScalar p) co using std::abs; using std::pow; - ArrayType logTdiag = m_A.diagonal().array().log(); res.coeffRef(0,0) = pow(m_A.coeff(0,0), p); for (Index i=1; i < m_A.cols(); ++i) { diff --git a/unsupported/doc/examples/CMakeLists.txt b/unsupported/doc/examples/CMakeLists.txt index 978f9afd8..c47646dfc 100644 --- a/unsupported/doc/examples/CMakeLists.txt +++ b/unsupported/doc/examples/CMakeLists.txt @@ -10,12 +10,10 @@ FOREACH(example_src ${examples_SRCS}) if(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO) target_link_libraries(example_${example} ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO}) endif() - GET_TARGET_PROPERTY(example_executable - example_${example} LOCATION) ADD_CUSTOM_COMMAND( TARGET example_${example} POST_BUILD - COMMAND ${example_executable} + COMMAND example_${example} ARGS >${CMAKE_CURRENT_BINARY_DIR}/${example}.out ) ADD_DEPENDENCIES(unsupported_examples example_${example}) diff --git a/unsupported/doc/snippets/CMakeLists.txt b/unsupported/doc/snippets/CMakeLists.txt index 4a4157933..f0c5cc2a8 100644 --- a/unsupported/doc/snippets/CMakeLists.txt +++ b/unsupported/doc/snippets/CMakeLists.txt @@ -14,12 +14,10 @@ FOREACH(snippet_src ${snippets_SRCS}) if(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO) target_link_libraries(${compile_snippet_target} ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO}) endif() - GET_TARGET_PROPERTY(compile_snippet_executable - ${compile_snippet_target} LOCATION) ADD_CUSTOM_COMMAND( TARGET ${compile_snippet_target} POST_BUILD - COMMAND ${compile_snippet_executable} + COMMAND ${compile_snippet_target} ARGS >${CMAKE_CURRENT_BINARY_DIR}/${snippet}.out ) ADD_DEPENDENCIES(unsupported_snippets ${compile_snippet_target}) diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt index a94a3b5e5..2e4cfdb2e 100644 --- a/unsupported/test/CMakeLists.txt +++ b/unsupported/test/CMakeLists.txt @@ -29,11 +29,7 @@ ei_add_test(NonLinearOptimization) ei_add_test(NumericalDiff) ei_add_test(autodiff) - -if (NOT CMAKE_CXX_COMPILER MATCHES "clang\\+\\+$") ei_add_test(BVH) -endif() - ei_add_test(matrix_exponential) ei_add_test(matrix_function) ei_add_test(matrix_power) @@ -73,8 +69,9 @@ if(NOT EIGEN_TEST_NO_OPENGL) find_package(GLUT) find_package(GLEW) if(OPENGL_FOUND AND GLUT_FOUND AND GLEW_FOUND) + include_directories(${OPENGL_INCLUDE_DIR} ${GLUT_INCLUDE_DIR} ${GLEW_INCLUDE_DIRS}) ei_add_property(EIGEN_TESTED_BACKENDS "OpenGL, ") - set(EIGEN_GL_LIB ${GLUT_LIBRARIES} ${GLEW_LIBRARIES}) + set(EIGEN_GL_LIB ${GLUT_LIBRARIES} ${GLEW_LIBRARIES} ${OPENGL_LIBRARIES}) ei_add_test(openglsupport "" "${EIGEN_GL_LIB}" ) else() ei_add_property(EIGEN_MISSING_BACKENDS "OpenGL, ") diff --git a/unsupported/test/mpreal/mpreal.h b/unsupported/test/mpreal/mpreal.h index ef0a6a9f0..dddda7dd9 100644 --- a/unsupported/test/mpreal/mpreal.h +++ b/unsupported/test/mpreal/mpreal.h @@ -1,59 +1,47 @@ /*
- Multi-precision floating point number class for C++.
+ MPFR C++: Multi-precision floating point number class for C++.
Based on MPFR library: http://mpfr.org
Project homepage: http://www.holoborodko.com/pavel/mpfr
Contact e-mail: pavel@holoborodko.com
- Copyright (c) 2008-2012 Pavel Holoborodko
+ Copyright (c) 2008-2014 Pavel Holoborodko
Contributors:
Dmitriy Gubanov, Konstantin Holoborodko, Brian Gladman,
Helmut Jarausch, Fokko Beekhof, Ulrich Mutze, Heinz van Saanen,
Pere Constans, Peter van Hoof, Gael Guennebaud, Tsai Chia Cheng,
- Alexei Zubanov, Jauhien Piatlicki, Victor Berger, John Westwood.
+ Alexei Zubanov, Jauhien Piatlicki, Victor Berger, John Westwood,
+ Petr Aleksandrov, Orion Poplawski, Charles Karney.
- ****************************************************************************
- This library is free software; you can redistribute it and/or
- modify it under the terms of the GNU Lesser General Public
- License as published by the Free Software Foundation; either
- version 2.1 of the License, or (at your option) any later version.
-
- This library is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- Lesser General Public License for more details.
+ Licensing:
+ (A) MPFR C++ is under GNU General Public License ("GPL").
+
+ (B) Non-free licenses may also be purchased from the author, for users who
+ do not want their programs protected by the GPL.
- You should have received a copy of the GNU Lesser General Public
- License along with this library; if not, write to the Free Software
- Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ The non-free licenses are for users that wish to use MPFR C++ in
+ their products but are unwilling to release their software
+ under the GPL (which would require them to release source code
+ and allow free redistribution).
- ****************************************************************************
- Redistribution and use in source and binary forms, with or without
- modification, are permitted provided that the following conditions
- are met:
-
- 1. Redistributions of source code must retain the above copyright
- notice, this list of conditions and the following disclaimer.
-
- 2. Redistributions in binary form must reproduce the above copyright
- notice, this list of conditions and the following disclaimer in the
- documentation and/or other materials provided with the distribution.
+ Such users can purchase an unlimited-use license from the author.
+ Contact us for more details.
- 3. The name of the author may not be used to endorse or promote products
- derived from this software without specific prior written permission.
-
- THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
- ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
- FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- SUCH DAMAGE.
+ GNU General Public License ("GPL") copyright permissions statement:
+ **************************************************************************
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef __MPREAL_H__
@@ -65,11 +53,21 @@ #include <stdexcept>
#include <cfloat>
#include <cmath>
+#include <cstring>
#include <limits>
// Options
#define MPREAL_HAVE_INT64_SUPPORT // Enable int64_t support if possible. Available only for MSVC 2010 & GCC.
#define MPREAL_HAVE_MSVC_DEBUGVIEW // Enable Debugger Visualizer for "Debug" builds in MSVC.
+#define MPREAL_HAVE_DYNAMIC_STD_NUMERIC_LIMITS // Enable extended std::numeric_limits<mpfr::mpreal> specialization.
+ // Meaning that "digits", "round_style" and similar members are defined as functions, not constants.
+ // See std::numeric_limits<mpfr::mpreal> at the end of the file for more information.
+
+// Library version
+#define MPREAL_VERSION_MAJOR 3
+#define MPREAL_VERSION_MINOR 5
+#define MPREAL_VERSION_PATCHLEVEL 9
+#define MPREAL_VERSION_STRING "3.5.9"
// Detect compiler using signatures from http://predef.sourceforge.net/
#if defined(__GNUC__) && defined(__INTEL_COMPILER)
@@ -82,6 +80,32 @@ #define IsInf(x) std::isinf(x) // GNU C/C++ (and/or other compilers), just hope for C99 conformance
#endif
+// A Clang feature extension to determine compiler features.
+#ifndef __has_feature
+ #define __has_feature(x) 0
+#endif
+
+// Detect support for r-value references (move semantic). Borrowed from Eigen.
+#if (__has_feature(cxx_rvalue_references) || \
+ defined(__GXX_EXPERIMENTAL_CXX0X__) || __cplusplus >= 201103L || \
+ (defined(_MSC_VER) && _MSC_VER >= 1600))
+
+ #define MPREAL_HAVE_MOVE_SUPPORT
+
+ // Use fields in mpfr_t structure to check if it was initialized / set dummy initialization
+ #define mpfr_is_initialized(x) (0 != (x)->_mpfr_d)
+ #define mpfr_set_uninitialized(x) ((x)->_mpfr_d = 0 )
+#endif
+
+// Detect support for explicit converters.
+#if (__has_feature(cxx_explicit_conversions) || \
+ defined(__GXX_EXPERIMENTAL_CXX0X__) || __cplusplus >= 201103L || \
+ (defined(_MSC_VER) && _MSC_VER >= 1800))
+
+ #define MPREAL_HAVE_EXPLICIT_CONVERTERS
+#endif
+
+// Detect available 64-bit capabilities
#if defined(MPREAL_HAVE_INT64_SUPPORT)
#define MPFR_USE_INTMAX_T // Should be defined before mpfr.h
@@ -97,7 +121,7 @@ #endif
#elif defined (__GNUC__) && defined(__linux__)
- #if defined(__amd64__) || defined(__amd64) || defined(__x86_64__) || defined(__x86_64) || defined(__ia64) || defined(__itanium__) || defined(_M_IA64)
+ #if defined(__amd64__) || defined(__amd64) || defined(__x86_64__) || defined(__x86_64) || defined(__ia64) || defined(__itanium__) || defined(_M_IA64) || defined (__PPC64__)
#undef MPREAL_HAVE_INT64_SUPPORT // Remove all shaman dances for x64 builds since
#undef MPFR_USE_INTMAX_T // GCC already supports x64 as of "long int" is 64-bit integer, nothing left to do
#else
@@ -111,7 +135,7 @@ #endif
#if defined(MPREAL_HAVE_MSVC_DEBUGVIEW) && defined(_MSC_VER) && defined(_DEBUG)
-#define MPREAL_MSVC_DEBUGVIEW_CODE DebugView = toString();
+ #define MPREAL_MSVC_DEBUGVIEW_CODE DebugView = toString();
#define MPREAL_MSVC_DEBUGVIEW_DATA std::string DebugView;
#else
#define MPREAL_MSVC_DEBUGVIEW_CODE
@@ -149,7 +173,6 @@ public: // Constructors && type conversions
mpreal();
mpreal(const mpreal& u);
- mpreal(const mpfr_t u);
mpreal(const mpf_t u);
mpreal(const mpz_t u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd());
mpreal(const mpq_t u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd());
@@ -159,6 +182,10 @@ public: mpreal(const unsigned int u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd());
mpreal(const long int u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd());
mpreal(const int u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd());
+
+ // Construct mpreal from mpfr_t structure.
+ // shared = true allows to avoid deep copy, so that mpreal and 'u' share the same data & pointers.
+ mpreal(const mpfr_t u, bool shared = false);
#if defined (MPREAL_HAVE_INT64_SUPPORT)
mpreal(const uint64_t u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd());
@@ -170,6 +197,11 @@ public: ~mpreal();
+#ifdef MPREAL_HAVE_MOVE_SUPPORT
+ mpreal& operator=(mpreal&& v);
+ mpreal(mpreal&& u);
+#endif
+
// Operations
// =
// +, -, *, /, ++, --, <<, >>
@@ -292,14 +324,34 @@ public: friend bool operator == (const mpreal& a, const double b);
// Type Conversion operators
+ bool toBool (mp_rnd_t mode = GMP_RNDZ) const;
long toLong (mp_rnd_t mode = GMP_RNDZ) const;
unsigned long toULong (mp_rnd_t mode = GMP_RNDZ) const;
+ float toFloat (mp_rnd_t mode = GMP_RNDN) const;
double toDouble (mp_rnd_t mode = GMP_RNDN) const;
long double toLDouble (mp_rnd_t mode = GMP_RNDN) const;
+#if defined (MPREAL_HAVE_EXPLICIT_CONVERTERS)
+ explicit operator bool () const { return toBool(); }
+ explicit operator int () const { return toLong(); }
+ explicit operator long () const { return toLong(); }
+ explicit operator long long () const { return toLong(); }
+ explicit operator unsigned () const { return toULong(); }
+ explicit operator unsigned long () const { return toULong(); }
+ explicit operator unsigned long long () const { return toULong(); }
+ explicit operator float () const { return toFloat(); }
+ explicit operator double () const { return toDouble(); }
+ explicit operator long double () const { return toLDouble(); }
+#endif
+
#if defined (MPREAL_HAVE_INT64_SUPPORT)
int64_t toInt64 (mp_rnd_t mode = GMP_RNDZ) const;
uint64_t toUInt64 (mp_rnd_t mode = GMP_RNDZ) const;
+
+ #if defined (MPREAL_HAVE_EXPLICIT_CONVERTERS)
+ explicit operator int64_t () const { return toInt64(); }
+ explicit operator uint64_t () const { return toUInt64(); }
+ #endif
#endif
// Get raw pointers so that mpreal can be directly used in raw mpfr_* functions
@@ -308,121 +360,125 @@ public: ::mpfr_srcptr mpfr_srcptr() const;
// Convert mpreal to string with n significant digits in base b
- // n = 0 -> convert with the maximum available digits
- std::string toString(int n = 0, int b = 10, mp_rnd_t mode = mpreal::get_default_rnd()) const;
+ // n = -1 -> convert with the maximum available digits
+ std::string toString(int n = -1, int b = 10, mp_rnd_t mode = mpreal::get_default_rnd()) const;
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
- std::string toString(const std::string& format) const;
+ std::string toString(const std::string& format) const;
#endif
+ std::ostream& output(std::ostream& os) const;
+
// Math Functions
- friend const mpreal sqr (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal sqrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal sqrt(const unsigned long int v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal cbrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal root(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal pow (const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal pow (const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal pow (const mpreal& a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal pow (const mpreal& a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal pow (const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal pow (const unsigned long int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal fabs(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
-
- friend const mpreal abs(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+ friend const mpreal sqr (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal sqrt(const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal sqrt(const unsigned long int v, mp_rnd_t rnd_mode);
+ friend const mpreal cbrt(const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal root(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode);
+ friend const mpreal pow (const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode);
+ friend const mpreal pow (const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode);
+ friend const mpreal pow (const mpreal& a, const unsigned long int b, mp_rnd_t rnd_mode);
+ friend const mpreal pow (const mpreal& a, const long int b, mp_rnd_t rnd_mode);
+ friend const mpreal pow (const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode);
+ friend const mpreal pow (const unsigned long int a, const unsigned long int b, mp_rnd_t rnd_mode);
+ friend const mpreal fabs(const mpreal& v, mp_rnd_t rnd_mode);
+
+ friend const mpreal abs(const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode);
+ friend inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode);
+ friend inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode);
+ friend inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode);
+ friend inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode);
friend int cmpabs(const mpreal& a,const mpreal& b);
- friend const mpreal log (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal log2 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal log10(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal exp (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal exp2 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal exp10(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal log1p(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal expm1(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
-
- friend const mpreal cos(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal sin(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal tan(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal sec(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal csc(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal cot(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
-
- friend const mpreal acos (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal asin (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal atan (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal acot (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal asec (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal acsc (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
-
- friend const mpreal cosh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal sinh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal tanh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal sech (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal csch (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal coth (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal acosh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal asinh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal atanh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal acoth (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal asech (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal acsch (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
-
- friend const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
-
- friend const mpreal fac_ui (unsigned long int v, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal eint (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
-
- friend const mpreal gamma (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal lngamma (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal lgamma (const mpreal& v, int *signp = 0, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal zeta (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal erf (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal erfc (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal besselj0 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal besselj1 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal besseljn (long n, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal bessely0 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal bessely1 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal besselyn (long n, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+ friend const mpreal log (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal log2 (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal log10(const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal exp (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal exp2 (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal exp10(const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal log1p(const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal expm1(const mpreal& v, mp_rnd_t rnd_mode);
+
+ friend const mpreal cos(const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal sin(const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal tan(const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal sec(const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal csc(const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal cot(const mpreal& v, mp_rnd_t rnd_mode);
+ friend int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode);
+
+ friend const mpreal acos (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal asin (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal atan (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode);
+ friend const mpreal acot (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal asec (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal acsc (const mpreal& v, mp_rnd_t rnd_mode);
+
+ friend const mpreal cosh (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal sinh (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal tanh (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal sech (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal csch (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal coth (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal acosh (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal asinh (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal atanh (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal acoth (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal asech (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal acsch (const mpreal& v, mp_rnd_t rnd_mode);
+
+ friend const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode);
+
+ friend const mpreal fac_ui (unsigned long int v, mp_prec_t prec, mp_rnd_t rnd_mode);
+ friend const mpreal eint (const mpreal& v, mp_rnd_t rnd_mode);
+
+ friend const mpreal gamma (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal lngamma (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal lgamma (const mpreal& v, int *signp, mp_rnd_t rnd_mode);
+ friend const mpreal zeta (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal erf (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal erfc (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal besselj0 (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal besselj1 (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal besseljn (long n, const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal bessely0 (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal bessely1 (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal besselyn (long n, const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode);
+ friend const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode);
+ friend const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode);
+ friend const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode);
friend int sgn(const mpreal& v); // returns -1 or +1
// MPFR 2.4.0 Specifics
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
- friend int sinh_cosh (mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal li2 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal rec_sqrt (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+ friend int sinh_cosh (mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal li2 (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode);
+ friend const mpreal rec_sqrt (const mpreal& v, mp_rnd_t rnd_mode);
// MATLAB's semantic equivalents
- friend const mpreal rem (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); // Remainder after division
- friend const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); // Modulus after division
+ friend const mpreal rem (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); // Remainder after division
+ friend const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); // Modulus after division
#endif
// MPFR 3.0.0 Specifics
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
- friend const mpreal digamma (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal ai (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal urandom (gmp_randstate_t& state, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); // use gmp_randinit_default() to init state, gmp_randclear() to clear
+ friend const mpreal digamma (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal ai (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal urandom (gmp_randstate_t& state, mp_rnd_t rnd_mode); // use gmp_randinit_default() to init state, gmp_randclear() to clear
+ friend const mpreal grandom (gmp_randstate_t& state, mp_rnd_t rnd_mode); // use gmp_randinit_default() to init state, gmp_randclear() to clear
+ friend const mpreal grandom (unsigned int seed);
#endif
// Uniformly distributed random number generation in [0,1] using
// Mersenne-Twister algorithm by default.
// Use parameter to setup seed, e.g.: random((unsigned)time(NULL))
// Check urandom() for more precise control.
- friend const mpreal random(unsigned int seed = 0);
-
+ friend const mpreal random(unsigned int seed);
+
// Exponent and mantissa manipulation
friend const mpreal frexp(const mpreal& v, mp_exp_t* exp);
friend const mpreal ldexp(const mpreal& v, mp_exp_t exp);
@@ -433,31 +489,31 @@ public: // Constants
// don't forget to call mpfr_free_cache() for every thread where you are using const-functions
- friend const mpreal const_log2 (mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal const_pi (mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal const_euler (mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal const_catalan (mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+ friend const mpreal const_log2 (mp_prec_t prec, mp_rnd_t rnd_mode);
+ friend const mpreal const_pi (mp_prec_t prec, mp_rnd_t rnd_mode);
+ friend const mpreal const_euler (mp_prec_t prec, mp_rnd_t rnd_mode);
+ friend const mpreal const_catalan (mp_prec_t prec, mp_rnd_t rnd_mode);
// returns +inf iff sign>=0 otherwise -inf
- friend const mpreal const_infinity(int sign = 1, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+ friend const mpreal const_infinity(int sign, mp_prec_t prec);
// Output/ Input
friend std::ostream& operator<<(std::ostream& os, const mpreal& v);
friend std::istream& operator>>(std::istream& is, mpreal& v);
// Integer Related Functions
- friend const mpreal rint (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+ friend const mpreal rint (const mpreal& v, mp_rnd_t rnd_mode);
friend const mpreal ceil (const mpreal& v);
friend const mpreal floor(const mpreal& v);
friend const mpreal round(const mpreal& v);
friend const mpreal trunc(const mpreal& v);
- friend const mpreal rint_ceil (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal rint_floor (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal rint_round (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal rint_trunc (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal frac (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal remainder ( const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+ friend const mpreal rint_ceil (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal rint_floor (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal rint_round (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal rint_trunc (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal frac (const mpreal& v, mp_rnd_t rnd_mode);
+ friend const mpreal remainder ( const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode);
+ friend const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode);
// Miscellaneous Functions
friend const mpreal nexttoward (const mpreal& x, const mpreal& y);
@@ -524,19 +580,22 @@ public: // Efficient swapping of two mpreal values - needed for std algorithms
friend void swap(mpreal& x, mpreal& y);
- friend const mpreal fmax(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
- friend const mpreal fmin(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+ friend const mpreal fmax(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode);
+ friend const mpreal fmin(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode);
private:
// Human friendly Debug Preview in Visual Studio.
// Put one of these lines:
//
- // mpfr::mpreal=<DebugView> ; Show value only
+ // mpfr::mpreal=<DebugView> ; Show value only
// mpfr::mpreal=<DebugView>, <mp[0]._mpfr_prec,u>bits ; Show value & precision
//
// at the beginning of
// [Visual Studio Installation Folder]\Common7\Packages\Debugger\autoexp.dat
MPREAL_MSVC_DEBUGVIEW_DATA
+
+ // "Smart" resources deallocation. Checks if instance initialized before deletion.
+ void clear(::mpfr_ptr);
};
//////////////////////////////////////////////////////////////////////////
@@ -551,64 +610,89 @@ public: // Default constructor: creates mp number and initializes it to 0.
inline mpreal::mpreal()
{
- mpfr_init2(mp,mpreal::get_default_prec());
- mpfr_set_ui(mp,0,mpreal::get_default_rnd());
+ mpfr_init2 (mpfr_ptr(), mpreal::get_default_prec());
+ mpfr_set_ui(mpfr_ptr(), 0, mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mpreal::mpreal(const mpreal& u)
{
- mpfr_init2(mp,mpfr_get_prec(u.mp));
- mpfr_set(mp,u.mp,mpreal::get_default_rnd());
+ mpfr_init2(mpfr_ptr(),mpfr_get_prec(u.mpfr_srcptr()));
+ mpfr_set (mpfr_ptr(),u.mpfr_srcptr(),mpreal::get_default_rnd());
+
+ MPREAL_MSVC_DEBUGVIEW_CODE;
+}
+
+#ifdef MPREAL_HAVE_MOVE_SUPPORT
+inline mpreal::mpreal(mpreal&& other)
+{
+ mpfr_set_uninitialized(mpfr_ptr()); // make sure "other" holds no pinter to actual data
+ mpfr_swap(mpfr_ptr(), other.mpfr_ptr());
+
+ MPREAL_MSVC_DEBUGVIEW_CODE;
+}
+
+inline mpreal& mpreal::operator=(mpreal&& other)
+{
+ mpfr_swap(mpfr_ptr(), other.mpfr_ptr());
MPREAL_MSVC_DEBUGVIEW_CODE;
+ return *this;
}
+#endif
-inline mpreal::mpreal(const mpfr_t u)
+inline mpreal::mpreal(const mpfr_t u, bool shared)
{
- mpfr_init2(mp,mpfr_get_prec(u));
- mpfr_set(mp,u,mpreal::get_default_rnd());
+ if(shared)
+ {
+ std::memcpy(mpfr_ptr(), u, sizeof(mpfr_t));
+ }
+ else
+ {
+ mpfr_init2(mpfr_ptr(), mpfr_get_prec(u));
+ mpfr_set (mpfr_ptr(), u, mpreal::get_default_rnd());
+ }
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mpreal::mpreal(const mpf_t u)
{
- mpfr_init2(mp,(mp_prec_t) mpf_get_prec(u)); // (gmp: mp_bitcnt_t) unsigned long -> long (mpfr: mp_prec_t)
- mpfr_set_f(mp,u,mpreal::get_default_rnd());
+ mpfr_init2(mpfr_ptr(),(mp_prec_t) mpf_get_prec(u)); // (gmp: mp_bitcnt_t) unsigned long -> long (mpfr: mp_prec_t)
+ mpfr_set_f(mpfr_ptr(),u,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mpreal::mpreal(const mpz_t u, mp_prec_t prec, mp_rnd_t mode)
{
- mpfr_init2(mp,prec);
- mpfr_set_z(mp,u,mode);
+ mpfr_init2(mpfr_ptr(), prec);
+ mpfr_set_z(mpfr_ptr(), u, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mpreal::mpreal(const mpq_t u, mp_prec_t prec, mp_rnd_t mode)
{
- mpfr_init2(mp,prec);
- mpfr_set_q(mp,u,mode);
+ mpfr_init2(mpfr_ptr(), prec);
+ mpfr_set_q(mpfr_ptr(), u, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode)
{
- mpfr_init2(mp, prec);
+ mpfr_init2(mpfr_ptr(), prec);
#if (MPREAL_DOUBLE_BITS_OVERFLOW > -1)
- if(fits_in_bits(u, MPREAL_DOUBLE_BITS_OVERFLOW))
- {
- mpfr_set_d(mp, u, mode);
- }else
- throw conversion_overflow();
+ if(fits_in_bits(u, MPREAL_DOUBLE_BITS_OVERFLOW))
+ {
+ mpfr_set_d(mpfr_ptr(), u, mode);
+ }else
+ throw conversion_overflow();
#else
- mpfr_set_d(mp, u, mode);
+ mpfr_set_d(mpfr_ptr(), u, mode);
#endif
MPREAL_MSVC_DEBUGVIEW_CODE;
@@ -616,40 +700,40 @@ inline mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode) inline mpreal::mpreal(const long double u, mp_prec_t prec, mp_rnd_t mode)
{
- mpfr_init2(mp,prec);
- mpfr_set_ld(mp,u,mode);
+ mpfr_init2 (mpfr_ptr(), prec);
+ mpfr_set_ld(mpfr_ptr(), u, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mpreal::mpreal(const unsigned long int u, mp_prec_t prec, mp_rnd_t mode)
{
- mpfr_init2(mp,prec);
- mpfr_set_ui(mp,u,mode);
+ mpfr_init2 (mpfr_ptr(), prec);
+ mpfr_set_ui(mpfr_ptr(), u, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mpreal::mpreal(const unsigned int u, mp_prec_t prec, mp_rnd_t mode)
{
- mpfr_init2(mp,prec);
- mpfr_set_ui(mp,u,mode);
+ mpfr_init2 (mpfr_ptr(), prec);
+ mpfr_set_ui(mpfr_ptr(), u, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mpreal::mpreal(const long int u, mp_prec_t prec, mp_rnd_t mode)
{
- mpfr_init2(mp,prec);
- mpfr_set_si(mp,u,mode);
+ mpfr_init2 (mpfr_ptr(), prec);
+ mpfr_set_si(mpfr_ptr(), u, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mpreal::mpreal(const int u, mp_prec_t prec, mp_rnd_t mode)
{
- mpfr_init2(mp,prec);
- mpfr_set_si(mp,u,mode);
+ mpfr_init2 (mpfr_ptr(), prec);
+ mpfr_set_si(mpfr_ptr(), u, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
@@ -657,16 +741,16 @@ inline mpreal::mpreal(const int u, mp_prec_t prec, mp_rnd_t mode) #if defined (MPREAL_HAVE_INT64_SUPPORT)
inline mpreal::mpreal(const uint64_t u, mp_prec_t prec, mp_rnd_t mode)
{
- mpfr_init2(mp,prec);
- mpfr_set_uj(mp, u, mode);
+ mpfr_init2 (mpfr_ptr(), prec);
+ mpfr_set_uj(mpfr_ptr(), u, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mpreal::mpreal(const int64_t u, mp_prec_t prec, mp_rnd_t mode)
{
- mpfr_init2(mp,prec);
- mpfr_set_sj(mp, u, mode);
+ mpfr_init2 (mpfr_ptr(), prec);
+ mpfr_set_sj(mpfr_ptr(), u, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
@@ -674,23 +758,31 @@ inline mpreal::mpreal(const int64_t u, mp_prec_t prec, mp_rnd_t mode) inline mpreal::mpreal(const char* s, mp_prec_t prec, int base, mp_rnd_t mode)
{
- mpfr_init2(mp, prec);
- mpfr_set_str(mp, s, base, mode);
+ mpfr_init2 (mpfr_ptr(), prec);
+ mpfr_set_str(mpfr_ptr(), s, base, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mpreal::mpreal(const std::string& s, mp_prec_t prec, int base, mp_rnd_t mode)
{
- mpfr_init2(mp, prec);
- mpfr_set_str(mp, s.c_str(), base, mode);
+ mpfr_init2 (mpfr_ptr(), prec);
+ mpfr_set_str(mpfr_ptr(), s.c_str(), base, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
+inline void mpreal::clear(::mpfr_ptr x)
+{
+#ifdef MPREAL_HAVE_MOVE_SUPPORT
+ if(mpfr_is_initialized(x))
+#endif
+ mpfr_clear(x);
+}
+
inline mpreal::~mpreal()
{
- mpfr_clear(mp);
+ clear(mpfr_ptr());
}
// internal namespace needed for template magic
@@ -761,6 +853,9 @@ const mpreal sqrt(const int v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); const mpreal sqrt(const long double v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
const mpreal sqrt(const double v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+// abs
+inline const mpreal abs(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd());
+
//////////////////////////////////////////////////////////////////////////
// pow
const mpreal pow(const mpreal& a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
@@ -813,6 +908,11 @@ const mpreal pow(const double a, const unsigned int b, mp_rnd_t rnd_mode = mprea const mpreal pow(const double a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
const mpreal pow(const double a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+
//////////////////////////////////////////////////////////////////////////
// Estimate machine epsilon for the given precision
// Returns smallest eps such that 1.0 + eps != 1.0
@@ -860,15 +960,15 @@ inline mpreal& mpreal::operator=(const mpreal& v) {
if (this != &v)
{
- mp_prec_t tp = mpfr_get_prec(mp);
- mp_prec_t vp = mpfr_get_prec(v.mp);
+ mp_prec_t tp = mpfr_get_prec( mpfr_srcptr());
+ mp_prec_t vp = mpfr_get_prec(v.mpfr_srcptr());
- if(tp != vp){
- mpfr_clear(mp);
- mpfr_init2(mp, vp);
- }
+ if(tp != vp){
+ clear(mpfr_ptr());
+ mpfr_init2(mpfr_ptr(), vp);
+ }
- mpfr_set(mp, v.mp, mpreal::get_default_rnd());
+ mpfr_set(mpfr_ptr(), v.mpfr_srcptr(), mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
}
@@ -877,7 +977,7 @@ inline mpreal& mpreal::operator=(const mpreal& v) inline mpreal& mpreal::operator=(const mpf_t v)
{
- mpfr_set_f(mp, v, mpreal::get_default_rnd());
+ mpfr_set_f(mpfr_ptr(), v, mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
@@ -885,7 +985,7 @@ inline mpreal& mpreal::operator=(const mpf_t v) inline mpreal& mpreal::operator=(const mpz_t v)
{
- mpfr_set_z(mp, v, mpreal::get_default_rnd());
+ mpfr_set_z(mpfr_ptr(), v, mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
@@ -893,7 +993,7 @@ inline mpreal& mpreal::operator=(const mpz_t v) inline mpreal& mpreal::operator=(const mpq_t v)
{
- mpfr_set_q(mp, v, mpreal::get_default_rnd());
+ mpfr_set_q(mpfr_ptr(), v, mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
@@ -901,7 +1001,7 @@ inline mpreal& mpreal::operator=(const mpq_t v) inline mpreal& mpreal::operator=(const long double v)
{
- mpfr_set_ld(mp, v, mpreal::get_default_rnd());
+ mpfr_set_ld(mpfr_ptr(), v, mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
@@ -910,22 +1010,22 @@ inline mpreal& mpreal::operator=(const long double v) inline mpreal& mpreal::operator=(const double v)
{
#if (MPREAL_DOUBLE_BITS_OVERFLOW > -1)
- if(fits_in_bits(v, MPREAL_DOUBLE_BITS_OVERFLOW))
- {
- mpfr_set_d(mp,v,mpreal::get_default_rnd());
- }else
- throw conversion_overflow();
+ if(fits_in_bits(v, MPREAL_DOUBLE_BITS_OVERFLOW))
+ {
+ mpfr_set_d(mpfr_ptr(),v,mpreal::get_default_rnd());
+ }else
+ throw conversion_overflow();
#else
- mpfr_set_d(mp,v,mpreal::get_default_rnd());
+ mpfr_set_d(mpfr_ptr(),v,mpreal::get_default_rnd());
#endif
- MPREAL_MSVC_DEBUGVIEW_CODE;
+ MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator=(const unsigned long int v)
{
- mpfr_set_ui(mp, v, mpreal::get_default_rnd());
+ mpfr_set_ui(mpfr_ptr(), v, mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
@@ -933,7 +1033,7 @@ inline mpreal& mpreal::operator=(const unsigned long int v) inline mpreal& mpreal::operator=(const unsigned int v)
{
- mpfr_set_ui(mp, v, mpreal::get_default_rnd());
+ mpfr_set_ui(mpfr_ptr(), v, mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
@@ -941,7 +1041,7 @@ inline mpreal& mpreal::operator=(const unsigned int v) inline mpreal& mpreal::operator=(const long int v)
{
- mpfr_set_si(mp, v, mpreal::get_default_rnd());
+ mpfr_set_si(mpfr_ptr(), v, mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
@@ -949,7 +1049,7 @@ inline mpreal& mpreal::operator=(const long int v) inline mpreal& mpreal::operator=(const int v)
{
- mpfr_set_si(mp, v, mpreal::get_default_rnd());
+ mpfr_set_si(mpfr_ptr(), v, mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
@@ -966,15 +1066,15 @@ inline mpreal& mpreal::operator=(const char* s) mpfr_t t;
- mpfr_init2(t, mpfr_get_prec(mp));
+ mpfr_init2(t, mpfr_get_prec(mpfr_srcptr()));
if(0 == mpfr_set_str(t, s, 10, mpreal::get_default_rnd()))
{
- mpfr_set(mp, t, mpreal::get_default_rnd());
+ mpfr_set(mpfr_ptr(), t, mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
}
- mpfr_clear(t);
+ clear(t);
return *this;
}
@@ -989,15 +1089,15 @@ inline mpreal& mpreal::operator=(const std::string& s) mpfr_t t;
- mpfr_init2(t, mpfr_get_prec(mp));
+ mpfr_init2(t, mpfr_get_prec(mpfr_srcptr()));
if(0 == mpfr_set_str(t, s.c_str(), 10, mpreal::get_default_rnd()))
{
- mpfr_set(mp, t, mpreal::get_default_rnd());
+ mpfr_set(mpfr_ptr(), t, mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
}
- mpfr_clear(t);
+ clear(t);
return *this;
}
@@ -1006,7 +1106,7 @@ inline mpreal& mpreal::operator=(const std::string& s) // + Addition
inline mpreal& mpreal::operator+=(const mpreal& v)
{
- mpfr_add(mp,mp,v.mp,mpreal::get_default_rnd());
+ mpfr_add(mpfr_ptr(), mpfr_srcptr(), v.mpfr_srcptr(), mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
@@ -1020,14 +1120,14 @@ inline mpreal& mpreal::operator+=(const mpf_t u) inline mpreal& mpreal::operator+=(const mpz_t u)
{
- mpfr_add_z(mp,mp,u,mpreal::get_default_rnd());
+ mpfr_add_z(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator+=(const mpq_t u)
{
- mpfr_add_q(mp,mp,u,mpreal::get_default_rnd());
+ mpfr_add_q(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
@@ -1042,7 +1142,7 @@ inline mpreal& mpreal::operator+= (const long double u) inline mpreal& mpreal::operator+= (const double u)
{
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
- mpfr_add_d(mp,mp,u,mpreal::get_default_rnd());
+ mpfr_add_d(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
#else
*this += mpreal(u);
#endif
@@ -1053,28 +1153,28 @@ inline mpreal& mpreal::operator+= (const double u) inline mpreal& mpreal::operator+=(const unsigned long int u)
{
- mpfr_add_ui(mp,mp,u,mpreal::get_default_rnd());
+ mpfr_add_ui(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator+=(const unsigned int u)
{
- mpfr_add_ui(mp,mp,u,mpreal::get_default_rnd());
+ mpfr_add_ui(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator+=(const long int u)
{
- mpfr_add_si(mp,mp,u,mpreal::get_default_rnd());
+ mpfr_add_si(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator+=(const int u)
{
- mpfr_add_si(mp,mp,u,mpreal::get_default_rnd());
+ mpfr_add_si(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
@@ -1094,9 +1194,9 @@ inline const mpreal mpreal::operator+()const { return mpreal(*this); } inline const mpreal operator+(const mpreal& a, const mpreal& b)
{
- mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));
- mpfr_add(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
- return c;
+ mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));
+ mpfr_add(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
+ return c;
}
inline mpreal& mpreal::operator++()
@@ -1127,21 +1227,21 @@ inline const mpreal mpreal::operator-- (int) // - Subtraction
inline mpreal& mpreal::operator-=(const mpreal& v)
{
- mpfr_sub(mp,mp,v.mp,mpreal::get_default_rnd());
+ mpfr_sub(mpfr_ptr(),mpfr_srcptr(),v.mpfr_srcptr(),mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator-=(const mpz_t v)
{
- mpfr_sub_z(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_sub_z(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator-=(const mpq_t v)
{
- mpfr_sub_q(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_sub_q(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
@@ -1156,7 +1256,7 @@ inline mpreal& mpreal::operator-=(const long double v) inline mpreal& mpreal::operator-=(const double v)
{
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
- mpfr_sub_d(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_sub_d(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
#else
*this -= mpreal(v);
#endif
@@ -1167,28 +1267,28 @@ inline mpreal& mpreal::operator-=(const double v) inline mpreal& mpreal::operator-=(const unsigned long int v)
{
- mpfr_sub_ui(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_sub_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator-=(const unsigned int v)
{
- mpfr_sub_ui(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_sub_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator-=(const long int v)
{
- mpfr_sub_si(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_sub_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator-=(const int v)
{
- mpfr_sub_si(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_sub_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
@@ -1196,15 +1296,15 @@ inline mpreal& mpreal::operator-=(const int v) inline const mpreal mpreal::operator-()const
{
mpreal u(*this);
- mpfr_neg(u.mp,u.mp,mpreal::get_default_rnd());
+ mpfr_neg(u.mpfr_ptr(),u.mpfr_srcptr(),mpreal::get_default_rnd());
return u;
}
inline const mpreal operator-(const mpreal& a, const mpreal& b)
{
- mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));
- mpfr_sub(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
- return c;
+ mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));
+ mpfr_sub(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
+ return c;
}
inline const mpreal operator-(const double b, const mpreal& a)
@@ -1252,21 +1352,21 @@ inline const mpreal operator-(const int b, const mpreal& a) // * Multiplication
inline mpreal& mpreal::operator*= (const mpreal& v)
{
- mpfr_mul(mp,mp,v.mp,mpreal::get_default_rnd());
+ mpfr_mul(mpfr_ptr(),mpfr_srcptr(),v.mpfr_srcptr(),mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator*=(const mpz_t v)
{
- mpfr_mul_z(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_mul_z(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator*=(const mpq_t v)
{
- mpfr_mul_q(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_mul_q(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
@@ -1281,7 +1381,7 @@ inline mpreal& mpreal::operator*=(const long double v) inline mpreal& mpreal::operator*=(const double v)
{
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
- mpfr_mul_d(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_mul_d(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
#else
*this *= mpreal(v);
#endif
@@ -1291,58 +1391,58 @@ inline mpreal& mpreal::operator*=(const double v) inline mpreal& mpreal::operator*=(const unsigned long int v)
{
- mpfr_mul_ui(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_mul_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator*=(const unsigned int v)
{
- mpfr_mul_ui(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_mul_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator*=(const long int v)
{
- mpfr_mul_si(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_mul_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator*=(const int v)
{
- mpfr_mul_si(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_mul_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline const mpreal operator*(const mpreal& a, const mpreal& b)
{
- mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));
- mpfr_mul(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
- return c;
+ mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));
+ mpfr_mul(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
+ return c;
}
//////////////////////////////////////////////////////////////////////////
// / Division
inline mpreal& mpreal::operator/=(const mpreal& v)
{
- mpfr_div(mp,mp,v.mp,mpreal::get_default_rnd());
+ mpfr_div(mpfr_ptr(),mpfr_srcptr(),v.mpfr_srcptr(),mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator/=(const mpz_t v)
{
- mpfr_div_z(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_div_z(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator/=(const mpq_t v)
{
- mpfr_div_q(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_div_q(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
@@ -1357,7 +1457,7 @@ inline mpreal& mpreal::operator/=(const long double v) inline mpreal& mpreal::operator/=(const double v)
{
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
- mpfr_div_d(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_div_d(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
#else
*this /= mpreal(v);
#endif
@@ -1367,72 +1467,72 @@ inline mpreal& mpreal::operator/=(const double v) inline mpreal& mpreal::operator/=(const unsigned long int v)
{
- mpfr_div_ui(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_div_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator/=(const unsigned int v)
{
- mpfr_div_ui(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_div_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator/=(const long int v)
{
- mpfr_div_si(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_div_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator/=(const int v)
{
- mpfr_div_si(mp,mp,v,mpreal::get_default_rnd());
+ mpfr_div_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline const mpreal operator/(const mpreal& a, const mpreal& b)
{
- mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));
- mpfr_div(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
- return c;
+ mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_srcptr()), mpfr_get_prec(b.mpfr_srcptr())));
+ mpfr_div(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
+ return c;
}
inline const mpreal operator/(const unsigned long int b, const mpreal& a)
{
- mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));
+ mpreal x(0, mpfr_get_prec(a.mpfr_srcptr()));
mpfr_ui_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());
return x;
}
inline const mpreal operator/(const unsigned int b, const mpreal& a)
{
- mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));
+ mpreal x(0, mpfr_get_prec(a.mpfr_srcptr()));
mpfr_ui_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());
return x;
}
inline const mpreal operator/(const long int b, const mpreal& a)
{
- mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));
- mpfr_si_div(x.mpfr_ptr(), b, a.mpfr_srcptr(),mpreal::get_default_rnd());
+ mpreal x(0, mpfr_get_prec(a.mpfr_srcptr()));
+ mpfr_si_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());
return x;
}
inline const mpreal operator/(const int b, const mpreal& a)
{
- mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));
- mpfr_si_div(x.mpfr_ptr(), b, a.mpfr_srcptr(),mpreal::get_default_rnd());
+ mpreal x(0, mpfr_get_prec(a.mpfr_srcptr()));
+ mpfr_si_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());
return x;
}
inline const mpreal operator/(const double b, const mpreal& a)
{
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
- mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));
- mpfr_d_div(x.mpfr_ptr(), b, a.mpfr_srcptr(),mpreal::get_default_rnd());
+ mpreal x(0, mpfr_get_prec(a.mpfr_srcptr()));
+ mpfr_d_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());
return x;
#else
mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));
@@ -1445,56 +1545,56 @@ inline const mpreal operator/(const double b, const mpreal& a) // Shifts operators - Multiplication/Division by power of 2
inline mpreal& mpreal::operator<<=(const unsigned long int u)
{
- mpfr_mul_2ui(mp,mp,u,mpreal::get_default_rnd());
+ mpfr_mul_2ui(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator<<=(const unsigned int u)
{
- mpfr_mul_2ui(mp,mp,static_cast<unsigned long int>(u),mpreal::get_default_rnd());
+ mpfr_mul_2ui(mpfr_ptr(),mpfr_srcptr(),static_cast<unsigned long int>(u),mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator<<=(const long int u)
{
- mpfr_mul_2si(mp,mp,u,mpreal::get_default_rnd());
+ mpfr_mul_2si(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator<<=(const int u)
{
- mpfr_mul_2si(mp,mp,static_cast<long int>(u),mpreal::get_default_rnd());
+ mpfr_mul_2si(mpfr_ptr(),mpfr_srcptr(),static_cast<long int>(u),mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator>>=(const unsigned long int u)
{
- mpfr_div_2ui(mp,mp,u,mpreal::get_default_rnd());
+ mpfr_div_2ui(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator>>=(const unsigned int u)
{
- mpfr_div_2ui(mp,mp,static_cast<unsigned long int>(u),mpreal::get_default_rnd());
+ mpfr_div_2ui(mpfr_ptr(),mpfr_srcptr(),static_cast<unsigned long int>(u),mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator>>=(const long int u)
{
- mpfr_div_2si(mp,mp,u,mpreal::get_default_rnd());
+ mpfr_div_2si(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator>>=(const int u)
{
- mpfr_div_2si(mp,mp,static_cast<long int>(u),mpreal::get_default_rnd());
+ mpfr_div_2si(mpfr_ptr(),mpfr_srcptr(),static_cast<long int>(u),mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
@@ -1543,7 +1643,7 @@ inline const mpreal operator>>(const mpreal& v, const int k) inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode)
{
mpreal x(v);
- mpfr_mul_2ui(x.mp,v.mp,k,rnd_mode);
+ mpfr_mul_2ui(x.mpfr_ptr(),v.mpfr_srcptr(),k,rnd_mode);
return x;
}
@@ -1551,61 +1651,63 @@ inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_m inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode)
{
mpreal x(v);
- mpfr_mul_2si(x.mp,v.mp,k,rnd_mode);
+ mpfr_mul_2si(x.mpfr_ptr(),v.mpfr_srcptr(),k,rnd_mode);
return x;
}
inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode)
{
mpreal x(v);
- mpfr_div_2ui(x.mp,v.mp,k,rnd_mode);
+ mpfr_div_2ui(x.mpfr_ptr(),v.mpfr_srcptr(),k,rnd_mode);
return x;
}
inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode)
{
mpreal x(v);
- mpfr_div_2si(x.mp,v.mp,k,rnd_mode);
+ mpfr_div_2si(x.mpfr_ptr(),v.mpfr_srcptr(),k,rnd_mode);
return x;
}
//////////////////////////////////////////////////////////////////////////
//Boolean operators
-inline bool operator > (const mpreal& a, const mpreal& b){ return (mpfr_greater_p(a.mp,b.mp) !=0); }
-inline bool operator >= (const mpreal& a, const mpreal& b){ return (mpfr_greaterequal_p(a.mp,b.mp) !=0); }
-inline bool operator < (const mpreal& a, const mpreal& b){ return (mpfr_less_p(a.mp,b.mp) !=0); }
-inline bool operator <= (const mpreal& a, const mpreal& b){ return (mpfr_lessequal_p(a.mp,b.mp) !=0); }
-inline bool operator == (const mpreal& a, const mpreal& b){ return (mpfr_equal_p(a.mp,b.mp) !=0); }
-inline bool operator != (const mpreal& a, const mpreal& b){ return (mpfr_lessgreater_p(a.mp,b.mp) !=0); }
-
-inline bool operator == (const mpreal& a, const unsigned long int b ){ return (mpfr_cmp_ui(a.mp,b) == 0); }
-inline bool operator == (const mpreal& a, const unsigned int b ){ return (mpfr_cmp_ui(a.mp,b) == 0); }
-inline bool operator == (const mpreal& a, const long int b ){ return (mpfr_cmp_si(a.mp,b) == 0); }
-inline bool operator == (const mpreal& a, const int b ){ return (mpfr_cmp_si(a.mp,b) == 0); }
-inline bool operator == (const mpreal& a, const long double b ){ return (mpfr_cmp_ld(a.mp,b) == 0); }
-inline bool operator == (const mpreal& a, const double b ){ return (mpfr_cmp_d(a.mp,b) == 0); }
-
-
-inline bool isnan (const mpreal& v){ return (mpfr_nan_p(v.mp) != 0); }
-inline bool isinf (const mpreal& v){ return (mpfr_inf_p(v.mp) != 0); }
-inline bool isfinite (const mpreal& v){ return (mpfr_number_p(v.mp) != 0); }
-inline bool iszero (const mpreal& v){ return (mpfr_zero_p(v.mp) != 0); }
-inline bool isint (const mpreal& v){ return (mpfr_integer_p(v.mp) != 0); }
+inline bool operator > (const mpreal& a, const mpreal& b){ return (mpfr_greater_p (a.mpfr_srcptr(),b.mpfr_srcptr()) !=0 ); }
+inline bool operator >= (const mpreal& a, const mpreal& b){ return (mpfr_greaterequal_p (a.mpfr_srcptr(),b.mpfr_srcptr()) !=0 ); }
+inline bool operator < (const mpreal& a, const mpreal& b){ return (mpfr_less_p (a.mpfr_srcptr(),b.mpfr_srcptr()) !=0 ); }
+inline bool operator <= (const mpreal& a, const mpreal& b){ return (mpfr_lessequal_p (a.mpfr_srcptr(),b.mpfr_srcptr()) !=0 ); }
+inline bool operator == (const mpreal& a, const mpreal& b){ return (mpfr_equal_p (a.mpfr_srcptr(),b.mpfr_srcptr()) !=0 ); }
+inline bool operator != (const mpreal& a, const mpreal& b){ return (mpfr_lessgreater_p (a.mpfr_srcptr(),b.mpfr_srcptr()) !=0 ); }
+
+inline bool operator == (const mpreal& a, const unsigned long int b ){ return (mpfr_cmp_ui(a.mpfr_srcptr(),b) == 0 ); }
+inline bool operator == (const mpreal& a, const unsigned int b ){ return (mpfr_cmp_ui(a.mpfr_srcptr(),b) == 0 ); }
+inline bool operator == (const mpreal& a, const long int b ){ return (mpfr_cmp_si(a.mpfr_srcptr(),b) == 0 ); }
+inline bool operator == (const mpreal& a, const int b ){ return (mpfr_cmp_si(a.mpfr_srcptr(),b) == 0 ); }
+inline bool operator == (const mpreal& a, const long double b ){ return (mpfr_cmp_ld(a.mpfr_srcptr(),b) == 0 ); }
+inline bool operator == (const mpreal& a, const double b ){ return (mpfr_cmp_d (a.mpfr_srcptr(),b) == 0 ); }
+
+
+inline bool isnan (const mpreal& op){ return (mpfr_nan_p (op.mpfr_srcptr()) != 0 ); }
+inline bool isinf (const mpreal& op){ return (mpfr_inf_p (op.mpfr_srcptr()) != 0 ); }
+inline bool isfinite (const mpreal& op){ return (mpfr_number_p (op.mpfr_srcptr()) != 0 ); }
+inline bool iszero (const mpreal& op){ return (mpfr_zero_p (op.mpfr_srcptr()) != 0 ); }
+inline bool isint (const mpreal& op){ return (mpfr_integer_p(op.mpfr_srcptr()) != 0 ); }
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
-inline bool isregular(const mpreal& v){ return (mpfr_regular_p(v.mp));}
+inline bool isregular(const mpreal& op){ return (mpfr_regular_p(op.mpfr_srcptr()));}
#endif
//////////////////////////////////////////////////////////////////////////
// Type Converters
-inline long mpreal::toLong (mp_rnd_t mode) const { return mpfr_get_si(mp, mode); }
-inline unsigned long mpreal::toULong (mp_rnd_t mode) const { return mpfr_get_ui(mp, mode); }
-inline double mpreal::toDouble (mp_rnd_t mode) const { return mpfr_get_d (mp, mode); }
-inline long double mpreal::toLDouble(mp_rnd_t mode) const { return mpfr_get_ld(mp, mode); }
+inline bool mpreal::toBool (mp_rnd_t /*mode*/) const { return mpfr_zero_p (mpfr_srcptr()) == 0; }
+inline long mpreal::toLong (mp_rnd_t mode) const { return mpfr_get_si (mpfr_srcptr(), mode); }
+inline unsigned long mpreal::toULong (mp_rnd_t mode) const { return mpfr_get_ui (mpfr_srcptr(), mode); }
+inline float mpreal::toFloat (mp_rnd_t mode) const { return mpfr_get_flt(mpfr_srcptr(), mode); }
+inline double mpreal::toDouble (mp_rnd_t mode) const { return mpfr_get_d (mpfr_srcptr(), mode); }
+inline long double mpreal::toLDouble(mp_rnd_t mode) const { return mpfr_get_ld (mpfr_srcptr(), mode); }
#if defined (MPREAL_HAVE_INT64_SUPPORT)
-inline int64_t mpreal::toInt64 (mp_rnd_t mode) const{ return mpfr_get_sj(mp, mode); }
-inline uint64_t mpreal::toUInt64(mp_rnd_t mode) const{ return mpfr_get_uj(mp, mode); }
+inline int64_t mpreal::toInt64 (mp_rnd_t mode) const{ return mpfr_get_sj(mpfr_srcptr(), mode); }
+inline uint64_t mpreal::toUInt64(mp_rnd_t mode) const{ return mpfr_get_uj(mpfr_srcptr(), mode); }
#endif
inline ::mpfr_ptr mpreal::mpfr_ptr() { return mp; }
@@ -1629,7 +1731,7 @@ inline std::string mpreal::toString(const std::string& format) const if( !format.empty() )
{
- if(!(mpfr_asprintf(&s,format.c_str(),mp) < 0))
+ if(!(mpfr_asprintf(&s, format.c_str(), mpfr_srcptr()) < 0))
{
out = std::string(s);
@@ -1644,20 +1746,19 @@ inline std::string mpreal::toString(const std::string& format) const inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const
{
+ // TODO: Add extended format specification (f, e, rounding mode) as it done in output operator
(void)b;
(void)mode;
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
- // Use MPFR native function for output
- char format[128];
- int digits;
+ std::ostringstream format;
- digits = n > 0 ? n : bits2digits(mpfr_get_prec(mp));
-
- sprintf(format,"%%.%dRNg",digits); // Default format
+ int digits = (n >= 0) ? n : bits2digits(mpfr_get_prec(mpfr_srcptr()));
+
+ format << "%." << digits << "RNg";
- return toString(std::string(format));
+ return toString(format.str());
#else
@@ -1675,8 +1776,8 @@ inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const if(mpfr_zero_p(mp)) return "0";
if(mpfr_nan_p(mp)) return "NaN";
- s = mpfr_get_str(NULL,&exp,b,0,mp,mode);
- ns = mpfr_get_str(NULL,&exp,b,n,mp,mode);
+ s = mpfr_get_str(NULL, &exp, b, 0, mp, mode);
+ ns = mpfr_get_str(NULL, &exp, b, (std::max)(0,n), mp, mode);
if(s!=NULL && ns!=NULL)
{
@@ -1761,17 +1862,43 @@ inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const //////////////////////////////////////////////////////////////////////////
// I/O
+inline std::ostream& mpreal::output(std::ostream& os) const
+{
+ std::ostringstream format;
+ const std::ios::fmtflags flags = os.flags();
+
+ format << ((flags & std::ios::showpos) ? "%+" : "%");
+ if (os.precision() >= 0)
+ format << '.' << os.precision() << "R*"
+ << ((flags & std::ios::floatfield) == std::ios::fixed ? 'f' :
+ (flags & std::ios::floatfield) == std::ios::scientific ? 'e' :
+ 'g');
+ else
+ format << "R*e";
+
+ char *s = NULL;
+ if(!(mpfr_asprintf(&s, format.str().c_str(),
+ mpfr::mpreal::get_default_rnd(),
+ mpfr_srcptr())
+ < 0))
+ {
+ os << std::string(s);
+ mpfr_free_str(s);
+ }
+ return os;
+}
+
inline std::ostream& operator<<(std::ostream& os, const mpreal& v)
{
- return os<<v.toString(static_cast<int>(os.precision()));
+ return v.output(os);
}
inline std::istream& operator>>(std::istream &is, mpreal& v)
{
- // ToDo, use cout::hexfloat and other flags to setup base
+ // TODO: use cout::hexfloat and other flags to setup base
std::string tmp;
is >> tmp;
- mpfr_set_str(v.mp, tmp.c_str(), 10, mpreal::get_default_rnd());
+ mpfr_set_str(v.mpfr_ptr(), tmp.c_str(), 10, mpreal::get_default_rnd());
return is;
}
@@ -1784,53 +1911,53 @@ inline mp_prec_t digits2bits(int d) {
const double LOG2_10 = 3.3219280948873624;
- return (mp_prec_t) std::ceil( d * LOG2_10 );
+ return mp_prec_t(std::ceil( d * LOG2_10 ));
}
inline int bits2digits(mp_prec_t b)
{
const double LOG10_2 = 0.30102999566398119;
- return (int) std::floor( b * LOG10_2 );
+ return int(std::floor( b * LOG10_2 ));
}
//////////////////////////////////////////////////////////////////////////
// Set/Get number properties
-inline int sgn(const mpreal& v)
+inline int sgn(const mpreal& op)
{
- int r = mpfr_signbit(v.mp);
- return (r>0?-1:1);
+ int r = mpfr_signbit(op.mpfr_srcptr());
+ return (r > 0? -1 : 1);
}
inline mpreal& mpreal::setSign(int sign, mp_rnd_t RoundingMode)
{
- mpfr_setsign(mp,mp,(sign < 0 ? 1 : 0),RoundingMode);
+ mpfr_setsign(mpfr_ptr(), mpfr_srcptr(), (sign < 0 ? 1 : 0), RoundingMode);
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline int mpreal::getPrecision() const
{
- return mpfr_get_prec(mp);
+ return int(mpfr_get_prec(mpfr_srcptr()));
}
inline mpreal& mpreal::setPrecision(int Precision, mp_rnd_t RoundingMode)
{
- mpfr_prec_round(mp, Precision, RoundingMode);
+ mpfr_prec_round(mpfr_ptr(), Precision, RoundingMode);
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::setInf(int sign)
{
- mpfr_set_inf(mp,sign);
+ mpfr_set_inf(mpfr_ptr(), sign);
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::setNan()
{
- mpfr_set_nan(mp);
+ mpfr_set_nan(mpfr_ptr());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
@@ -1839,9 +1966,9 @@ inline mpreal& mpreal::setZero(int sign) {
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
- mpfr_set_zero(mp, sign);
+ mpfr_set_zero(mpfr_ptr(), sign);
#else
- mpfr_set_si(mp, 0, (mpfr_get_default_rounding_mode)());
+ mpfr_set_si(mpfr_ptr(), 0, (mpfr_get_default_rounding_mode)());
setSign(sign);
#endif
@@ -1851,23 +1978,23 @@ inline mpreal& mpreal::setZero(int sign) inline mp_prec_t mpreal::get_prec() const
{
- return mpfr_get_prec(mp);
+ return mpfr_get_prec(mpfr_srcptr());
}
inline void mpreal::set_prec(mp_prec_t prec, mp_rnd_t rnd_mode)
{
- mpfr_prec_round(mp,prec,rnd_mode);
+ mpfr_prec_round(mpfr_ptr(),prec,rnd_mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mp_exp_t mpreal::get_exp ()
{
- return mpfr_get_exp(mp);
+ return mpfr_get_exp(mpfr_srcptr());
}
inline int mpreal::set_exp (mp_exp_t e)
{
- int x = mpfr_set_exp(mp, e);
+ int x = mpfr_set_exp(mpfr_ptr(), e);
MPREAL_MSVC_DEBUGVIEW_CODE;
return x;
}
@@ -1885,7 +2012,7 @@ inline const mpreal ldexp(const mpreal& v, mp_exp_t exp) mpreal x(v);
// rounding is not important since we just increasing the exponent
- mpfr_mul_2si(x.mp,x.mp,exp,mpreal::get_default_rnd());
+ mpfr_mul_2si(x.mpfr_ptr(), x.mpfr_srcptr(), exp, mpreal::get_default_rnd());
return x;
}
@@ -1900,9 +2027,9 @@ inline mpreal machine_epsilon(const mpreal& x) /* the smallest eps such that x + eps != x */
if( x < 0)
{
- return nextabove(-x)+x;
+ return nextabove(-x) + x;
}else{
- return nextabove(x)-x;
+ return nextabove( x) - x;
}
}
@@ -1922,37 +2049,37 @@ inline mpreal maxval(mp_prec_t prec) inline bool isEqualUlps(const mpreal& a, const mpreal& b, int maxUlps)
{
- return abs(a - b) <= machine_epsilon((max)(abs(a), abs(b))) * maxUlps;
+ return abs(a - b) <= machine_epsilon((max)(abs(a), abs(b))) * maxUlps;
}
inline bool isEqualFuzzy(const mpreal& a, const mpreal& b, const mpreal& eps)
{
- return abs(a - b) <= (min)(abs(a), abs(b)) * eps;
+ return abs(a - b) <= eps;
}
inline bool isEqualFuzzy(const mpreal& a, const mpreal& b)
{
- return isEqualFuzzy(a, b, machine_epsilon((min)(abs(a), abs(b))));
+ return isEqualFuzzy(a, b, machine_epsilon((max)(1, (min)(abs(a), abs(b)))));
}
inline const mpreal modf(const mpreal& v, mpreal& n)
{
- mpreal frac(v);
+ mpreal f(v);
// rounding is not important since we are using the same number
- mpfr_frac(frac.mp,frac.mp,mpreal::get_default_rnd());
- mpfr_trunc(n.mp,v.mp);
- return frac;
+ mpfr_frac (f.mpfr_ptr(),f.mpfr_srcptr(),mpreal::get_default_rnd());
+ mpfr_trunc(n.mpfr_ptr(),v.mpfr_srcptr());
+ return f;
}
inline int mpreal::check_range (int t, mp_rnd_t rnd_mode)
{
- return mpfr_check_range(mp,t,rnd_mode);
+ return mpfr_check_range(mpfr_ptr(),t,rnd_mode);
}
inline int mpreal::subnormalize (int t,mp_rnd_t rnd_mode)
{
- int r = mpfr_subnormalize(mp,t,rnd_mode);
+ int r = mpfr_subnormalize(mpfr_ptr(),t,rnd_mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
return r;
}
@@ -2005,8 +2132,11 @@ inline mp_exp_t mpreal::get_emax_max (void) mpfr_##f(y.mpfr_ptr(), x.mpfr_srcptr(), r); \
return y;
-inline const mpreal sqr (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(sqr ); }
-inline const mpreal sqrt (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(sqrt); }
+inline const mpreal sqr (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd())
+{ MPREAL_UNARY_MATH_FUNCTION_BODY(sqr ); }
+
+inline const mpreal sqrt (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd())
+{ MPREAL_UNARY_MATH_FUNCTION_BODY(sqrt); }
inline const mpreal sqrt(const unsigned long int x, mp_rnd_t r)
{
@@ -2032,14 +2162,14 @@ inline const mpreal sqrt(const int v, mp_rnd_t rnd_mode) else return mpreal().setNan(); // NaN
}
-inline const mpreal root(const mpreal& x, unsigned long int k, mp_rnd_t r)
+inline const mpreal root(const mpreal& x, unsigned long int k, mp_rnd_t r = mpreal::get_default_rnd())
{
mpreal y(0, mpfr_get_prec(x.mpfr_srcptr()));
mpfr_root(y.mpfr_ptr(), x.mpfr_srcptr(), k, r);
return y;
}
-inline const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t r)
+inline const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t r = mpreal::get_default_rnd())
{
mpreal y(0, mpfr_get_prec(a.mpfr_srcptr()));
mpfr_dim(y.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), r);
@@ -2048,161 +2178,130 @@ inline const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t r) inline int cmpabs(const mpreal& a,const mpreal& b)
{
- return mpfr_cmpabs(a.mp,b.mp);
+ return mpfr_cmpabs(a.mpfr_ptr(), b.mpfr_srcptr());
}
-inline int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode)
+inline int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
- return mpfr_sin_cos(s.mp,c.mp,v.mp,rnd_mode);
+ return mpfr_sin_cos(s.mpfr_ptr(), c.mpfr_ptr(), v.mpfr_srcptr(), rnd_mode);
}
inline const mpreal sqrt (const long double v, mp_rnd_t rnd_mode) { return sqrt(mpreal(v),rnd_mode); }
inline const mpreal sqrt (const double v, mp_rnd_t rnd_mode) { return sqrt(mpreal(v),rnd_mode); }
-inline const mpreal cbrt (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(cbrt ); }
-inline const mpreal fabs (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(abs ); }
-inline const mpreal abs (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(abs ); }
-inline const mpreal log (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(log ); }
-inline const mpreal log2 (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(log2 ); }
-inline const mpreal log10 (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(log10); }
-inline const mpreal exp (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(exp ); }
-inline const mpreal exp2 (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(exp2 ); }
-inline const mpreal exp10 (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(exp10); }
-inline const mpreal cos (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(cos ); }
-inline const mpreal sin (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(sin ); }
-inline const mpreal tan (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(tan ); }
-inline const mpreal sec (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(sec ); }
-inline const mpreal csc (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(csc ); }
-inline const mpreal cot (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(cot ); }
-inline const mpreal acos (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(acos); }
-inline const mpreal asin (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(asin); }
-inline const mpreal atan (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(atan); }
-
-inline const mpreal acot (const mpreal& v, mp_rnd_t r) { return atan (1/v, r); }
-inline const mpreal asec (const mpreal& v, mp_rnd_t r) { return acos (1/v, r); }
-inline const mpreal acsc (const mpreal& v, mp_rnd_t r) { return asin (1/v, r); }
-inline const mpreal acoth (const mpreal& v, mp_rnd_t r) { return atanh(1/v, r); }
-inline const mpreal asech (const mpreal& v, mp_rnd_t r) { return acosh(1/v, r); }
-inline const mpreal acsch (const mpreal& v, mp_rnd_t r) { return asinh(1/v, r); }
-
-inline const mpreal cosh (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(cosh ); }
-inline const mpreal sinh (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(sinh ); }
-inline const mpreal tanh (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(tanh ); }
-inline const mpreal sech (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(sech ); }
-inline const mpreal csch (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(csch ); }
-inline const mpreal coth (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(coth ); }
-inline const mpreal acosh (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(acosh); }
-inline const mpreal asinh (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(asinh); }
-inline const mpreal atanh (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(atanh); }
-
-inline const mpreal log1p (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(log1p ); }
-inline const mpreal expm1 (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(expm1 ); }
-inline const mpreal eint (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(eint ); }
-inline const mpreal gamma (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(gamma ); }
-inline const mpreal lngamma (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(lngamma); }
-inline const mpreal zeta (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(zeta ); }
-inline const mpreal erf (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(erf ); }
-inline const mpreal erfc (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(erfc ); }
-inline const mpreal besselj0(const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(j0 ); }
-inline const mpreal besselj1(const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(j1 ); }
-inline const mpreal bessely0(const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(y0 ); }
-inline const mpreal bessely1(const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(y1 ); }
-
-inline const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode)
-{
- mpreal a;
- mp_prec_t yp, xp;
-
- yp = y.get_prec();
- xp = x.get_prec();
-
- a.set_prec(yp>xp?yp:xp);
-
- mpfr_atan2(a.mp, y.mp, x.mp, rnd_mode);
-
+inline const mpreal cbrt (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(cbrt ); }
+inline const mpreal fabs (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(abs ); }
+inline const mpreal abs (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(abs ); }
+inline const mpreal log (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(log ); }
+inline const mpreal log2 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(log2 ); }
+inline const mpreal log10 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(log10); }
+inline const mpreal exp (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(exp ); }
+inline const mpreal exp2 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(exp2 ); }
+inline const mpreal exp10 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(exp10); }
+inline const mpreal cos (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(cos ); }
+inline const mpreal sin (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(sin ); }
+inline const mpreal tan (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(tan ); }
+inline const mpreal sec (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(sec ); }
+inline const mpreal csc (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(csc ); }
+inline const mpreal cot (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(cot ); }
+inline const mpreal acos (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(acos ); }
+inline const mpreal asin (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(asin ); }
+inline const mpreal atan (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(atan ); }
+
+inline const mpreal acot (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return atan (1/v, r); }
+inline const mpreal asec (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return acos (1/v, r); }
+inline const mpreal acsc (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return asin (1/v, r); }
+inline const mpreal acoth (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return atanh(1/v, r); }
+inline const mpreal asech (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return acosh(1/v, r); }
+inline const mpreal acsch (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return asinh(1/v, r); }
+
+inline const mpreal cosh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(cosh ); }
+inline const mpreal sinh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(sinh ); }
+inline const mpreal tanh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(tanh ); }
+inline const mpreal sech (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(sech ); }
+inline const mpreal csch (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(csch ); }
+inline const mpreal coth (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(coth ); }
+inline const mpreal acosh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(acosh); }
+inline const mpreal asinh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(asinh); }
+inline const mpreal atanh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(atanh); }
+
+inline const mpreal log1p (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(log1p ); }
+inline const mpreal expm1 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(expm1 ); }
+inline const mpreal eint (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(eint ); }
+inline const mpreal gamma (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(gamma ); }
+inline const mpreal lngamma (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(lngamma); }
+inline const mpreal zeta (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(zeta ); }
+inline const mpreal erf (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(erf ); }
+inline const mpreal erfc (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(erfc ); }
+inline const mpreal besselj0(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(j0 ); }
+inline const mpreal besselj1(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(j1 ); }
+inline const mpreal bessely0(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(y0 ); }
+inline const mpreal bessely1(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(y1 ); }
+
+inline const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
+{
+ mpreal a(0,(std::max)(y.getPrecision(), x.getPrecision()));
+ mpfr_atan2(a.mpfr_ptr(), y.mpfr_srcptr(), x.mpfr_srcptr(), rnd_mode);
return a;
}
-inline const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
+inline const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
- mpreal a;
- mp_prec_t yp, xp;
-
- yp = y.get_prec();
- xp = x.get_prec();
-
- a.set_prec(yp>xp?yp:xp);
-
- mpfr_hypot(a.mp, x.mp, y.mp, rnd_mode);
-
+ mpreal a(0,(std::max)(y.getPrecision(), x.getPrecision()));
+ mpfr_hypot(a.mpfr_ptr(), x.mpfr_srcptr(), y.mpfr_srcptr(), rnd_mode);
return a;
}
-inline const mpreal remainder (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
+inline const mpreal remainder (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
- mpreal a;
- mp_prec_t yp, xp;
-
- yp = y.get_prec();
- xp = x.get_prec();
-
- a.set_prec(yp>xp?yp:xp);
-
- mpfr_remainder(a.mp, x.mp, y.mp, rnd_mode);
-
+ mpreal a(0,(std::max)(y.getPrecision(), x.getPrecision()));
+ mpfr_remainder(a.mpfr_ptr(), x.mpfr_srcptr(), y.mpfr_srcptr(), rnd_mode);
return a;
}
-inline const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
+inline const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
- mpreal a;
- mp_prec_t yp, xp;
-
- yp = y.get_prec();
- xp = x.get_prec();
-
- a.set_prec(yp>xp?yp:xp);
-
- mpfr_remquo(a.mp,q, x.mp, y.mp, rnd_mode);
-
+ mpreal a(0,(std::max)(y.getPrecision(), x.getPrecision()));
+ mpfr_remquo(a.mpfr_ptr(),q, x.mpfr_srcptr(), y.mpfr_srcptr(), rnd_mode);
return a;
}
-inline const mpreal fac_ui (unsigned long int v, mp_prec_t prec, mp_rnd_t rnd_mode)
+inline const mpreal fac_ui (unsigned long int v, mp_prec_t prec = mpreal::get_default_prec(),
+ mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal x(0, prec);
- mpfr_fac_ui(x.mp,v,rnd_mode);
+ mpfr_fac_ui(x.mpfr_ptr(),v,rnd_mode);
return x;
}
-inline const mpreal lgamma (const mpreal& v, int *signp, mp_rnd_t rnd_mode)
+inline const mpreal lgamma (const mpreal& v, int *signp = 0, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal x(v);
int tsignp;
- if(signp) mpfr_lgamma(x.mp,signp,v.mp,rnd_mode);
- else mpfr_lgamma(x.mp,&tsignp,v.mp,rnd_mode);
+ if(signp) mpfr_lgamma(x.mpfr_ptr(), signp,v.mpfr_srcptr(),rnd_mode);
+ else mpfr_lgamma(x.mpfr_ptr(),&tsignp,v.mpfr_srcptr(),rnd_mode);
return x;
}
-inline const mpreal besseljn (long n, const mpreal& x, mp_rnd_t r)
+inline const mpreal besseljn (long n, const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd())
{
- mpreal y(0, mpfr_get_prec(x.mpfr_srcptr()));
+ mpreal y(0, x.getPrecision());
mpfr_jn(y.mpfr_ptr(), n, x.mpfr_srcptr(), r);
return y;
}
-inline const mpreal besselyn (long n, const mpreal& x, mp_rnd_t r)
+inline const mpreal besselyn (long n, const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd())
{
- mpreal y(0, mpfr_get_prec(x.mpfr_srcptr()));
+ mpreal y(0, x.getPrecision());
mpfr_yn(y.mpfr_ptr(), n, x.mpfr_srcptr(), r);
return y;
}
-inline const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode)
+inline const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal a;
mp_prec_t p1, p2, p3;
@@ -2217,7 +2316,7 @@ inline const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, m return a;
}
-inline const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode)
+inline const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal a;
mp_prec_t p1, p2, p3;
@@ -2232,7 +2331,7 @@ inline const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, m return a;
}
-inline const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode)
+inline const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal a;
mp_prec_t p1, p2;
@@ -2247,7 +2346,7 @@ inline const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode) return a;
}
-inline const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode)
+inline const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal x;
mpfr_ptr* t;
@@ -2264,23 +2363,23 @@ inline const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_m // MPFR 2.4.0 Specifics
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
-inline int sinh_cosh(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode)
+inline int sinh_cosh(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
return mpfr_sinh_cosh(s.mp,c.mp,v.mp,rnd_mode);
}
-inline const mpreal li2 (const mpreal& x, mp_rnd_t r)
+inline const mpreal li2 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd())
{
MPREAL_UNARY_MATH_FUNCTION_BODY(li2);
}
-inline const mpreal rem (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
+inline const mpreal rem (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
/* R = rem(X,Y) if Y != 0, returns X - n * Y where n = trunc(X/Y). */
return fmod(x, y, rnd_mode);
}
-inline const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
+inline const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
(void)rnd_mode;
@@ -2305,7 +2404,7 @@ inline const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode) return m;
}
-inline const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
+inline const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal a;
mp_prec_t yp, xp;
@@ -2320,7 +2419,7 @@ inline const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode) return a;
}
-inline const mpreal rec_sqrt(const mpreal& v, mp_rnd_t rnd_mode)
+inline const mpreal rec_sqrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal x(v);
mpfr_rec_sqrt(x.mp,v.mp,rnd_mode);
@@ -2331,41 +2430,41 @@ inline const mpreal rec_sqrt(const mpreal& v, mp_rnd_t rnd_mode) //////////////////////////////////////////////////////////////////////////
// MPFR 3.0.0 Specifics
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
-inline const mpreal digamma (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(digamma); }
-inline const mpreal ai (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(ai); }
+inline const mpreal digamma (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(digamma); }
+inline const mpreal ai (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(ai); }
#endif // MPFR 3.0.0 Specifics
//////////////////////////////////////////////////////////////////////////
// Constants
-inline const mpreal const_log2 (mp_prec_t p, mp_rnd_t r)
+inline const mpreal const_log2 (mp_prec_t p = mpreal::get_default_prec(), mp_rnd_t r = mpreal::get_default_rnd())
{
mpreal x(0, p);
mpfr_const_log2(x.mpfr_ptr(), r);
return x;
}
-inline const mpreal const_pi (mp_prec_t p, mp_rnd_t r)
+inline const mpreal const_pi (mp_prec_t p = mpreal::get_default_prec(), mp_rnd_t r = mpreal::get_default_rnd())
{
mpreal x(0, p);
mpfr_const_pi(x.mpfr_ptr(), r);
return x;
}
-inline const mpreal const_euler (mp_prec_t p, mp_rnd_t r)
+inline const mpreal const_euler (mp_prec_t p = mpreal::get_default_prec(), mp_rnd_t r = mpreal::get_default_rnd())
{
mpreal x(0, p);
mpfr_const_euler(x.mpfr_ptr(), r);
return x;
}
-inline const mpreal const_catalan (mp_prec_t p, mp_rnd_t r)
+inline const mpreal const_catalan (mp_prec_t p = mpreal::get_default_prec(), mp_rnd_t r = mpreal::get_default_rnd())
{
mpreal x(0, p);
mpfr_const_catalan(x.mpfr_ptr(), r);
return x;
}
-inline const mpreal const_infinity (int sign, mp_prec_t p, mp_rnd_t /*r*/)
+inline const mpreal const_infinity (int sign = 1, mp_prec_t p = mpreal::get_default_prec())
{
mpreal x(0, p);
mpfr_set_inf(x.mpfr_ptr(), sign);
@@ -2402,12 +2501,12 @@ inline const mpreal trunc(const mpreal& v) return x;
}
-inline const mpreal rint (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint ); }
-inline const mpreal rint_ceil (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint_ceil ); }
-inline const mpreal rint_floor (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint_floor); }
-inline const mpreal rint_round (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint_round); }
-inline const mpreal rint_trunc (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint_trunc); }
-inline const mpreal frac (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(frac ); }
+inline const mpreal rint (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint ); }
+inline const mpreal rint_ceil (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint_ceil ); }
+inline const mpreal rint_floor (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint_floor); }
+inline const mpreal rint_round (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint_round); }
+inline const mpreal rint_trunc (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint_trunc); }
+inline const mpreal frac (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(frac ); }
//////////////////////////////////////////////////////////////////////////
// Miscellaneous Functions
@@ -2415,14 +2514,14 @@ inline void swap (mpreal& a, mpreal& b) { mpfr_swap(a.mp,b inline const mpreal (max)(const mpreal& x, const mpreal& y){ return (x>y?x:y); }
inline const mpreal (min)(const mpreal& x, const mpreal& y){ return (x<y?x:y); }
-inline const mpreal fmax(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
+inline const mpreal fmax(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal a;
mpfr_max(a.mp,x.mp,y.mp,rnd_mode);
return a;
}
-inline const mpreal fmin(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
+inline const mpreal fmin(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal a;
mpfr_min(a.mp,x.mp,y.mp,rnd_mode);
@@ -2457,14 +2556,22 @@ inline const mpreal urandomb (gmp_randstate_t& state) return x;
}
-#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
+#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,1,0))
// use gmp_randinit_default() to init state, gmp_randclear() to clear
-inline const mpreal urandom (gmp_randstate_t& state, mp_rnd_t rnd_mode)
+inline const mpreal urandom (gmp_randstate_t& state, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal x;
mpfr_urandom(x.mp,state,rnd_mode);
return x;
}
+
+inline const mpreal grandom (gmp_randstate_t& state, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
+{
+ mpreal x;
+ mpfr_grandom(x.mp, NULL, state, rnd_mode);
+ return x;
+}
+
#endif
#if (MPFR_VERSION <= MPFR_VERSION_NUM(2,4,2))
@@ -2480,7 +2587,7 @@ inline const mpreal random2 (mp_size_t size, mp_exp_t exp) // a = random(seed); <- initialization & first random number generation
// a = random(); <- next random numbers generation
// seed != 0
-inline const mpreal random(unsigned int seed)
+inline const mpreal random(unsigned int seed = 0)
{
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
@@ -2504,6 +2611,25 @@ inline const mpreal random(unsigned int seed) }
+#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
+inline const mpreal grandom(unsigned int seed = 0)
+{
+ static gmp_randstate_t state;
+ static bool isFirstTime = true;
+
+ if(isFirstTime)
+ {
+ gmp_randinit_default(state);
+ gmp_randseed_ui(state,0);
+ isFirstTime = false;
+ }
+
+ if(seed != 0) gmp_randseed_ui(state,seed);
+
+ return mpfr::grandom(state);
+}
+#endif
+
//////////////////////////////////////////////////////////////////////////
// Set/Get global properties
inline void mpreal::set_default_prec(mp_prec_t prec)
@@ -2523,21 +2649,21 @@ inline bool mpreal::fits_in_bits(double x, int n) return IsInf(x) || (std::modf ( std::ldexp ( std::frexp ( x, &i ), n ), &t ) == 0.0);
}
-inline const mpreal pow(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode)
+inline const mpreal pow(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal x(a);
mpfr_pow(x.mp,x.mp,b.mp,rnd_mode);
return x;
}
-inline const mpreal pow(const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode)
+inline const mpreal pow(const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal x(a);
mpfr_pow_z(x.mp,x.mp,b,rnd_mode);
return x;
}
-inline const mpreal pow(const mpreal& a, const unsigned long int b, mp_rnd_t rnd_mode)
+inline const mpreal pow(const mpreal& a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal x(a);
mpfr_pow_ui(x.mp,x.mp,b,rnd_mode);
@@ -2549,7 +2675,7 @@ inline const mpreal pow(const mpreal& a, const unsigned int b, mp_rnd_t rnd_mode return pow(a,static_cast<unsigned long int>(b),rnd_mode);
}
-inline const mpreal pow(const mpreal& a, const long int b, mp_rnd_t rnd_mode)
+inline const mpreal pow(const mpreal& a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal x(a);
mpfr_pow_si(x.mp,x.mp,b,rnd_mode);
@@ -2571,7 +2697,7 @@ inline const mpreal pow(const mpreal& a, const double b, mp_rnd_t rnd_mode) return pow(a,mpreal(b),rnd_mode);
}
-inline const mpreal pow(const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode)
+inline const mpreal pow(const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal x(a);
mpfr_ui_pow(x.mp,a,b.mp,rnd_mode);
@@ -2586,13 +2712,13 @@ inline const mpreal pow(const unsigned int a, const mpreal& b, mp_rnd_t rnd_mode inline const mpreal pow(const long int a, const mpreal& b, mp_rnd_t rnd_mode)
{
if (a>=0) return pow(static_cast<unsigned long int>(a),b,rnd_mode);
- else return pow(mpreal(a),b,rnd_mode);
+ else return pow(mpreal(a),b,rnd_mode);
}
inline const mpreal pow(const int a, const mpreal& b, mp_rnd_t rnd_mode)
{
if (a>=0) return pow(static_cast<unsigned long int>(a),b,rnd_mode);
- else return pow(mpreal(a),b,rnd_mode);
+ else return pow(mpreal(a),b,rnd_mode);
}
inline const mpreal pow(const long double a, const mpreal& b, mp_rnd_t rnd_mode)
@@ -2621,13 +2747,13 @@ inline const mpreal pow(const unsigned long int a, const unsigned int b, mp_rnd_ inline const mpreal pow(const unsigned long int a, const long int b, mp_rnd_t rnd_mode)
{
if(b>0) return pow(a,static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui
- else return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow
+ else return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow
}
inline const mpreal pow(const unsigned long int a, const int b, mp_rnd_t rnd_mode)
{
if(b>0) return pow(a,static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui
- else return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow
+ else return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow
}
inline const mpreal pow(const unsigned long int a, const long double b, mp_rnd_t rnd_mode)
@@ -2824,7 +2950,7 @@ inline const mpreal pow(const double a, const int b, mp_rnd_t rnd_mode) // Non-throwing swap C++ idiom: http://en.wikibooks.org/wiki/More_C%2B%2B_Idioms/Non-throwing_swap
namespace std
{
- // only allowed to extend namespace std with specializations
+ // we are allowed to extend namespace std with specializations only
template <>
inline void swap(mpfr::mpreal& x, mpfr::mpreal& y)
{
@@ -2852,20 +2978,6 @@ namespace std static const bool tinyness_before = true;
static const float_denorm_style has_denorm = denorm_absent;
-
- inline static float_round_style round_style()
- {
- mp_rnd_t r = mpfr::mpreal::get_default_rnd();
-
- switch (r)
- {
- case MPFR_RNDN: return round_to_nearest;
- case MPFR_RNDZ: return round_toward_zero;
- case MPFR_RNDU: return round_toward_infinity;
- case MPFR_RNDD: return round_toward_neg_infinity;
- default: return round_indeterminate;
- }
- }
inline static mpfr::mpreal (min) (mp_prec_t precision = mpfr::mpreal::get_default_prec()) { return mpfr::minval(precision); }
inline static mpfr::mpreal (max) (mp_prec_t precision = mpfr::mpreal::get_default_prec()) { return mpfr::maxval(precision); }
@@ -2873,7 +2985,7 @@ namespace std // Returns smallest eps such that 1 + eps != 1 (classic machine epsilon)
inline static mpfr::mpreal epsilon(mp_prec_t precision = mpfr::mpreal::get_default_prec()) { return mpfr::machine_epsilon(precision); }
-
+
// Returns smallest eps such that x + eps != x (relative machine epsilon)
inline static mpfr::mpreal epsilon(const mpfr::mpreal& x) { return mpfr::machine_epsilon(x); }
@@ -2881,7 +2993,7 @@ namespace std {
mp_rnd_t r = mpfr::mpreal::get_default_rnd();
- if(r == MPFR_RNDN) return mpfr::mpreal(0.5, precision);
+ if(r == GMP_RNDN) return mpfr::mpreal(0.5, precision);
else return mpfr::mpreal(1.0, precision);
}
@@ -2896,10 +3008,30 @@ namespace std MPREAL_PERMISSIVE_EXPR static const int min_exponent10 = (int) (MPFR_EMIN_DEFAULT * 0.3010299956639811);
MPREAL_PERMISSIVE_EXPR static const int max_exponent10 = (int) (MPFR_EMAX_DEFAULT * 0.3010299956639811);
- // Should be constant according to standard, but 'digits' depends on precision in MPFR
+#ifdef MPREAL_HAVE_DYNAMIC_STD_NUMERIC_LIMITS
+
+ // Following members should be constant according to standard, but they can be variable in MPFR
+ // So we define them as functions here.
+ //
+ // This is preferable way for std::numeric_limits<mpfr::mpreal> specialization.
+ // But it is incompatible with standard std::numeric_limits and might not work with other libraries, e.g. boost.
+ // See below for compatible implementation.
+ inline static float_round_style round_style()
+ {
+ mp_rnd_t r = mpfr::mpreal::get_default_rnd();
+
+ switch (r)
+ {
+ case GMP_RNDN: return round_to_nearest;
+ case GMP_RNDZ: return round_toward_zero;
+ case GMP_RNDU: return round_toward_infinity;
+ case GMP_RNDD: return round_toward_neg_infinity;
+ default: return round_indeterminate;
+ }
+ }
- inline static int digits() { return mpfr::mpreal::get_default_prec(); }
- inline static int digits(const mpfr::mpreal& x) { return x.getPrecision(); }
+ inline static int digits() { return int(mpfr::mpreal::get_default_prec()); }
+ inline static int digits(const mpfr::mpreal& x) { return x.getPrecision(); }
inline static int digits10(mp_prec_t precision = mpfr::mpreal::get_default_prec())
{
@@ -2915,6 +3047,25 @@ namespace std {
return digits10(precision);
}
+#else
+ // Digits and round_style are NOT constants when it comes to mpreal.
+ // If possible, please use functions digits() and round_style() defined above.
+ //
+ // These (default) values are preserved for compatibility with existing libraries, e.g. boost.
+ // Change them accordingly to your application.
+ //
+ // For example, if you use 256 bits of precision uniformly in your program, then:
+ // digits = 256
+ // digits10 = 77
+ // max_digits10 = 78
+ //
+ // Approximate formula for decimal digits is: digits10 = floor(log10(2) * digits). See bits2digits() for more details.
+
+ static const std::float_round_style round_style = round_to_nearest;
+ static const int digits = 53;
+ static const int digits10 = 15;
+ static const int max_digits10 = 16;
+#endif
};
}
diff --git a/unsupported/test/polynomialsolver.cpp b/unsupported/test/polynomialsolver.cpp index 13f92169e..de79f1538 100644 --- a/unsupported/test/polynomialsolver.cpp +++ b/unsupported/test/polynomialsolver.cpp @@ -104,9 +104,7 @@ void evalSolverSugarFunction( const POLYNOMIAL& pols, const ROOTS& roots, const // 1) the roots found are correct // 2) the roots have distinct moduli - typedef typename POLYNOMIAL::Scalar Scalar; typedef typename REAL_ROOTS::Scalar Real; - typedef PolynomialSolver<Scalar, Deg > PolynomialSolverType; //Test realRoots std::vector< Real > calc_realRoots; |