add a support module for MPFR C++ with basic unit testing

1 files changed, 160 insertions, 0 deletions

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