page@swan.ulowell.edu (Bob Page) (12/31/88)
Submitted-by: kim@uts.amdahl.com (Kim E. DeVaughn)
Posting-number: Volume 2, Issue 111
Archive-name: applications/complex.1
Here are some complex number functions. There are versions for both
x,y and polar representations. You know, wonderful things like +, -,
*, /, sin, cos, tan, exp, ln and power. Sorry, no hyperbolic complex
functions...
Enjoy!
# This is a shell archive.
# Remove everything above and including the cut line.
# Then run the rest of the file through sh.
#----cut here-----cut here-----cut here-----cut here----#
#!/bin/sh
# shar: Shell Archiver
# Run the following text with /bin/sh to create:
# complex.c
# This archive created: Fri Dec 30 13:04:13 1988
cat << \SHAR_EOF > complex.c
/*
The implementations of these functions are hereby released to the public
domain. I think it would be silly to try to sell this because:
1. limited market.
2. anyone with half a brain can do the same after looking in a frosh
level calculus textbook.
If you have find any mistakes, have a better implementation or
have any questions contact me at
GEnie - R.DEVENEZIA
or
225B Hillcrest Manor Ct
Utica, NY 13501
-------------- complex math routines ---------------
Basic notation:
z = x + iy one mapping (x, y)
= r7[ cos(theta) + i sin(theta) ] infinite mappings (r, 2pi7k7theta)
= r7cis(theta);
where cis(theta) = cos(theta) + i sin(theta)
i7x
= r7e infinite mappings (r, 2pi7k7theta)
Result Methodology:
Since most compilers can't return a structure as a function result
this implementation maintains a 'result pool'. All functions
return pointers to an entry in the result pool (set to 100 currently).
Why do it this way? Well, I want to nest functions within functions,
and it ain't handy (hell, it ain't possible) when the result is passed
back through the argument list (you need all sorts of intermediate
variables, uggh).
For example, I can now write an expression like
zprime = cMult (cMult (z,z), cExp (z));
voila! z' = z27exp(z)
Dangers: Any variable that is to be a result of these functions
must be a pointer to allocated memory.
Examples:
If predefined variables of type complex are going to be used
the addresses of these variables must be passed to the functions.
i.e.
complex z1, z2; / * complex variable * /
complex *result; / * pointer to complex variable * /
real modulus;
z1.R = 0;
z1.i = 2; / * 2i * /
z2.R = 1;
z2.i = 2; / * 1 + 2i * /
modulus = cMod (&z1); / * argument is a pointer * /
result = cSin (&z2); / * argument is a pointer, as is result * /
/ * result is a pointer to an entry in the result pool * /
modulus = cMod (cMult (&z1, &z2)); / * nested call * /
If predefined variables of type *complex are going to be used
you must first allocate memory and pass the pointer returned
[i.e. the variable] to the functions.
i.e.
complex *z1, *z2; / * pointers to complex variable * /
complex *result; / * pointer to complex variable * /
real modulus;
z1 = malloc (sizeof(complex));
z1->R = 0;
z1->i = 2; / * 2i * /
z2->R = 1;
z2->i = 2; / * 1 + 2i * /
modulus = cMod (z1); / * argument is a pointer * /
result = cSin (z2); / * argument is a pointer, as is result * /
/ * result is a pointer to an entry in the result pool * /
modulus = cMod (cMult (z1, z2));
and be sure to free up any memory you allocate using Amiga memory
allocation functions.
Personally I prefer the first method.
available functions:
c prefix requires complex arguments
p prefix requires polar arguments
cNextResult - maintains result pool, don't play with it
pNextResult - maintains result pool, don't play with it
cRe - Real part
pRe - Real part
cIm - Imaginary part
pIm - Imaginary part
cMod - Modulus (radius)
pMod - Modulus (radius)
cArg - Argument (angle)
pArg - Argument (angle)
cTOp - convert complex to polar
pTOc - convert polar to complex
cConj - complex conjugate
pConj - complex conjugate
MakeComplex - return a complex* from two reals
MakePolar - return a polar* from two reals
cAdd - addition
cSub - subtraction
cMult - multiplication
cDiv - division, check for /0
cExp - exponentiation
cLn - natural logarithm, check for /0
cSin - sine
cCos - cosine
cTan - tangent, check for /0
cIPower - integer power
cPower - complex to complex power
pAdd - addition
pSub - subtraction
pMult - multiplication
pDiv - division, check for /0
pExp - exponentiation
pLn - natural logarithm, check for /0
pSin - sine
pCos - cosine
pTan - tangent, check for /0
pIPower - integer power
pPower - complex to complex power
*/
#ifndef sqrt
#include <math.h>
#endif
#define TRUE 1
#define FALSE 0
#define pi 3.1415926
#define cPoolSize 100 /* don't put those freakin' commas */
#define pPoolSize 100 /* after a define!! */
typedef struct {
float R;
float i;
} complex;
typedef struct {
float r;
float theta;
} polar;
static complex cResultPool [cPoolSize];
static polar pResultPool [pPoolSize];
static complex *ComplexResult;
static polar *PolarResult;
static int ComplexIndex = 0; /* index into result pool */
static int PolarIndex = 0; /* index into result pool */
int divisionByZero = FALSE;
void
cNextResult () /* arrrgh! remember, no commas here! */
{
ComplexResult = &cResultPool [ComplexIndex++];
if (ComplexIndex == cPoolSize)
ComplexIndex = 0;
}
/* cNextResult */
void
pNextResult ()
{
PolarResult = &pResultPool [PolarIndex++];
if (PolarIndex == pPoolSize)
PolarIndex = 0;
}
/* pNextResult */
float
cRe (z) /* real part */
complex *z;
{
return (z->R);
}
/* cRe */
float
pRe (z) /* real part */
polar *z;
{
return (z->r * cos (z->theta));
}
/* pRe */
float
cIm (z) /* imaginary part */
complex *z;
{
return (z->i);
}
/* cIm */
float
pIm (z) /* imaginary part */
polar *z;
{
return (z->r * sin (z->theta));
}
/* pIm */
float
cMod (z) /* modulus (radius) */
complex *z;
{
return (sqrt ( z->R * z->R + z->i * z->i ));
}
/* cMod */
float
pMod (z) /* modulus (radius) */
polar *z;
{
return (z->r);
}
/* pMod */
float
cArg (z) /* angle */
complex *z;
/*
arg z is a one to infinite mapping. (infinite number of solutions)
arg z is a one to one mapping (one soln) using the principal determination.
ie.
a + bi = (R,theta)
where theta = principal_theta 1 27pi7k; where k is an integer
and theta in (-pi,pi] is the principal determination
*/
{
float temp;
if (z->R == 0.0) /* no radius */
if (z->i > 0.0) return ( pi / 2.0);
else
if (z->i < 0.0) return (-pi / 2.0);
else /* origin */
return (0.0);
else
/* Manx's atan returns the arc tangent in the range -pi/2 to pi/2 */
return (atan (z->i / z->R));
}
/* cArg */
float
pArg (z) /* angle */
polar *z;
{
return (z->theta);
}
/* pArg */
polar *
cTOp (z) /* convert complex to polar */
complex *z;
{
pNextResult;
PolarResult->r = cMod (z);
PolarResult->theta = cArg (z);
return (PolarResult);
}
/* cTOp */
complex *
pTOc (z) /* convert polar to complex */
polar *z;
{
cNextResult;
ComplexResult->R = pRe (z);
ComplexResult->i = pIm (z);
return (ComplexResult);
}
/* pTOc */
complex *
cConj (z) /* complex conjugate */
complex *z;
{
cNextResult;
ComplexResult->R = z->R;
ComplexResult->i = -z->i;
return (ComplexResult);
}
/* cConj */
polar *
pConj (z) /* complex conjugate (polar form) */
polar *z;
{
pNextResult;
PolarResult->r = z->r;
PolarResult->theta = - (z->theta);
return (PolarResult);
}
/* pConj */
complex *
MakeComplex (R, i)
float R, i;
{
cNextResult;
ComplexResult->R = R;
ComplexResult->i = i;
return (ComplexResult);
}
/* MakeComplex */
polar *
MakePolar (r, theta)
float r, theta;
{
pNextResult;
PolarResult->r = r;
PolarResult->theta = theta;
return (PolarResult);
}
/* MakePolar */
/*-------------------------------------------------------------------------*/
/* Complex */
complex *
cAdd (z1, z2) /* addition, z1+z2 */
complex *z1, *z2;
{
cNextResult;
ComplexResult->R = z1->R + z2->R;
ComplexResult->i = z1->i + z2->i;
return (ComplexResult);
}
/* cAdd */
complex *
cSub (z1, z2) /* subtraction, z1-z2 */
complex *z1, *z2;
{
cNextResult;
ComplexResult->R = z1->R - z2->R;
ComplexResult->i = z1->i - z2->i;
return (ComplexResult);
}
/* cSub */
complex * /* multiplication, z1*z2 */
cMult (z1, z2)
complex *z1, *z2;
{
cNextResult;
ComplexResult->R = z1->R * z2->R - z1->i * z2->i;
ComplexResult->i = z1->R * z2->i + z1->i * z2->R;
return (ComplexResult);
}
/* cMult */
complex *
cDiv (z1, z2) /* division, z1/z2 */
complex *z1, *z2;
{
float z2_modulusSquared;
cNextResult;
z2_modulusSquared = z2->R * z2->R + z2->i * z2->i;
if (z2_modulusSquared == 0.0) {
divisionByZero = TRUE;
ComplexResult->R = 0.0;
ComplexResult->i = 0.0;
}
else {
divisionByZero = FALSE;
ComplexResult->R = (z1->i * z2->i + z1->R * z2->R)
/ z2_modulusSquared;
ComplexResult->i = (z1->i * z2->R - z1->R * z2->i)
/ z2_modulusSquared;
}
return (ComplexResult);
}
/* cDiv */
complex *
cExp (z) /* exponentiation */
complex *z;
/*
infinity
z -- n
e == 1 + \ z
/ ----
-- n!
n=1
OR the method used (by the following equality):
z x + iy x iy x x
e == e == e 7 e == e * cos(y) + i * e * sin(y)
*/
{
float eToTheX;
cNextResult;
eToTheX = exp (z->R);
ComplexResult->R = eToTheX * cos (z->i);
ComplexResult->i = eToTheX * sin (z->i);
return (ComplexResult);
}
/* cExp */
complex *
cLn (z) /* natural logarithm */
complex *z;
/*
ln z is a one to infinite mapping. (infinite number of solutions)
ln z is a one to one mapping (one soln) using the principal determination.
complex real
ln z = ln(|z|) + i * arg(z)
*/
{
cNextResult;
if (cMod (z) == 0.0) {
divisionByZero = TRUE;
ComplexResult->R = 0.0;
ComplexResult->i = 0.0;
}
else {
divisionByZero = FALSE;
ComplexResult->R = log ( cMod (z) );
ComplexResult->i = cArg (z); /* principal determination */
}
return (ComplexResult);
}
/* cLn */
complex *
cSin (z) /* sine */
complex *z;
/*
infinity
-- k
sin z == \ (-1) 2k+1
/ ---------- z
-- (2k + 1)!
k=0
OR the method used (by the following equality):
1 iz -iz
sin z == --- [ e - e ]
27i
also,
== sin(x)7cosh(y) + i cos(x)7sinh(y)
where cosh is the hyperbolic cosine and
sinh is the hyperbolic sine.
*/
{
return (cMult (MakeComplex (0.0, -0.5), /* 1/27i */
cSub (cExp (cMult (MakeComplex (0.0, 1.0), z)),
cExp (cMult (MakeComplex (0.0,-1.0), z))
)
)
);
}
/* cSin */
complex *
cCos (z) /* cosine */
complex *z;
/*
infinity
-- k
cos z == \ (-1) 2k
/ ------ z
-- (2k)!
k=0
OR the method used (by the following equality):
1 iz -iz
cos z == --- [ e + e ]
2
also,
== cos(x)7cosh(y) - i sin(x)7sinh(y)
where cosh is the hyperbolic cosine and
sinh is the hyperbolic sine.
*/
{
return (cMult (MakeComplex (0.5, 0.0),
cAdd (cExp (cMult (MakeComplex (0.0, 1.0), z)),
cExp (cMult (MakeComplex (0.0,-1.0), z))
)
)
);
}
/* cCos */
complex *
cTan (z) /* tangent */
complex *z;
/*
sin z
tan z = -------
cos z
*/
{
if (cMod (z) == 0.0) {
divisionByZero = TRUE;
cNextResult;
ComplexResult->R = 0.0;
ComplexResult->i = 0.0;
return (ComplexResult);
}
else {
divisionByZero = FALSE;
return (cDiv (cSin (z), cCos (z)));
}
}
/* cTan */
complex *
cIPower (z, n) /* complex number to an integer power */
complex *z;
int n;
/*
n n
z = |z| 7 cis (n7theta)
*/
{
float modulusToTheN, ntheta;
cNextResult;
modulusToTheN = pow (cMod (z), n); /* Manx's power function */
ntheta = n * cArg (z);
ComplexResult->R = modulusToTheN * cos (ntheta);
ComplexResult->i = modulusToTheN * sin (ntheta);
return (ComplexResult);
}
/* cIPower */
complex *
cPower (z1, z2) /* complex number raised to a complex number */
complex *z1, *z2;
/*
z2 z2*log(z1)
z1 = exp
z2
Note: There are an infinite number of points, (r,theta), which equal z1 .
There is only one point if we take the principal determinations of
log and exp.
*/
{
return (cExp (cMult (z2, cLn (z1))));
}
/* cPower */
/*-------------------------------------------------------------------------*/
/* Polar */
polar *
pAdd (z1, z2) /* addition */
polar *z1, *z2;
{
return (cTOp (MakeComplex (pRe (z1) + pRe (z2),
pIm (z1) + pIm (z2))
)
);
}
/* pAdd */
polar *
pSub (z1, z2) /* subtraction */
polar *z1, *z2;
{
return (cTOp (MakeComplex (pRe (z1) - pRe (z2),
pIm (z1) - pIm (z2))
)
);
}
/* pSub */
polar *
pMult (z1, z2) /* multiplication */
polar *z1, *z2;
{
float theta;
pNextResult;
PolarResult->r = z1->r * z2->r;
theta = z1->theta + z2->theta;
/*
principal determination ?
if ((theta <= -pi) OR (theta > pi))
theta = theta + 2*pi* INT(theta/2/pi)
}
*/
PolarResult->theta = theta;
return (PolarResult);
}
/* pMult */
polar *
pDiv (z1, z2) /* division */
polar *z1, *z2;
{
pNextResult;
if (z2->r != 0.0) {
divisionByZero = FALSE;
PolarResult->r = z1->r / z2->r;
}
else {
divisionByZero = TRUE;
PolarResult->r = 0.0;
}
PolarResult->theta = z1->theta - z2->theta;
return (PolarResult);
}
/* pDiv */
polar *
pExp (z) /* exponentiation */
polar *z;
{
return (cTOp (cExp (pTOc (z))));
}
/* pExp */
polar *
pLn (z) /* natural logarithm */
polar *z;
{
return (cTOp (cLn (pTOc (z))));
}
/* pLn */
polar *
pSin (z) /* sine */
polar *z;
{
return (cTOp (cSin (pTOc (z))));
}
/* pSin */
polar *
pCos (z) /* cosine */
polar *z;
{
return (cTOp (cCos (pTOc (z))));
}
/* pCos */
polar *
pTan (z) /* tangent */
polar *z;
{
return (cTOp (cTan (pTOc (z))));
}
/* pTan */
polar *
pIPower (z, n) /* polar number to an integer power */
polar *z;
int n;
{
pNextResult;
PolarResult->r = n * z->r;
PolarResult->theta = n * z->theta;
return (PolarResult);
}
/* pIPower */
polar *
pPower (z1, z2) /* polar number raised to a polar number */
polar *z1, *z2;
{
return (cTOp (cPower (pTOc (z1), pTOc (z2))));
}
/* pPower */
SHAR_EOF
# End of shell archive
exit 0
--
Bob Page, U of Lowell CS Dept. page@swan.ulowell.edu ulowell!page
Have five nice days.