vcglib/eigenlib/unsupported/Eigen/MPRealSupport

149 lines
4.4 KiB
Plaintext

// 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