Home | History | Annotate | Download | only in Eigen 
      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 (a] holoborodko.com>
      5 // Copyright (C) 2010 Konstantin Holoborodko <konstantin (a] holoborodko.com>
      6 // Copyright (C) 2010 Gael Guennebaud <gael.guennebaud (a] 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 
     22 namespace 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>
     43 using namespace mpfr;
     44 using namespace Eigen;
     45 int 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 
     95 namespace 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