typedef 
	     	struct 
		{
	double real;
	double imag;
		} 
complex;

/* handy but not necessary: */
complex zero = {0.0,0.0};
complex one  = {1.0,0.0};

/* Some complex functions. All have  c_  in
front, just as a reminder that they are not
for real numbers but for complex. They of course
work for real numbers, if you just use 0.0 for the
imaginary part.  */

complex c_add(alpha,beta)
complex alpha, beta;
{
   complex z;
   z.real = alpha.real + beta.real;
   z.imag = alpha.imag + beta.imag;
   return (z);
}

complex c_subtract(alpha,beta)
complex alpha, beta;
{
   complex z;
   z.real = alpha.real - beta.real;
   z.imag = alpha.imag - beta.imag;
   return (z);
}

complex c_mult(alpha,beta)
complex alpha, beta;
{
   complex z;
   z.real = alpha.real * beta.real - alpha.imag * beta.imag;
   z.imag = alpha.real * beta.imag + alpha.imag * beta.real;
   return (z);
}

complex c_square(alpha)
complex alpha;
{
   complex z;
   z = c_mult(alpha,alpha);
   return (z);
}

/* This division subroutine is set up to give z/0.0 = 0.0. Thus
it should never crashe, but of course it could lead to the wrong
answer if you accidentally divide by zero. That's why the warning
is there on stderr. */

complex c_divide(alpha,beta)
complex alpha, beta;
{
   complex z;
   double R, modulus();
   R = modulus(beta);
   if (R==0.0)
   {
	z.real = 0.0;
	z.imag = 0.0;
	fprintf(stderr,"\n\nWARNING: DIVISION BY ZERO in c_divide\n\n");
   }
   else
   {
	z.real = (alpha.real*beta.real + alpha.imag*beta.imag)/(R*R);
	z.imag = (alpha.imag*beta.real - alpha.real*beta.imag)/(R*R);
   }
   return (z);
}

/* The next routine is included as an illustration of what might
be done. The exact way to do it will depend on the context. */

complex_print(alpha)
complex alpha;
{
   printf("%5.3f + %5.3fi",alpha.real,alpha.imag);
}

/* Another word for absolute value is modulus: */

double modulus(alpha)
complex alpha;
{
   double z;
   z = sqrt(alpha.real*alpha.real + alpha.imag*alpha.imag);
   return(z);
}

/* Warning: a non-zero complex number has two square roots; this
is only one of them. There is an easy way to get the other one from
this one. */

complex square_root(alpha)
complex alpha;
{
   complex z;
   if (modulus(alpha)==0.0) return (alpha);
   else
   {
	double angle, R;
	angle = atan2(alpha.imag,alpha.real);
	angle = angle/2.0;
	R = modulus(alpha);
	R = exp(log(R)/2);
	z.real = R * cos(angle); z.imag = R * sin(angle);
	return(z);
   }
}

/* Warning: a non-zero complex number has three cube roots; this
is only one of them. There is an easy way to get the other two from
this one. Note also, that it should be easy to modify this for the
n-th root, if that were needed in some context.  */

complex cube_root(alpha)
complex alpha;
{
   complex z;
   if (modulus(alpha)==0.0) return (alpha);
   else
   {
	double angle, R;
	angle = atan2(alpha.imag,alpha.real);
	angle = angle/3.0;
	R = modulus(alpha);
	R = exp(log(R)/3);
	z.real = R * cos(angle); z.imag = R * sin(angle);
	return(z);
   }
}
complex neg = {-1.0,0.0};