diff options
Diffstat (limited to 'unsupported/Eigen/MPRealSupport')
-rw-r--r-- | unsupported/Eigen/MPRealSupport | 148 |
1 files changed, 148 insertions, 0 deletions
diff --git a/unsupported/Eigen/MPRealSupport b/unsupported/Eigen/MPRealSupport new file mode 100644 index 000000000..3895623fe --- /dev/null +++ b/unsupported/Eigen/MPRealSupport @@ -0,0 +1,148 @@ +// This file is part of a joint effort between Eigen, a lightweight C++ template library +// for linear algebra, and MPFR C++, a C++ interface to MPFR library (http://www.holoborodko.com/pavel/) +// +// Copyright (C) 2010 Pavel Holoborodko <pavel@holoborodko.com> +// Copyright (C) 2010 Konstantin Holoborodko <konstantin@holoborodko.com> +// Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@inria.fr> +// +// 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 +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. +// +// Contributors: +// Brian Gladman, Helmut Jarausch, Fokko Beekhof, Ulrich Mutze, Heinz van Saanen, Pere Constans + +#ifndef EIGEN_MPREALSUPPORT_MODULE_H +#define EIGEN_MPREALSUPPORT_MODULE_H + +#include <ctime> +#include <mpreal.h> +#include <Eigen/Core> + +namespace Eigen { + + /** \ingroup Unsupported_modules + * \defgroup MPRealSupport_Module MPFRC++ Support module + * + * \code + * #include <Eigen/MPRealSupport> + * \endcode + * + * This module provides support for multi precision floating point numbers + * via the <a href="http://www.holoborodko.com/pavel/mpfr">MPFR C++</a> + * library which itself is built upon <a href="http://www.mpfr.org/">MPFR</a>/<a href="http://gmplib.org/">GMP</a>. + * + * You can find a copy of MPFR C++ that is known to be compatible in the unsupported/test/mpreal folder. + * + * Here is an example: + * +\code +#include <iostream> +#include <Eigen/MPRealSupport> +#include <Eigen/LU> +using namespace mpfr; +using namespace Eigen; +int main() +{ + // set precision to 256 bits (double has only 53 bits) + mpreal::set_default_prec(256); + // Declare matrix and vector types with multi-precision scalar type + typedef Matrix<mpreal,Dynamic,Dynamic> MatrixXmp; + typedef Matrix<mpreal,Dynamic,1> VectorXmp; + + MatrixXmp A = MatrixXmp::Random(100,100); + VectorXmp b = VectorXmp::Random(100); + + // Solve Ax=b using LU + VectorXmp x = A.lu().solve(b); + std::cout << "relative error: " << (A*x - b).norm() / b.norm() << std::endl; + return 0; +} +\endcode + * + */ + + template<> struct NumTraits<mpfr::mpreal> + : GenericNumTraits<mpfr::mpreal> + { + enum { + IsInteger = 0, + IsSigned = 1, + IsComplex = 0, + RequireInitialization = 1, + ReadCost = 10, + AddCost = 10, + MulCost = 40 + }; + + typedef mpfr::mpreal Real; + typedef mpfr::mpreal NonInteger; + + inline static mpfr::mpreal highest() { return mpfr::mpreal_max(mpfr::mpreal::get_default_prec()); } + inline static mpfr::mpreal lowest() { return -mpfr::mpreal_max(mpfr::mpreal::get_default_prec()); } + + inline static Real epsilon() + { + return mpfr::machine_epsilon(mpfr::mpreal::get_default_prec()); + } + inline static Real dummy_precision() + { + unsigned int weak_prec = ((mpfr::mpreal::get_default_prec()-1)*90)/100; + return mpfr::machine_epsilon(weak_prec); + } + }; + +namespace internal { + + template<> mpfr::mpreal random<mpfr::mpreal>() + { +#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0)) + static gmp_randstate_t state; + static bool isFirstTime = true; + + if(isFirstTime) + { + gmp_randinit_default(state); + gmp_randseed_ui(state,(unsigned)time(NULL)); + isFirstTime = false; + } + + return mpfr::urandom(state)*2-1; +#else + return mpfr::mpreal(random<double>()); +#endif + } + + template<> mpfr::mpreal random<mpfr::mpreal>(const mpfr::mpreal& a, const mpfr::mpreal& b) + { + return a + (b-a) * random<mpfr::mpreal>(); + } + + bool isMuchSmallerThan(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& prec) + { + return mpfr::abs(a) <= mpfr::abs(b) * prec; + } + + inline bool isApprox(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& prec) + { + return mpfr::abs(a - b) <= (mpfr::min)(mpfr::abs(a), mpfr::abs(b)) * prec; + } + + inline bool isApproxOrLessThan(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& prec) + { + return a <= b || isApprox(a, b, prec); + } + + template<> inline long double cast<mpfr::mpreal,long double>(const mpfr::mpreal& x) + { return x.toLDouble(); } + template<> inline double cast<mpfr::mpreal,double>(const mpfr::mpreal& x) + { return x.toDouble(); } + template<> inline long cast<mpfr::mpreal,long>(const mpfr::mpreal& x) + { return x.toLong(); } + template<> inline int cast<mpfr::mpreal,int>(const mpfr::mpreal& x) + { return int(x.toLong()); } + +} // end namespace internal +} + +#endif // EIGEN_MPREALSUPPORT_MODULE_H |