c++ - boost::multiprecision::mpfr_float を使用した Eigen3 動的マトリックス

Question

行列を作成し、Eigen3 ライブラリを使用して使用したいと思います。数値型は Boost.Multiprecision の mpfr_float ラッパーです。行列をうまく作成できますが、行列の追加を除いて、私が試したすべての操作で失敗します。2 つの恒等行列を乗算するだけで、ガベージ結果が生成されます。

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 行おきに生成します。 :

  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

行列のサイズが奇数の場合、下部に 2 行の nan が生成されます...

NumTraitsこれは、デフォルトの精度、定義する構造体の詳細、または定義するかどうかにも依存していないようです。から継承できるGenericTraits<mpfr_float>かどうか。RequireInitialization = 1、またはと言え0ます。NaN を取得します。連立方程式を解くために LU 反転を試みると、返される行列は完全に NaN になります。行列のサイズが 1x1 の場合、行列の乗算から NaN が 1 つ得られます。さまざまな静的関数を変更しても影響はありません。

最も奇妙な部分は、mpfr_float を実部と虚部の基礎となる型としてカスタムの複雑なクラス (データ損失の理由から std::complex ではない) を定義すると、関数行列が得られることです。

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 / などに何ができますか?

Answer 1

この奇妙なnan動作は、boost/multiprecision/mpfr.hpp ファイルのバグが原因であり、Boost バージョン 1.56、1.57、およびおそらくそれ以前で発生しました。代入演算子は自己代入をテストしなかったため、数値は に設定されましたnan。

これが 1 行おきに発生した理由と、行列のサイズが奇数の場合の最後の行で発生した理由は、まだ不明です。ただし、ここでGitHub のコミットのように、自己割り当てのテストを追加すると、問題が解決します。メンテナからの公式パッチが近日中にリリースされますが、1.58 にはならない予定です。バグを見つけてくれた Boost ユーザーのメーリング リストの John に感謝します。

Boost のインストールでこのバグを修正したい場合boost/multiprecision/mpfr.hppは、参照代入演算子の内容 (リターンを含まない) を でラップするだけif(this != &o){...}です。

私が持っている 1 つの長引く質問は、Eigen の実装に関するものです。自己割り当ての原因は何ですか? それは一般的に問題ですか？