Динамическая матрица Eigen3 с boost :: multiprecision :: mpfr_float

Я хотел бы сделать матрицы и использовать их, используя библиотеку Eigen3, с моим типом чисел Boost.Multiprecision — оберткой mpfr_float. Я могу сделать матрицы просто отлично, но работать с ними не удастся для всего, что я пробовал, кроме добавления матрицы. Простое умножение двух одинаковых матриц приводит к мусорным результатам!

Вот MWE:

#include <eigen3/Eigen/Dense>
#include <eigen3/Eigen/LU>
#include <boost/multiprecision/mpfr.hpp>
#include <iostream>

namespace Eigen{

using boost::multiprecision::mpfr_float;
template<> struct NumTraits<boost::multiprecision::mpfr_float>
{

typedef boost::multiprecision::mpfr_float Real;
typedef boost::multiprecision::mpfr_float NonInteger;
typedef boost::multiprecision::mpfr_float Nested;
enum {
IsComplex = 0,
IsInteger = 0,
IsSigned = 1,
RequireInitialization = 1,
ReadCost = 20, //these have no impact
AddCost = 30,
MulCost = 40
};inline static Real highest() { // these seem to have no impact
return (mpfr_float(1) - epsilon()) * pow(mpfr_float(2),mpfr_get_emax()-1);
}

inline static Real lowest() {
return -highest();
}

inline static Real dummy_precision(){
return pow(mpfr_float(10),-int(mpfr_float::default_precision()-3));
}

inline static Real epsilon(){
return pow(mpfr_float(10),-int( mpfr_float::default_precision()));
}
//http://www.manpagez.com/info/mpfr/mpfr-2.3.2/mpfr_31.php
};
} // namespace eigenint main()
{
int size = 10;

typedef Eigen::Matrix<boost::multiprecision::mpfr_float, Eigen::Dynamic, Eigen::Dynamic> mp_matrix;
mp_matrix A = mp_matrix::Identity(size, size);
std::cout << A * A << std::endl; // produces nan's every other row!!!

return 0;
}

Он производит матрицу тождественности просто отлично, но на моей машине, используя последние версии (и другие) зависимостей для этого кода, распространяемые на домашнем пиве (Boost 1.57, Eigen 3.2.4), моя программа генерирует все остальные строки NaN в матрице :

  1   0   0   0   0   0   0   0   0   0
nan nan nan nan nan nan nan nan nan nan
0   0   1   0   0   0   0   0   0   0
nan nan nan nan nan nan nan nan nan nan
0   0   0   0   1   0   0   0   0   0
nan nan nan nan nan nan nan nan nan nan
0   0   0   0   0   0   1   0   0   0
nan nan nan nan nan nan nan nan nan nan
0   0   0   0   0   0   0   0   1   0
nan nan nan nan nan nan nan nan nan nan

Нечетные размеры матриц дают два ряда в нижней части …

Кажется, это не зависит от точности по умолчанию или подробностей NumTraits struct я определяю, или даже определяю ли я одну. Я могу унаследовать от GenericTraits<mpfr_float>, или нет; я могу сказать RequireInitialization = 1, или же 0, Я получаю NaN’s. Если я пытаюсь инвертировать LU для решения системы, возвращаемая матрица полностью NaN. Если размер матрицы 1×1, я даже получаю один NaN из умножения матрицы. Изменение различных статических функций также не оказывает влияния.

Я чувствую, что самая странная часть в том, что если я определяю пользовательский комплексный класс (не std :: complex, по причинам потери данных), с типом mpfr_float в качестве базового типа для реальной и мнимой частей, я действительно получаю функциональные матрицы.

редактировать: вот NumTraits сложного типа:

/**
\brief this templated struct permits us to use the Float type in Eigen matrices.
*/
template<> struct NumTraits<mynamespace::complex> : NumTraits<boost::multiprecision::mpfr_float> // permits to get the epsilon, dummy_precision, lowest, highest functions
{
typedef boost::multiprecision::mpfr_float Real;
typedef boost::multiprecision::mpfr_float NonInteger;
typedef mynamespace::complex Nested;
enum {
IsComplex = 1,
IsInteger = 0,
IsSigned = 1,
RequireInitialization = 1, // yes, require initialization, otherwise get crashes
ReadCost = 2 * NumTraits<Real>::ReadCost,
AddCost = 2 * NumTraits<Real>::AddCost,
MulCost = 4 * NumTraits<Real>::MulCost + 2 * NumTraits<Real>::AddCost
};
};

Вот сложный класс, который я написал:

#include <boost/multiprecision/mpfr.hpp>
#include <boost/multiprecision/random.hpp>#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>
#include <boost/serialization/split_member.hpp>

#include <eigen3/Eigen/Core>

#include <assert.h>
namespace mynamespace {
using boost::multiprecision::mpfr_float;

class complex {

private:

mpfr_float real_, imag_; ///< the real and imaginary parts of the complex number// let the boost serialization library have access to the private members of this class.
friend class boost::serialization::access;

template<class Archive>
void save(Archive & ar, const unsigned int version) const {
// note, version is always the latest when saving
ar & real_;
ar & imag_;
}/**
\brief load method for archiving a bertini::complex
*/
template<class Archive>
void load(Archive & ar, const unsigned int version) {
ar & real_;
ar & imag_;
}

BOOST_SERIALIZATION_SPLIT_MEMBER()public:

complex():real_(), imag_(){}

complex(double re) : real_(re), imag_("0.0"){}

complex(const mpfr_float & re) : real_(re), imag_("0.0"){}

complex(const std::string & re) : real_(re), imag_("0.0"){}

complex(const mpfr_float & re, const mpfr_float & im) : real_(re), imag_(im) {}

complex(double re, double im) : real_(re), imag_(im) {}

complex(const std::string & re, const std::string & im) : real_(re), imag_(im) {}

complex(const mpfr_float & re, const std::string & im) : real_(re), imag_(im) {}

complex(const std::string & re, const mpfr_float & im) : real_(re), imag_(im) {}

complex(complex&& other) : complex() {
swap(*this, other);
}

complex(const complex & other) : real_(other.real_), imag_(other.imag_) {}

friend void swap(complex& first, complex& second)  {
using std::swap;
swap(first.real_,second.real_);
swap(first.imag_,second.imag_);
}

complex& operator=(complex other) {
swap(*this, other);
return *this;
}

mpfr_float real() const {return real_;}

mpfr_float imag() const {return imag_;}

void real(const mpfr_float & new_real){real_ = new_real;}

void imag(const mpfr_float & new_imag){imag_ = new_imag;}

void real(const std::string & new_real){real_ = mpfr_float(new_real);}

void imag(const std::string & new_imag){imag_ = mpfr_float(new_imag);}complex& operator+=(const complex & rhs) {
real_+=rhs.real_;
imag_+=rhs.imag_;
return *this;
}

complex& operator-=(const complex & rhs) {
real_-=rhs.real_;
imag_-=rhs.imag_;
return *this;
}

complex& operator*=(const complex & rhs) {
mpfr_float a = real_*rhs.real_ - imag_*rhs.imag_; // cache the real part of the result
imag_ = real_*rhs.imag_ + imag_*rhs.real_;
real_ = a;
return *this;
}

complex& operator/=(const complex & rhs) {
mpfr_float d = rhs.abs2();
mpfr_float a = real_*rhs.real_ + imag_*rhs.imag_; // cache the numerator of the real part of the result
imag_ = imag_*rhs.real_ - real_*rhs.imag_/d;
real_ = a/d;

return *this;
}

complex operator-() const
return complex(-real(), -imag());
}

mpfr_float abs2() const {
return pow(real(),2)+pow(imag(),2);
}

mpfr_float abs() const {
return sqrt(abs2());
}

mpfr_float arg() const {
return boost::multiprecision::atan2(imag(),real());
}

mpfr_float norm() const {
return abs2();
}

complex conj() const {
return complex(real(), -imag());
}

void precision(unsigned int prec) {
real_.precision(prec);
imag_.precision(prec);
}

unsigned int precision() const {
assert(real_.precision()==imag_.precision());
return real_.precision();
}

friend std::ostream& operator<<(std::ostream& out, const complex & z) {
out << "(" << z.real() << "," << z.imag() << ")";
return out;
}

friend std::istream& operator>>(std::istream& in, complex & z) {
std::string gotten;
in >> gotten;

if (gotten[0]=='(') {
if (*(gotten.end()-1)!=')') {
in.setstate(std::ios::failbit);
z.real("NaN");
z.imag("NaN");
return in;
}
else{
// try to find a comma in the string.
size_t comma_pos = gotten.find(",");

// if the second character, have no numbers in the real part.
// if the second to last character, have no numbers in the imag part.

if (comma_pos!=std::string::npos){
if (comma_pos==1 || comma_pos==gotten.size()-2) {
in.setstate(std::ios::failbit);
z.real("NaN");
z.imag("NaN");
return in;
}
else{
z.real(gotten.substr(1, comma_pos-1));
z.imag(gotten.substr(comma_pos+1, gotten.size()-2 - (comma_pos)));
return in;
}
}
// did not find a comma
else{
z.real(gotten.substr(1,gotten.size()-2));
z.imag("0.0");
return in;
}

}
}
else{
z.real(gotten);
z.imag("0.0");
return in;
}
}
}; // end declaration of the mynamespace::complex number class

inline complex operator+(complex lhs, const complex & rhs){
lhs += rhs;
return lhs;
}

inline complex operator+(complex lhs, const mpfr_float & rhs) {
lhs.real(lhs.real()+rhs);
return lhs;
}

inline complex operator+(const mpfr_float & lhs, complex rhs) {
return rhs+lhs;
}

inline complex operator-(complex lhs, const complex & rhs){
lhs -= rhs;
return lhs;
}

inline complex operator-(complex lhs, const mpfr_float & rhs) {
lhs.real(lhs.real()-rhs);
return lhs;
}

inline complex operator-(const mpfr_float & lhs, complex rhs) {
rhs.real(lhs - rhs.real());
return rhs;
}

inline complex operator*(complex lhs, const complex & rhs){
lhs *= rhs;
return lhs;
}

inline complex operator*(complex lhs, const mpfr_float & rhs) {
lhs.real(lhs.real()*rhs);
lhs.imag(lhs.imag()*rhs);
return lhs;
}

inline complex operator*(const mpfr_float & lhs, complex rhs) {
return rhs*lhs; // it commutes!
}

inline complex operator/(complex lhs, const complex & rhs){
lhs /= rhs;
return lhs;
}

inline complex operator/(complex lhs, const mpfr_float & rhs) {
lhs.real(lhs.real()/rhs);
lhs.imag(lhs.imag()/rhs);
return lhs;
}

inline complex operator/(const mpfr_float & lhs, const complex & rhs) {
mpfr_float d = rhs.abs2();
return complex(lhs*rhs.real()/d, -lhs*rhs.imag()/d);
}

inline mpfr_float real(const complex & z) {
return z.real();
}

inline mpfr_float imag(const complex & z) {
return z.imag();
}

inline complex conj(const complex & z) {
return z.conj();
}

inline mpfr_float abs2(const complex & z) {
return z.abs2();
}

inline mpfr_float abs(const complex & z) {
return boost::multiprecision::sqrt(abs2(z));
}

inline mpfr_float arg(const complex & z) {
return boost::multiprecision::atan2(z.imag(),z.real());
}

inline complex inverse(const complex & z) {
mpfr_float d = z.abs2();

return complex(z.real()/d, -z.imag()/d);
}

inline complex square(const complex & z) {
return complex(z.real()*z.real() - z.imag()*z.imag(), mpfr_float("2.0")*z.real()*z.imag());
}

inline complex pow(const complex & z, int power) {
if (power < 0) {
return pow(inverse(z), -power);
}
else if (power==0)
return complex("1.0","0.0");
else if(power==1)
return z;
else if(power==2)
return z*z;
else if(power==3)
return z*z*z;
else {
unsigned int p(power);
complex result("1.0","0.0"), z_to_the_current_power_of_two = z;
// have copy of p in memory, can freely modify it.
do {
if ( (p & 1) == 1 ) { // get the lowest bit of the number
result *= z_to_the_current_power_of_two;
}
z_to_the_current_power_of_two *= z_to_the_current_power_of_two; // square z_to_the_current_power_of_two
} while (p  >>= 1);

return result;
}
}

inline complex polar(const mpfr_float & rho, const mpfr_float & theta) {
return complex(rho*cos(theta), rho*sin(theta));
}

inline complex sqrt(const complex & z) {
return polar(sqrt(abs(z)), arg(z)/2);
}
} // re: namespace

Что я делаю неправильно? Что я могу сделать с Eigen / NumTraits / etc, чтобы заставить матричные операции работать правильно?