Динамическая матрица 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 eigen


int 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, которую я определяю, или даже от того, определяю ли я ее. Могу я унаследовать от GenericTraits<mpfr_float> или нет; Я могу сказать RequireInitialization = 1 или 0. Я получаю NaN. Если я попытаюсь инвертировать LU для решения системы, возвращенная матрица будет полностью NaN. Если размер матрицы равен 1x1, я даже получаю один NaN от умножения матриц. Изменение различных статических функций также не влияет.

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

edit: Вот 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 / и т. Д., Чтобы матричные операции работали правильно?