/*pcomplex.c - data structure for, and operations on, complex numbers*/

/*THIS VERSION IS WRONG.  WHY?*/

#include "pcomplex.h"

complex *ComplexMultiply(a, b)
complex *a, *b;
  {
  complex product_data, *product;
  product = &product_data;
  Re(product) = Re(a)*Re(b) - Im(a)*Im(b);	/*notice that because fields*/
  Im(product) = Re(a)*Im(b) + Im(a)*Re(b);	/*are accessed through macros,*/
  return product;				/*we don't need to change much*/
  }						/*from the previous version*/

/* ... */

/*THE FOLLOWING VERSION IS BETTER, BUT CAN STILL BREAK UNDER A SPECIFIC
  CALLING CIRCUMSTANCE.  WHAT IS THAT CIRCUMSTANCE?  HOW DO YOU FIX IT, OR IS
  IT WORTH BOTHERING WITH?*/

void ComplexMultiply(a, b, product)
complex *a, *b, *product;
  {
  Re(product) = Re(a)*Re(b) - Im(a)*Im(b);	/*notice that because fields*/
  Im(product) = Re(a)*Im(b) + Im(a)*Re(b);	/*are accessed through macros,*/
  }						/*we don't need to change much*/
  						/*from the previous version*/

void ComplexDivide(a, b, quotient)
complex *a, *b, *quotient;
  {
  complex bconjugate_data, *bconjugate;
  double magnitude_sq;
  bconjugate = &bconjugate_data;
  magnitude_sq = Re(b)*Re(b) + Im(b)*Im(b);
  Re(bconjugate) = Re(b)/magnitude_sq;
  Im(bconjugate) = -Im(b)/magnitude_sq;	/*construct normalised conjugate*/
  quotient = ComplexMultiply(a, bconjugate);
  }

void ComplexAdd(a, b, sum)
complex *a, *b, *sum;
  {
  Re(sum) = Re(a)+Re(b);
  Im(sum) = Im(a)+Im(b);
  }

complex ComplexSubtract(a, b, difference)
complex *a, *b, *difference;
  {
  Re(difference) = Re(a)-Re(b);
  Im(difference) = Im(a)-Im(b);
  }