complex math | GameDev.net

Afaik · 2002-12-29T00:53:00+00:00

Complex math as in imaginary numbers. Does MSVC 6++ support complex math? If not, where can I find a C/C++ complex math library?

#include <complex>
#include <iostream>

int main() {
std::complex<double> z1, z2;
std::cout << "z1 -> ";
std::cin >> z1;
std::cout << "z2 -> ";
std::cin >> z2;

std::complex<double> z3 =
z1 * std::exp(z2);
std::cout << "Result: " << z3
<< std::endl;
}

Here is a complex library from my CAS. You can rewrite it as a template, is you wish...

    
/*
 *
 *	Math Studio NG, (C)2002 David Sehnal
 *	
 *	Complex Arithmetics, Machine Precision (double)
 *
 */
#pragma once

// compare 2 machine precision numbers
inline CN(const double &a, const double &b)
{
       return fabs(a-b) < /* epsilon */ 1e-15;
}

class Complex
{
public:
	double re, im;

	BOOL IsReal() const { return CN(im, 0.0); }

	operator double() const { return re; }
	operator int() const { return (int)re; }
	operator float() const { return (float)re; }

	void operator =(const Complex &r) { re = r.re, im = r.im; }
	void operator =(const double r) { re = r, im = 0.0; }
	void operator =(const int r) { re = (double)r, im = 0.0; }

	Complex operator +(const Complex &r) const { return Complex(re + r.re, im + r.im); }
	Complex operator -(const Complex &r) const { return Complex(re - r.re, im - r.im); }
	inline Complex operator *(const Complex &r) const;
	inline Complex operator /(const Complex &r) const;

	Complex operator +(const double r) const { return Complex(re + r, im); }
	Complex operator -(const double r) const { return Complex(re - r, im); }
	inline Complex operator *(const double r) const { return Complex(re * r, im * r); }
	inline Complex operator /(const double  r) const { return Complex(re / r, im / r); };

	Complex operator +(const int r) const { return Complex(re + r, im); }
	Complex operator -(const int r) const { return Complex(re - r, im); }
	inline Complex operator *(const int r) const { return Complex(re * r, im * r); }
	inline Complex operator /(const int  r) const { return Complex(re / r, im / r); }

	Complex operator -() { return Complex(-re, -im); } 
	
	void operator +=(const Complex &r) { re += r.re, im += r.im; }
	void operator -=(const Complex &r) { re -= r.re, im -= r.im; }
	void operator *=(const Complex &r) { Complex tmp = *this * r; re = tmp.re, im = tmp.im; }
	void operator /=(const Complex &r) { Complex tmp = *this / r; re = tmp.re, im = tmp.im; }

	void operator +=(const double r) { re += r, im += r; }
	void operator -=(const double r) { re -= r, im -= r; }
	void operator *=(const double r) { re *= r, im *= r; }
	void operator /=(const double r) { re /= r, im /= r; }

	void operator +=(const int r) { re += r, im += r; }
	void operator -=(const int r) { re -= r, im -= r; }
	void operator *=(const int r) { re *= r, im *= r; }
	void operator /=(const int r) { re /= r, im /= r; }

	BOOL operator ==(const Complex &r) const { return CN(re, r.re) && CN(im, r.im); }
	BOOL operator !=(const Complex &r) const { return !CN(re, r.re) || !CN(im, r.im); }
	BOOL operator <(const Complex &r) const { return re < r.re; }
	BOOL operator >(const Complex &r) const { return re > r.re; }
	BOOL operator <=(const Complex &r) const { return re <= r.re; }
	BOOL operator >=(const Complex &r) const { return re >= r.re; }

	Complex(double _re = 0.0, double _im = 0.0) : re(_re), im(_im) {}
	~Complex() { }
};

inline double Re(const Complex &c) { return c.re; }
inline double Im(const Complex &c) { return c.im; }

inline BOOL CN(const Complex &v1, const Complex &v2)
{
	return fabs(v1.re - v2.re) <= 1e-15 && fabs(v1.im - v2.im) <= 1e-15;
}
inline BOOL IsComplex(const Complex &) { return TRUE; }

inline int sign(const Complex &val)
{ 
	if (val == Complex()) return 0;
	if (val.re > 0.0 || (CN(val.re, 0.0) && val.im > 0.0)) return 1;
	return -1;
}
inline Complex frac(Complex val) { return val.re - floor(val.re); }

Complex Complex::operator *(const Complex &r) const
{
	if (IsReal() && r.IsReal())
	{
		return re * r.re;
	}

	double ac = re * r.re;
	double bd = im * r.im;
	
	Complex ret;

	ret.re = ac - bd;
	ret.im = (re + im) * (r.re + r.im) - ac - bd;

	return ret;
}

Complex Complex::operator /(const Complex &r) const
{
	if (IsReal() && r.IsReal())
	{
		return re / r.re;
	}


	Complex ret;

	if (abs(r.re) > abs(r.im))
	{
		double dc = r.im/r.re;
		double l = double(1.)/(r.re + r.im*dc);

		ret.re = (re + im*dc)*l;
		ret.im = (im - re*dc)*l;

		return ret;
	}

	double cd = r.re/r.im;
	double l = double(1.)/(r.im + r.re*cd);

	ret.re = (im + re*cd)*l;
	ret.im = (im*cd - re)*l;

	return ret;
}

inline Complex operator +(const double l, const Complex &r) 
{ 
	return r + l; 
}

inline Complex operator -(const double l, const Complex &r) 
{ 
	return Complex(l - r.re, r.im); 
}

inline Complex operator *(const double l, const Complex &r)
{
	return r * l;
}

Complex operator /(const double _l, const Complex &r)
{
	if (r.IsReal())
	{
		return _l / r.re;
	}

	Complex ret;

	if (abs(r.re) > abs(r.im))
	{
		double dc = r.im/r.re;
		double l = double(1.)/(r.re + r.im*dc);

		ret.re = _l*l;
		ret.im = -_l*dc*l;

		return ret;
	}

	double cd = r.re/r.im;
	double l = double(1.)/(r.im + r.re*cd);

	ret.re = _l*cd*l;
	ret.im = -_l*l;

	return ret;
}

inline Complex operator +(const int l, const Complex &r) 
{ 
	return r + l; 
}

inline Complex operator -(const int l, const Complex &r) 
{ 
	return Complex(l - r.re, r.im); 
}

inline Complex operator *(const int l, const Complex &r)
{
	return r * l;
}

Complex operator /(const int __l, const Complex &r)
{
	double _l = (double) __l;

	if (r.IsReal())
	{
		return _l / r.re;
	}

	Complex ret;

	if (abs(r.re) > abs(r.im))
	{
		double dc = r.im/r.re;
		double l = double(1.)/(r.re + r.im*dc);

		ret.re = _l*l;
		ret.im = -_l*dc*l;

		return ret;
	}

	double cd = r.re/r.im;
	double l = double(1.)/(r.im + r.re*cd);

	ret.re = _l*cd*l;
	ret.im = -_l*l;

	return ret;
}

double abs(const Complex &x)
{
	if (x.IsReal()) return abs(x.re);

	
	double ret;

	if (abs(x.re) >= abs(x.im))
	{
		double t = x.im/x.re;

		ret = abs(x.re) * sqrt(1. + t*t);

		return ret;
	}	

	double t = (x.re/x.im);

	ret = abs(x.im) * sqrt(1. + t*t);

	return ret;
}

inline double arg(const Complex &x)
{
	if (x.IsReal())
		return (x.re >= 0.0) ? 0.0 : PI;

	return atan2(x.im, x.re);
}

Complex sqrt(const Complex &x) 
{
	if (x.IsReal())
	{
		if (x.re >= 0.0) return sqrt(x.re);
		else return Complex(0., sqrt(x.re));
	}

	double aa = abs(x.re), ab = abs(x.im);

	double w;
	
	if (aa >= ab)
	{
		w = sqrt(aa * (1. + sqrt(1. + ((x.im * x.im)/ (x.re * x.re))))/2.);
	}
	else
	{
		double adb = x.re/x.im;
		w = sqrt(ab * (adb + sqrt(1. + (adb * adb)))/2.);
	}

	if (CN(w, 0.)) return 0;

	else if (x.re >= 0.)
	{
		return Complex(w, x.im / (2. * w));
	}
	else if (x.re < 0. && x.im >= 0.)
	{
		return Complex(ab / (2. * w), w);
	}
	else
	{
		return Complex(ab / (2. * w), -w);
	}
}


Complex log(const Complex &x)
{
	if (x.IsReal() && x.re >= 0.0) return log(x.re);

	return Complex(log(abs(x)), arg(x));
}

Complex log10(const Complex &x)
{
	if (x.IsReal() && x.re >= 0.0) return log10(x.re);

	return Complex(log10(abs(x)), arg(x));
}

Complex exp(const Complex &x)
{
	if (x.IsReal()) return exp(x.re);

	double l = exp(x.re);

	return Complex(l * cos(x.im), l * sin(x.im));
}


Complex pow(const Complex &x, const Complex &y)
{
	if (x.IsReal() && y.IsReal() && x.re >= 0.0) return pow(x.re, y.re);

	return exp(log(x) * y);
}

---------
Dave

[edited by - Dave007 on December 30, 2002 1:23:23 PM]

The STL and the STL math libs have complex numbers, it''s std::complex as stated previously. You don''t have to know how they work, and you don''t need to worry about templates because the default template argument provided gives you std::complex< typename T = double >, so you get double precision complexes by default (i.e. std::complex == std::complex< double > ). All the maths functions which apply to complex numbers have overloads for complex types, I believe. I don''t think there''s a sqrt func which takes a real and returns a complex though (because overloads can''t differ only by return type). But you can cast a real as a std::complex and then you''ll get a compleex returned.

Templates are good for this because you can specify the precision required without having to rewrite code using copy & paste when you decide you have too much (or not enough) precision.

"Most people think, great God will come from the sky, take away everything, and make everybody feel high" - Bob Marley

Topic Locked