
//	Module:		Complex
//	Version:	2.00  28-Oct-1989
//	Language:	C++ 2.0;  Environ: Any;  Compilers: Zortech C++ 2.01
//	Purpose:	Provides the class "Complex" for C++ programs.  The
//			majority of the class is implemented inline for
//			efficiency.  Only the division, power, and i/o methods
//			are actual functions.
//	Written by:	Scott Robert Ladd, 705 West Virginia, Gunnison CO 81230
//			BBS (303)641-6438; FidoNet 1:104/708

#include <math.h>
#include <stdlib.h>
#include <stream.hpp>
#include "complex.hpp"

static void	defaultHandler ();
void ( * Complex::errorHandler) () = defaultHandler;

static
void
defaultHandler (
) {
	cout << "\aERROR in complex object: DIVIDE BY ZERO\n";
	exit ( 1);
}

Complex
operator/ (
	const Complex &	c1,
	const Complex &	c2
) {
	Complex		result;
	double		den;
	den = norm ( c2);
	if ( den != 0.0) {
		result.re = ( c1.re * c2.re + c1.im * c2.im) / den;
		result.im = ( c1.im * c2.re - c1.re * c2.im) / den;
	} else
		Complex::errorHandler ();
	return result;
}

Complex
Complex::operator/= (
	const Complex &	c
) {
	double		den;
	den = norm ( c);
	if ( den != 0.0) {
		double		holdr = re;
		re = ( re * c.re + im * c.im) / den;
		im = ( im * c.re - holdr * c.im) / den;
	} else
		Complex::errorHandler ();
	return *this;
}

Complex
pow (
	const Complex &	base,
	const Complex &	power
) {
	Complex		result;
	if ( power.re == 0.0 && power.im == 0.0) {
		result.re = 1.0;
		result.im = 0.0;
	} else
		if ( base.re != 0.0 || base.im != 0.0)
			result = exp ( log ( base) * power);
		else
			Complex::errorHandler ();
	return result;
}

Complex
sqrt (
	const Complex &	c
) {
	Complex		result;
	double		r, i, ratio, w;
	if ( c.re != 0.0 || c.im != 0.0) {
		r = ( c.re < 0.0 ? -c.re : c.re);
		i = ( c.im < 0.0 ? -c.im : c.im);
		if ( r > i) {
			ratio = i / r;
			w = sqrt ( r) * sqrt ( 0.5 * ( 1.0 +
				sqrt ( 1.0 + ratio * ratio)));
		} else {
			ratio = r / i;
			w = sqrt ( i) * sqrt ( 0.5 * ( ratio +
				sqrt ( 1.0 + ratio * ratio)));
		}
		if ( c.re > 0.0) {
			result.re = w;
			result.im = c.im / ( 2.0 * w);
		} else {
			result.im = ( c.im > 0.0 ? w : -w);
			result.re = c.im / ( 2.0 * result.im);
		}
	} else
		Complex::errorHandler ();
	return result;
}

ostream &
operator<< (
	ostream &	os,
	const Complex &	c
) {
	os << form ( "(%+1g%+1gi)", c.re, c.im);
	return os;
}

istream &
operator>> (
	istream &	is,
	Complex &	c
) {
	char		ch;
	c.re = 0.0;
	c.im = 0.0;
	is >> ch;
	if ( ch == '(') {
		is >> c.re >> ch;
		if ( ch == ',')
			is >> c.im >> ch;
		if ( ch != ')')
			is.clear ( _bad);
	} else {
		is.putback ( ch);
		is >> c.re;
	}
	return is;
}

