MPRealSupport

Go to the documentation of this file.

// 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-2012 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/.
 
#ifndef EIGEN_MPREALSUPPORT_MODULE_H
#define EIGEN_MPREALSUPPORT_MODULE_H
 
#include "../../Eigen/Core"
#include <mpreal.h>
 
namespace Eigen {
  
  template<> struct NumTraits<mpfr::mpreal>
    : GenericNumTraits<mpfr::mpreal>
  {
    enum {
      IsInteger = 0,
      IsSigned = 1,
      IsComplex = 0,
      RequireInitialization = 1,
      ReadCost = HugeCost,
      AddCost  = HugeCost,
      MulCost  = HugeCost
    };
 
    typedef mpfr::mpreal Real;
    typedef mpfr::mpreal NonInteger;
    
    static inline Real highest  (long Precision = mpfr::mpreal::get_default_prec()) { return  mpfr::maxval(Precision); }
    static inline Real lowest   (long Precision = mpfr::mpreal::get_default_prec()) { return -mpfr::maxval(Precision); }
 
    // Constants
    static inline Real Pi      (long Precision = mpfr::mpreal::get_default_prec())  { return mpfr::const_pi(Precision);        }
    static inline Real Euler   (long Precision = mpfr::mpreal::get_default_prec())  { return mpfr::const_euler(Precision);     }
    static inline Real Log2    (long Precision = mpfr::mpreal::get_default_prec())  { return mpfr::const_log2(Precision);      }
    static inline Real Catalan (long Precision = mpfr::mpreal::get_default_prec())  { return mpfr::const_catalan(Precision);   }
 
    static inline Real epsilon (long Precision = mpfr::mpreal::get_default_prec())  { return mpfr::machine_epsilon(Precision); }
    static inline Real epsilon (const Real& x)                                      { return mpfr::machine_epsilon(x); }
 
#ifdef MPREAL_HAVE_DYNAMIC_STD_NUMERIC_LIMITS
    static inline int digits10 (long Precision = mpfr::mpreal::get_default_prec())  { return std::numeric_limits<Real>::digits10(Precision); }
    static inline int digits10 (const Real& x)                                      { return std::numeric_limits<Real>::digits10(x); }
    
    static inline int digits ()               { return std::numeric_limits<Real>::digits(); }
    static inline int digits (const Real& x)  { return std::numeric_limits<Real>::digits(x); }
#endif
 
    static inline Real dummy_precision()
    {
      mpfr_prec_t weak_prec = ((mpfr::mpreal::get_default_prec()-1) * 90) / 100;
      return mpfr::machine_epsilon(weak_prec);
    }
  };
 
  namespace internal {
 
  template<> inline mpfr::mpreal random<mpfr::mpreal>()
  {
    return mpfr::random();
  }
 
  template<> inline mpfr::mpreal random<mpfr::mpreal>(const mpfr::mpreal& a, const mpfr::mpreal& b)
  {
    return a + (b-a) * random<mpfr::mpreal>();
  }
 
  inline bool isMuchSmallerThan(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& eps)
  {
    return mpfr::abs(a) <= mpfr::abs(b) * eps;
  }
 
  inline bool isApprox(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& eps)
  {
    return mpfr::isEqualFuzzy(a,b,eps);
  }
 
  inline bool isApproxOrLessThan(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& eps)
  {
    return a <= b || mpfr::isEqualFuzzy(a,b,eps);
  }
 
  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()); }
 
  // Specialize GEBP kernel and traits for mpreal (no need for peeling, nor complicated stuff)
  // This also permits to directly call mpfr's routines and avoid many temporaries produced by mpreal
    template<>
    class gebp_traits<mpfr::mpreal, mpfr::mpreal, false, false>
    {
    public:
      typedef mpfr::mpreal ResScalar;
      enum {
        Vectorizable = false,
        LhsPacketSize = 1,
        RhsPacketSize = 1,
        ResPacketSize = 1,
        NumberOfRegisters = 1,
        nr = 1,
        mr = 1,
        LhsProgress = 1,
        RhsProgress = 1
      };
      typedef ResScalar LhsPacket;
      typedef ResScalar RhsPacket;
      typedef ResScalar ResPacket;
      typedef LhsPacket LhsPacket4Packing;
      
    };
 
 
 
    template<typename Index, typename DataMapper, bool ConjugateLhs, bool ConjugateRhs>
    struct gebp_kernel<mpfr::mpreal,mpfr::mpreal,Index,DataMapper,1,1,ConjugateLhs,ConjugateRhs>
    {
      typedef mpfr::mpreal mpreal;
 
      EIGEN_DONT_INLINE
      void operator()(const DataMapper& res, const mpreal* blockA, const mpreal* blockB, 
                      Index rows, Index depth, Index cols, const mpreal& alpha,
                      Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0)
      {
        if(rows==0 || cols==0 || depth==0)
          return;
 
        mpreal  acc1(0,mpfr_get_prec(blockA[0].mpfr_srcptr())),
                tmp (0,mpfr_get_prec(blockA[0].mpfr_srcptr()));
 
        if(strideA==-1) strideA = depth;
        if(strideB==-1) strideB = depth;
 
        for(Index i=0; i<rows; ++i)
        {
          for(Index j=0; j<cols; ++j)
          {
            const mpreal *A = blockA + i*strideA + offsetA;
            const mpreal *B = blockB + j*strideB + offsetB;
            
            acc1 = 0;
            for(Index k=0; k<depth; k++)
            {
              mpfr_mul(tmp.mpfr_ptr(), A[k].mpfr_srcptr(), B[k].mpfr_srcptr(), mpreal::get_default_rnd());
              mpfr_add(acc1.mpfr_ptr(), acc1.mpfr_ptr(), tmp.mpfr_ptr(),  mpreal::get_default_rnd());
            }
            
            mpfr_mul(acc1.mpfr_ptr(), acc1.mpfr_srcptr(), alpha.mpfr_srcptr(), mpreal::get_default_rnd());
            mpfr_add(res(i,j).mpfr_ptr(), res(i,j).mpfr_srcptr(), acc1.mpfr_srcptr(),  mpreal::get_default_rnd());
          }
        }
      }
    };
  } // end namespace internal
}
 
#endif // EIGEN_MPREALSUPPORT_MODULE_H