多倍長ライブラリMPIRで複素数std::complexを扱う方法(C++)

はてな匿名ダイアリー

2010-05-10

■多倍長ライブラリMPIRで複素数 std::complexを扱う方法(C++)

C++プログラムで多倍長な複素数を使いたかったのでメモ．

C++において多倍長整数，多倍長浮動小数点を扱うライブラリとしてGMPが有名だ．

Visual C++を使っているのなら，GMPと互換のあるMPIRが導入しやすくて良い．

以下，VC10 (Visual Studio 2010)環境での話．

このMPIRはどうもC++の複素数クラス std::complexと相性が悪いようだ．

以下のコードがエラーとなってしまう．

#include <complex&gt; 
#include <iostream&gt; 
#include <iomanip&gt; 
#include <mpirxx.h&gt; 

int main()
{
	std::complex<mpf_class&gt; a(1.0, 2.0);
	std::complex<mpf_class&gt; b(0.0, 1.0);
	std::complex<mpf_class&gt; c1 = a + b;
	std::complex<mpf_class&gt; c2 = a - b;
	std::complex<mpf_class&gt; c3 = a * b;
//	std::complex<mpf_class&gt; c4 = a / b;	// error
	
	std::cout << "a  =" << a  << std::endl
	          << "b  =" << b  << std::endl
	          << "c1 =" << c1 << std::endl
	          << "c2 =" << c2 << std::endl
	          << "c3 =" << c3 << std::endl;
//	std::cout << "c4 =" << c4 << std::endl;
	return 0;
}

「operator/=」から呼び出される「_Div」内部でエラーが出る．

そこでテンプレートの特殊化をする．

#pragma once
#include <mpirxx.h&gt;
#if defined(_MSC_VER) &amp;amp;&amp;amp; (_MSC_VER == 1500)      /* VC9 (Visual Studio 2008) */
  #pragma comment(lib, "C:\\lib\\MPIR\\vc9\\mpirxx.lib")
  #pragma comment(lib, "C:\\lib\\MPIR\\vc9\\mpir.lib")
#else if defined(_MSC_VER) &amp;amp;&amp;amp; (_MSC_VER == 1600) /* VC10 (Visual Studio 2010) */
  #pragma comment(lib, "C:\\lib\\MPIR\\vc10\\mpirxx.lib")
  #pragma comment(lib, "C:\\lib\\MPIR\\vc10\\mpir.lib")
#endif

#include<complex&gt;
namespace std {
	template<&gt;
	complex<mpf_class&gt;&amp;amp; complex<mpf_class&gt;::operator/=(const complex<mpf_class&gt;&amp;amp; _Right)
	{	// divide by other complex   //this-&gt;_Div(_Right);
		mpf_class _Rightreal = (mpf_class)_Right.real();
		mpf_class _Rightimag = (mpf_class)_Right.imag();

		mpf_class bunbo = _Rightreal * _Rightreal + _Rightimag * _Rightimag;
		mpf_class re = ( this-&gt;real() * _Rightreal + this-&gt;imag() * _Rightimag ) / bunbo;
		mpf_class im = ( this-&gt;imag() * _Rightreal - this-&gt;real() * _Rightimag ) / bunbo;

		this-&gt;real(re);
		this-&gt;imag(im);
		return (*this);
	}
}