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

/**
  * \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/mpfr">MPFR C++</a>
  * library which itself is built upon <a href="http://www.mpfr.org/">MPFR</a>/<a href="http://gmplib.org/">GMP</a>.
  *
  * You can find a copy of MPFR C++ that is known to be compatible in the unsupported/test/mpreal folder.
  *
  * Here is an example:
  *
\code
#include <iostream>
#include <Eigen/MPRealSupport>
#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,
    RequireInitialization = 1,
    ReadCost = 10,
    AddCost = 10,
    MulCost = 40
  };

  typedef mpfr::mpreal Real;
  typedef mpfr::mpreal NonInteger;

  inline static Real highest(long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::maxval(Precision); }
  inline static Real lowest(long Precision = mpfr::mpreal::get_default_prec()) { return -mpfr::maxval(Precision); }

  // Constants
  inline static Real Pi(long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_pi(Precision); }
  inline static Real Euler(long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_euler(Precision); }
  inline static Real Log2(long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_log2(Precision); }
  inline static Real Catalan(long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_catalan(Precision); }

  inline static Real epsilon(long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::machine_epsilon(Precision); }
  inline static Real epsilon(const Real &x) { return mpfr::machine_epsilon(x); }

  inline static Real dummy_precision() {
    unsigned int 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 {
    nr = 2, // must be 2 for proper packing...
    mr = 1,
    WorkSpaceFactor = nr,
    LhsProgress = 1,
    RhsProgress = 1
  };
};

template<typename Index, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
struct gebp_kernel<mpfr::mpreal, mpfr::mpreal, Index, mr, nr, ConjugateLhs, ConjugateRhs> {
  typedef mpfr::mpreal mpreal;

  EIGEN_DONT_INLINE
  void operator()(mpreal *res,
                  Index resStride,
                  const mpreal *blockA,
                  const mpreal *blockB,
                  Index rows,
                  Index depth,
                  Index cols,
                  mpreal alpha,
                  Index strideA = -1,
                  Index strideB = -1,
                  Index offsetA = 0,
                  Index offsetB = 0,
                  mpreal * /*unpackedB*/ = 0) {
    mpreal acc1, acc2, tmp;

    if (strideA == -1) strideA = depth;
    if (strideB == -1) strideB = depth;

    for (Index j = 0; j < cols; j += nr) {
      Index actual_nr = (std::min<Index>)(nr, cols - j);
      mpreal *C1 = res + j * resStride;
      mpreal *C2 = res + (j + 1) * resStride;
      for (Index i = 0; i < rows; i++) {
        mpreal *B = const_cast<mpreal *>(blockB) + j * strideB + offsetB * actual_nr;
        mpreal *A = const_cast<mpreal *>(blockA) + i * strideA + offsetA;
        acc1 = 0;
        acc2 = 0;
        for (Index k = 0; k < depth; k++) {
          mpfr_mul(tmp.mpfr_ptr(), A[k].mpfr_ptr(), B[0].mpfr_ptr(), mpreal::get_default_rnd());
          mpfr_add(acc1.mpfr_ptr(), acc1.mpfr_ptr(), tmp.mpfr_ptr(), mpreal::get_default_rnd());

          if (actual_nr == 2) {
            mpfr_mul(tmp.mpfr_ptr(), A[k].mpfr_ptr(), B[1].mpfr_ptr(), mpreal::get_default_rnd());
            mpfr_add(acc2.mpfr_ptr(), acc2.mpfr_ptr(), tmp.mpfr_ptr(), mpreal::get_default_rnd());
          }

          B += actual_nr;
        }

        mpfr_mul(acc1.mpfr_ptr(), acc1.mpfr_ptr(), alpha.mpfr_ptr(), mpreal::get_default_rnd());
        mpfr_add(C1[i].mpfr_ptr(), C1[i].mpfr_ptr(), acc1.mpfr_ptr(), mpreal::get_default_rnd());

        if (actual_nr == 2) {
          mpfr_mul(acc2.mpfr_ptr(), acc2.mpfr_ptr(), alpha.mpfr_ptr(), mpreal::get_default_rnd());
          mpfr_add(C2[i].mpfr_ptr(), C2[i].mpfr_ptr(), acc2.mpfr_ptr(), mpreal::get_default_rnd());
        }
      }
    }
  }
};

} // end namespace internal
}

#endif // EIGEN_MPREALSUPPORT_MODULE_H
