diff --git a/Eigen/src/Core/DenseStorage.h b/Eigen/src/Core/DenseStorage.h index 3c02a1025..a8bb8a624 100644 --- a/Eigen/src/Core/DenseStorage.h +++ b/Eigen/src/Core/DenseStorage.h @@ -420,7 +420,7 @@ template class DenseStorage(m_data, m_rows*m_cols); - if (size) + if (size>0) // >0 and not simply !=0 to let the compiler knows that size cannot be negative m_data = internal::conditional_aligned_new_auto(size); else m_data = 0; @@ -497,7 +497,7 @@ template class DenseStorage(m_data, _Rows*m_cols); - if (size) + if (size>0) // >0 and not simply !=0 to let the compiler knows that size cannot be negative m_data = internal::conditional_aligned_new_auto(size); else m_data = 0; @@ -573,7 +573,7 @@ template class DenseStorage(m_data, _Cols*m_rows); - if (size) + if (size>0) // >0 and not simply !=0 to let the compiler knows that size cannot be negative m_data = internal::conditional_aligned_new_auto(size); else m_data = 0; diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h index 72116e144..34dd15d85 100644 --- a/Eigen/src/Core/MathFunctions.h +++ b/Eigen/src/Core/MathFunctions.h @@ -1217,7 +1217,8 @@ inline int log2(int x) /** \returns the square root of \a x. * - * It is essentially equivalent to \code using std::sqrt; return sqrt(x); \endcode, + * It is essentially equivalent to + * \code using std::sqrt; return sqrt(x); \endcode * but slightly faster for float/double and some compilers (e.g., gcc), thanks to * specializations when SSE is enabled. * diff --git a/Eigen/src/Core/arch/SSE/PacketMath.h b/Eigen/src/Core/arch/SSE/PacketMath.h index 2e815c0c5..99d55d5e9 100755 --- a/Eigen/src/Core/arch/SSE/PacketMath.h +++ b/Eigen/src/Core/arch/SSE/PacketMath.h @@ -28,7 +28,7 @@ namespace internal { #endif #endif -#if (defined EIGEN_VECTORIZE_AVX) && (EIGEN_COMP_GNUC_STRICT || EIGEN_COMP_MINGW) && (__GXX_ABI_VERSION < 1004) +#if ((defined EIGEN_VECTORIZE_AVX) && (EIGEN_COMP_GNUC_STRICT || EIGEN_COMP_MINGW) && (__GXX_ABI_VERSION < 1004)) || EIGEN_OS_QNX // With GCC's default ABI version, a __m128 or __m256 are the same types and therefore we cannot // have overloads for both types without linking error. // One solution is to increase ABI version using -fabi-version=4 (or greater). diff --git a/Eigen/src/Core/util/ConfigureVectorization.h b/Eigen/src/Core/util/ConfigureVectorization.h index e75c7d89e..a2743624e 100644 --- a/Eigen/src/Core/util/ConfigureVectorization.h +++ b/Eigen/src/Core/util/ConfigureVectorization.h @@ -379,10 +379,12 @@ #include #endif -#if defined(EIGEN_HIP_DEVICE_COMPILE) - +#if defined(EIGEN_HIPCC) #define EIGEN_VECTORIZE_GPU #include +#endif + +#if defined(EIGEN_HIP_DEVICE_COMPILE) #define EIGEN_HAS_HIP_FP16 #include diff --git a/Eigen/src/Core/util/IntegralConstant.h b/Eigen/src/Core/util/IntegralConstant.h index bf99cd8ab..caeea232d 100644 --- a/Eigen/src/Core/util/IntegralConstant.h +++ b/Eigen/src/Core/util/IntegralConstant.h @@ -55,7 +55,9 @@ public: operator int() const { return value; } FixedInt() {} FixedInt( VariableAndFixedInt other) { - EIGEN_ONLY_USED_FOR_DEBUG(other); + #ifndef EIGEN_INTERNAL_DEBUGGING + EIGEN_UNUSED_VARIABLE(other); + #endif eigen_internal_assert(int(other)==N); } diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h index 9dd2e0252..6664770f3 100644 --- a/Eigen/src/Core/util/Memory.h +++ b/Eigen/src/Core/util/Memory.h @@ -96,10 +96,16 @@ inline void throw_std_bad_alloc() /** \internal Like malloc, but the returned pointer is guaranteed to be 16-byte aligned. * Fast, but wastes 16 additional bytes of memory. Does not throw any exception. */ -inline void* handmade_aligned_malloc(std::size_t size, std::size_t alignment = EIGEN_DEFAULT_ALIGN_BYTES) +EIGEN_DEVICE_FUNC inline void* handmade_aligned_malloc(std::size_t size, std::size_t alignment = EIGEN_DEFAULT_ALIGN_BYTES) { eigen_assert(alignment >= sizeof(void*) && (alignment & -alignment) == alignment && "Alignment must be at least sizeof(void*) and a power of 2"); + +#if defined(EIGEN_HIP_DEVICE_COMPILE) + void *original = ::malloc(size+alignment); +#else void *original = std::malloc(size+alignment); +#endif + if (original == 0) return 0; void *aligned = reinterpret_cast((reinterpret_cast(original) & ~(std::size_t(alignment-1))) + alignment); *(reinterpret_cast(aligned) - 1) = original; @@ -107,9 +113,15 @@ inline void* handmade_aligned_malloc(std::size_t size, std::size_t alignment = E } /** \internal Frees memory allocated with handmade_aligned_malloc */ -inline void handmade_aligned_free(void *ptr) +EIGEN_DEVICE_FUNC inline void handmade_aligned_free(void *ptr) { - if (ptr) std::free(*(reinterpret_cast(ptr) - 1)); + if (ptr) { +#if defined(EIGEN_HIP_DEVICE_COMPILE) + ::free(*(reinterpret_cast(ptr) - 1)); +#else + std::free(*(reinterpret_cast(ptr) - 1)); +#endif + } } /** \internal @@ -872,6 +884,15 @@ public: ~aligned_allocator() {} + #if EIGEN_COMP_GNUC_STRICT && EIGEN_GNUC_AT_LEAST(7,0) + // In gcc std::allocator::max_size() is bugged making gcc triggers a warning: + // eigen/Eigen/src/Core/util/Memory.h:189:12: warning: argument 1 value '18446744073709551612' exceeds maximum object size 9223372036854775807 + // See https://gcc.gnu.org/bugzilla/show_bug.cgi?id=87544 + size_type max_size() const { + return (std::numeric_limits::max)()/sizeof(T); + } + #endif + pointer allocate(size_type num, const void* /*hint*/ = 0) { internal::check_size_for_overflow(num); diff --git a/Eigen/src/plugins/ReshapedMethods.h b/Eigen/src/plugins/ReshapedMethods.h index 2466d24cb..538636ab5 100644 --- a/Eigen/src/plugins/ReshapedMethods.h +++ b/Eigen/src/plugins/ReshapedMethods.h @@ -105,7 +105,7 @@ EIGEN_DEVICE_FUNC inline Reshaped::value, internal::get_compiletime_reshape_size::value, - Order==AutoOrder?Flags&RowMajorBit:Order> + (Order==AutoOrder?Flags&RowMajorBit:Order)> reshaped(NRowsType nRows, NColsType nCols) EIGEN_RESHAPED_METHOD_CONST { return Reshaped EIGEN_DEVICE_FUNC -inline Reshaped +inline Reshaped reshaped() EIGEN_RESHAPED_METHOD_CONST { EIGEN_STATIC_ASSERT(Order==RowMajor || Order==ColMajor || Order==AutoOrder, INVALID_TEMPLATE_PARAMETER); diff --git a/test/indexed_view.cpp b/test/indexed_view.cpp index 9a284cf0a..8ede612d1 100644 --- a/test/indexed_view.cpp +++ b/test/indexed_view.cpp @@ -15,6 +15,14 @@ #ifdef EIGEN_TEST_PART_3 // Make sure we also check c++98 max implementation #define EIGEN_MAX_CPP_VER 03 + +// We need to disable this warning when compiling with c++11 while limiting Eigen to c++98 +// Ideally we would rather configure the compiler to build in c++98 mode but this needs +// to be done at the CMakeLists.txt level. +#if defined(__GNUC__) && (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 8)) + #pragma GCC diagnostic ignored "-Wdeprecated" +#endif + #endif #include diff --git a/test/ref.cpp b/test/ref.cpp index a94ea9677..250135bdb 100644 --- a/test/ref.cpp +++ b/test/ref.cpp @@ -255,8 +255,8 @@ void test_ref_overloads() void test_ref_fixed_size_assert() { - Vector4f v4; - VectorXf vx(10); + Vector4f v4 = Vector4f::Random(); + VectorXf vx = VectorXf::Random(10); VERIFY_RAISES_STATIC_ASSERT( Ref y = v4; (void)y; ); VERIFY_RAISES_STATIC_ASSERT( Ref y = vx.head<4>(); (void)y; ); VERIFY_RAISES_STATIC_ASSERT( Ref y = v4; (void)y; ); diff --git a/test/stddeque.cpp b/test/stddeque.cpp index 0a6aa2572..ea85ea968 100644 --- a/test/stddeque.cpp +++ b/test/stddeque.cpp @@ -18,7 +18,7 @@ void check_stddeque_matrix(const MatrixType& m) Index rows = m.rows(); Index cols = m.cols(); MatrixType x = MatrixType::Random(rows,cols), y = MatrixType::Random(rows,cols); - std::deque > v(10, MatrixType(rows,cols)), w(20, y); + std::deque > v(10, MatrixType::Zero(rows,cols)), w(20, y); v.front() = x; w.front() = w.back(); VERIFY_IS_APPROX(w.front(), w.back()); @@ -33,7 +33,7 @@ void check_stddeque_matrix(const MatrixType& m) ++wi; } - v.resize(21); + v.resize(21,MatrixType::Zero(rows,cols)); v.back() = x; VERIFY_IS_APPROX(v.back(), x); v.resize(22,y); @@ -46,8 +46,8 @@ template void check_stddeque_transform(const TransformType&) { typedef typename TransformType::MatrixType MatrixType; - TransformType x(MatrixType::Random()), y(MatrixType::Random()); - std::deque > v(10), w(20, y); + TransformType x(MatrixType::Random()), y(MatrixType::Random()), ti=TransformType::Identity(); + std::deque > v(10,ti), w(20, y); v.front() = x; w.front() = w.back(); VERIFY_IS_APPROX(w.front(), w.back()); @@ -62,7 +62,7 @@ void check_stddeque_transform(const TransformType&) ++wi; } - v.resize(21); + v.resize(21,ti); v.back() = x; VERIFY_IS_APPROX(v.back(), x); v.resize(22,y); @@ -75,8 +75,8 @@ template void check_stddeque_quaternion(const QuaternionType&) { typedef typename QuaternionType::Coefficients Coefficients; - QuaternionType x(Coefficients::Random()), y(Coefficients::Random()); - std::deque > v(10), w(20, y); + QuaternionType x(Coefficients::Random()), y(Coefficients::Random()), qi=QuaternionType::Identity(); + std::deque > v(10,qi), w(20, y); v.front() = x; w.front() = w.back(); VERIFY_IS_APPROX(w.front(), w.back()); @@ -91,7 +91,7 @@ void check_stddeque_quaternion(const QuaternionType&) ++wi; } - v.resize(21); + v.resize(21,qi); v.back() = x; VERIFY_IS_APPROX(v.back(), x); v.resize(22,y); diff --git a/test/stddeque_overload.cpp b/test/stddeque_overload.cpp index eebe93d81..0f59f0695 100644 --- a/test/stddeque_overload.cpp +++ b/test/stddeque_overload.cpp @@ -31,7 +31,7 @@ void check_stddeque_matrix(const MatrixType& m) Index rows = m.rows(); Index cols = m.cols(); MatrixType x = MatrixType::Random(rows,cols), y = MatrixType::Random(rows,cols); - std::deque v(10, MatrixType(rows,cols)), w(20, y); + std::deque v(10, MatrixType::Zero(rows,cols)), w(20, y); v[5] = x; w[6] = v[5]; VERIFY_IS_APPROX(w[6], v[5]); @@ -64,8 +64,8 @@ template void check_stddeque_transform(const TransformType&) { typedef typename TransformType::MatrixType MatrixType; - TransformType x(MatrixType::Random()), y(MatrixType::Random()); - std::deque v(10), w(20, y); + TransformType x(MatrixType::Random()), y(MatrixType::Random()), ti=TransformType::Identity(); + std::deque v(10,ti), w(20, y); v[5] = x; w[6] = v[5]; VERIFY_IS_APPROX(w[6], v[5]); @@ -75,7 +75,7 @@ void check_stddeque_transform(const TransformType&) VERIFY_IS_APPROX(w[i], v[i]); } - v.resize(21); + v.resize(21,ti); v[20] = x; VERIFY_IS_APPROX(v[20], x); v.resize(22,y); @@ -98,8 +98,8 @@ template void check_stddeque_quaternion(const QuaternionType&) { typedef typename QuaternionType::Coefficients Coefficients; - QuaternionType x(Coefficients::Random()), y(Coefficients::Random()); - std::deque v(10), w(20, y); + QuaternionType x(Coefficients::Random()), y(Coefficients::Random()), qi=QuaternionType::Identity(); + std::deque v(10,qi), w(20, y); v[5] = x; w[6] = v[5]; VERIFY_IS_APPROX(w[6], v[5]); @@ -109,7 +109,7 @@ void check_stddeque_quaternion(const QuaternionType&) VERIFY_IS_APPROX(w[i], v[i]); } - v.resize(21); + v.resize(21,qi); v[20] = x; VERIFY_IS_APPROX(v[20], x); v.resize(22,y); diff --git a/test/stdlist.cpp b/test/stdlist.cpp index c9e5b7286..1af9e6ecb 100644 --- a/test/stdlist.cpp +++ b/test/stdlist.cpp @@ -18,7 +18,7 @@ void check_stdlist_matrix(const MatrixType& m) Index rows = m.rows(); Index cols = m.cols(); MatrixType x = MatrixType::Random(rows,cols), y = MatrixType::Random(rows,cols); - std::list > v(10, MatrixType(rows,cols)), w(20, y); + std::list > v(10, MatrixType::Zero(rows,cols)), w(20, y); v.front() = x; w.front() = w.back(); VERIFY_IS_APPROX(w.front(), w.back()); @@ -33,7 +33,7 @@ void check_stdlist_matrix(const MatrixType& m) ++wi; } - v.resize(21); + v.resize(21, MatrixType::Zero(rows,cols)); v.back() = x; VERIFY_IS_APPROX(v.back(), x); v.resize(22,y); @@ -46,8 +46,8 @@ template void check_stdlist_transform(const TransformType&) { typedef typename TransformType::MatrixType MatrixType; - TransformType x(MatrixType::Random()), y(MatrixType::Random()); - std::list > v(10), w(20, y); + TransformType x(MatrixType::Random()), y(MatrixType::Random()), ti=TransformType::Identity(); + std::list > v(10,ti), w(20, y); v.front() = x; w.front() = w.back(); VERIFY_IS_APPROX(w.front(), w.back()); @@ -62,7 +62,7 @@ void check_stdlist_transform(const TransformType&) ++wi; } - v.resize(21); + v.resize(21, ti); v.back() = x; VERIFY_IS_APPROX(v.back(), x); v.resize(22,y); @@ -75,8 +75,8 @@ template void check_stdlist_quaternion(const QuaternionType&) { typedef typename QuaternionType::Coefficients Coefficients; - QuaternionType x(Coefficients::Random()), y(Coefficients::Random()); - std::list > v(10), w(20, y); + QuaternionType x(Coefficients::Random()), y(Coefficients::Random()), qi=QuaternionType::Identity(); + std::list > v(10,qi), w(20, y); v.front() = x; w.front() = w.back(); VERIFY_IS_APPROX(w.front(), w.back()); @@ -91,7 +91,7 @@ void check_stdlist_quaternion(const QuaternionType&) ++wi; } - v.resize(21); + v.resize(21,qi); v.back() = x; VERIFY_IS_APPROX(v.back(), x); v.resize(22,y); diff --git a/test/stdlist_overload.cpp b/test/stdlist_overload.cpp index 93b6fc9ed..a78516e24 100644 --- a/test/stdlist_overload.cpp +++ b/test/stdlist_overload.cpp @@ -47,7 +47,7 @@ void check_stdlist_matrix(const MatrixType& m) Index rows = m.rows(); Index cols = m.cols(); MatrixType x = MatrixType::Random(rows,cols), y = MatrixType::Random(rows,cols); - std::list v(10, MatrixType(rows,cols)), w(20, y); + std::list v(10, MatrixType::Zero(rows,cols)), w(20, y); typename std::list::iterator itv = get(v, 5); typename std::list::iterator itw = get(w, 6); *itv = x; @@ -86,8 +86,8 @@ template void check_stdlist_transform(const TransformType&) { typedef typename TransformType::MatrixType MatrixType; - TransformType x(MatrixType::Random()), y(MatrixType::Random()); - std::list v(10), w(20, y); + TransformType x(MatrixType::Random()), y(MatrixType::Random()), ti=TransformType::Identity(); + std::list v(10,ti), w(20, y); typename std::list::iterator itv = get(v, 5); typename std::list::iterator itw = get(w, 6); *itv = x; @@ -103,7 +103,7 @@ void check_stdlist_transform(const TransformType&) ++itw; } - v.resize(21); + v.resize(21, ti); set(v, 20, x); VERIFY_IS_APPROX(*get(v, 20), x); v.resize(22,y); @@ -126,8 +126,8 @@ template void check_stdlist_quaternion(const QuaternionType&) { typedef typename QuaternionType::Coefficients Coefficients; - QuaternionType x(Coefficients::Random()), y(Coefficients::Random()); - std::list v(10), w(20, y); + QuaternionType x(Coefficients::Random()), y(Coefficients::Random()), qi=QuaternionType::Identity(); + std::list v(10,qi), w(20, y); typename std::list::iterator itv = get(v, 5); typename std::list::iterator itw = get(w, 6); *itv = x; @@ -143,7 +143,7 @@ void check_stdlist_quaternion(const QuaternionType&) ++itw; } - v.resize(21); + v.resize(21,qi); set(v, 20, x); VERIFY_IS_APPROX(*get(v, 20), x); v.resize(22,y); diff --git a/test/stdvector.cpp b/test/stdvector.cpp index e2b7bd061..18de240c6 100644 --- a/test/stdvector.cpp +++ b/test/stdvector.cpp @@ -17,7 +17,7 @@ void check_stdvector_matrix(const MatrixType& m) Index rows = m.rows(); Index cols = m.cols(); MatrixType x = MatrixType::Random(rows,cols), y = MatrixType::Random(rows,cols); - std::vector > v(10, MatrixType(rows,cols)), w(20, y); + std::vector > v(10, MatrixType::Zero(rows,cols)), w(20, y); v[5] = x; w[6] = v[5]; VERIFY_IS_APPROX(w[6], v[5]); @@ -86,8 +86,8 @@ template void check_stdvector_quaternion(const QuaternionType&) { typedef typename QuaternionType::Coefficients Coefficients; - QuaternionType x(Coefficients::Random()), y(Coefficients::Random()); - std::vector > v(10), w(20, y); + QuaternionType x(Coefficients::Random()), y(Coefficients::Random()), qi=QuaternionType::Identity(); + std::vector > v(10,qi), w(20, y); v[5] = x; w[6] = v[5]; VERIFY_IS_APPROX(w[6], v[5]); @@ -117,6 +117,16 @@ void check_stdvector_quaternion(const QuaternionType&) } } +// the code below triggered an invalid warning with gcc >= 7 +// eigen/Eigen/src/Core/util/Memory.h:189:12: warning: argument 1 value '18446744073709551612' exceeds maximum object size 9223372036854775807 +// This has been reported to gcc there: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=87544 +void std_vector_gcc_warning() +{ + typedef Eigen::Vector3f T; + std::vector > v; + v.push_back(T()); +} + EIGEN_DECLARE_TEST(stdvector) { // some non vectorizable fixed sizes diff --git a/test/stdvector_overload.cpp b/test/stdvector_overload.cpp index 5c042a64c..da04f8a84 100644 --- a/test/stdvector_overload.cpp +++ b/test/stdvector_overload.cpp @@ -31,7 +31,7 @@ void check_stdvector_matrix(const MatrixType& m) Index rows = m.rows(); Index cols = m.cols(); MatrixType x = MatrixType::Random(rows,cols), y = MatrixType::Random(rows,cols); - std::vector v(10, MatrixType(rows,cols)), w(20, y); + std::vector v(10, MatrixType::Zero(rows,cols)), w(20, y); v[5] = x; w[6] = v[5]; VERIFY_IS_APPROX(w[6], v[5]); @@ -100,8 +100,8 @@ template void check_stdvector_quaternion(const QuaternionType&) { typedef typename QuaternionType::Coefficients Coefficients; - QuaternionType x(Coefficients::Random()), y(Coefficients::Random()); - std::vector v(10), w(20, y); + QuaternionType x(Coefficients::Random()), y(Coefficients::Random()), qi=QuaternionType::Identity(); + std::vector v(10,qi), w(20, y); v[5] = x; w[6] = v[5]; VERIFY_IS_APPROX(w[6], v[5]); diff --git a/unsupported/Eigen/CXX11/src/Tensor/TensorContraction.h b/unsupported/Eigen/CXX11/src/Tensor/TensorContraction.h index 8f66e587a..61a4e1a3a 100644 --- a/unsupported/Eigen/CXX11/src/Tensor/TensorContraction.h +++ b/unsupported/Eigen/CXX11/src/Tensor/TensorContraction.h @@ -186,21 +186,21 @@ struct TensorContractionKernel { /*ConjugateLhs*/ false, /*ConjugateRhs*/ false> GebpKernel; - EIGEN_DONT_INLINE + EIGEN_DEVICE_FUNC EIGEN_DONT_INLINE static void packLhs(LhsScalar* lhsBlock, const typename LhsMapper::SubMapper& data_mapper, const StorageIndex depth, const StorageIndex rows) { LhsPacker()(lhsBlock, data_mapper, depth, rows, /*stride*/ 0, /*offset*/ 0); } - EIGEN_DONT_INLINE + EIGEN_DEVICE_FUNC EIGEN_DONT_INLINE static void packRhs(RhsScalar* rhsBlock, const typename RhsMapper::SubMapper& data_mapper, const StorageIndex depth, const StorageIndex cols) { RhsPacker()(rhsBlock, data_mapper, depth, cols); } - EIGEN_DONT_INLINE + EIGEN_DEVICE_FUNC EIGEN_DONT_INLINE static void invoke(const OutputMapper& output_mapper, const LhsScalar* lhsBlock, const RhsScalar* rhsBlock, const StorageIndex rows, const StorageIndex depth, @@ -667,8 +667,8 @@ struct TensorContractionEvaluatorBase this->m_device.memset(buffer, 0, m * n * sizeof(Scalar)); this->template evalGemmPartial(buffer, - 0, k, 1); + rhs_inner_dim_reordered, + Alignment, true>(buffer, 0, k, 1); } template + template EIGEN_DEVICE_FUNC void evalGemmPartial(Scalar* buffer, Index k_start, Index k_end, int num_threads) const { eigen_assert(k_end >= k_start && k_start >= 0 && k_end <= this->m_k_size); // columns in slice on left side, rows on right side diff --git a/unsupported/Eigen/CXX11/src/Tensor/TensorContractionThreadPool.h b/unsupported/Eigen/CXX11/src/Tensor/TensorContractionThreadPool.h index 42c6ae528..24ba3e431 100644 --- a/unsupported/Eigen/CXX11/src/Tensor/TensorContractionThreadPool.h +++ b/unsupported/Eigen/CXX11/src/Tensor/TensorContractionThreadPool.h @@ -794,7 +794,7 @@ struct TensorEvaluator(k, block_size); // we use 'result' for the first block's partial result. MaxSizeVector block_buffers(num_blocks - 1); - Barrier barrier(num_blocks); + Barrier barrier(internal::convert_index(num_blocks)); auto process_block = [=, &barrier](Scalar* buf, Index begin, Index end) { ::memset(buf, 0, m * n * sizeof(Scalar)); TENSOR_CONTRACTION_DISPATCH( diff --git a/unsupported/Eigen/CXX11/src/Tensor/TensorMorphing.h b/unsupported/Eigen/CXX11/src/Tensor/TensorMorphing.h index 213379dbd..69a59906b 100644 --- a/unsupported/Eigen/CXX11/src/Tensor/TensorMorphing.h +++ b/unsupported/Eigen/CXX11/src/Tensor/TensorMorphing.h @@ -195,6 +195,14 @@ struct TensorEvaluator, Device> m_impl.getResourceRequirements(resources); } + // required in block(OutputTensorBlock* output_block) const + // For C++03 compatibility this must be defined outside the method + struct BlockIteratorState { + Index stride; + Index span; + Index size; + Index count; + }; // TODO(andydavis) Reduce the overhead of this function. EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void block( OutputTensorBlock* output_block) const { @@ -219,12 +227,6 @@ struct TensorEvaluator, Device> } // Initialize output block iterator state. - struct BlockIteratorState { - Index stride; - Index span; - Index size; - Index count; - }; array block_iter_state; for (Index i = 0; i < NumOutputDims; ++i) { diff --git a/unsupported/Eigen/CXX11/src/Tensor/TensorReduction.h b/unsupported/Eigen/CXX11/src/Tensor/TensorReduction.h index 949764f3a..2c69e4fd4 100644 --- a/unsupported/Eigen/CXX11/src/Tensor/TensorReduction.h +++ b/unsupported/Eigen/CXX11/src/Tensor/TensorReduction.h @@ -218,6 +218,7 @@ struct InnerMostDimReducer { } }; +#if !defined(EIGEN_HIPCC) template struct InnerMostDimReducer { static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename Self::CoeffReturnType @@ -257,7 +258,8 @@ struct InnerMostDimReducer { } } }; - +#endif + template struct InnerMostDimPreserver { static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void reduce(const Self&, typename Self::Index, Op&, typename Self::PacketReturnType*) { diff --git a/unsupported/Eigen/CXX11/src/Tensor/TensorReductionGpu.h b/unsupported/Eigen/CXX11/src/Tensor/TensorReductionGpu.h index 88940e6e6..375c570b3 100644 --- a/unsupported/Eigen/CXX11/src/Tensor/TensorReductionGpu.h +++ b/unsupported/Eigen/CXX11/src/Tensor/TensorReductionGpu.h @@ -292,7 +292,7 @@ __global__ void FullReductionKernelHalfFloat(Reducer reducer, const Self input, } template -__global__ void ReductionCleanupKernelHalfFloat(Op& reducer, half* output, half2* scratch) { +__global__ void ReductionCleanupKernelHalfFloat(Op reducer, half* output, half2* scratch) { eigen_assert(threadIdx.x == 1); half tmp = __low2half(*scratch); reducer.reduce(__high2half(*scratch), &tmp); diff --git a/unsupported/Eigen/CXX11/src/Tensor/TensorScan.h b/unsupported/Eigen/CXX11/src/Tensor/TensorScan.h index 641366d9d..64f10d0a4 100644 --- a/unsupported/Eigen/CXX11/src/Tensor/TensorScan.h +++ b/unsupported/Eigen/CXX11/src/Tensor/TensorScan.h @@ -124,7 +124,11 @@ struct TensorEvaluator, Device> { m_stride = m_stride * dims[i]; } } else { - for (int i = NumDims - 1; i > op.axis(); --i) { + // dims can only be indexed through unsigned integers, + // so let's use an unsigned type to let the compiler knows. + // This prevents stupid warnings: ""'*((void*)(& evaluator)+64)[18446744073709551615]' may be used uninitialized in this function" + unsigned int axis = internal::convert_index(op.axis()); + for (unsigned int i = NumDims - 1; i > axis; --i) { m_stride = m_stride * dims[i]; } } diff --git a/unsupported/test/cxx11_tensor_reduction.cpp b/unsupported/test/cxx11_tensor_reduction.cpp index 9458f4e4e..c0d4c786f 100644 --- a/unsupported/test/cxx11_tensor_reduction.cpp +++ b/unsupported/test/cxx11_tensor_reduction.cpp @@ -225,11 +225,11 @@ static void test_simple_reductions() { Tensor ints(10); std::iota(ints.data(), ints.data() + ints.dimension(0), 0); - TensorFixedSize > all; - all = ints.all(); - VERIFY(!all()); - all = (ints >= ints.constant(0)).all(); - VERIFY(all()); + TensorFixedSize > all_; + all_ = ints.all(); + VERIFY(!all_()); + all_ = (ints >= ints.constant(0)).all(); + VERIFY(all_()); TensorFixedSize > any; any = (ints > ints.constant(10)).any(); diff --git a/unsupported/test/mpreal/mpreal.h b/unsupported/test/mpreal/mpreal.h index 8404f1ff8..d68d7a78f 100644 --- a/unsupported/test/mpreal/mpreal.h +++ b/unsupported/test/mpreal/mpreal.h @@ -1,34 +1,34 @@ /* - MPFR C++: 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-2015 Pavel Holoborodko + Copyright (c) 2008-2016 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, + 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, Petr Aleksandrov, Orion Poplawski, Charles Karney, Arash Partow, - Rodney James, Jorge Leitao. + Rodney James, Jorge Leitao, Jerome Benoit. 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 + + (B) Non-free licenses may also be purchased from the author, for users who do not want their programs protected by the GPL. - 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 + 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). Such users can purchase an unlimited-use license from the author. Contact us for more details. - + GNU General Public License ("GPL") copyright permissions statement: ************************************************************************** This program is free software: you can redistribute it and/or modify @@ -58,6 +58,7 @@ #include #include #include +#include // Options #define MPREAL_HAVE_MSVC_DEBUGVIEW // Enable Debugger Visualizer for "Debug" builds in MSVC. @@ -68,16 +69,18 @@ // Library version #define MPREAL_VERSION_MAJOR 3 #define MPREAL_VERSION_MINOR 6 -#define MPREAL_VERSION_PATCHLEVEL 2 -#define MPREAL_VERSION_STRING "3.6.2" +#define MPREAL_VERSION_PATCHLEVEL 5 +#define MPREAL_VERSION_STRING "3.6.5" // Detect compiler using signatures from http://predef.sourceforge.net/ -#if defined(__GNUC__) - #define IsInf(x) (isinf)(x) // GNU C++/Intel ICC compiler on Linux -#elif defined(_MSC_VER) // Microsoft Visual C++ - #define IsInf(x) (!_finite(x)) +#if defined(__GNUC__) && defined(__INTEL_COMPILER) + #define IsInf(x) isinf EIGEN_NOT_A_MACRO (x) // Intel ICC compiler on Linux + +#elif defined(_MSC_VER) // Microsoft Visual C++ + #define IsInf(x) (!_finite(x)) + #else - #define IsInf(x) (std::isinf)(x) // GNU C/C++ (and/or other compilers), just hope for C99 conformance + #define IsInf(x) std::isinf EIGEN_NOT_A_MACRO (x) // GNU C/C++ (and/or other compilers), just hope for C99 conformance #endif // A Clang feature extension to determine compiler features. @@ -85,22 +88,26 @@ #define __has_feature(x) 0 #endif -// Detect support for r-value references (move semantic). Borrowed from Eigen. +// Detect support for r-value references (move semantic). +// Move semantic should be enabled with great care in multi-threading environments, +// especially if MPFR uses custom memory allocators. +// Everything should be thread-safe and support passing ownership over thread boundary. #if (__has_feature(cxx_rvalue_references) || \ defined(__GXX_EXPERIMENTAL_CXX0X__) || __cplusplus >= 201103L || \ - (defined(_MSC_VER) && _MSC_VER >= 1600)) + (defined(_MSC_VER) && _MSC_VER >= 1600) && !defined(MPREAL_DISABLE_MOVE_SEMANTIC)) #define MPREAL_HAVE_MOVE_SUPPORT - // Use fields in mpfr_t structure to check if it was initialized / set dummy initialization + // 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. +// Detect support for explicit converters. #if (__has_feature(cxx_explicit_conversions) || \ - (defined(__GXX_EXPERIMENTAL_CXX0X__) && __GNUC_MINOR__ >= 5) || __cplusplus >= 201103L || \ - (defined(_MSC_VER) && _MSC_VER >= 1800)) + (defined(__GXX_EXPERIMENTAL_CXX0X__) && __GNUC_MINOR >= 5) || __cplusplus >= 201103L || \ + (defined(_MSC_VER) && _MSC_VER >= 1800) || \ + (defined(__INTEL_COMPILER) && __INTEL_COMPILER >= 1300)) #define MPREAL_HAVE_EXPLICIT_CONVERTERS #endif @@ -111,8 +118,8 @@ #define MPREAL_MSVC_DEBUGVIEW_CODE DebugView = toString(); #define MPREAL_MSVC_DEBUGVIEW_DATA std::string DebugView; #else - #define MPREAL_MSVC_DEBUGVIEW_CODE - #define MPREAL_MSVC_DEBUGVIEW_DATA + #define MPREAL_MSVC_DEBUGVIEW_CODE + #define MPREAL_MSVC_DEBUGVIEW_DATA #endif #include @@ -122,12 +129,12 @@ #endif // Less important options -#define MPREAL_DOUBLE_BITS_OVERFLOW -1 // Triggers overflow exception during conversion to double if mpreal +#define MPREAL_DOUBLE_BITS_OVERFLOW -1 // Triggers overflow exception during conversion to double if mpreal // cannot fit in MPREAL_DOUBLE_BITS_OVERFLOW bits // = -1 disables overflow checks (default) // Fast replacement for mpfr_set_zero(x, +1): -// (a) uses low-level data members, might not be compatible with new versions of MPFR +// (a) uses low-level data members, might not be forward compatible // (b) sign is not set, add (x)->_mpfr_sign = 1; #define mpfr_set_zero_fast(x) ((x)->_mpfr_exp = __MPFR_EXP_ZERO) @@ -142,19 +149,19 @@ namespace mpfr { class mpreal { private: mpfr_t mp; - + public: - + // Get default rounding mode & precision inline static mp_rnd_t get_default_rnd() { return (mp_rnd_t)(mpfr_get_default_rounding_mode()); } - inline static mp_prec_t get_default_prec() { return mpfr_get_default_prec(); } + inline static mp_prec_t get_default_prec() { return (mpfr_get_default_prec)(); } // Constructors && type conversions mpreal(); mpreal(const mpreal& 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()); + 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()); mpreal(const double u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); mpreal(const long double u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); mpreal(const unsigned long long int u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); @@ -163,15 +170,15 @@ 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); + // 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); mpreal(const char* s, mp_prec_t prec = mpreal::get_default_prec(), int base = 10, mp_rnd_t mode = mpreal::get_default_rnd()); mpreal(const std::string& s, mp_prec_t prec = mpreal::get_default_prec(), int base = 10, mp_rnd_t mode = mpreal::get_default_rnd()); - ~mpreal(); + ~mpreal(); #ifdef MPREAL_HAVE_MOVE_SUPPORT mpreal& operator=(mpreal&& v); @@ -180,7 +187,7 @@ public: // Operations // = - // +, -, *, /, ++, --, <<, >> + // +, -, *, /, ++, --, <<, >> // *=, +=, -=, /=, // <, >, ==, <=, >= @@ -190,7 +197,7 @@ public: mpreal& operator=(const mpz_t v); mpreal& operator=(const mpq_t v); mpreal& operator=(const long double v); - mpreal& operator=(const double v); + mpreal& operator=(const double v); mpreal& operator=(const unsigned long int v); mpreal& operator=(const unsigned long long int v); mpreal& operator=(const long long int v); @@ -224,7 +231,7 @@ public: const mpreal operator+() const; mpreal& operator++ (); - const mpreal operator++ (int); + const mpreal operator++ (int); // - mpreal& operator-=(const mpreal& v); @@ -242,7 +249,7 @@ public: friend const mpreal operator-(const long int b, const mpreal& a); friend const mpreal operator-(const int b, const mpreal& a); friend const mpreal operator-(const double b, const mpreal& a); - mpreal& operator-- (); + mpreal& operator-- (); const mpreal operator-- (int); // * @@ -255,7 +262,7 @@ public: mpreal& operator*=(const unsigned int v); mpreal& operator*=(const long int v); mpreal& operator*=(const int v); - + // / mpreal& operator/=(const mpreal& v); mpreal& operator/=(const mpz_t v); @@ -343,17 +350,19 @@ public: 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); friend const mpreal log2 (const mpreal& v, mp_rnd_t rnd_mode); friend const mpreal logb (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 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 nextpow2(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); @@ -395,17 +404,17 @@ public: 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 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 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[], const unsigned long int n, int& status, mp_rnd_t rnd_mode); - friend int sgn(const mpreal& v); // returns -1 or +1 + friend int sgn (const mpreal& v); // MPFR 2.4.0 Specifics #if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) @@ -438,7 +447,7 @@ public: // Splits mpreal value into fractional and integer parts. // Returns fractional part and stores integer part in n. - friend const mpreal modf(const mpreal& v, mpreal& n); + friend const mpreal modf(const mpreal& v, mpreal& n); // Constants // don't forget to call mpfr_free_cache() for every thread where you are using const-functions @@ -467,14 +476,14 @@ public: 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); friend const mpreal nextabove (const mpreal& x); friend const mpreal nextbelow (const mpreal& x); // use gmp_randinit_default() to init state, gmp_randclear() to clear - friend const mpreal urandomb (gmp_randstate_t& state); + friend const mpreal urandomb (gmp_randstate_t& state); // MPFR < 2.4.2 Specifics #if (MPFR_VERSION <= MPFR_VERSION_NUM(2,4,2)) @@ -482,7 +491,7 @@ public: #endif // Instance Checkers - friend bool (isnan) (const mpreal& v); + friend bool isnan EIGEN_NOT_A_MACRO (const mpreal& v); friend bool (isinf) (const mpreal& v); friend bool (isfinite) (const mpreal& v); @@ -501,15 +510,15 @@ public: // Aliases for get_prec(), set_prec() - needed for compatibility with std::complex interface inline mpreal& setPrecision(int Precision, mp_rnd_t RoundingMode = get_default_rnd()); inline int getPrecision() const; - + // Set mpreal to +/- inf, NaN, +/-0 - mpreal& setInf (int Sign = +1); + mpreal& setInf (int Sign = +1); mpreal& setNan (); mpreal& setZero (int Sign = +1); mpreal& setSign (int Sign, mp_rnd_t RoundingMode = get_default_rnd()); //Exponent - mp_exp_t get_exp(); + mp_exp_t get_exp() const; int set_exp(mp_exp_t e); int check_range (int t, mp_rnd_t rnd_mode = get_default_rnd()); int subnormalize (int t, mp_rnd_t rnd_mode = get_default_rnd()); @@ -532,7 +541,7 @@ 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); friend const mpreal fmin(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); @@ -542,7 +551,7 @@ private: // // mpfr::mpreal= ; Show value only // mpfr::mpreal=, bits ; Show value & precision - // + // // at the beginning of // [Visual Studio Installation Folder]\Common7\Packages\Debugger\autoexp.dat MPREAL_MSVC_DEBUGVIEW_DATA @@ -561,15 +570,15 @@ public: ////////////////////////////////////////////////////////////////////////// // Constructors & converters // Default constructor: creates mp number and initializes it to 0. -inline mpreal::mpreal() -{ - mpfr_init2(mpfr_ptr(), mpreal::get_default_prec()); +inline mpreal::mpreal() +{ + mpfr_init2(mpfr_ptr(), mpreal::get_default_prec()); mpfr_set_zero_fast(mpfr_ptr()); MPREAL_MSVC_DEBUGVIEW_CODE; } -inline mpreal::mpreal(const mpreal& u) +inline mpreal::mpreal(const mpreal& u) { mpfr_init2(mpfr_ptr(),mpfr_get_prec(u.mpfr_srcptr())); mpfr_set (mpfr_ptr(),u.mpfr_srcptr(),mpreal::get_default_rnd()); @@ -580,7 +589,7 @@ inline mpreal::mpreal(const mpreal& u) #ifdef MPREAL_HAVE_MOVE_SUPPORT inline mpreal::mpreal(mpreal&& other) { - mpfr_set_uninitialized(mpfr_ptr()); // make sure "other" holds no pointer to actual data + mpfr_set_uninitialized(mpfr_ptr()); // make sure "other" holds null-pointer (in uninitialized state) mpfr_swap(mpfr_ptr(), other.mpfr_ptr()); MPREAL_MSVC_DEBUGVIEW_CODE; @@ -588,9 +597,11 @@ inline mpreal::mpreal(mpreal&& other) inline mpreal& mpreal::operator=(mpreal&& other) { - mpfr_swap(mpfr_ptr(), other.mpfr_ptr()); - - MPREAL_MSVC_DEBUGVIEW_CODE; + if (this != &other) + { + mpfr_swap(mpfr_ptr(), other.mpfr_ptr()); // destructor for "other" will be called just afterwards + MPREAL_MSVC_DEBUGVIEW_CODE; + } return *this; } #endif @@ -639,20 +650,20 @@ inline mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode) mpfr_init2(mpfr_ptr(), prec); #if (MPREAL_DOUBLE_BITS_OVERFLOW > -1) - if(fits_in_bits(u, MPREAL_DOUBLE_BITS_OVERFLOW)) - { - mpfr_set_d(mpfr_ptr(), 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(mpfr_ptr(), u, mode); + mpfr_set_d(mpfr_ptr(), u, mode); #endif MPREAL_MSVC_DEBUGVIEW_CODE; } inline mpreal::mpreal(const long double u, mp_prec_t prec, mp_rnd_t mode) -{ +{ mpfr_init2 (mpfr_ptr(), prec); mpfr_set_ld(mpfr_ptr(), u, mode); @@ -660,7 +671,7 @@ inline mpreal::mpreal(const long double u, mp_prec_t prec, mp_rnd_t mode) } inline mpreal::mpreal(const unsigned long long int u, mp_prec_t prec, mp_rnd_t mode) -{ +{ mpfr_init2 (mpfr_ptr(), prec); mpfr_set_uj(mpfr_ptr(), u, mode); @@ -668,7 +679,7 @@ inline mpreal::mpreal(const unsigned long long int u, mp_prec_t prec, mp_rnd_t m } inline mpreal::mpreal(const long long int u, mp_prec_t prec, mp_rnd_t mode) -{ +{ mpfr_init2 (mpfr_ptr(), prec); mpfr_set_sj(mpfr_ptr(), u, mode); @@ -676,7 +687,7 @@ inline mpreal::mpreal(const long long int u, mp_prec_t prec, mp_rnd_t mode) } inline mpreal::mpreal(const unsigned long int u, mp_prec_t prec, mp_rnd_t mode) -{ +{ mpfr_init2 (mpfr_ptr(), prec); mpfr_set_ui(mpfr_ptr(), u, mode); @@ -684,7 +695,7 @@ inline mpreal::mpreal(const unsigned long int u, mp_prec_t prec, mp_rnd_t mode) } inline mpreal::mpreal(const unsigned int u, mp_prec_t prec, mp_rnd_t mode) -{ +{ mpfr_init2 (mpfr_ptr(), prec); mpfr_set_ui(mpfr_ptr(), u, mode); @@ -692,7 +703,7 @@ inline mpreal::mpreal(const unsigned int u, mp_prec_t prec, mp_rnd_t mode) } inline mpreal::mpreal(const long int u, mp_prec_t prec, mp_rnd_t mode) -{ +{ mpfr_init2 (mpfr_ptr(), prec); mpfr_set_si(mpfr_ptr(), u, mode); @@ -700,7 +711,7 @@ inline mpreal::mpreal(const long int u, mp_prec_t prec, mp_rnd_t mode) } inline mpreal::mpreal(const int u, mp_prec_t prec, mp_rnd_t mode) -{ +{ mpfr_init2 (mpfr_ptr(), prec); mpfr_set_si(mpfr_ptr(), u, mode); @@ -710,7 +721,7 @@ inline mpreal::mpreal(const int 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 (mpfr_ptr(), prec); - mpfr_set_str(mpfr_ptr(), s, base, mode); + mpfr_set_str(mpfr_ptr(), s, base, mode); MPREAL_MSVC_DEBUGVIEW_CODE; } @@ -718,7 +729,7 @@ inline mpreal::mpreal(const char* s, mp_prec_t prec, int base, mp_rnd_t mode) inline mpreal::mpreal(const std::string& s, mp_prec_t prec, int base, mp_rnd_t mode) { mpfr_init2 (mpfr_ptr(), prec); - mpfr_set_str(mpfr_ptr(), s.c_str(), base, mode); + mpfr_set_str(mpfr_ptr(), s.c_str(), base, mode); MPREAL_MSVC_DEBUGVIEW_CODE; } @@ -726,15 +737,15 @@ inline mpreal::mpreal(const std::string& s, mp_prec_t prec, int base, mp_rnd_t m inline void mpreal::clear(::mpfr_ptr x) { #ifdef MPREAL_HAVE_MOVE_SUPPORT - if(mpfr_is_initialized(x)) + if(mpfr_is_initialized(x)) #endif mpfr_clear(x); } -inline mpreal::~mpreal() -{ +inline mpreal::~mpreal() +{ clear(mpfr_ptr()); -} +} // internal namespace needed for template magic namespace internal{ @@ -742,55 +753,55 @@ namespace internal{ // Use SFINAE to restrict arithmetic operations instantiation only for numeric types // This is needed for smooth integration with libraries based on expression templates, like Eigen. // TODO: Do the same for boolean operators. - template struct result_type {}; - - template <> struct result_type {typedef mpreal type;}; - template <> struct result_type {typedef mpreal type;}; - template <> struct result_type {typedef mpreal type;}; - template <> struct result_type {typedef mpreal type;}; - template <> struct result_type {typedef mpreal type;}; - template <> struct result_type {typedef mpreal type;}; - template <> struct result_type {typedef mpreal type;}; - template <> struct result_type {typedef mpreal type;}; - template <> struct result_type {typedef mpreal type;}; - template <> struct result_type {typedef mpreal type;}; - template <> struct result_type {typedef mpreal type;}; + template struct result_type {}; + + template <> struct result_type {typedef mpreal type;}; + template <> struct result_type {typedef mpreal type;}; + template <> struct result_type {typedef mpreal type;}; + template <> struct result_type {typedef mpreal type;}; + template <> struct result_type {typedef mpreal type;}; + template <> struct result_type {typedef mpreal type;}; + template <> struct result_type {typedef mpreal type;}; + template <> struct result_type {typedef mpreal type;}; + template <> struct result_type {typedef mpreal type;}; + template <> struct result_type {typedef mpreal type;}; + template <> struct result_type {typedef mpreal type;}; } // + Addition -template -inline const typename internal::result_type::type +template +inline const typename internal::result_type::type operator+(const mpreal& lhs, const Rhs& rhs){ return mpreal(lhs) += rhs; } -template -inline const typename internal::result_type::type - operator+(const Lhs& lhs, const mpreal& rhs){ return mpreal(rhs) += lhs; } +template +inline const typename internal::result_type::type + operator+(const Lhs& lhs, const mpreal& rhs){ return mpreal(rhs) += lhs; } // - Subtraction -template -inline const typename internal::result_type::type +template +inline const typename internal::result_type::type operator-(const mpreal& lhs, const Rhs& rhs){ return mpreal(lhs) -= rhs; } -template -inline const typename internal::result_type::type +template +inline const typename internal::result_type::type operator-(const Lhs& lhs, const mpreal& rhs){ return mpreal(lhs) -= rhs; } // * Multiplication -template -inline const typename internal::result_type::type +template +inline const typename internal::result_type::type operator*(const mpreal& lhs, const Rhs& rhs){ return mpreal(lhs) *= rhs; } -template -inline const typename internal::result_type::type - operator*(const Lhs& lhs, const mpreal& rhs){ return mpreal(rhs) *= lhs; } +template +inline const typename internal::result_type::type + operator*(const Lhs& lhs, const mpreal& rhs){ return mpreal(rhs) *= lhs; } // / Division -template -inline const typename internal::result_type::type +template +inline const typename internal::result_type::type operator/(const mpreal& lhs, const Rhs& rhs){ return mpreal(lhs) /= rhs; } -template -inline const typename internal::result_type::type +template +inline const typename internal::result_type::type operator/(const Lhs& lhs, const mpreal& rhs){ return mpreal(lhs) /= rhs; } ////////////////////////////////////////////////////////////////////////// @@ -840,17 +851,17 @@ const mpreal pow(const long int a, const double b, mp_rnd_t rnd_mode = mpreal::g const mpreal pow(const int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); const mpreal pow(const int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); const mpreal pow(const int a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const int a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); +const mpreal pow(const int a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); const mpreal pow(const int a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const int a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); +const mpreal pow(const int a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const long double a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); +const mpreal pow(const long double a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); const mpreal pow(const long double a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); const mpreal pow(const long double a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); const mpreal pow(const long double a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); const mpreal pow(const long double a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const double a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); +const mpreal pow(const double a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); const mpreal pow(const double a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); const mpreal pow(const double a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); const mpreal pow(const double a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); @@ -867,9 +878,9 @@ inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpr inline mpreal machine_epsilon(mp_prec_t prec = mpreal::get_default_prec()); // Returns smallest eps such that x + eps != x (relative machine epsilon) -inline mpreal machine_epsilon(const mpreal& x); +inline mpreal machine_epsilon(const mpreal& x); -// Gives max & min values for the required precision, +// Gives max & min values for the required precision, // minval is 'safe' meaning 1 / minval does not overflow // maxval is 'safe' meaning 1 / maxval does not underflow inline mpreal minval(mp_prec_t prec = mpreal::get_default_prec()); @@ -882,7 +893,7 @@ inline bool isEqualFuzzy(const mpreal& a, const mpreal& b, const mpreal& eps); inline bool isEqualFuzzy(const mpreal& a, const mpreal& b); // 'Bitwise' equality check -// maxUlps - a and b can be apart by maxUlps binary numbers. +// maxUlps - a and b can be apart by maxUlps binary numbers. inline bool isEqualUlps(const mpreal& a, const mpreal& b, int maxUlps); ////////////////////////////////////////////////////////////////////////// @@ -908,13 +919,13 @@ inline mpreal& mpreal::operator=(const mpreal& v) { if (this != &v) { - mp_prec_t tp = mpfr_get_prec( mpfr_srcptr()); - mp_prec_t vp = mpfr_get_prec(v.mpfr_srcptr()); + mp_prec_t tp = mpfr_get_prec( mpfr_srcptr()); + mp_prec_t vp = mpfr_get_prec(v.mpfr_srcptr()); - if(tp != vp){ - clear(mpfr_ptr()); - mpfr_init2(mpfr_ptr(), vp); - } + if(tp != vp){ + clear(mpfr_ptr()); + mpfr_init2(mpfr_ptr(), vp); + } mpfr_set(mpfr_ptr(), v.mpfr_srcptr(), mpreal::get_default_rnd()); @@ -926,7 +937,7 @@ inline mpreal& mpreal::operator=(const mpreal& v) inline mpreal& mpreal::operator=(const mpf_t v) { mpfr_set_f(mpfr_ptr(), v, mpreal::get_default_rnd()); - + MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } @@ -934,7 +945,7 @@ inline mpreal& mpreal::operator=(const mpf_t v) inline mpreal& mpreal::operator=(const mpz_t v) { mpfr_set_z(mpfr_ptr(), v, mpreal::get_default_rnd()); - + MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } @@ -947,73 +958,73 @@ inline mpreal& mpreal::operator=(const mpq_t v) return *this; } -inline mpreal& mpreal::operator=(const long double v) -{ +inline mpreal& mpreal::operator=(const long double v) +{ mpfr_set_ld(mpfr_ptr(), v, mpreal::get_default_rnd()); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } -inline mpreal& mpreal::operator=(const 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(mpfr_ptr(),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(mpfr_ptr(),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(mpfr_ptr(), v, mpreal::get_default_rnd()); +inline mpreal& mpreal::operator=(const unsigned long int v) +{ + mpfr_set_ui(mpfr_ptr(), v, mpreal::get_default_rnd()); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } -inline mpreal& mpreal::operator=(const unsigned int v) -{ - mpfr_set_ui(mpfr_ptr(), v, mpreal::get_default_rnd()); +inline mpreal& mpreal::operator=(const unsigned int v) +{ + mpfr_set_ui(mpfr_ptr(), v, mpreal::get_default_rnd()); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } -inline mpreal& mpreal::operator=(const unsigned long long int v) -{ - mpfr_set_uj(mpfr_ptr(), v, mpreal::get_default_rnd()); +inline mpreal& mpreal::operator=(const unsigned long long int v) +{ + mpfr_set_uj(mpfr_ptr(), v, mpreal::get_default_rnd()); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } -inline mpreal& mpreal::operator=(const long long int v) -{ - mpfr_set_sj(mpfr_ptr(), v, mpreal::get_default_rnd()); +inline mpreal& mpreal::operator=(const long long int v) +{ + mpfr_set_sj(mpfr_ptr(), v, mpreal::get_default_rnd()); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } -inline mpreal& mpreal::operator=(const long int v) -{ - mpfr_set_si(mpfr_ptr(), v, mpreal::get_default_rnd()); +inline mpreal& mpreal::operator=(const long int v) +{ + mpfr_set_si(mpfr_ptr(), v, mpreal::get_default_rnd()); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } inline mpreal& mpreal::operator=(const int v) -{ - mpfr_set_si(mpfr_ptr(), v, mpreal::get_default_rnd()); +{ + mpfr_set_si(mpfr_ptr(), v, mpreal::get_default_rnd()); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; @@ -1034,7 +1045,7 @@ inline mpreal& mpreal::operator=(const char* s) if(0 == mpfr_set_str(t, s, 10, mpreal::get_default_rnd())) { - mpfr_set(mpfr_ptr(), t, mpreal::get_default_rnd()); + mpfr_set(mpfr_ptr(), t, mpreal::get_default_rnd()); MPREAL_MSVC_DEBUGVIEW_CODE; } @@ -1057,7 +1068,7 @@ inline mpreal& mpreal::operator=(const std::string& s) if(0 == mpfr_set_str(t, s.c_str(), 10, mpreal::get_default_rnd())) { - mpfr_set(mpfr_ptr(), t, mpreal::get_default_rnd()); + mpfr_set(mpfr_ptr(), t, mpreal::get_default_rnd()); MPREAL_MSVC_DEBUGVIEW_CODE; } @@ -1065,7 +1076,7 @@ inline mpreal& mpreal::operator=(const std::string& s) return *this; } -template +template inline mpreal& mpreal::operator= (const std::complex& z) { return *this = z.real(); @@ -1103,9 +1114,9 @@ inline mpreal& mpreal::operator+=(const mpq_t u) inline mpreal& mpreal::operator+= (const long double u) { - *this += mpreal(u); + *this += mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; + return *this; } inline mpreal& mpreal::operator+= (const double u) @@ -1161,12 +1172,12 @@ 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++() +inline mpreal& mpreal::operator++() { return *this += 1; } @@ -1178,7 +1189,7 @@ inline const mpreal mpreal::operator++ (int) return x; } -inline mpreal& mpreal::operator--() +inline mpreal& mpreal::operator--() { return *this -= 1; } @@ -1215,9 +1226,9 @@ inline mpreal& mpreal::operator-=(const mpq_t v) inline mpreal& mpreal::operator-=(const long double v) { - *this -= mpreal(v); + *this -= mpreal(v); MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; + return *this; } inline mpreal& mpreal::operator-=(const double v) @@ -1225,7 +1236,7 @@ inline mpreal& mpreal::operator-=(const double v) #if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) mpfr_sub_d(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); #else - *this -= mpreal(v); + *this -= mpreal(v); #endif MPREAL_MSVC_DEBUGVIEW_CODE; @@ -1269,9 +1280,9 @@ inline const mpreal mpreal::operator-()const 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) @@ -1340,9 +1351,9 @@ inline mpreal& mpreal::operator*=(const mpq_t v) inline mpreal& mpreal::operator*=(const long double v) { - *this *= mpreal(v); + *this *= mpreal(v); MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; + return *this; } inline mpreal& mpreal::operator*=(const double v) @@ -1350,7 +1361,7 @@ inline mpreal& mpreal::operator*=(const double v) #if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) mpfr_mul_d(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); #else - *this *= mpreal(v); + *this *= mpreal(v); #endif MPREAL_MSVC_DEBUGVIEW_CODE; return *this; @@ -1386,9 +1397,9 @@ inline mpreal& mpreal::operator*=(const int v) 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; } ////////////////////////////////////////////////////////////////////////// @@ -1418,7 +1429,7 @@ inline mpreal& mpreal::operator/=(const long double v) { *this /= mpreal(v); MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; + return *this; } inline mpreal& mpreal::operator/=(const double v) @@ -1426,7 +1437,7 @@ inline mpreal& mpreal::operator/=(const double v) #if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) mpfr_div_d(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); #else - *this /= mpreal(v); + *this /= mpreal(v); #endif MPREAL_MSVC_DEBUGVIEW_CODE; return *this; @@ -1462,9 +1473,9 @@ inline mpreal& mpreal::operator/=(const int v) inline const mpreal operator/(const mpreal& a, const mpreal& b) { - 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; + 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) @@ -1639,9 +1650,9 @@ inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode) ////////////////////////////////////////////////////////////////////////// //Relational operators -// WARNING: +// WARNING: // -// Please note that following checks for double-NaN are guaranteed to work only in IEEE math mode: +// Please note that following checks for double-NaN are guaranteed to work only in IEEE math mode: // // isnan(b) = (b != b) // isnan(b) = !(b == b) (we use in code below) @@ -1659,7 +1670,7 @@ inline bool operator > (const mpreal& a, const double b ){ return ! 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 unsigned long int b ){ return !isnan EIGEN_NOT_A_MACRO (a) && (mpfr_cmp_ui(a.mpfr_srcptr(),b) >= 0 ); } -// inline bool operator >= (const mpreal& a, const unsigned int b ){ return !isnan EIGEN_NOT_A_MACRO (isnan()a) && (mpfr_cmp_ui(a.mpfr_srcptr(),b) >= 0 ); } +inline bool operator >= (const mpreal& a, const unsigned int b ){ return !isnan EIGEN_NOT_A_MACRO (a) && (mpfr_cmp_ui(a.mpfr_srcptr(),b) >= 0 ); } inline bool operator >= (const mpreal& a, const long int b ){ return !isnan EIGEN_NOT_A_MACRO (a) && (mpfr_cmp_si(a.mpfr_srcptr(),b) >= 0 ); } inline bool operator >= (const mpreal& a, const int b ){ return !isnan EIGEN_NOT_A_MACRO (a) && (mpfr_cmp_si(a.mpfr_srcptr(),b) >= 0 ); } inline bool operator >= (const mpreal& a, const long double b ){ return !isnan EIGEN_NOT_A_MACRO (a) && (b == b) && (mpfr_cmp_ld(a.mpfr_srcptr(),b) >= 0 ); } @@ -1697,7 +1708,7 @@ inline bool operator != (const mpreal& a, const int b ){ return ! inline bool operator != (const mpreal& a, const long double b ){ return !(a == b); } inline bool operator != (const mpreal& a, const double b ){ return !(a == b); } -inline bool (isnan) (const mpreal& op){ return (mpfr_nan_p (op.mpfr_srcptr()) != 0 ); } +inline bool isnan EIGEN_NOT_A_MACRO (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 ); } @@ -1705,7 +1716,7 @@ inline bool isint (const mpreal& op){ return (mpfr_integer_p(op.mpfr_srcpt #if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0)) inline bool isregular(const mpreal& op){ return (mpfr_regular_p(op.mpfr_srcptr()));} -#endif +#endif ////////////////////////////////////////////////////////////////////////// // Type Converters @@ -1762,21 +1773,21 @@ inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const std::ostringstream format; - int digits = (n >= 0) ? n : 1 + bits2digits(mpfr_get_prec(mpfr_srcptr())); - + int digits = (n >= 0) ? n : 2 + bits2digits(mpfr_get_prec(mpfr_srcptr())); + format << "%." << digits << "RNg"; return toString(format.str()); #else - char *s, *ns = NULL; + char *s, *ns = NULL; size_t slen, nslen; mp_exp_t exp; std::string out; if(mpfr_inf_p(mp)) - { + { if(mpfr_sgn(mp)>0) return "+Inf"; else return "-Inf"; } @@ -1791,7 +1802,7 @@ inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const { slen = strlen(s); nslen = strlen(ns); - if(nslen<=slen) + if(nslen<=slen) { mpfr_free_str(s); s = ns; @@ -1808,7 +1819,7 @@ inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const { // Remove zeros starting from right end char* ptr = s+slen-1; - while (*ptr=='0' && ptr>s+exp) ptr--; + while (*ptr=='0' && ptr>s+exp) ptr--; if(ptr==s+exp) out = std::string(s,exp+1); else out = std::string(s,exp+1)+'.'+std::string(s+exp+1,ptr-(s+exp+1)+1); @@ -1819,7 +1830,7 @@ inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const { // Remove zeros starting from right end char* ptr = s+slen-1; - while (*ptr=='0' && ptr>s+exp-1) ptr--; + while (*ptr=='0' && ptr>s+exp-1) ptr--; if(ptr==s+exp-1) out = std::string(s,exp); else out = std::string(s,exp)+'.'+std::string(s+exp,ptr-(s+exp)+1); @@ -1832,7 +1843,7 @@ inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const { // Remove zeros starting from right end char* ptr = s+slen-1; - while (*ptr=='0' && ptr>s+1) ptr--; + while (*ptr=='0' && ptr>s+1) ptr--; if(ptr==s+1) out = std::string(s,2); else out = std::string(s,2)+'.'+std::string(s+2,ptr-(s+2)+1); @@ -1843,7 +1854,7 @@ inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const { // Remove zeros starting from right end char* ptr = s+slen-1; - while (*ptr=='0' && ptr>s) ptr--; + while (*ptr=='0' && ptr>s) ptr--; if(ptr==s) out = std::string(s,1); else out = std::string(s,1)+'.'+std::string(s+1,ptr-(s+1)+1); @@ -1870,7 +1881,7 @@ 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 +inline std::ostream& mpreal::output(std::ostream& os) const { std::ostringstream format; const std::ios::fmtflags flags = os.flags(); @@ -1931,14 +1942,9 @@ inline int bits2digits(mp_prec_t b) ////////////////////////////////////////////////////////////////////////// // Set/Get number properties -inline int sgn(const mpreal& op) -{ - return mpfr_sgn(op.mpfr_srcptr()); -} - inline mpreal& mpreal::setSign(int sign, mp_rnd_t RoundingMode) { - mpfr_setsign(mpfr_ptr(), mpfr_srcptr(), (sign < 0 ? 1 : 0), RoundingMode); + mpfr_setsign(mpfr_ptr(), mpfr_srcptr(), sign < 0, RoundingMode); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } @@ -1955,14 +1961,14 @@ inline mpreal& mpreal::setPrecision(int Precision, mp_rnd_t RoundingMode) return *this; } -inline mpreal& mpreal::setInf(int sign) -{ +inline mpreal& mpreal::setInf(int sign) +{ mpfr_set_inf(mpfr_ptr(), sign); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; -} +} -inline mpreal& mpreal::setNan() +inline mpreal& mpreal::setNan() { mpfr_set_nan(mpfr_ptr()); MPREAL_MSVC_DEBUGVIEW_CODE; @@ -1976,7 +1982,7 @@ inline mpreal& mpreal::setZero(int sign) #else mpfr_set_si(mpfr_ptr(), 0, (mpfr_get_default_rounding_mode)()); setSign(sign); -#endif +#endif MPREAL_MSVC_DEBUGVIEW_CODE; return *this; @@ -1993,7 +1999,7 @@ inline void mpreal::set_prec(mp_prec_t prec, mp_rnd_t rnd_mode) MPREAL_MSVC_DEBUGVIEW_CODE; } -inline mp_exp_t mpreal::get_exp () +inline mp_exp_t mpreal::get_exp () const { return mpfr_get_exp(mpfr_srcptr()); } @@ -2022,7 +2028,7 @@ inline const mpreal ldexp(const mpreal& v, mp_exp_t exp) mpreal x(v); // rounding is not important since we are just increasing the exponent (= exact operation) - mpfr_mul_2si(x.mpfr_ptr(), x.mpfr_srcptr(), exp, mpreal::get_default_rnd()); + mpfr_mul_2si(x.mpfr_ptr(), x.mpfr_srcptr(), exp, mpreal::get_default_rnd()); return x; } @@ -2038,7 +2044,7 @@ inline mpreal machine_epsilon(mp_prec_t prec) } inline mpreal machine_epsilon(const mpreal& x) -{ +{ /* the smallest eps such that x + eps != x */ if( x < 0) { @@ -2059,7 +2065,7 @@ inline mpreal minval(mp_prec_t prec) inline mpreal maxval(mp_prec_t prec) { /* max = (1 - eps) * 2^emax, eps is machine epsilon */ - return (mpreal(1, prec) - machine_epsilon(prec)) << mpreal::get_emax(); + return (mpreal(1, prec) - machine_epsilon(prec)) << mpreal::get_emax(); } inline bool isEqualUlps(const mpreal& a, const mpreal& b, int maxUlps) @@ -2091,12 +2097,18 @@ inline bool signbit(const mpreal& x) return mpfr_signbit(x.mpfr_srcptr()); } +inline mpreal& setsignbit(mpreal& x, bool minus, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) +{ + mpfr_setsign(x.mpfr_ptr(), x.mpfr_srcptr(), minus, rnd_mode); + return x; +} + inline const mpreal modf(const mpreal& v, mpreal& n) { mpreal f(v); // rounding is not important since we are using the same number - mpfr_frac (f.mpfr_ptr(),f.mpfr_srcptr(),mpreal::get_default_rnd()); + mpfr_frac (f.mpfr_ptr(),f.mpfr_srcptr(),mpreal::get_default_rnd()); mpfr_trunc(n.mpfr_ptr(),v.mpfr_srcptr()); return f; } @@ -2159,7 +2171,7 @@ inline mp_exp_t mpreal::get_emax_max (void) #define MPREAL_UNARY_MATH_FUNCTION_BODY(f) \ mpreal y(0, mpfr_get_prec(x.mpfr_srcptr())); \ mpfr_##f(y.mpfr_ptr(), x.mpfr_srcptr(), r); \ - return y; + return y; inline const mpreal sqr (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(sqr ); } @@ -2182,7 +2194,7 @@ inline const mpreal sqrt(const unsigned int v, mp_rnd_t rnd_mode) inline const mpreal sqrt(const long int v, mp_rnd_t rnd_mode) { if (v>=0) return sqrt(static_cast(v),rnd_mode); - else return mpreal().setNan(); // NaN + else return mpreal().setNan(); // NaN } inline const mpreal sqrt(const int v, mp_rnd_t rnd_mode) @@ -2193,9 +2205,9 @@ inline const mpreal sqrt(const int v, mp_rnd_t rnd_mode) 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; + mpreal y(0, mpfr_get_prec(x.mpfr_srcptr())); + mpfr_rootn_ui(y.mpfr_ptr(), x.mpfr_srcptr(), k, r); + return y; } inline const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t r = mpreal::get_default_rnd()) @@ -2270,6 +2282,16 @@ inline const mpreal besselj1(const mpreal& x, mp_rnd_t r = mpreal::get_default_r 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 nextpow2(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) +{ + mpreal y(0, x.getPrecision()); + + if(!iszero(x)) + y = ceil(log2(abs(x,r),r)); + + return y; +} + 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())); @@ -2284,8 +2306,46 @@ inline const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = return a; } -inline const mpreal remainder (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) +inline const mpreal hypot(const mpreal& a, const mpreal& b, const mpreal& c) { + if(isnan EIGEN_NOT_A_MACRO (a) || isnan EIGEN_NOT_A_MACRO (b) || isnan EIGEN_NOT_A_MACRO(c)) return mpreal().setNan(); + else + { + mpreal absa = abs(a), absb = abs(b), absc = abs(c); + mpreal w = (std::max)(absa, (std::max)(absb, absc)); + mpreal r; + + if (!iszero(w)) + { + mpreal iw = 1/w; + r = w * sqrt(sqr(absa*iw) + sqr(absb*iw) + sqr(absc*iw)); + } + + return r; + } +} + +inline const mpreal hypot(const mpreal& a, const mpreal& b, const mpreal& c, const mpreal& d) +{ + if(isnan EIGEN_NOT_A_MACRO (a) || isnan EIGEN_NOT_A_MACRO (b) || isnan EIGEN_NOT_A_MACRO (c) || isnan EIGEN_NOT_A_MACRO (d)) return mpreal().setNan(); + else + { + mpreal absa = abs(a), absb = abs(b), absc = abs(c), absd = abs(d); + mpreal w = (std::max)(absa, (std::max)(absb, (std::max)(absc, absd))); + mpreal r; + + if (!iszero(w)) + { + mpreal iw = 1/w; + r = w * sqrt(sqr(absa*iw) + sqr(absb*iw) + sqr(absc*iw) + sqr(absd*iw)); + } + + return r; + } +} + +inline const mpreal remainder (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) +{ mpreal a(0,(std::max)(y.getPrecision(), x.getPrecision())); mpfr_remainder(a.mpfr_ptr(), x.mpfr_srcptr(), y.mpfr_srcptr(), rnd_mode); return a; @@ -2299,7 +2359,7 @@ inline const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t } 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()) + mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { mpreal x(0, prec); mpfr_fac_ui(x.mpfr_ptr(),v,rnd_mode); @@ -2338,9 +2398,9 @@ inline const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, m mpreal a; mp_prec_t p1, p2, p3; - p1 = v1.get_prec(); - p2 = v2.get_prec(); - p3 = v3.get_prec(); + p1 = v1.get_prec(); + p2 = v2.get_prec(); + p3 = v3.get_prec(); a.set_prec(p3>p2?(p3>p1?p3:p1):(p2>p1?p2:p1)); @@ -2353,9 +2413,9 @@ inline const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, m mpreal a; mp_prec_t p1, p2, p3; - p1 = v1.get_prec(); - p2 = v2.get_prec(); - p3 = v3.get_prec(); + p1 = v1.get_prec(); + p2 = v2.get_prec(); + p3 = v3.get_prec(); a.set_prec(p3>p2?(p3>p1?p3:p1):(p2>p1?p2:p1)); @@ -2368,8 +2428,8 @@ inline const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode = mpreal a; mp_prec_t p1, p2; - p1 = v1.get_prec(); - p2 = v2.get_prec(); + p1 = v1.get_prec(); + p2 = v2.get_prec(); a.set_prec(p1>p2?p1:p2); @@ -2387,7 +2447,7 @@ inline const mpreal sum (const mpreal tab[], const unsigned long int n, int& sta mpreal x; status = mpfr_sum(x.mpfr_ptr(), (mpfr_ptr*)p, n, mode); - + delete [] p; return x; } @@ -2401,9 +2461,9 @@ inline int sinh_cosh(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = return mpfr_sinh_cosh(s.mp,c.mp,v.mp,rnd_mode); } -inline const mpreal li2 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) -{ - MPREAL_UNARY_MATH_FUNCTION_BODY(li2); +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 = mpreal::get_default_rnd()) @@ -2415,16 +2475,16 @@ inline const mpreal rem (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = m inline const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { (void)rnd_mode; - - /* + + /* m = mod(x,y) if y != 0, returns x - n*y where n = floor(x/y) The following are true by convention: - mod(x,0) is x - mod(x,x) is 0 - - mod(x,y) for x != y and y != 0 has the same sign as y. - + - mod(x,y) for x != y and y != 0 has the same sign as y. + */ if(iszero(y)) return x; @@ -2432,9 +2492,7 @@ inline const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = m mpreal m = x - floor(x / y) * y; - m.setSign(sgn(y)); // make sure result has the same sign as Y - - return m; + return copysign(abs(m),y); // make sure result has the same sign as Y } inline const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) @@ -2442,8 +2500,8 @@ inline const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal a; mp_prec_t yp, xp; - yp = y.get_prec(); - xp = x.get_prec(); + yp = y.get_prec(); + xp = x.get_prec(); a.set_prec(yp>xp?yp:xp); @@ -2543,7 +2601,16 @@ inline const mpreal frac (const mpreal& x, mp_rnd_t r = mpreal::get_defaul ////////////////////////////////////////////////////////////////////////// // Miscellaneous Functions -inline void swap (mpreal& a, mpreal& b) { mpfr_swap(a.mp,b.mp); } +inline int sgn(const mpreal& op) +{ + // Please note, this is classic signum function which ignores sign of zero. + // Use signbit if you need sign of zero. + return mpfr_sgn(op.mpfr_srcptr()); +} + +////////////////////////////////////////////////////////////////////////// +// Miscellaneous Functions +inline void swap (mpreal& a, mpreal& b) { mpfr_swap(a.mpfr_ptr(),b.mpfr_ptr()); } 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= MPFR_VERSION_NUM(3,1,0)) +#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,1,0) && MPFR_VERSION < MPFR_VERSION_NUM(4,0,0)) +// TODO: +// Use mpfr_nrandom since mpfr_grandom is deprecated +#if defined(_MSC_VER) +#pragma warning( push ) +#pragma warning( disable : 1478) +#endif inline const mpreal grandom (gmp_randstate_t& state, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { mpreal x; mpfr_grandom(x.mpfr_ptr(), NULL, state, rnd_mode); return x; } +#if defined(_MSC_VER) +#pragma warning( pop ) +#endif inline const mpreal grandom(unsigned int seed = 0) { @@ -2664,17 +2740,17 @@ inline const mpreal grandom(unsigned int seed = 0) ////////////////////////////////////////////////////////////////////////// // Set/Get global properties inline void mpreal::set_default_prec(mp_prec_t prec) -{ - mpfr_set_default_prec(prec); +{ + mpfr_set_default_prec(prec); } inline void mpreal::set_default_rnd(mp_rnd_t rnd_mode) -{ - mpfr_set_default_rounding_mode(rnd_mode); +{ + mpfr_set_default_rounding_mode(rnd_mode); } inline bool mpreal::fits_in_bits(double x, int n) -{ +{ int i; double t; return IsInf(x) || (std::modf ( std::ldexp ( std::frexp ( x, &i ), n ), &t ) == 0.0); @@ -2924,7 +3000,7 @@ inline const mpreal pow(const int a, const double b, mp_rnd_t rnd_mode) else return pow(mpreal(a),mpreal(b),rnd_mode); //mpfr_pow } -// pow long double +// pow long double inline const mpreal pow(const long double a, const long double b, mp_rnd_t rnd_mode) { return pow(mpreal(a),mpreal(b),rnd_mode); @@ -2981,11 +3057,11 @@ 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 { - // we are allowed to extend namespace std with specializations only + // we are allowed to extend namespace std with specializations only template <> - inline void swap(mpfr::mpreal& x, mpfr::mpreal& y) - { - return mpfr::swap(x, y); + inline void swap(mpfr::mpreal& x, mpfr::mpreal& y) + { + return mpfr::swap(x, y); } template<> @@ -2996,7 +3072,7 @@ namespace std static const bool is_signed = true; static const bool is_integer = false; static const bool is_exact = false; - static const int radix = 2; + static const int radix = 2; static const bool has_infinity = true; static const bool has_quiet_NaN = true; @@ -3016,7 +3092,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); } @@ -3024,8 +3100,8 @@ namespace std { mp_rnd_t r = mpfr::mpreal::get_default_rnd(); - if(r == GMP_RNDN) return mpfr::mpreal(0.5, precision); - else return mpfr::mpreal(1.0, precision); + if(r == GMP_RNDN) return mpfr::mpreal(0.5, precision); + else return mpfr::mpreal(1.0, precision); } inline static const mpfr::mpreal infinity() { return mpfr::const_infinity(); } @@ -3036,17 +3112,17 @@ namespace std // Please note, exponent range is not fixed in MPFR static const int min_exponent = MPFR_EMIN_DEFAULT; static const int max_exponent = MPFR_EMAX_DEFAULT; - 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); + 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); #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. + // So we define them as functions here. // // This is preferable way for std::numeric_limits 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. + // 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(); @@ -3054,9 +3130,9 @@ namespace std 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; + 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; } } @@ -3083,13 +3159,13 @@ namespace std // 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. + // Change them accordingly to your application. // // For example, if you use 256 bits of precision uniformly in your program, then: // digits = 256 - // digits10 = 77 + // 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;