mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-06-04 18:54:00 +08:00
add a support module for MPFR C++ with basic unit testing
This commit is contained in:
Gael Guennebaud
parent
bfbe61454e
commit
87e89fea4e
21
cmake/FindGMP.cmake
Normal file
21
cmake/FindGMP.cmake
Normal file
@ -0,0 +1,21 @@
|
|||||||
|
# Try to find the GNU Multiple Precision Arithmetic Library (GMP)
|
||||||
|
# See http://gmplib.org/
|
||||||
|
|
||||||
|
if (GMP_INCLUDES AND GMP_LIBRARIES)
|
||||||
|
set(GMP_FIND_QUIETLY TRUE)
|
||||||
|
endif (GMP_INCLUDES AND GMP_LIBRARIES)
|
||||||
|
|
||||||
|
find_path(GMP_INCLUDES
|
||||||
|
NAMES
|
||||||
|
gmp.h
|
||||||
|
PATHS
|
||||||
|
$ENV{GMPDIR}
|
||||||
|
${INCLUDE_INSTALL_DIR}
|
||||||
|
)
|
||||||
|
|
||||||
|
find_library(GMP_LIBRARIES gmp PATHS $ENV{GMPDIR} ${LIB_INSTALL_DIR})
|
||||||
|
|
||||||
|
include(FindPackageHandleStandardArgs)
|
||||||
|
find_package_handle_standard_args(GMP DEFAULT_MSG
|
||||||
|
GMP_INCLUDES GMP_LIBRARIES)
|
||||||
|
mark_as_advanced(GMP_INCLUDES GMP_LIBRARIES)
|
||||||
21
cmake/FindMPFR.cmake
Normal file
21
cmake/FindMPFR.cmake
Normal file
@ -0,0 +1,21 @@
|
|||||||
|
# Try to find the MPFR library
|
||||||
|
# See http://www.mpfr.org/
|
||||||
|
|
||||||
|
if (MPFR_INCLUDES AND MPFR_LIBRARIES)
|
||||||
|
set(MPFR_FIND_QUIETLY TRUE)
|
||||||
|
endif (MPFR_INCLUDES AND MPFR_LIBRARIES)
|
||||||
|
|
||||||
|
find_path(MPFR_INCLUDES
|
||||||
|
NAMES
|
||||||
|
mpfr.h
|
||||||
|
PATHS
|
||||||
|
$ENV{GMPDIR}
|
||||||
|
${INCLUDE_INSTALL_DIR}
|
||||||
|
)
|
||||||
|
|
||||||
|
find_library(MPFR_LIBRARIES mpfr PATHS $ENV{GMPDIR} ${LIB_INSTALL_DIR})
|
||||||
|
|
||||||
|
include(FindPackageHandleStandardArgs)
|
||||||
|
find_package_handle_standard_args(MPFR DEFAULT_MSG
|
||||||
|
MPFR_INCLUDES MPFR_LIBRARIES)
|
||||||
|
mark_as_advanced(MPFR_INCLUDES MPFR_LIBRARIES)
|
||||||
@ -61,4 +61,7 @@ namespace Eigen {
|
|||||||
/** \ingroup Unsupported_modules
|
/** \ingroup Unsupported_modules
|
||||||
* \defgroup SparseExtra_Module */
|
* \defgroup SparseExtra_Module */
|
||||||
|
|
||||||
|
/** \ingroup Unsupported_modules
|
||||||
|
* \defgroup MpfrcxxSupport_Module */
|
||||||
|
|
||||||
} // namespace Eigen
|
} // namespace Eigen
|
||||||
|
|||||||
160
unsupported/Eigen/MPRealSupport
Normal file
160
unsupported/Eigen/MPRealSupport
Normal file
@ -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 <pavel@holoborodko.com>
|
||||||
|
// Copyright (C) 2010 Konstantin Holoborodko <konstantin@holoborodko.com>
|
||||||
|
// Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@inria.fr>
|
||||||
|
//
|
||||||
|
// 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 <http://www.gnu.org/licenses/>.
|
||||||
|
//
|
||||||
|
// 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 <mpreal.h>
|
||||||
|
#include <Eigen/Core>
|
||||||
|
|
||||||
|
namespace Eigen {
|
||||||
|
|
||||||
|
/** \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/?page_id=12">MPFR C++</a>
|
||||||
|
* library which itself is built upon <a href="http://www.mpfr.org/">MPFR</a>/<a href="http://gmplib.org/">GMP</a>.
|
||||||
|
*
|
||||||
|
* Here is an example:
|
||||||
|
*
|
||||||
|
\code
|
||||||
|
#include <iostream>
|
||||||
|
#include <Eigen/Mpfrc++Support>
|
||||||
|
#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,
|
||||||
|
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<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(ei_random<double>());
|
||||||
|
#endif
|
||||||
|
}
|
||||||
|
|
||||||
|
template<> mpfr::mpreal ei_random<mpfr::mpreal>(const mpfr::mpreal& a, const mpfr::mpreal& b)
|
||||||
|
{
|
||||||
|
return a + (b-a) * ei_random<mpfr::mpreal>();
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
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
|
||||||
@ -2,6 +2,8 @@ FILE(GLOB examples_SRCS "*.cpp")
|
|||||||
|
|
||||||
ADD_CUSTOM_TARGET(unsupported_examples)
|
ADD_CUSTOM_TARGET(unsupported_examples)
|
||||||
|
|
||||||
|
INCLUDE_DIRECTORIES(../../../unsupported ../../../unsupported/test)
|
||||||
|
|
||||||
FOREACH(example_src ${examples_SRCS})
|
FOREACH(example_src ${examples_SRCS})
|
||||||
GET_FILENAME_COMPONENT(example ${example_src} NAME_WE)
|
GET_FILENAME_COMPONENT(example ${example_src} NAME_WE)
|
||||||
ADD_EXECUTABLE(example_${example} ${example_src})
|
ADD_EXECUTABLE(example_${example} ${example_src})
|
||||||
|
|||||||
@ -55,10 +55,10 @@ include_directories(../../test ../../unsupported ../../Eigen)
|
|||||||
|
|
||||||
if(ADOLC_FOUND)
|
if(ADOLC_FOUND)
|
||||||
include_directories(${ADOLC_INCLUDES})
|
include_directories(${ADOLC_INCLUDES})
|
||||||
ei_add_property(EIGEN_TESTED_BACKENDS "Adolc")
|
ei_add_property(EIGEN_TESTED_BACKENDS "Adolc, ")
|
||||||
ei_add_test(forward_adolc " " ${ADOLC_LIBRARIES})
|
ei_add_test(forward_adolc " " ${ADOLC_LIBRARIES})
|
||||||
else(ADOLC_FOUND)
|
else(ADOLC_FOUND)
|
||||||
ei_add_property(EIGEN_MISSING_BACKENDS "Adolc")
|
ei_add_property(EIGEN_MISSING_BACKENDS "Adolc, ")
|
||||||
endif(ADOLC_FOUND)
|
endif(ADOLC_FOUND)
|
||||||
|
|
||||||
ei_add_test(NonLinearOptimization)
|
ei_add_test(NonLinearOptimization)
|
||||||
@ -70,6 +70,15 @@ ei_add_test(matrix_function)
|
|||||||
ei_add_test(alignedvector3)
|
ei_add_test(alignedvector3)
|
||||||
ei_add_test(FFT)
|
ei_add_test(FFT)
|
||||||
|
|
||||||
|
find_package(MPFR)
|
||||||
|
if(MPFR_FOUND)
|
||||||
|
include_directories(${MPFR_INCLUDES})
|
||||||
|
ei_add_property(EIGEN_TESTED_BACKENDS "MPFR C++, ")
|
||||||
|
ei_add_test(mpreal_support " " ${MPFR_LIBRARIES} )
|
||||||
|
else()
|
||||||
|
ei_add_property(EIGEN_MISSING_BACKENDS "MPFR C++, ")
|
||||||
|
endif()
|
||||||
|
|
||||||
ei_add_test(sparse_llt " " "${SPARSE_LIBS}")
|
ei_add_test(sparse_llt " " "${SPARSE_LIBS}")
|
||||||
ei_add_test(sparse_ldlt " " "${SPARSE_LIBS}")
|
ei_add_test(sparse_ldlt " " "${SPARSE_LIBS}")
|
||||||
ei_add_test(sparse_lu " " "${SPARSE_LIBS}")
|
ei_add_test(sparse_lu " " "${SPARSE_LIBS}")
|
||||||
|
|||||||
408
unsupported/test/mpreal.cpp
Normal file
408
unsupported/test/mpreal.cpp
Normal file
@ -0,0 +1,408 @@
|
|||||||
|
/*
|
||||||
|
Multi-precision real number class. C++ wrapper fo MPFR library.
|
||||||
|
Project homepage: http://www.holoborodko.com/pavel/
|
||||||
|
Contact e-mail: pavel@holoborodko.com
|
||||||
|
|
||||||
|
Copyright (c) 2008-2010 Pavel Holoborodko
|
||||||
|
|
||||||
|
This library 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 2.1 of the License, or (at your option) any later version.
|
||||||
|
|
||||||
|
This library 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 for more details.
|
||||||
|
|
||||||
|
You should have received a copy of the GNU Lesser General Public
|
||||||
|
License along with this library; if not, write to the Free Software
|
||||||
|
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
|
||||||
|
|
||||||
|
Contributors:
|
||||||
|
Brian Gladman, Helmut Jarausch, Fokko Beekhof, Ulrich Mutze,
|
||||||
|
Heinz van Saanen, Pere Constans, Dmitriy Gubanov
|
||||||
|
*/
|
||||||
|
|
||||||
|
#include <cstring>
|
||||||
|
#include "mpreal.h"
|
||||||
|
|
||||||
|
using std::ws;
|
||||||
|
using std::cerr;
|
||||||
|
using std::endl;
|
||||||
|
using std::string;
|
||||||
|
using std::ostream;
|
||||||
|
using std::istream;
|
||||||
|
|
||||||
|
namespace mpfr{
|
||||||
|
|
||||||
|
mp_rnd_t mpreal::default_rnd = mpfr_get_default_rounding_mode();
|
||||||
|
mp_prec_t mpreal::default_prec = mpfr_get_default_prec();
|
||||||
|
int mpreal::default_base = 10;
|
||||||
|
int mpreal::double_bits = -1;
|
||||||
|
|
||||||
|
// Default constructor: creates mp number and initializes it to 0.
|
||||||
|
mpreal::mpreal()
|
||||||
|
{
|
||||||
|
mpfr_init2(mp,default_prec);
|
||||||
|
mpfr_set_ui(mp,0,default_rnd);
|
||||||
|
}
|
||||||
|
|
||||||
|
mpreal::mpreal(const mpreal& u)
|
||||||
|
{
|
||||||
|
mpfr_init2(mp,mpfr_get_prec(u.mp));
|
||||||
|
mpfr_set(mp,u.mp,default_rnd);
|
||||||
|
}
|
||||||
|
|
||||||
|
mpreal::mpreal(const mpfr_t u)
|
||||||
|
{
|
||||||
|
mpfr_init2(mp,mpfr_get_prec(u));
|
||||||
|
mpfr_set(mp,u,default_rnd);
|
||||||
|
}
|
||||||
|
|
||||||
|
mpreal::mpreal(const mpf_t u)
|
||||||
|
{
|
||||||
|
mpfr_init2(mp,mpf_get_prec(u));
|
||||||
|
mpfr_set_f(mp,u,default_rnd);
|
||||||
|
}
|
||||||
|
|
||||||
|
mpreal::mpreal(const mpz_t u, mp_prec_t prec, mp_rnd_t mode)
|
||||||
|
{
|
||||||
|
mpfr_init2(mp,prec);
|
||||||
|
mpfr_set_z(mp,u,mode);
|
||||||
|
}
|
||||||
|
|
||||||
|
mpreal::mpreal(const mpq_t u, mp_prec_t prec, mp_rnd_t mode)
|
||||||
|
{
|
||||||
|
mpfr_init2(mp,prec);
|
||||||
|
mpfr_set_q(mp,u,mode);
|
||||||
|
}
|
||||||
|
|
||||||
|
mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode)
|
||||||
|
{
|
||||||
|
if(double_bits == -1 || fits_in_bits(u, double_bits))
|
||||||
|
{
|
||||||
|
mpfr_init2(mp,prec);
|
||||||
|
mpfr_set_d(mp,u,mode);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
throw conversion_overflow();
|
||||||
|
}
|
||||||
|
|
||||||
|
mpreal::mpreal(const long double u, mp_prec_t prec, mp_rnd_t mode)
|
||||||
|
{
|
||||||
|
mpfr_init2(mp,prec);
|
||||||
|
mpfr_set_ld(mp,u,mode);
|
||||||
|
}
|
||||||
|
|
||||||
|
mpreal::mpreal(const unsigned long int u, mp_prec_t prec, mp_rnd_t mode)
|
||||||
|
{
|
||||||
|
mpfr_init2(mp,prec);
|
||||||
|
mpfr_set_ui(mp,u,mode);
|
||||||
|
}
|
||||||
|
|
||||||
|
mpreal::mpreal(const unsigned int u, mp_prec_t prec, mp_rnd_t mode)
|
||||||
|
{
|
||||||
|
mpfr_init2(mp,prec);
|
||||||
|
mpfr_set_ui(mp,u,mode);
|
||||||
|
}
|
||||||
|
|
||||||
|
mpreal::mpreal(const long int u, mp_prec_t prec, mp_rnd_t mode)
|
||||||
|
{
|
||||||
|
mpfr_init2(mp,prec);
|
||||||
|
mpfr_set_si(mp,u,mode);
|
||||||
|
}
|
||||||
|
|
||||||
|
mpreal::mpreal(const int u, mp_prec_t prec, mp_rnd_t mode)
|
||||||
|
{
|
||||||
|
mpfr_init2(mp,prec);
|
||||||
|
mpfr_set_si(mp,u,mode);
|
||||||
|
}
|
||||||
|
|
||||||
|
mpreal::mpreal(const char* s, mp_prec_t prec, int base, mp_rnd_t mode)
|
||||||
|
{
|
||||||
|
mpfr_init2(mp,prec);
|
||||||
|
mpfr_set_str(mp, s, base, mode);
|
||||||
|
}
|
||||||
|
|
||||||
|
mpreal::~mpreal()
|
||||||
|
{
|
||||||
|
mpfr_clear(mp);
|
||||||
|
}
|
||||||
|
|
||||||
|
// Operators - Assignment
|
||||||
|
mpreal& mpreal::operator=(const char* s)
|
||||||
|
{
|
||||||
|
mpfr_t t;
|
||||||
|
if(0==mpfr_init_set_str(t,s,default_base,default_rnd))
|
||||||
|
{
|
||||||
|
mpfr_set(mp,t,mpreal::default_rnd);
|
||||||
|
mpfr_clear(t);
|
||||||
|
}else{
|
||||||
|
mpfr_clear(t);
|
||||||
|
// cerr<<"fail to convert string"<<endl;
|
||||||
|
}
|
||||||
|
|
||||||
|
return *this;
|
||||||
|
}
|
||||||
|
|
||||||
|
const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode)
|
||||||
|
{
|
||||||
|
mpreal a;
|
||||||
|
mp_prec_t p1, p2, p3;
|
||||||
|
|
||||||
|
p1 = v1.get_prec();
|
||||||
|
p2 = v2.get_prec();
|
||||||
|
p3 = v3.get_prec();
|
||||||
|
|
||||||
|
a.set_prec(p3>p2?(p3>p1?p3:p1):(p2>p1?p2:p1));
|
||||||
|
|
||||||
|
mpfr_fma(a.mp,v1.mp,v2.mp,v3.mp,rnd_mode);
|
||||||
|
return a;
|
||||||
|
}
|
||||||
|
|
||||||
|
const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode)
|
||||||
|
{
|
||||||
|
mpreal a;
|
||||||
|
mp_prec_t p1, p2, p3;
|
||||||
|
|
||||||
|
p1 = v1.get_prec();
|
||||||
|
p2 = v2.get_prec();
|
||||||
|
p3 = v3.get_prec();
|
||||||
|
|
||||||
|
a.set_prec(p3>p2?(p3>p1?p3:p1):(p2>p1?p2:p1));
|
||||||
|
|
||||||
|
mpfr_fms(a.mp,v1.mp,v2.mp,v3.mp,rnd_mode);
|
||||||
|
return a;
|
||||||
|
}
|
||||||
|
|
||||||
|
const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode)
|
||||||
|
{
|
||||||
|
mpreal a;
|
||||||
|
mp_prec_t p1, p2;
|
||||||
|
|
||||||
|
p1 = v1.get_prec();
|
||||||
|
p2 = v2.get_prec();
|
||||||
|
|
||||||
|
a.set_prec(p1>p2?p1:p2);
|
||||||
|
|
||||||
|
mpfr_agm(a.mp, v1.mp, v2.mp, rnd_mode);
|
||||||
|
|
||||||
|
return a;
|
||||||
|
}
|
||||||
|
|
||||||
|
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);
|
||||||
|
|
||||||
|
return a;
|
||||||
|
}
|
||||||
|
|
||||||
|
const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode)
|
||||||
|
{
|
||||||
|
mpreal x;
|
||||||
|
mpfr_ptr* t;
|
||||||
|
unsigned long int i;
|
||||||
|
|
||||||
|
t = new mpfr_ptr[n];
|
||||||
|
for (i=0;i<n;i++) t[i] = (mpfr_ptr)tab[i].mp;
|
||||||
|
mpfr_sum(x.mp,t,n,rnd_mode);
|
||||||
|
delete[] t;
|
||||||
|
return x;
|
||||||
|
}
|
||||||
|
|
||||||
|
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);
|
||||||
|
|
||||||
|
return a;
|
||||||
|
}
|
||||||
|
|
||||||
|
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);
|
||||||
|
|
||||||
|
return a;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <class T>
|
||||||
|
std::string to_string(T t, std::ios_base & (*f)(std::ios_base&))
|
||||||
|
{
|
||||||
|
std::ostringstream oss;
|
||||||
|
oss << f << t;
|
||||||
|
return oss.str();
|
||||||
|
}
|
||||||
|
|
||||||
|
mpreal::operator std::string() const
|
||||||
|
{
|
||||||
|
return to_string();
|
||||||
|
}
|
||||||
|
|
||||||
|
string mpreal::to_string(size_t n, int b, mp_rnd_t mode) const
|
||||||
|
{
|
||||||
|
char *s, *ns = NULL;
|
||||||
|
size_t slen, nslen;
|
||||||
|
mp_exp_t exp;
|
||||||
|
string out;
|
||||||
|
|
||||||
|
if(mpfr_inf_p(mp))
|
||||||
|
{
|
||||||
|
if(mpfr_sgn(mp)>0) return "+@Inf@";
|
||||||
|
else return "-@Inf@";
|
||||||
|
}
|
||||||
|
|
||||||
|
if(mpfr_zero_p(mp)) return "0";
|
||||||
|
if(mpfr_nan_p(mp)) return "@NaN@";
|
||||||
|
|
||||||
|
|
||||||
|
s = mpfr_get_str(NULL,&exp,b,0,mp,mode);
|
||||||
|
ns = mpfr_get_str(NULL,&exp,b,n,mp,mode);
|
||||||
|
|
||||||
|
if(s!=NULL && ns!=NULL)
|
||||||
|
{
|
||||||
|
slen = strlen(s);
|
||||||
|
nslen = strlen(ns);
|
||||||
|
if(nslen<=slen)
|
||||||
|
{
|
||||||
|
mpfr_free_str(s);
|
||||||
|
s = ns;
|
||||||
|
slen = nslen;
|
||||||
|
}
|
||||||
|
else {
|
||||||
|
mpfr_free_str(ns);
|
||||||
|
}
|
||||||
|
|
||||||
|
// Make human eye-friendly formatting if possible
|
||||||
|
if (exp>0 && static_cast<size_t>(exp)<slen)
|
||||||
|
{
|
||||||
|
if(s[0]=='-')
|
||||||
|
{
|
||||||
|
// Remove zeros starting from right end
|
||||||
|
char* ptr = s+slen-1;
|
||||||
|
while (*ptr=='0' && ptr>s+exp) ptr--;
|
||||||
|
|
||||||
|
if(ptr==s+exp) out = string(s,exp+1);
|
||||||
|
else out = string(s,exp+1)+'.'+string(s+exp+1,ptr-(s+exp+1)+1);
|
||||||
|
|
||||||
|
//out = string(s,exp+1)+'.'+string(s+exp+1);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
// Remove zeros starting from right end
|
||||||
|
char* ptr = s+slen-1;
|
||||||
|
while (*ptr=='0' && ptr>s+exp-1) ptr--;
|
||||||
|
|
||||||
|
if(ptr==s+exp-1) out = string(s,exp);
|
||||||
|
else out = string(s,exp)+'.'+string(s+exp,ptr-(s+exp)+1);
|
||||||
|
|
||||||
|
//out = string(s,exp)+'.'+string(s+exp);
|
||||||
|
}
|
||||||
|
|
||||||
|
}else{ // exp<0 || exp>slen
|
||||||
|
if(s[0]=='-')
|
||||||
|
{
|
||||||
|
// Remove zeros starting from right end
|
||||||
|
char* ptr = s+slen-1;
|
||||||
|
while (*ptr=='0' && ptr>s+1) ptr--;
|
||||||
|
|
||||||
|
if(ptr==s+1) out = string(s,2);
|
||||||
|
else out = string(s,2)+'.'+string(s+2,ptr-(s+2)+1);
|
||||||
|
|
||||||
|
//out = string(s,2)+'.'+string(s+2);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
// Remove zeros starting from right end
|
||||||
|
char* ptr = s+slen-1;
|
||||||
|
while (*ptr=='0' && ptr>s) ptr--;
|
||||||
|
|
||||||
|
if(ptr==s) out = string(s,1);
|
||||||
|
else out = string(s,1)+'.'+string(s+1,ptr-(s+1)+1);
|
||||||
|
|
||||||
|
//out = string(s,1)+'.'+string(s+1);
|
||||||
|
}
|
||||||
|
|
||||||
|
// Make final string
|
||||||
|
if(--exp)
|
||||||
|
{
|
||||||
|
if(exp>0) out += "e+"+mpfr::to_string<mp_exp_t>(exp,std::dec);
|
||||||
|
else out += "e"+mpfr::to_string<mp_exp_t>(exp,std::dec);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
mpfr_free_str(s);
|
||||||
|
return out;
|
||||||
|
}else{
|
||||||
|
return "conversion error!";
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
//////////////////////////////////////////////////////////////////////////
|
||||||
|
// I/O
|
||||||
|
ostream& operator<<(ostream& os, const mpreal& v)
|
||||||
|
{
|
||||||
|
return os<<v.to_string(os.precision());
|
||||||
|
}
|
||||||
|
|
||||||
|
istream& operator>>(istream &is, mpreal& v)
|
||||||
|
{
|
||||||
|
char c;
|
||||||
|
string s = "";
|
||||||
|
mpfr_t t;
|
||||||
|
|
||||||
|
if(is.good())
|
||||||
|
{
|
||||||
|
is>>ws;
|
||||||
|
while ((c = is.get())!=EOF)
|
||||||
|
{
|
||||||
|
if(c ==' ' || c == '\t' || c == '\n' || c == '\r')
|
||||||
|
{
|
||||||
|
is.putback(c);
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
s += c;
|
||||||
|
}
|
||||||
|
|
||||||
|
if(s.size() != 0)
|
||||||
|
{
|
||||||
|
// Protect current value from alternation in case of input error
|
||||||
|
// so some error handling(roll back) procedure can be used
|
||||||
|
if(0==mpfr_init_set_str(t,s.c_str(),mpreal::default_base,mpreal::default_rnd))
|
||||||
|
{
|
||||||
|
mpfr_set(v.mp,t,mpreal::default_rnd);
|
||||||
|
mpfr_clear(t);
|
||||||
|
|
||||||
|
}else{
|
||||||
|
mpfr_clear(t);
|
||||||
|
cerr<<"error reading from istream"<<endl;
|
||||||
|
// throw an exception
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return is;
|
||||||
|
}
|
||||||
|
}
|
||||||
3122
unsupported/test/mpreal.h
Normal file
3122
unsupported/test/mpreal.h
Normal file
File diff suppressed because it is too large
Load Diff
42
unsupported/test/mpreal_support.cpp
Normal file
42
unsupported/test/mpreal_support.cpp
Normal file
@ -0,0 +1,42 @@
|
|||||||
|
#include "main.h"
|
||||||
|
#include <Eigen/MPRealSupport>
|
||||||
|
|
||||||
|
using namespace mpfr;
|
||||||
|
using namespace std;
|
||||||
|
using namespace Eigen;
|
||||||
|
|
||||||
|
void test_mpreal_support()
|
||||||
|
{
|
||||||
|
// set precision to 256 bits (double has only 53 bits)
|
||||||
|
mpreal::set_default_prec(256);
|
||||||
|
typedef Matrix<mpreal,Eigen::Dynamic,Eigen::Dynamic> MatrixXmp;
|
||||||
|
|
||||||
|
std::cerr << "epsilon = " << NumTraits<mpreal>::epsilon() << "\n";
|
||||||
|
std::cerr << "dummy_precision = " << NumTraits<mpreal>::dummy_precision() << "\n";
|
||||||
|
std::cerr << "highest = " << NumTraits<mpreal>::highest() << "\n";
|
||||||
|
std::cerr << "lowest = " << NumTraits<mpreal>::lowest() << "\n";
|
||||||
|
|
||||||
|
for(int i = 0; i < g_repeat; i++) {
|
||||||
|
int s = ei_random<int>(1,100);
|
||||||
|
MatrixXmp A = MatrixXmp::Random(s,s);
|
||||||
|
MatrixXmp B = MatrixXmp::Random(s,s);
|
||||||
|
MatrixXmp S = A.adjoint() * A;
|
||||||
|
MatrixXmp X;
|
||||||
|
|
||||||
|
// Cholesky
|
||||||
|
X = S.selfadjointView<Lower>().llt().solve(B);
|
||||||
|
VERIFY_IS_APPROX((S.selfadjointView<Lower>()*X).eval(),B);
|
||||||
|
|
||||||
|
// partial LU
|
||||||
|
X = A.lu().solve(B);
|
||||||
|
VERIFY_IS_APPROX((A*X).eval(),B);
|
||||||
|
|
||||||
|
// symmetric eigenvalues
|
||||||
|
SelfAdjointEigenSolver<MatrixXmp> eig(S);
|
||||||
|
VERIFY_IS_EQUAL(eig.info(), Success);
|
||||||
|
VERIFY_IS_APPROX((S.selfadjointView<Lower>() * eig.eigenvectors()),
|
||||||
|
eig.eigenvectors() * eig.eigenvalues().asDiagonal());
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
#include "mpreal.cpp"
|
||||||
Loading…
x
Reference in New Issue
Block a user