• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1// This file is part of a joint effort between Eigen, a lightweight C++ template library
2// for linear algebra, and MPFR C++, a C++ interface to MPFR library (http://www.holoborodko.com/pavel/)
3//
4// Copyright (C) 2010 Pavel Holoborodko <pavel@holoborodko.com>
5// Copyright (C) 2010 Konstantin Holoborodko <konstantin@holoborodko.com>
6// Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@inria.fr>
7//
8// This Source Code Form is subject to the terms of the Mozilla
9// Public License v. 2.0. If a copy of the MPL was not distributed
10// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
11//
12// Contributors:
13// Brian Gladman, Helmut Jarausch, Fokko Beekhof, Ulrich Mutze, Heinz van Saanen, Pere Constans
14
15#ifndef EIGEN_MPREALSUPPORT_MODULE_H
16#define EIGEN_MPREALSUPPORT_MODULE_H
17
18#include <ctime>
19#include <mpreal.h>
20#include <Eigen/Core>
21
22namespace Eigen {
23
24  /** \ingroup Unsupported_modules
25    * \defgroup MPRealSupport_Module MPFRC++ Support module
26    *
27    * \code
28    * #include <Eigen/MPRealSupport>
29    * \endcode
30    *
31    * This module provides support for multi precision floating point numbers
32    * via the <a href="http://www.holoborodko.com/pavel/mpfr">MPFR C++</a>
33    * library which itself is built upon <a href="http://www.mpfr.org/">MPFR</a>/<a href="http://gmplib.org/">GMP</a>.
34    *
35    * You can find a copy of MPFR C++ that is known to be compatible in the unsupported/test/mpreal folder.
36    *
37    * Here is an example:
38    *
39\code
40#include <iostream>
41#include <Eigen/MPRealSupport>
42#include <Eigen/LU>
43using namespace mpfr;
44using namespace Eigen;
45int main()
46{
47  // set precision to 256 bits (double has only 53 bits)
48  mpreal::set_default_prec(256);
49  // Declare matrix and vector types with multi-precision scalar type
50  typedef Matrix<mpreal,Dynamic,Dynamic>  MatrixXmp;
51  typedef Matrix<mpreal,Dynamic,1>        VectorXmp;
52
53  MatrixXmp A = MatrixXmp::Random(100,100);
54  VectorXmp b = VectorXmp::Random(100);
55
56  // Solve Ax=b using LU
57  VectorXmp x = A.lu().solve(b);
58  std::cout << "relative error: " << (A*x - b).norm() / b.norm() << std::endl;
59  return 0;
60}
61\endcode
62    *
63    */
64
65  template<> struct NumTraits<mpfr::mpreal>
66    : GenericNumTraits<mpfr::mpreal>
67  {
68    enum {
69      IsInteger = 0,
70      IsSigned = 1,
71      IsComplex = 0,
72      RequireInitialization = 1,
73      ReadCost = 10,
74      AddCost = 10,
75      MulCost = 40
76    };
77
78    typedef mpfr::mpreal Real;
79    typedef mpfr::mpreal NonInteger;
80
81    inline static mpfr::mpreal highest() { return  mpfr::mpreal_max(mpfr::mpreal::get_default_prec()); }
82    inline static mpfr::mpreal lowest()  { return -mpfr::mpreal_max(mpfr::mpreal::get_default_prec()); }
83
84    inline static Real epsilon()
85    {
86      return mpfr::machine_epsilon(mpfr::mpreal::get_default_prec());
87    }
88    inline static Real dummy_precision()
89    {
90      unsigned int weak_prec = ((mpfr::mpreal::get_default_prec()-1)*90)/100;
91      return mpfr::machine_epsilon(weak_prec);
92    }
93  };
94
95namespace internal {
96
97  template<> mpfr::mpreal random<mpfr::mpreal>()
98  {
99#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
100    static gmp_randstate_t state;
101    static bool isFirstTime = true;
102
103    if(isFirstTime)
104    {
105      gmp_randinit_default(state);
106      gmp_randseed_ui(state,(unsigned)time(NULL));
107      isFirstTime = false;
108    }
109
110    return mpfr::urandom(state)*2-1;
111#else
112    return mpfr::mpreal(random<double>());
113#endif
114  }
115
116  template<> mpfr::mpreal random<mpfr::mpreal>(const mpfr::mpreal& a, const mpfr::mpreal& b)
117  {
118    return a + (b-a) * random<mpfr::mpreal>();
119  }
120
121  bool isMuchSmallerThan(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& prec)
122  {
123    return mpfr::abs(a) <= mpfr::abs(b) * prec;
124  }
125
126  inline bool isApprox(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& prec)
127  {
128    return mpfr::abs(a - b) <= (mpfr::min)(mpfr::abs(a), mpfr::abs(b)) * prec;
129  }
130
131  inline bool isApproxOrLessThan(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& prec)
132  {
133    return a <= b || isApprox(a, b, prec);
134  }
135
136  template<> inline long double cast<mpfr::mpreal,long double>(const mpfr::mpreal& x)
137  { return x.toLDouble(); }
138  template<> inline double cast<mpfr::mpreal,double>(const mpfr::mpreal& x)
139  { return x.toDouble(); }
140  template<> inline long cast<mpfr::mpreal,long>(const mpfr::mpreal& x)
141  { return x.toLong(); }
142  template<> inline int cast<mpfr::mpreal,int>(const mpfr::mpreal& x)
143  { return int(x.toLong()); }
144
145} // end namespace internal
146}
147
148#endif // EIGEN_MPREALSUPPORT_MODULE_H
149