/* complex.c:	Complex arithmetic routines.
 *
 * Written by Daniel Barrett.  100% PUBLIC DOMAIN. */

#include "decl.h"

/* Complex arithmetic functions Add(), Subtract(), Multiply() and Divide()
 * perform as their names indicate.  They perform their operation on their
 * first 2 arguments, and return the result in the third argument. */

Add(a, b, sum)
complex a, b, *sum;
{
	sum->n[REAL] = a.n[REAL] + b.n[REAL];
	sum->n[IMAG] = a.n[IMAG] + b.n[IMAG];
}

	
Subtract(a, b, diff)
complex a, b, *diff;
{
	diff->n[REAL] = a.n[REAL] - b.n[REAL];
	diff->n[IMAG] = a.n[IMAG] - b.n[IMAG];
}

	
Multiply(a, b, prod)
complex a, b, *prod;
{
	prod->n[REAL] = (a.n[REAL] * b.n[REAL]) - (a.n[IMAG] * b.n[IMAG]);
	prod->n[IMAG] = (a.n[REAL] * b.n[IMAG]) + (a.n[IMAG] * b.n[REAL]);
}

	
Divide(a, b, quot)
complex a, b, *quot;
{
	double denom;

	denom = SQR(b.n[REAL]) + SQR(b.n[IMAG]);
	if (denom == 0.0)
		printf("Divide by zero error!\n"), exit(20);

	quot->n[REAL] = ((a.n[REAL] * b.n[REAL]) + (a.n[IMAG] * b.n[IMAG]))
			/ denom;
	quot->n[IMAG] = ((a.n[IMAG] * b.n[REAL]) - (a.n[REAL] * b.n[IMAG]))
			/ denom;
}

	
SynthDivide(poly, comp, stop)
/* Computes the synthetic division of the polynomial poly[] by (z-comp), 
 * where "z" is assumed the complex variable in our polynomial. */
complex poly[], comp;
int stop;
{
	int i;
	complex prod;

	for (i=1; i<=stop; i++) {
		Multiply(poly[i-1], comp, &prod);
		Add(poly[i], prod, &poly[i]);
	}
}

	
Iterate(poly, zOld, n, zNew)
/* One iteration of Newton's method.  "zOld" is the old guess of the value
 * of a root of poly[].  "zNew" is the newly calculated guess. */
complex poly[], zOld;
int n;
complex *zNew;
{
	complex quotient;
	Divide(poly[n], poly[n-1], &quotient);
	Subtract(zOld, quotient, zNew);
}

	
Guess(z)
/* A random complex number, 50% chance of being negative. */
complex *z;
{
#ifdef UNIX
#define	ran	drand48
#endif
	double ran();

	z->n[REAL] = ran();
	z->n[IMAG] = ran();
	z->n[REAL] = (ran() < 0.50) ? z->n[REAL] : -(z->n[REAL]);
	z->n[IMAG] = (ran() < 0.50) ? z->n[IMAG] : -(z->n[IMAG]);
}

	
BOOL ErrorTooBig(y, z, eps)
/* Compute the Euclidean distance between y and z in the complex plane.
 * This is just our good old friend, the "distance formula".  (Add the
 * squares of y and z, and take the square root.)  Is the distance less
 * than epsilon?  If so, the error is allowable; else, it's too big. */
complex y, z;
double eps;
{
	complex diff;
	double sqrt();

	Subtract(y, z, &diff);
	return(	sqrt(SQR(diff.n[REAL]) + SQR(diff.n[IMAG])) < eps
		? FALSE
		: TRUE);
}