From 87e89fea4e05f2f516bf6efb3986a8ee9aee6de6 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Thu, 15 Jul 2010 16:29:04 +0200 Subject: add a support module for MPFR C++ with basic unit testing --- unsupported/Eigen/MPRealSupport | 160 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 160 insertions(+) create mode 100644 unsupported/Eigen/MPRealSupport (limited to 'unsupported/Eigen/MPRealSupport') diff --git a/unsupported/Eigen/MPRealSupport b/unsupported/Eigen/MPRealSupport new file mode 100644 index 000000000..abc719ba7 --- /dev/null +++ b/unsupported/Eigen/MPRealSupport @@ -0,0 +1,160 @@ +// 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 +// Copyright (C) 2010 Konstantin Holoborodko +// Copyright (C) 2010 Gael Guennebaud +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// this library. If not, see . +// +// 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 +#include + +namespace Eigen { + + /** \defgroup MPRealSupport_Module MPFRC++ Support module + * + * \code + * #include + * \endcode + * + * This module provides support for multi precision floating point numbers + * via the MPFR C++ + * library which itself is built upon MPFR/GMP. + * + * Here is an example: + * +\code +#include +#include +#include +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 MatrixXmp; + typedef Matrix 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 + : GenericNumTraits + { + enum { + IsInteger = 0, + IsSigned = 1, + IsComplex = 0, + 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); + } + }; + + template<> mpfr::mpreal ei_random() + { +#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(ei_random()); +#endif + } + + template<> mpfr::mpreal ei_random(const mpfr::mpreal& a, const mpfr::mpreal& b) + { + return a + (b-a) * ei_random(); + } +} + +namespace mpfr { + + inline const mpreal& ei_conj(const mpreal& x) { return x; } + inline const mpreal& ei_real(const mpreal& x) { return x; } + inline mpreal ei_imag(const mpreal&) { return 0.0; } + inline mpreal ei_abs(const mpreal& x) { return fabs(x); } + inline mpreal ei_abs2(const mpreal& x) { return x*x; } + inline mpreal ei_sqrt(const mpreal& x) { return sqrt(x); } + inline mpreal ei_exp(const mpreal& x) { return exp(x); } + inline mpreal ei_log(const mpreal& x) { return log(x); } + inline mpreal ei_sin(const mpreal& x) { return sin(x); } + inline mpreal ei_cos(const mpreal& x) { return cos(x); } + inline mpreal ei_pow(const mpreal& x, mpreal& y) { return pow(x, y); } + + bool ei_isMuchSmallerThan(const mpreal& a, const mpreal& b, const mpreal& prec) + { + return ei_abs(a) <= abs(b) * prec; + } + + inline bool ei_isApprox(const mpreal& a, const mpreal& b, const mpreal& prec) + { + return ei_abs(a - b) <= min(abs(a), abs(b)) * prec; + } + + inline bool ei_isApproxOrLessThan(const mpreal& a, const mpreal& b, const mpreal& prec) + { + return a <= b || ei_isApprox(a, b, prec); + } +} + +#endif // EIGEN_MPREALSUPPORT_MODULE_H -- cgit v1.2.3