GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-09-12 09:23:12 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merged in advanpix/eigen-mp-devs (pull request PR-32)

Browse Source 
Fixes for SparseMatrix to support non-POD scalar types
This commit is contained in:
Gael Guennebaud 2013-09-03 22:05:14 +02:00
commit ed78a76161
5 changed files with 174 additions and 89 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

29
Eigen/src/Core/util/Memory.h
View File

@ -6,6 +6,7 @@
// Copyright (C) 2009 Kenneth Riddile <kfriddile@yahoo.com> // Copyright (C) 2009 Kenneth Riddile <kfriddile@yahoo.com>
// Copyright (C) 2010 Hauke Heibel <hauke.heibel@gmail.com> // Copyright (C) 2010 Hauke Heibel <hauke.heibel@gmail.com>
// Copyright (C) 2010 Thomas Capricelli <orzel@freehackers.org> // Copyright (C) 2010 Thomas Capricelli <orzel@freehackers.org>
// Copyright (C) 2013 Pavel Holoborodko <pavel@holoborodko.com>
// //
// This Source Code Form is subject to the terms of the Mozilla // 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 // Public License v. 2.0. If a copy of the MPL was not distributed
@ -514,6 +515,34 @@ template<typename T> struct smart_copy_helper<T,false> {
{ std::copy(start, end, target); } { std::copy(start, end, target); }
}; };
// intelligent memmove. falls back to std::memmove for POD types, uses std::copy otherwise.
template<typename T, bool UseMemmove> struct smart_memmove_helper;
template<typename T> void smart_memmove(const T* start, const T* end, T* target)
{
smart_memmove_helper<T,!NumTraits<T>::RequireInitialization>::run(start, end, target);
}
template<typename T> struct smart_memmove_helper<T,true> {
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<typename T> struct smart_memmove_helper<T,false> {
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) *** *** Implementation of runtime stack allocation (falling back to malloc) ***

4
Eigen/src/SparseCore/CompressedStorage.h
View File

@ -51,8 +51,8 @@ class CompressedStorage
CompressedStorage& operator=(const CompressedStorage& other) CompressedStorage& operator=(const CompressedStorage& other)
{ {
resize(other.size()); resize(other.size());
memcpy(m_values, other.m_values, m_size * sizeof(Scalar)); internal::smart_copy(other.m_values, other.m_values + m_size, m_values);
memcpy(m_indices, other.m_indices, m_size * sizeof(Index)); internal::smart_copy(other.m_indices, other.m_indices + m_size, m_indices);
return *this; return *this;
} }

27
Eigen/src/SparseCore/SparseBlock.h
View File

@ -67,12 +67,13 @@ public:
Index m_outerStart; Index m_outerStart;
const internal::variable_if_dynamic<Index, OuterSize> m_outerSize; const internal::variable_if_dynamic<Index, OuterSize> m_outerSize;
public:
EIGEN_INHERIT_ASSIGNMENT_OPERATORS(BlockImpl) EIGEN_INHERIT_ASSIGNMENT_OPERATORS(BlockImpl)
}; };
/*************************************************************************** /***************************************************************************
* specialisation for SparseMatrix * specialization for SparseMatrix
***************************************************************************/ ***************************************************************************/
template<typename _Scalar, int _Options, typename _Index, int BlockRows, int BlockCols> template<typename _Scalar, int _Options, typename _Index, int BlockRows, int BlockCols>
@ -125,7 +126,7 @@ public:
{ {
typedef typename internal::remove_all<typename SparseMatrixType::Nested>::type _NestedMatrixType; typedef typename internal::remove_all<typename SparseMatrixType::Nested>::type _NestedMatrixType;
_NestedMatrixType& matrix = const_cast<_NestedMatrixType&>(m_matrix);; _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. // 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 // 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 // 2 - let's check whether there is enough allocated memory
Index nnz = tmp.nonZeros(); Index nnz = tmp.nonZeros();
Index start = m_outerStart==0 ? 0 : matrix.outerIndexPtr()[m_outerStart]; // starting position of the current block 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 block_size = end - start; // available room in the current block
Index tail_size = m_matrix.outerIndexPtr()[m_matrix.outerSize()] - end; Index tail_size = m_matrix.outerIndexPtr()[m_matrix.outerSize()] - end;
@ -147,14 +148,14 @@ public:
// realloc manually to reduce copies // realloc manually to reduce copies
typename SparseMatrixType::Storage newdata(m_matrix.data().allocatedSize() - block_size + nnz); typename SparseMatrixType::Storage newdata(m_matrix.data().allocatedSize() - block_size + nnz);
std::memcpy(&newdata.value(0), &m_matrix.data().value(0), start*sizeof(Scalar)); internal::smart_copy(&m_matrix.data().value(0), &m_matrix.data().value(0) + start, &newdata.value(0));
std::memcpy(&newdata.index(0), &m_matrix.data().index(0), start*sizeof(Index)); 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)); internal::smart_copy(&tmp.data().value(0), &tmp.data().value(0) + nnz, &newdata.value(start));
std::memcpy(&newdata.index(start), &tmp.data().index(0), nnz*sizeof(Index)); 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)); internal::smart_copy(&matrix.data().value(end), &matrix.data().value(end) + tail_size, &newdata.value(start+nnz));
std::memcpy(&newdata.index(start+nnz), &matrix.data().index(end), tail_size*sizeof(Index)); 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); 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 // no need to realloc, simply copy the tail at its respective position and insert tmp
matrix.data().resize(start + nnz + tail_size); matrix.data().resize(start + nnz + tail_size);
std::memmove(&matrix.data().value(start+nnz), &matrix.data().value(end), tail_size*sizeof(Scalar)); internal::smart_memmove(&matrix.data().value(end), &matrix.data().value(end) + tail_size, &matrix.data().value(start + nnz));
std::memmove(&matrix.data().index(start+nnz), &matrix.data().index(end), tail_size*sizeof(Index)); 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)); internal::smart_copy(&tmp.data().value(0), &tmp.data().value(0) + nnz, &matrix.value(start));
std::memcpy(&matrix.data().index(start), &tmp.data().index(0), nnz*sizeof(Index)); internal::smart_copy(&tmp.data().index(0), &tmp.data().index(0) + nnz, &matrix.index(start));
} }
// update innerNonZeros // update innerNonZeros

2
Eigen/src/SparseCore/SparseMatrix.h
View File

@ -711,7 +711,7 @@ class SparseMatrix
initAssignment(other); initAssignment(other);
if(other.isCompressed()) 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; m_data = other.m_data;
} }
else else

193
unsupported/test/mpreal/mpreal.h
View File

@ -5,7 +5,7 @@
Project homepage: http://www.holoborodko.com/pavel/mpfr Project homepage: http://www.holoborodko.com/pavel/mpfr
Contact e-mail: pavel@holoborodko.com Contact e-mail: pavel@holoborodko.com
Copyright (c) 2008-2012 Pavel Holoborodko Copyright (c) 2008-2013 Pavel Holoborodko
Contributors: Contributors:
Dmitriy Gubanov, Konstantin Holoborodko, Brian Gladman, Dmitriy Gubanov, Konstantin Holoborodko, Brian Gladman,
@ -65,6 +65,7 @@
#include <stdexcept> #include <stdexcept>
#include <cfloat> #include <cfloat>
#include <cmath> #include <cmath>
#include <cstring>
#include <limits> #include <limits>
// Options // Options
@ -82,6 +83,20 @@
#define IsInf(x) std::isinf(x) // GNU C/C++ (and/or other compilers), just hope for C99 conformance #define IsInf(x) std::isinf(x) // GNU C/C++ (and/or other compilers), just hope for C99 conformance
#endif #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) #if defined(MPREAL_HAVE_INT64_SUPPORT)
#define MPFR_USE_INTMAX_T // Should be defined before mpfr.h #define MPFR_USE_INTMAX_T // Should be defined before mpfr.h
@ -111,7 +126,7 @@
#endif #endif
#if defined(MPREAL_HAVE_MSVC_DEBUGVIEW) && defined(_MSC_VER) && defined(_DEBUG) #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; #define MPREAL_MSVC_DEBUGVIEW_DATA std::string DebugView;
#else #else
#define MPREAL_MSVC_DEBUGVIEW_CODE #define MPREAL_MSVC_DEBUGVIEW_CODE
@ -149,7 +164,6 @@ public:
// Constructors && type conversions // Constructors && type conversions
mpreal(); mpreal();
mpreal(const mpreal& u); mpreal(const mpreal& u);
mpreal(const mpfr_t u);
mpreal(const mpf_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 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 mpq_t u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd());
@ -160,6 +174,10 @@ public:
mpreal(const long 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()); 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) #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()); mpreal(const uint64_t u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd());
mpreal(const int64_t u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); mpreal(const int64_t u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd());
@ -170,6 +188,11 @@ public:
~mpreal(); ~mpreal();
#ifdef MPREAL_HAVE_MOVE_SUPPORT
mpreal& operator=(mpreal&& v);
mpreal(mpreal&& u);
#endif
// Operations // 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 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 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 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 #endif
// Uniformly distributed random number generation in [0,1] using // Uniformly distributed random number generation in [0,1] using
@ -537,6 +562,9 @@ private:
// at the beginning of // at the beginning of
// [Visual Studio Installation Folder]\Common7\Packages\Debugger\autoexp.dat // [Visual Studio Installation Folder]\Common7\Packages\Debugger\autoexp.dat
MPREAL_MSVC_DEBUGVIEW_DATA 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; 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_ptr()->_mpfr_d = 0; // make sure "other" holds no pinter to actual data
mpfr_set(mp,u,mpreal::get_default_rnd());
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; 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; 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() inline mpreal::~mpreal()
{ {
mpfr_clear(mp); clear(mpfr_ptr());
} }
// internal namespace needed for template magic // 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); mp_prec_t vp = mpfr_get_prec(v.mp);
if(tp != vp){ if(tp != vp){
mpfr_clear(mp); clear(mpfr_ptr());
mpfr_init2(mp, vp); mpfr_init2(mpfr_ptr(), vp);
} }
mpfr_set(mp, v.mp, mpreal::get_default_rnd()); mpfr_set(mp, v.mp, mpreal::get_default_rnd());
@ -974,7 +1034,7 @@ inline mpreal& mpreal::operator=(const char* s)
MPREAL_MSVC_DEBUGVIEW_CODE; MPREAL_MSVC_DEBUGVIEW_CODE;
} }
mpfr_clear(t); clear(t);
return *this; return *this;
} }
@ -997,7 +1057,7 @@ inline mpreal& mpreal::operator=(const std::string& s)
MPREAL_MSVC_DEBUGVIEW_CODE; MPREAL_MSVC_DEBUGVIEW_CODE;
} }
mpfr_clear(t); clear(t);
return *this; return *this;
} }
@ -1202,7 +1262,7 @@ inline const mpreal mpreal::operator-()const
inline const mpreal operator-(const mpreal& a, const mpreal& b) 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()); mpfr_sub(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
return c; return c;
} }
@ -1319,7 +1379,7 @@ inline mpreal& mpreal::operator*=(const int v)
inline const mpreal operator*(const mpreal& a, const mpreal& b) 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()); mpfr_mul(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
return c; return c;
} }
@ -1395,7 +1455,7 @@ inline mpreal& mpreal::operator/=(const int v)
inline const mpreal operator/(const mpreal& a, const mpreal& b) 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()); mpfr_div(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
return c; return c;
} }
@ -1927,12 +1987,12 @@ inline bool isEqualUlps(const mpreal& a, const mpreal& b, int maxUlps)
inline bool isEqualFuzzy(const mpreal& a, const mpreal& b, const mpreal& eps) 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) 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) 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) inline int cmpabs(const mpreal& a,const mpreal& b)
{ {
return mpfr_cmpabs(a.mp,b.mp); return mpfr_cmpabs(a.mpfr_ptr(), b.mpfr_srcptr());
} }
inline int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode) inline int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode)
{ {
return mpfr_sin_cos(s.mp,c.mp,v.mp,rnd_mode); return mpfr_sin_cos(s.mpfr_ptr(), c.mpfr_ptr(), v.mpfr_srcptr(), rnd_mode);
} }
inline const mpreal sqrt (const long double v, mp_rnd_t rnd_mode) { return sqrt(mpreal(v),rnd_mode); } inline const mpreal sqrt (const 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 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 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 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 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 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 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 acot (const mpreal& v, mp_rnd_t r) { return atan (1/v, r); }
inline const mpreal asec (const mpreal& v, mp_rnd_t r) { return acos (1/v, r); } inline const mpreal 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) inline const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode)
{ {
mpreal a; mpreal a(0,(std::max)(y.getPrecision(), x.getPrecision()));
mp_prec_t yp, xp; mpfr_atan2(a.mpfr_ptr(), y.mpfr_srcptr(), x.mpfr_srcptr(), rnd_mode);
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);
return a; return a;
} }
inline const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode) inline const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
{ {
mpreal a; mpreal a(0,(std::max)(y.getPrecision(), x.getPrecision()));
mp_prec_t yp, xp; mpfr_hypot(a.mpfr_ptr(), x.mpfr_srcptr(), y.mpfr_srcptr(), rnd_mode);
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);
return a; return a;
} }
inline const mpreal remainder (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode) inline const mpreal remainder (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
{ {
mpreal a; mpreal a(0,(std::max)(y.getPrecision(), x.getPrecision()));
mp_prec_t yp, xp; mpfr_remainder(a.mpfr_ptr(), x.mpfr_srcptr(), y.mpfr_srcptr(), rnd_mode);
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);
return a; return a;
} }
inline const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode) inline const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
{ {
mpreal a; mpreal a(0,(std::max)(y.getPrecision(), x.getPrecision()));
mp_prec_t yp, xp; mpfr_remquo(a.mpfr_ptr(),q, x.mpfr_srcptr(), y.mpfr_srcptr(), rnd_mode);
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);
return a; return a;
} }
inline const mpreal fac_ui (unsigned long int v, mp_prec_t prec, mp_rnd_t rnd_mode) inline const mpreal fac_ui (unsigned long int v, mp_prec_t prec, mp_rnd_t rnd_mode)
{ {
mpreal x(0, prec); mpreal x(0, prec);
mpfr_fac_ui(x.mp,v,rnd_mode); mpfr_fac_ui(x.mpfr_ptr(),v,rnd_mode);
return x; return x;
} }
@ -2181,7 +2209,7 @@ inline const mpreal lgamma (const mpreal& v, int *signp, mp_rnd_t rnd_mode)
mpreal x(v); mpreal x(v);
int tsignp; 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); else mpfr_lgamma(x.mp,&tsignp,v.mp,rnd_mode);
return x; return x;
@ -2365,7 +2393,7 @@ inline const mpreal const_catalan (mp_prec_t p, mp_rnd_t r)
return x; 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); mpreal x(0, p);
mpfr_set_inf(x.mpfr_ptr(), sign); 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); mpfr_urandom(x.mp,state,rnd_mode);
return x; 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 #endif
#if (MPFR_VERSION <= MPFR_VERSION_NUM(2,4,2)) #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 // Set/Get global properties
inline void mpreal::set_default_prec(mp_prec_t prec) inline void mpreal::set_default_prec(mp_prec_t prec)
@ -2859,10 +2914,10 @@ namespace std
switch (r) switch (r)
{ {
case MPFR_RNDN: return round_to_nearest; case GMP_RNDN: return round_to_nearest;
case MPFR_RNDZ: return round_toward_zero; case GMP_RNDZ: return round_toward_zero;
case MPFR_RNDU: return round_toward_infinity; case GMP_RNDU: return round_toward_infinity;
case MPFR_RNDD: return round_toward_neg_infinity; case GMP_RNDD: return round_toward_neg_infinity;
default: return round_indeterminate; default: return round_indeterminate;
} }
} }
@ -2881,7 +2936,7 @@ namespace std
{ {
mp_rnd_t r = mpfr::mpreal::get_default_rnd(); 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); else return mpfr::mpreal(1.0, precision);
} }