complex.c

该程序是用vc开发的对动态数组进行管理的DLL

/* Copyright (c) Colorado School of Mines, 2003.*/
/* All rights reserved.                       */

/*********************** self documentation **********************/
/*****************************************************************************
COMPLEX - Functions to manipulate fcomplex numbers

cadd	add two fcomplex numbers
csub	subtract two fcomplex numbers
cmul	multiply two fcomplex numbers
cdiv	divide two fcomplex numbers
cmplx	make a fcomplex number from two real numbers
conjg	fcomplex conjugate of a fcomplex number 
cneg	negate a fcomplex number
cinv	invert a fcomplex number
csqrt	fcomplex square root of a fcomplex number
cexp	fcomplex exponential of a fcomplex number
crmul	multiply a fcomplex number by a real number 
rcabs	real magnitude of a fcomplex number

******************************************************************************
Structure:
typedef struct _complexStruct {  fcomplex number
	float r,i;
} fcomplex;

******************************************************************************
Function Prototypes:
fcomplex cadd (fcomplex a, fcomplex b);
fcomplex csub (fcomplex a, fcomplex b);
fcomplex cmul (fcomplex a, fcomplex b);
fcomplex cdiv (fcomplex a, fcomplex b);
float rcabs (fcomplex z);
fcomplex cmplx (float re, float im);
fcomplex conjg (fcomplex z);
fcomplex cneg (fcomplex z);
fcomplex cinv (fcomplex z);
fcomplex csqrt (fcomplex z);
fcomplex cexp (fcomplex z);
fcomplex crmul (fcomplex a, float x);

******************************************************************************
Notes:
The function "rcabs" was originally called "fcabs". This produced
a collision on some systems so a new name was chosen.

******************************************************************************
Reference:
Adapted from Press et al, 1988, Numerical Recipes in C (Appendix E).

******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 06/02/89
Modified:  Dave Hale, Colorado School of Mines, 04/26/90
	Added function cinv().
*****************************************************************************/
/**************** end self doc ********************************/

#include "cwp.h"

fcomplex cadd(fcomplex a, fcomplex b)
{
	fcomplex c;
	c.r = a.r+b.r;
	c.i = a.i+b.i;
	return c;
}

fcomplex csub(fcomplex a, fcomplex b)
{
	fcomplex c;
	c.r = a.r-b.r;
	c.i = a.i-b.i;
	return c;
}

fcomplex cmul(fcomplex a, fcomplex b)
{
	fcomplex c;
	c.r = a.r*b.r-a.i*b.i;
	c.i = a.i*b.r+a.r*b.i;
	return c;
}

fcomplex cdiv(fcomplex a, fcomplex b)
{
	fcomplex c;
	float r,den;
	if (fabs(b.r)>=fabs(b.i)) {
		r = b.i/b.r;
		den = b.r+r*b.i;
		c.r = (a.r+r*a.i)/den;
		c.i = (a.i-r*a.r)/den;
	} else {
		r = b.r/b.i;
		den = b.i+r*b.r;
		c.r = (a.r*r+a.i)/den;
		c.i = (a.i*r-a.r)/den;
	}
	return c;
}

fcomplex cmplx(float re, float im)
{
	fcomplex c;
	c.r = re;
	c.i = im;
	return c;
}

fcomplex conjg(fcomplex z)
{
	fcomplex c;
	c.r = z.r;
	c.i = -z.i;
	return c;
}

fcomplex cneg(fcomplex z)
{
	fcomplex c;
	c.r = -z.r;
	c.i = -z.i;
	return c;
}

fcomplex cinv(fcomplex z)
{
	fcomplex c;
	float s;
	s = 1.0f/(z.r*z.r+z.i*z.i);
	c.r = z.r*s;
	c.i = -z.i*s;
	return c;
}

fcomplex csqrt(fcomplex z)
{
	fcomplex c;
	float x,y,w,r;
	if (z.r==0.0 && z.i==0.0) {
		c.r = c.i = 0.0;
		return c;
	} else {
		x = (float)fabs(z.r);
		y = (float)fabs(z.i);
		if (x>=y) {
			r = y/x;
			w = (float)(sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r))));
		} else {
			r = x/y;
			w = (float)(sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r))));
		}
		if (z.r>=0.0) {
			c.r = w;
			c.i = z.i/(2.0f*w);
		} else {
			c.i = (z.i>=0.0) ? w : -w;
			c.r = z.i/(2.0f*c.i);
		}
		return c;
	}
}

fcomplex cexp(fcomplex z)
{
	float a;
	fcomplex c;
	a = (float)exp(z.r);
	c.r = (float)(a*cos(z.i));
	c.i = (float)(a*sin(z.i));
	return c;
}

fcomplex crmul(fcomplex a, float x)
{
	fcomplex c;
	c.r = x*a.r;
	c.i = x*a.i;
	return c;
}

float rcabs(fcomplex z)
{
	float x,y,ans,temp;
	x = (float)fabs(z.r);
	y = (float)fabs(z.i);
	if (x==0.0)
		ans = y;
	else if (y==0.0)
		ans = x;
	else if (x>y) {
		temp = y/x;
		ans = (float)(x*sqrt(1.0+temp*temp));
	} else {
		temp =x/y;
		ans = (float)(y*sqrt(1.0+temp*temp));
	}
	return ans;
}