eigen.h

#pragma once
 
#include <gmpxx.h>
 
#include <Eigen/Core>
#include <Eigen/LU>
 
namespace Eigen {
template <>
struct NumTraits<mpq_class> : GenericNumTraits<mpq_class> {
  typedef mpq_class Real;
  typedef mpq_class NonInteger;
  typedef mpq_class Nested;
 
  static inline Real epsilon() { return 0; }
  static inline Real dummy_precision() { return 0; }
  static inline int digits10() { return 0; }
 
  enum {
    IsInteger = 0,
    IsSigned = 1,
    IsComplex = 0,
    RequireInitialization = 1,
    ReadCost = 6,
    AddCost = 150,
    MulCost = 100
  };
};
 
namespace internal {
 
template <>
struct cast_impl<mpq_class, double> {
  static double run(const mpq_class& x) { return x.get_d(); }
};
#if 0
template <>
struct scalar_score_coeff_op<mpq_class> {
  struct result_type {
    std::size_t len;
    explicit result_type(const int i = 0) : len(i) {}  // Eigen uses Score(0) and Score()
    explicit result_type(const mpq_class& q) : len(mpz_size(q.get_num_mpz_t()) + mpz_size(q.get_den_mpz_t()) - 1) {}
    std::strong_ordering operator<=>(const result_type& y) const {
      if (len == y.len) return std::strong_ordering::equal;
      return (len < y.len) ? std::strong_ordering::less : std::strong_ordering::greater;
    }
  };
  result_type operator()(const mpq_class& x) const { return result_type{x}; }
};
#endif
}  // namespace internal
}  // namespace Eigen
 
namespace delpi {
using EMatrix = Eigen::MatrixX<mpq_class>;
using EVector = Eigen::VectorX<mpq_class>;
using DMatrix = Eigen::MatrixXd;
using DVector = Eigen::VectorXd;
template <class T>
using Matrix = Eigen::MatrixX<T>;
template <class T>
using Vector = Eigen::VectorX<T>;
template <class T>
using VectorView = Eigen::Map<Eigen::VectorX<T>>;
template <class T>
using MatrixView = Eigen::Map<Eigen::MatrixX<T>>;
using Index = Eigen::Index;
using VectorD = Vector<double>;
using VectorQ = Vector<mpq_class>;
using MatrixD = Matrix<double>;
using MatrixQ = Matrix<mpq_class>;
 
}  // namespace delpi
 
#ifdef DELPI_INCLUDE_FMT
 
#include "delpi/util/logging.h"
 
OSTREAM_FORMATTER(delpi::Matrix<mpq_class>)
OSTREAM_FORMATTER(delpi::Matrix<double>)
 
#endif