Homework Assignment 5 Supplement

Here is the source code to complex.h and complex.c.

complex.h

/*
 *  =================================================================== 
 *  Copyright (c) 1992-1997,                                                 
 *  Mark Austin.
 *  University of Maryland, MD 20742                                   
 *                                                                     
 *  ------------------------------------------------------------------- 
 *  Header file for complex number data structure...
 *  ------------------------------------------------------------------- 
 *                                                                     
 *  Written by: M.A. Austin                      July, 1992 - Jan, 1993
 *  ===================================================================
 */

#ifndef COMPLEX_H
#define COMPLEX_H

/* Com Number Data Structure */

typedef struct {
        double      dReal;
        double dImaginary;
} COMPLEX;

#ifdef  __STDC__

COMPLEX *comAlloc( void );
void     comFree( COMPLEX * );
COMPLEX *comNegate( COMPLEX * );
double   comAbs( COMPLEX * );
COMPLEX *comAdd( COMPLEX * , COMPLEX * );
COMPLEX *comSub( COMPLEX * , COMPLEX * );
COMPLEX *comMult( COMPLEX * , COMPLEX * );
COMPLEX *comDiv( COMPLEX * , COMPLEX * );
COMPLEX *comScale( double, COMPLEX * );
COMPLEX *comConjugate( COMPLEX * );
COMPLEX *comSqrt( COMPLEX * );

#else  /* case not STDC */

COMPLEX *comAlloc();
void     comFree();
COMPLEX *comNegate();
double   comAbs();
COMPLEX *comAdd();
COMPLEX *comSub();
COMPLEX *comMult();
COMPLEX *comDiv();
COMPLEX *comScale();
COMPLEX *comConjugate();
COMPLEX *comSqrt();

#endif /* end not STDC */
#endif /* not COMPLEX_H */

complex.c

/*
 *  =================================================================== 
 *  complex.c : Library of functions for Complex Number Arithmetic.
 *                                                                     
 *  Written by: M.A. Austin                                 March, 1993
 *  ===================================================================
 */

#include 
#include 

#include "miscellaneous.h"
#include "complex.h"

/* #define DEBUG */

/*
 *  =================================
 *  Allocate and Free Complex Numbers
 *  =================================
 */

COMPLEX *comAlloc() {
   return((COMPLEX *) safeMalloc( sizeof(COMPLEX), __FILE__, __LINE__ ));
}

#ifdef __STDC__
void comFree(COMPLEX *spComplex)
#else  /* case not STDC */
void comFree(spComplex)
COMPLEX *spComplex;
#endif /* end case not STDC */
{
   free((char *) spComplex);
}


/*
 *  =============================
 *  Operations on Complex Numbers 
 *  =============================
 */

#ifdef __STDC__
COMPLEX *comNegate(COMPLEX *spA )
#else  /* case not STDC */
COMPLEX *comNegate( *spA )
COMPLEX *spA;
#endif /* end case not STDC */
{
COMPLEX *spComplex = comAlloc();

   spComplex->dReal      = -spA->dReal;
   spComplex->dImaginary = -spA->dImaginary;

   return (spComplex);
}

#ifdef __STDC__
COMPLEX *comAdd( COMPLEX *spComplexA, COMPLEX *spComplexB )
#else  /* case not STDC */
COMPLEX *comAdd( spComplexA, spComplexB)
COMPLEX *spComplexA, *spComplexB;
#endif /* end case not STDC */
{
COMPLEX *spComplex = comAlloc();

   spComplex->dReal      = spComplexA->dReal      + spComplexB->dReal;
   spComplex->dImaginary = spComplexA->dImaginary + spComplexB->dImaginary;

   return (spComplex);
}

#ifdef __STDC__
COMPLEX *comSub( COMPLEX *spComplexA, COMPLEX *spComplexB )
#else  /* case not STDC */
COMPLEX *comSub( spComplexA, spComplexB)
COMPLEX *spComplexA, *spComplexB;
#endif /* end case not STDC */
{
COMPLEX *spComplex = comAlloc();

   spComplex->dReal      = spComplexA->dReal      - spComplexB->dReal;
   spComplex->dImaginary = spComplexA->dImaginary - spComplexB->dImaginary;

   return (spComplex);
}

#ifdef __STDC__
COMPLEX *comMult( COMPLEX *spComplexA, COMPLEX *spComplexB )
#else  /* start case not STDC */
COMPLEX *comMult( spComplexA, spComplexB )
COMPLEX *spComplexA, *spComplexB;
#endif /* end case not STDC */
{
COMPLEX *spComplex = comAlloc();

   spComplex->dReal      = spComplexA->dReal*spComplexB->dReal -
                           spComplexA->dImaginary*spComplexB->dImaginary;
   spComplex->dImaginary = spComplexA->dImaginary*spComplexB->dReal +
                           spComplexA->dReal*spComplexB->dImaginary;

   return (spComplex);
}

#ifdef __STDC__
COMPLEX *comDiv( COMPLEX *spComplexA, COMPLEX *spComplexB)
#else  /* case not STDC */
COMPLEX *comDiv( spComplexA, spComplexB)
COMPLEX *spComplexA, *spComplexB;
#endif /* end case not STDC */
{
COMPLEX *spComplex = comAlloc();
float  fR, fDen;

   if(fabs(spComplexB->dReal) >= fabs(spComplexB->dImaginary)) {
      fR   = spComplexB->dImaginary/spComplexB->dReal;
      fDen = spComplexB->dReal + fR*spComplexB->dImaginary;
      spComplex->dReal      = (spComplexA->dReal + fR*spComplexA->dImaginary)/fDen;
      spComplex->dImaginary = (spComplexA->dImaginary - fR*spComplexA->dReal)/fDen;
   } else {
      fR   = spComplexB->dReal/spComplexB->dImaginary;
      fDen = spComplexB->dImaginary + fR*spComplexB->dReal;
      spComplex->dReal      = (fR*spComplexA->dReal + spComplexA->dImaginary)/fDen;
      spComplex->dImaginary = (fR*spComplexA->dImaginary - spComplexA->dReal)/fDen;
   }

   return (spComplex);
}

#ifdef __STDC__
COMPLEX *comScale( double dFactor, COMPLEX *spComplexA )
#else  /* case not STDC */
COMPLEX *comScale( dFactor, spComplexA)
double      dFactor;
COMPLEX *spComplexA;
#endif /* end case not STDC */
{
COMPLEX *spComplex = comAlloc();

   spComplex->dReal      = dFactor*spComplexA->dReal;
   spComplex->dImaginary = dFactor*spComplexA->dImaginary;

   return (spComplex);
}

#ifdef __STDC__
COMPLEX *comConjugate( COMPLEX *spComplexA)
#else  /* case not STDC */
COMPLEX *comConjugate( spComplexA )
COMPLEX *spComplexA;
#endif /* end case not STDC */
{
COMPLEX *spComplex = comAlloc();

   spComplex->dReal      =  spComplexA->dReal;
   spComplex->dImaginary = -spComplexA->dImaginary;

   return (spComplex);
}

#ifdef __STDC__
double
comAbs( COMPLEX *spComplexA )
#else  /* case not STDC */
double comAbs( spComplexA )
COMPLEX *spComplexA;
#endif /* end case not STDC */
{
double dX, dY, dTemp;
double dAnswer;

   dX = fabs(spComplexA->dReal);
   dY = fabs(spComplexA->dImaginary);
   if(dX == 0)
      dAnswer = dY;
   else if(dY == 0)
      dAnswer = dX;
   else if(dX > dY) {
     dTemp = dY/dX;
     dAnswer = dX*sqrt((double) (1.0 + dTemp*dTemp));
   } else {
     dTemp = dX/dY;
     dAnswer = dY*sqrt((double) (1.0 + dTemp*dTemp));
   } 

   return (dAnswer);
}

#ifdef __STDC__
COMPLEX *comSqrt( COMPLEX *spComplexA )
#else  /* case not STDC */
COMPLEX *comSqrt( spComplexA )
COMPLEX *spComplexA;
#endif /* end case not STDC */
{
COMPLEX *spComplex = comAlloc();
double dX, dY, dW, dR;

    if((spComplexA->dReal == 0) && (spComplexA->dImaginary == 0.0)) {
        spComplex->dReal      = 0.0;
        spComplex->dImaginary = 0.0;
    } else {
        dX = fabs(spComplexA->dReal);
        dY = fabs(spComplexA->dImaginary);

        if(dX >= dY) {
           dR = dY/dX;
           dW = sqrt(dX)*sqrt(0.5*(1.0 + sqrt(1+dR*dR)));
        } else {
           dR = dX/dY;
           dW = sqrt(dY)*sqrt(0.5*(dR + sqrt(1+dR*dR)));
        }

        if(spComplexA->dReal >= 0.0) {
           spComplex->dReal = dW;
           spComplex->dImaginary = spComplexA->dImaginary/(2.0*dW);
        } else {
           spComplex->dImaginary = (spComplexA->dImaginary > 0.0) ? dW : -dW;
           spComplex->dReal = spComplexA->dImaginary/(2.0*spComplex->dImaginary);
        }
    }

    return (spComplex);
}