C++で複素数クラスを作成する

programming C++

STLのcomplexクラスについて

C++で複素数を扱うにはSTLで用意されているcomplexクラスを利用します。complexクラスはテンプレートで定義されているので、型を指定して扱うことができます。例えば、double型の複素数を扱いたい場合は、complexと記述します。
さらに、複素関数も下のように様々なものが用意されています。

#include <complex>
#include <iostream>
using namespace std;
int main(){
	const complex<double> i( 0.0, 1.0 );
	complex<double> z1( -1.0, 1.0 );
	complex<double> z2( 2.0, -2.0 );
	complex<double> w = polar( 1.0, 3.141592 / 2.0 );//極座標
	//複素数の基本関数
	cout << z1.real() << " , "<< z1.imag() << "i"<<endl;
	cout << arg(z1) << endl;
	cout << norm(z1) << endl;
	cout << conj(z1)<<endl;
	//複素数の演算
	cout << z1 + 1.0 << endl;
	cout << 1.0 + z1 << endl;
	cout << z1 * i << endl;
	cout << (z1 == z2) << endl;
	cout << (z1 != z2) << endl;
	//複素関数
	cout << exp(w) << endl;
	cout << sin(w) << endl;
	cout << cos(w) << endl;
	cout << sinh(w) << endl;
	cout << cosh(w) << endl;
	cout << log(w) << endl;
	return 0;
}

typedefを使ってcomplexなどを予め定義しておくと、多少ですがクラス名を短く書くことができます。

typedef complex<double> dcomplex;
typedef complex<int> icomplex;

dcomplex z(1.0,2.0);
z = dcomplex(0.0,1.0);

本格的に利用する場合は、既に用意されているcomplexクラスを使うことで事足りると思いますが、
今回は勉強のために自分で改めて"Complex"クラスを作成してみることにしてみます。

Complexクラスの作成

複素数の型はfloat型にしてみます。メンバ変数は実部と虚部の2つです。それぞれprivateメンバにして、外部から呼び出すためにこれらの値へのアクセッサーを用意することにします。

class Complex{
private:
	float re;
	float im;
public:
	// constructor
	Complex():re(0),im(0){}
	Complex(float re,float im):re(re),im(im){}
	Complex(float re):re(re),im(0){}
	// accessor
	float real()const{	return re;	}
	float imag()const{	return im;	}

代入と+,-演算子

代入は非constなインスタンスのみ、+や-の演算子はconstなインスタンスに対しても適用できます。

	// 代入
	Complex& Complex::operator=(const Complex& z){
		this->re=z.re;		this->im=z.im;
		return *this;
	}
	// +Complex,-Complex
	Complex operator+() const{ return *this; }
	Complex operator-() const{ return Complex(-re,-im); }

単項四則演算子の定義

"+=, -=, *=, /="など自身を書き換えるタイプの演算子をメンバ関数として定義します。
注意するべき点はそれぞれに対して、相手が複素数の場合と、実数の場合の2つのパターンを考えなければならないことです。

	//Complex+=Complex
	Complex& operator+=(const Complex& z){
		this->re += z.re;	this->im += z.im;
		return *this;
	}
	//Complex+=real float
	Complex& operator+=(const float& k){
		this->re += k;
		return *this;
	}
	//Complex-=Complex
	Complex& operator-=(const Complex& z){
		this->re -= z.re;	this->im -= z.im;
		return *this;
	}
	//Complex-=real float
	Complex& operator-=(const float& k){
		this->re -= k;
		return *this;
	}
	//Complex*=Complex
	Complex& operator*=(const Complex& z){
		float re_tmp = this->re*z.re - this->im*z.im;
		float im_tmp = this->re*z.im + this->im*z.re;
		this->re = re_tmp;	this->im = im_tmp;
		return *this;
	}
	//Complex*=real float
	Complex& operator*=(const float k){
		this->re *= k;this->im *= k;return *this;
	}
	//Complex/=Complex
	Complex& operator/=(const Complex& z){
		float norm_z = z.re*z.re+z.im*z.im;
		float re_tmp = (this->re * z.re + this->im * z.im) / norm_z;
		float im_tmp = (this->im * z.re - this->re * z.im) / norm_z;
		this->re = re_tmp;	this->im = im_tmp;
		return *this;
	}
	//Complex/=real float
	Complex& operator/=(const float& k){
		this->re /= k;	this->im /= k;
		return *this;
	}

二項演算子関数の定義

"+, -, *, /"の演算子をグローバル関数として定義します。
ちなみに複素数x,yの和差積商の演算は下のようになります。型として複素数クラスだけ定義する場合はどうしても商の計算が重くなりますね。

z = x + y
\mathbf{Re}[z] = \mathbf{Re}[x] + \mathbf{Re}[y] , \mathbf{Im}[z] = \mathbf{Im}[x] + \mathbf{Im}[y]
z = x - y
\mathbf{Re}[z] = \mathbf{Re}[x] - \mathbf{Re}[y] , \mathbf{Im}[z] = \mathbf{Im}[x] - \mathbf{Im}[y]
z = x * y
\mathbf{Re}[z] = \mathbf{Re}[x] * \mathbf{Re}[y] - \mathbf{Im}[x] * \mathbf{Im}[y]
\mathbf{Im}[z] = \mathbf{Re}[x] * \mathbf{Im}[y] + \mathbf{Im}[x] * \mathbf{Re}[y]
z = x / y
\mathbf{Re}[z] = \frac{\mathbf{Re}[x] * \mathbf{Re}[y] + \mathbf{Im}[x] * \mathbf{Im}[y]}{|x|^2+|y|^2}
\mathbf{Im}[z] = \frac{\mathbf{Im}[x] * \mathbf{Re}[y] - \mathbf{Re}[x] * \mathbf{Im}[y]}{|x|^2+|y|^2}
こちらも注意点として、演算する相手を考えなければなりません。例えば"+"の演算では、複素数+複素数,複素数+実数,実数+複素数の3パターンを考慮する必要があります。実部と虚部を計算して、直接インスタンスを返却しています。

//Complex+Complex
inline Complex operator+(const Complex& x,const Complex& y){
	return Complex(x.real() + y.real(),x.imag() + y.imag());
}
//Complex+real
inline Complex operator+(const Complex& x,const float& y){
	return Complex(x.real() + y, x.imag());
}
//real+Complex
inline Complex operator+(const float& x,const Complex& y){
	return Complex(x + y.real(), y.imag());
};
//Complex-Complex
inline Complex operator-(const Complex& x,const Complex& y){
	return Complex(x.real() - y.real(), x.imag() - y.imag());
}
//Complex-real
inline Complex operator-(const Complex& x,const float& y){
	return Complex(x.real() - y, x.imag());
}
//real-Complex 
inline Complex operator-(const float& x,const Complex& y){
	return Complex(x - y.real(), -y.imag());
}
//Complex*Complex
inline Complex operator*(const Complex& x,const Complex& y){
	return Complex (
		x.real() * y.real() - x.imag() * y.imag()	,\
		x.real() * y.imag() + x.imag() * y.real()	);
}
//Complex*float
inline Complex operator*(const Complex& x,float y){
	return Complex (
		x.real() * y	,\
		x.imag() * y	);
}
//float*Complex
inline Complex operator*(float x,const Complex& y){
	return Complex (
		x * y.real()	,\
		x * y.imag()	);
}
//Complex/Complex
inline Complex operator/(const Complex& x,const Complex& y){
	float norm_y = y.real()*y.real()+y.imag()*y.imag();
	return Complex (
		(x.real() * y.real() + x.imag() * y.imag()) / norm_y	,\
		(x.imag() * y.real() - x.real() * y.imag()) / norm_y	);
}
//Complex/float
inline Complex operator/(const Complex& x,float y){
	return Complex (
		x.real() / y	,\
		x.imag() / y	);
}
//float/Complex
inline Complex operator/(float x,const Complex& y){
	float norm_y = y.real() * y.real() + y.imag() * y.imag();
	return Complex (
		x * y.real() / norm_y		,\
		-x * y.imag() / norm_y	);
}

比較演算子

実部も虚部も等しいときのみ複素数は等しいと定義します。この演算子も同様に相手が実数の場合も考慮する必要があります。
比較演算子の定義で、特に浮動小数点で計算をすすめるような場合には誤差を考慮した比較というものを実装することもできます。
例えば、1のつもりでも実際には誤差が生じるために0.99994であるような場合、0.0001という範囲の誤差以下の場合は同じであるとみなすような比較演算子を定義することもできるでしょう。

inline bool operator==(const Complex& x,const Complex& y){
	return ( x.real() == y.real() && x.imag() == y.imag() );
}
inline bool operator==(const Complex& x,const float& y){
	return ( x.real() == y && x.imag() == 0 );
}
inline bool operator==(const float& x,const Complex& y){
	return ( x == y.real() && y.imag() == 0 );
}
inline bool operator!=(const Complex& x,const Complex& y){
	return ( x.real() != y.real() || x.imag() != y.imag() );
}
inline bool operator!=(const Complex& x,const float& y){
	return ( x.real() != y || x.imag() != 0 );
}
inline bool operator!=(const float& x,const Complex& y){
	return ( x != y.real() || y.imag() != 0 );
}

複素数に関する基本的な関数

複素数クラスを使って最低限の計算するために必要な関数を追加してみます。例えば長さ、偏角、共役複素数がそのような関数です。
STLではこのような関数はで定義されています。例えば自前のクラスを引数とするabs関数を、予め用意されたabs関数にさらにオーバーロードして定義することは可能です。
(ちなみにSTLのcomplexには極座標関数polarが定義されているので、返却値のみ異なるような同じ名前の極座標関数polarを作成することはできません。)

//絶対値
inline float abs(const Complex& z){
	return sqrt( z.real() * z.real() + z.imag() * z.imag() );
}
//ノルム
inline float norm(const Complex& z){
	return z.real() * z.real() + z.imag() * z.imag();
}
//偏角
inline float arg(const Complex& z){
	return atan2( z.imag(), z.real() );
}
//共役複素数
inline Complex conj(const Complex& z){
	return Complex( z.real(), -z.imag() );
}
//極座標による複素数の生成
inline Complex Polar( const float r, const float theta )
{
	return Complex(r*cos(theta),r*sin(theta));
}

これで複素数クラスはひとまず完成です。
STLでは代表的な複素関数が全て提供されていますので、ここでもいくつかの関数をSTLにならって実装してみます。

複素関数

inline Complex exp(const Complex& z)
{
	exp(z.real()) + Complex(0,1.0) * exp(z.imag());
}
inline Complex sinh(const Complex& z)
{
	(exp(z) - exp(-z)) / 2.0f;
}
inline Complex cosh(const Complex& z)
{
	(exp(z) + exp(-z)) / 2.0f;
}
inline Complex tanh(const Complex& z)
{
	(exp(z) - exp(-z)) / (exp(z) + exp(-z));
}


inline Complex sin(const Complex& z)
{
	const Complex i(0.0f,1.0f);
	- i * ( exp( i*z ) + exp( -i*z )) / 2.0f;
}
inline Complex cos(const Complex& z)
{
	const Complex i(0.0f,1.0f);
	(exp( i*z ) + exp( -i*z )) / 2.0f;
}
inline Complex tan(const Complex& z)
{
	const Complex i(0.0f,1.0f);
	- i * (exp( i*z ) - exp( -i*z )) / (exp( i*z ) + exp( -i*z ));
}

inline Complex log(const Complex& z){
//NANやInfinityなどの場合分けが必要
	log(z.abs()) + Complex(0,1.0) * z.arg();
}

iの使い方

STLのcomplexクラスや自作のComplexクラスを使って、数学の表記と同じように"i"を直接コードに書くことができると便利です。しかしプログラム言語の仕様として"***i"と書くと自動的に複素数型として認識してくれるような機能がなければ、ソースの中の全ての部分で一貫して"i"を虚数として使えるようにするのは難しいと思います。
"i"はfor文の中で利用したりしますので、字数が少なくて紛らわしいことが本質的な原因でしょう。
スコープが局所的で良いのなら簡単な方法として、iを利用する直前に定義してしまうのが手っ取り早いと思います。

	const Complex i(0.0f,1.0f);

ただし、forでいつも通りi,jを使っていて、forの中で虚数のiを使うような場合があると間違いの原因になりそうです。

	z *= i;
	for(int i=0;i<3;i++){
		z *= i;
	}