From a147500dee85092fbfbf755123553f412a218ec2 Mon Sep 17 00:00:00 2001 From: Pavel Holoborodko Date: Sun, 25 Aug 2013 18:00:28 +0900 Subject: [PATCH 1/5] Added smart_memmove with support of non-POD scalars (e.g. needed in SparseBlock.h). --- Eigen/src/Core/util/Memory.h | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h index 451535a0c..3a357e079 100644 --- a/Eigen/src/Core/util/Memory.h +++ b/Eigen/src/Core/util/Memory.h @@ -6,6 +6,7 @@ // Copyright (C) 2009 Kenneth Riddile // Copyright (C) 2010 Hauke Heibel // Copyright (C) 2010 Thomas Capricelli +// Copyright (C) 2013 Pavel Holoborodko // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed @@ -514,6 +515,34 @@ template struct smart_copy_helper { { std::copy(start, end, target); } }; +// intelligent memmove. falls back to std::memmove for POD types, uses std::copy otherwise. +template struct smart_memmove_helper; + +template void smart_memmove(const T* start, const T* end, T* target) +{ + smart_memmove_helper::RequireInitialization>::run(start, end, target); +} + +template struct smart_memmove_helper { + static inline void run(const T* start, const T* end, T* target) + { std::memmove(target, start, std::ptrdiff_t(end)-std::ptrdiff_t(start)); } +}; + +template struct smart_memmove_helper { + static inline void run(const T* start, const T* end, T* target) + { + if (uintptr_t(target) < uintptr_t(start)) + { + std::copy(start, end, target); + } + else + { + std::ptrdiff_t count = (std::ptrdiff_t(end)-std::ptrdiff_t(start)) / sizeof(T); + std::copy_backward(start, end, target + count); + } + } +}; + /***************************************************************************** *** Implementation of runtime stack allocation (falling back to malloc) *** From 1472f4bc61192e4ffd06c6baf419747d7cc4c90b Mon Sep 17 00:00:00 2001 From: Pavel Holoborodko Date: Sun, 25 Aug 2013 18:02:07 +0900 Subject: [PATCH 2/5] Fixed bug #647 by using smart_copy instead of bitwise memcpy. --- Eigen/src/SparseCore/CompressedStorage.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Eigen/src/SparseCore/CompressedStorage.h b/Eigen/src/SparseCore/CompressedStorage.h index 3321fab4a..10d516a66 100644 --- a/Eigen/src/SparseCore/CompressedStorage.h +++ b/Eigen/src/SparseCore/CompressedStorage.h @@ -51,8 +51,8 @@ class CompressedStorage CompressedStorage& operator=(const CompressedStorage& other) { resize(other.size()); - memcpy(m_values, other.m_values, m_size * sizeof(Scalar)); - memcpy(m_indices, other.m_indices, m_size * sizeof(Index)); + internal::smart_copy(other.m_values, other.m_values + m_size, m_values); + internal::smart_copy(other.m_indices, other.m_indices + m_size, m_indices); return *this; } From e6462c2ce3e2d871892bae13f795d56afc7ec109 Mon Sep 17 00:00:00 2001 From: Pavel Holoborodko Date: Sun, 25 Aug 2013 18:03:49 +0900 Subject: [PATCH 3/5] Switched to smart_copy to support non-trivial scalar types --- Eigen/src/SparseCore/SparseMatrix.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Eigen/src/SparseCore/SparseMatrix.h b/Eigen/src/SparseCore/SparseMatrix.h index 970c6be71..4dad50562 100644 --- a/Eigen/src/SparseCore/SparseMatrix.h +++ b/Eigen/src/SparseCore/SparseMatrix.h @@ -711,7 +711,7 @@ class SparseMatrix initAssignment(other); if(other.isCompressed()) { - memcpy(m_outerIndex, other.m_outerIndex, (m_outerSize+1)*sizeof(Index)); + internal::smart_copy(other.m_outerIndex, other.m_outerIndex + m_outerSize + 1, m_outerIndex); m_data = other.m_data; } else From 41321e436606fa9399cf9b89586b04c2ba95c23e Mon Sep 17 00:00:00 2001 From: Pavel Holoborodko Date: Sun, 25 Aug 2013 18:12:15 +0900 Subject: [PATCH 4/5] Replaced memcpy & memmove to smart_* alternatives for non-POD scalar types --- Eigen/src/SparseCore/SparseBlock.h | 29 +++++++++++++++-------------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/Eigen/src/SparseCore/SparseBlock.h b/Eigen/src/SparseCore/SparseBlock.h index b6a2d232c..efd44284f 100644 --- a/Eigen/src/SparseCore/SparseBlock.h +++ b/Eigen/src/SparseCore/SparseBlock.h @@ -66,13 +66,14 @@ public: typename XprType::Nested m_matrix; Index m_outerStart; const internal::variable_if_dynamic m_outerSize; - + + public: EIGEN_INHERIT_ASSIGNMENT_OPERATORS(BlockImpl) }; /*************************************************************************** -* specialisation for SparseMatrix +* specialization for SparseMatrix ***************************************************************************/ template @@ -125,7 +126,7 @@ public: { typedef typename internal::remove_all::type _NestedMatrixType; _NestedMatrixType& matrix = const_cast<_NestedMatrixType&>(m_matrix);; - // This assignement is slow if this vector set is not empty + // This assignment is slow if this vector set is not empty // and/or it is not at the end of the nonzeros of the underlying matrix. // 1 - eval to a temporary to avoid transposition and/or aliasing issues @@ -134,7 +135,7 @@ public: // 2 - let's check whether there is enough allocated memory Index nnz = tmp.nonZeros(); Index start = m_outerStart==0 ? 0 : matrix.outerIndexPtr()[m_outerStart]; // starting position of the current block - Index end = m_matrix.outerIndexPtr()[m_outerStart+m_outerSize.value()]; // ending posiiton of the current block + Index end = m_matrix.outerIndexPtr()[m_outerStart+m_outerSize.value()]; // ending position of the current block Index block_size = end - start; // available room in the current block Index tail_size = m_matrix.outerIndexPtr()[m_matrix.outerSize()] - end; @@ -147,14 +148,14 @@ public: // realloc manually to reduce copies typename SparseMatrixType::Storage newdata(m_matrix.data().allocatedSize() - block_size + nnz); - std::memcpy(&newdata.value(0), &m_matrix.data().value(0), start*sizeof(Scalar)); - std::memcpy(&newdata.index(0), &m_matrix.data().index(0), start*sizeof(Index)); + internal::smart_copy(&m_matrix.data().value(0), &m_matrix.data().value(0) + start, &newdata.value(0)); + internal::smart_copy(&m_matrix.data().index(0), &m_matrix.data().index(0) + start, &newdata.index(0)); - std::memcpy(&newdata.value(start), &tmp.data().value(0), nnz*sizeof(Scalar)); - std::memcpy(&newdata.index(start), &tmp.data().index(0), nnz*sizeof(Index)); + internal::smart_copy(&tmp.data().value(0), &tmp.data().value(0) + nnz, &newdata.value(start)); + internal::smart_copy(&tmp.data().index(0), &tmp.data().index(0) + nnz, &newdata.index(start)); - std::memcpy(&newdata.value(start+nnz), &matrix.data().value(end), tail_size*sizeof(Scalar)); - std::memcpy(&newdata.index(start+nnz), &matrix.data().index(end), tail_size*sizeof(Index)); + internal::smart_copy(&matrix.data().value(end), &matrix.data().value(end) + tail_size, &newdata.value(start+nnz)); + internal::smart_copy(&matrix.data().index(end), &matrix.data().index(end) + tail_size, &newdata.index(start+nnz)); newdata.resize(m_matrix.outerIndexPtr()[m_matrix.outerSize()] - block_size + nnz); @@ -165,11 +166,11 @@ public: // no need to realloc, simply copy the tail at its respective position and insert tmp matrix.data().resize(start + nnz + tail_size); - std::memmove(&matrix.data().value(start+nnz), &matrix.data().value(end), tail_size*sizeof(Scalar)); - std::memmove(&matrix.data().index(start+nnz), &matrix.data().index(end), tail_size*sizeof(Index)); + internal::smart_memmove(&matrix.data().value(end), &matrix.data().value(end) + tail_size, &matrix.data().value(start + nnz)); + internal::smart_memmove(&matrix.data().index(end), &matrix.data().index(end) + tail_size, &matrix.data().index(start + nnz)); - std::memcpy(&matrix.data().value(start), &tmp.data().value(0), nnz*sizeof(Scalar)); - std::memcpy(&matrix.data().index(start), &tmp.data().index(0), nnz*sizeof(Index)); + internal::smart_copy(&tmp.data().value(0), &tmp.data().value(0) + nnz, &matrix.value(start)); + internal::smart_copy(&tmp.data().index(0), &tmp.data().index(0) + nnz, &matrix.index(start)); } // update innerNonZeros From d2c4f4ab21bd39485cb2acf3658aac1e2628d031 Mon Sep 17 00:00:00 2001 From: Pavel Holoborodko Date: Mon, 26 Aug 2013 00:22:18 +0900 Subject: [PATCH 5/5] Updated mpfr::mpreal. Move semantic support, RVO, other new features --- unsupported/test/mpreal/mpreal.h | 199 ++++++++++++++++++++----------- 1 file changed, 127 insertions(+), 72 deletions(-) diff --git a/unsupported/test/mpreal/mpreal.h b/unsupported/test/mpreal/mpreal.h index ef0a6a9f0..75bda7db7 100644 --- a/unsupported/test/mpreal/mpreal.h +++ b/unsupported/test/mpreal/mpreal.h @@ -5,7 +5,7 @@ Project homepage: http://www.holoborodko.com/pavel/mpfr Contact e-mail: pavel@holoborodko.com - Copyright (c) 2008-2012 Pavel Holoborodko + Copyright (c) 2008-2013 Pavel Holoborodko Contributors: Dmitriy Gubanov, Konstantin Holoborodko, Brian Gladman, @@ -65,6 +65,7 @@ #include #include #include +#include #include // Options @@ -82,6 +83,20 @@ #define IsInf(x) std::isinf(x) // GNU C/C++ (and/or other compilers), just hope for C99 conformance #endif +// Detect support for r-value references (move semantic). Borrowed from Eigen. +// A Clang feature extension to determine compiler features. +// We use it to determine 'cxx_rvalue_references' +#ifndef __has_feature +# define __has_feature(x) 0 +#endif + +#if (__has_feature(cxx_rvalue_references) || \ + defined(__GXX_EXPERIMENTAL_CXX0X__) || \ + (defined(_MSC_VER) && _MSC_VER >= 1600)) +#define MPREAL_HAVE_MOVE_SUPPORT +#endif + +// Detect available 64-bit capabilities #if defined(MPREAL_HAVE_INT64_SUPPORT) #define MPFR_USE_INTMAX_T // Should be defined before mpfr.h @@ -111,7 +126,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 +164,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 +173,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' shared 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 +188,11 @@ public: ~mpreal(); +#ifdef MPREAL_HAVE_MOVE_SUPPORT + mpreal& operator=(mpreal&& v); + mpreal(mpreal&& u); +#endif + // Operations // = // +, -, *, /, ++, --, <<, >> @@ -415,6 +438,8 @@ public: friend const mpreal digamma (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); friend const mpreal ai (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); friend const mpreal urandom (gmp_randstate_t& state, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); // use gmp_randinit_default() to init state, gmp_randclear() to clear + friend const mpreal grandom (gmp_randstate_t& state, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); // use gmp_randinit_default() to init state, gmp_randclear() to clear + friend const mpreal grandom (unsigned int seed = 0); #endif // Uniformly distributed random number generation in [0,1] using @@ -422,7 +447,7 @@ public: // 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); - + // 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); @@ -531,12 +556,15 @@ private: // Human friendly Debug Preview in Visual Studio. // Put one of these lines: // - // mpfr::mpreal= ; Show value only + // 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 + + // "Smart" resources deallocation. Checks if instance initialized before deletion. + void clear(::mpfr_ptr); }; ////////////////////////////////////////////////////////////////////////// @@ -565,10 +593,36 @@ inline mpreal::mpreal(const mpreal& u) MPREAL_MSVC_DEBUGVIEW_CODE; } -inline mpreal::mpreal(const mpfr_t u) +#ifdef MPREAL_HAVE_MOVE_SUPPORT +inline mpreal::mpreal(mpreal&& other) { - mpfr_init2(mp,mpfr_get_prec(u)); - mpfr_set(mp,u,mpreal::get_default_rnd()); + mpfr_ptr()->_mpfr_d = 0; // 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, bool shared) +{ + if(shared) + { + std::memcpy(mp, u, sizeof(mpfr_t)); + } + else + { + mpfr_init2(mp, mpfr_get_prec(u)); + mpfr_set (mp, u, mpreal::get_default_rnd()); + } MPREAL_MSVC_DEBUGVIEW_CODE; } @@ -688,9 +742,15 @@ inline mpreal::mpreal(const std::string& s, mp_prec_t prec, int base, mp_rnd_t m MPREAL_MSVC_DEBUGVIEW_CODE; } +inline void mpreal::clear(::mpfr_ptr x) +{ + // Use undocumented field in mpfr_t structure to check if it was initialized + if(0 != x->_mpfr_d) mpfr_clear(x); +} + inline mpreal::~mpreal() { - mpfr_clear(mp); + clear(mpfr_ptr()); } // internal namespace needed for template magic @@ -864,8 +924,8 @@ inline mpreal& mpreal::operator=(const mpreal& v) mp_prec_t vp = mpfr_get_prec(v.mp); if(tp != vp){ - mpfr_clear(mp); - mpfr_init2(mp, vp); + clear(mpfr_ptr()); + mpfr_init2(mpfr_ptr(), vp); } mpfr_set(mp, v.mp, mpreal::get_default_rnd()); @@ -974,7 +1034,7 @@ inline mpreal& mpreal::operator=(const char* s) MPREAL_MSVC_DEBUGVIEW_CODE; } - mpfr_clear(t); + clear(t); return *this; } @@ -997,7 +1057,7 @@ inline mpreal& mpreal::operator=(const std::string& s) MPREAL_MSVC_DEBUGVIEW_CODE; } - mpfr_clear(t); + clear(t); return *this; } @@ -1202,7 +1262,7 @@ 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()))); + 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; } @@ -1319,7 +1379,7 @@ 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()))); + 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; } @@ -1395,7 +1455,7 @@ 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()))); + 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; } @@ -1922,17 +1982,17 @@ 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) @@ -2048,12 +2108,12 @@ 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) { - 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); } @@ -2074,9 +2134,9 @@ inline const mpreal tan (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FU 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 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); } @@ -2110,68 +2170,36 @@ inline const mpreal bessely1(const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_ 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); - + 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) { - 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) { - 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) { - 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) { mpreal x(0, prec); - mpfr_fac_ui(x.mp,v,rnd_mode); + mpfr_fac_ui(x.mpfr_ptr(),v,rnd_mode); return x; } @@ -2181,7 +2209,7 @@ inline const mpreal lgamma (const mpreal& v, int *signp, mp_rnd_t rnd_mode) mpreal x(v); int tsignp; - if(signp) mpfr_lgamma(x.mp,signp,v.mp,rnd_mode); + if(signp) mpfr_lgamma(x.mp, signp,v.mp,rnd_mode); else mpfr_lgamma(x.mp,&tsignp,v.mp,rnd_mode); return x; @@ -2365,7 +2393,7 @@ inline const mpreal const_catalan (mp_prec_t p, mp_rnd_t 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, mp_prec_t p, mp_rnd_t r) { mpreal x(0, p); mpfr_set_inf(x.mpfr_ptr(), sign); @@ -2465,6 +2493,14 @@ inline const mpreal urandom (gmp_randstate_t& state, mp_rnd_t rnd_mode) mpfr_urandom(x.mp,state,rnd_mode); return x; } + +inline const mpreal grandom (gmp_randstate_t& state, mp_rnd_t rnd_mode) +{ + mpreal x; + mpfr_grandom(x.mp, NULL, state, rnd_mode); + return x; +} + #endif #if (MPFR_VERSION <= MPFR_VERSION_NUM(2,4,2)) @@ -2504,6 +2540,25 @@ inline const mpreal random(unsigned int seed) } +#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0)) +inline const mpreal grandom(unsigned int seed) +{ + 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) @@ -2859,10 +2914,10 @@ namespace std 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; + 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; } } @@ -2881,7 +2936,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); }