rs@uunet.UU.NET (Rich Salz) (07/16/87)
Submitted-by: gwyn@brl.arpa (Doug Gwyn)
Posting-Number: Volume 10, Issue 56
Archive-name: complex-lib
[ Here is a package to do complex arithmetic; blame me for the Makefile.
--r$ ]
#! /bin/sh
# This is a shell archive. Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file". To overwrite existing
# files, type "sh file -c". You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g.. If this archive is complete, you
# will see the following message at the end:
# "End of shell archive."
# Contents: Makefile complex.3 complex.h cx_test.c cxadd.c cxampl.c
# cxconj.c cxcons.c cxcopy.c cxdiv.c cxmul.c cxphas.c cxphsr.c
# cxscal.c cxsqrt.c cxsub.c
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f Makefile -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"Makefile\"
else
echo shar: Extracting \"Makefile\" \(395 characters\)
sed "s/^X//" >Makefile <<'END_OF_Makefile'
XALL=complex.3 complex.h libcomplex.a
XOBJS=\
X cxadd.o cxampl.o cxconj.o cxcons.o cxcopy.o cxdiv.o cxmul.o \
X cxphas.o cxphsr.o cxscal.o cxsqrt.o cxsub.o
X
Xall: $(ALL)
X
Xinstall: $(ALL)
X @echo install $(ALL) according to local convention.
X
Xcx_test: cx_test.c libcomplex.a
X $(CC) $(CFLAGS) -o cx_test cx_test.c libcomplex.a
X
Xlibcomplex.a: $(OBJS)
X ar r libcomplex.a $(OBJS)
X
X$(OBJS): complex.h
END_OF_Makefile
if test 395 -ne `wc -c <Makefile`; then
echo shar: \"Makefile\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f complex.3 -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"complex.3\"
else
echo shar: Extracting \"complex.3\" \(5830 characters\)
sed "s/^X//" >complex.3 <<'END_OF_complex.3'
X'\" e
X.TH COMPLEX 3V LOCAL
X'\" last edit: 86/02/03 D A Gwyn
X'\" SCCS ID: @(#)complex.3 1.2 (modified for public version)
X.EQ
Xdelim @@
X.EN
X.SH NAME
Xcomplex \- complex arithmetic operations
X.SH SYNOPSIS
X.B
X#include <complex.h> /* assuming appropriate cc \-I option */
X.br
X/* All the following functions are declared in this header file. */
X.P
X.B complex *CxAdd(ap,bp);
X.br
X.B complex *ap, *bp;
X.P
X.B complex *CxSub(ap,bp);
X.br
X.B complex *ap, *bp;
X.P
X.B complex *CxMul(ap,bp);
X.br
X.B complex *ap, *bp;
X.P
X.B complex *CxDiv(ap,bp);
X.br
X.B complex *ap, *bp;
X.P
X.B complex *CxSqrt(cp);
X.br
X.B complex *cp;
X.P
X.B complex *CxScal(cp,\^s);
X.br
X.B complex *cp;
X.br
X.B double s;
X.P
X.B complex *CxNeg(cp);
X.br
X.B complex *cp;
X.P
X.B complex *CxConj(cp);
X.br
X.B complex *cp;
X.P
X.B complex *CxCopy(ap,bp);
X.br
X.B complex *ap, *bp;
X.P
X.B complex *CxCons(cp,\^r,\^i);
X.br
X.B complex *cp;
X.br
X.B double r, i;
X.P
X.B complex *CxPhsr(cp,m,p);
X.br
X.B complex *cp;
X.br
X.B double m, p;
X.P
X.B double CxReal(cp);
X.br
X.B complex *cp;
X.P
X.B double CxImag(cp);
X.br
X.B complex *cp;
X.P
X.B double CxAmpl(cp);
X.br
X.B complex *cp;
X.P
X.B double CxPhas(cp);
X.br
X.B complex *cp;
X.P
X.B complex *CxAllo(\ );
X.P
X.B void CxFree(cp);
X.br
X.B complex *cp;
X.SH DESCRIPTION
XThese routines perform arithmetic
Xand other useful operations on complex numbers.
XAn appropriate data structure
X.B complex
Xis defined in the header file;
Xall access to
X.B complex
Xdata should be
X.I via
Xthese predefined functions.
X(See
X.SM HINTS
Xfor further information.)
X.P
XIn the following descriptions,
Xthe names
X.IR a ,
X.IR b ,
Xand
X.I c
Xrepresent the
X.B complex
Xdata addressed by the corresponding pointers
X.IR ap ,
X.IR bp ,
Xand
X.IR cp .
X.P
X.I CxAdd\^
Xadds
X.I b
Xto
X.I a
Xand returns a pointer to the result.
X.P
X.I CxSub
Xsubtracts
X.I b
Xfrom
X.I a
Xand returns a pointer to the result.
X.P
X.I CxMul\^
Xmultiplies
X.I a
Xby
X.I b
Xand returns a pointer to the result.
X.P
X.I CxDiv
Xdivides
X.I a
Xby
X.I b
Xand returns a pointer to the result.
XThe divisor must not be precisely zero.
X.P
X.I CxSqrt
Xreplaces
X.I c
Xby the ``principal value'' of its square root
X(one having a non-negative imaginary part)
Xand returns a pointer to the result.
X.P
X.I CxScal\^
Xmultiplies
X.I c
Xby the scalar
X.I s
Xand returns a pointer to the result.
X.P
X.I CxNeg
Xnegates
X.I c
Xand returns a pointer to the result.
X.P
X.I CxConj
Xconjugates
X.I c
Xand returns a pointer to the result.
X.P
X.I CxCopy
Xassigns the value of
X.I b
Xto
X.I a
Xand returns a pointer to the result.
X.P
X.I CxCons
Xconstructs the complex number
X.I c
Xfrom its real and imaginary parts
X.I r
Xand
X.IR i ,
Xrespectively,
Xand returns a pointer to the result.
X.P
X.I CxPhsr
Xconstructs the complex number
X.I c
Xfrom its ``phasor'' amplitude and phase (given in radians)
X.I m
Xand
X.IR p ,
Xrespectively,
Xand returns a pointer to the result.
X.P
X.I CxReal\^
Xreturns the real part of the complex number
X.IR c .
X.P
X.I CxImag
Xreturns the imaginary part of the complex number
X.IR c .
X.P
X.I CxAmpl\^
Xreturns the amplitude of the complex number
X.IR c .
X.P
X.I CxPhas
Xreturns the phase of the complex number
X.IR c ,
Xas radians in the range @(- pi , pi ]@.
X.P
X.I CxAllo
Xallocates storage for a
X.B complex
Xdatum; it returns
X.SM
X.B NULL
X(defined as 0 in
X.BR <stdio.h> )
Xif not enough storage is available.
X.P
X.I CxFree
Xreleases storage previously allocated by
X.IR CxAllo .
XThe contents of such storage must not be used afterward.
X.SH HINTS
XThe
X.B complex
Xdata type consists of real and imaginary components;
X.I CxReal\^
Xand
X.I CxImag
Xare actually macros that access these components directly.
XThis allows addresses of the components to be taken,
Xas in the following \s-1EXAMPLE\s0.
X.P
XThe complex functions are designed to be nested;
Xsee the following \s-1EXAMPLE\s0.
XFor this reason,
Xmany of them modify the contents of their first parameter.
X.I CxCopy
Xcan be used to create a ``working copy'' of
X.B complex
Xdata that would otherwise be modified.
X.P
XThe square-root function is inherently double-valued;
Xin most applications, both roots should receive equal consideration.
XThe second root is the negative of the ``principal value''.
X.bp
X.SH EXAMPLE
XThe following program is compiled by the command
X.br
X $ \fIcc \|\-I/usr/local/include \|example.c \|/usr/local/lib/libcomplex.a \|\-lm\fP
X.br
XIt reads in two complex vectors,
Xthen computes and prints their inner product.
X.sp
X.P
X #include <stdio.h>
X.br
X #include <complex.h>
X.sp
X main( argc, argv )
X.br
X int argc;
X.br
X char *argv[\|];
X.br
X {
X.br
X int n; /* # elements in each array */
X.br
X int i; /* indexes arrays */
X.br
X complex a[10], b[10]; /* input vectors */
X.br
X complex s; /* accumulates scalar product */
X.br
X complex *c = CxAllo(\|); /* holds cross-term */
X.sp
X if ( c == NULL )
X.br
X {
X.br
X (void)fprintf( stderr, ``not enough memory\en'' );
X.br
X return 1;
X.br
X }
X.br
X (void)printf( ``\enenter number of elements: '' );
X.br
X (void)scanf( `` %d'', &n );
X.br
X /* (There really should be some input validation here.) */
X.br
X (void) printf( ``\enenter real, imaginary pairs for first array:\en'' );
X.br
X for ( i = 0; i < n; ++i )
X.br
X (void)scanf( `` %lg %lg'', &CxReal( &a[i] ), &CxImag( &a[i] ) );
X.br
X (void)printf( ``\enenter real, imaginary pairs for second array:\en'' );
X.br
X for ( i = 0; i < n; ++i )
X.br
X (void)scanf( `` %lg %lg'', &CxReal( &b[i] ), &CxImag( &b[i] ) );
X.br
X (void)CxCons( &s, 0.0, 0.0 ); /* initialize accumulator */
X.br
X for ( i = 0; i < n; ++i )
X.br
X (void)CxAdd( &s, CxMul( &a[i], CxConj( CxCopy( c, &b[i] ) ) ) );
X.br
X (void)printf( ``\enproduct is (%g,%g)\en'', CxReal( &s ), CxImag( &s ) );
X.br
X CxFree( c );
X.br
X return 0;
X.br
X }
X.SH FILES
X/usr/local/include/complex.h header file containing definitions
X.br
X/usr/local/lib/libcomplex.a complex run-time support library
X.SH AUTHORS
XDouglas A. Gwyn, BRL/VLD-VMB
X.br
XJeff Hanes, BRL/VLD-VMB (original version of
X.IR CxSqrt\^ )
END_OF_complex.3
if test 5830 -ne `wc -c <complex.3`; then
echo shar: \"complex.3\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f complex.h -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"complex.h\"
else
echo shar: Extracting \"complex.h\" \(966 characters\)
sed "s/^X//" >complex.h <<'END_OF_complex.h'
X/*
X <complex.h> -- definitions for complex arithmetic routines
X
X last edit: 86/01/04 D A Gwyn
X
X SCCS ID: @(#)complex.h 1.1 (modified for public version)
X*/
X
X/* "complex number" data type: */
X
Xtypedef struct
X {
X double re; /* real part */
X double im; /* imaginary part */
X } complex;
X
X/* "The future is now": */
X
X#ifdef __STDC__ /* X3J11 */
X#define _CxGenPtr void * /* generic pointer type */
X#else /* K&R */
X#define _CxGenPtr char * /* generic pointer type */
X#endif
X
X/* functions that are correctly done as macros: */
X
X#define CxAllo() ((complex *)malloc( sizeof (complex) ))
X#define CxFree( cp ) free( (_CxGenPtr)(cp) )
X#define CxNeg( cp ) CxScal( cp, -1.0 )
X#define CxReal( cp ) (cp)->re
X#define CxImag( cp ) (cp)->im
X
Xextern void free();
Xextern _CxGenPtr malloc();
X
X/* library functions: */
X
Xextern double CxAmpl(), CxPhas();
Xextern complex *CxAdd(), *CxConj(), *CxCons(), *CxCopy(), *CxDiv(),
X *CxMul(), *CxPhsr(), *CxScal(), *CxSqrt(), *CxSub();
END_OF_complex.h
if test 966 -ne `wc -c <complex.h`; then
echo shar: \"complex.h\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f cx_test.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"cx_test.c\"
else
echo shar: Extracting \"cx_test.c\" \(4250 characters\)
sed "s/^X//" >cx_test.c <<'END_OF_cx_test.c'
X/*
X ctest -- complex arithmetic test
X
X last edit: 86/01/04 D A Gwyn
X
X SCCS ID: @(#)cx_test.c 1.1 (modified for public version)
X*/
X
X#include <stdio.h>
X#include <math.h>
X
X#include <complex.h>
X
X#define DEGRAD 57.2957795130823208767981548141051703324054724665642
X /* degrees per radian */
X#define Abs( x ) ((x) < 0 ? -(x) : (x))
X#define Max( a, b ) ((a) > (b) ? (a) : (b))
X
Xextern void exit();
X
X#define Printf (void)printf
X
X#define TOL 1.0e-10 /* tolerance for checks */
X
Xstatic int errs = 0; /* tally errors */
X
Xstatic void CCheck(), RCheck();
Xstatic double RelDif();
X
X
X/*ARGSUSED*/
Xmain( argc, argv )
X int argc;
X char *argv[];
X {
X complex a, *bp, *cp;
X
X /* CxAllo test */
X bp = CxAllo();
X if ( bp == NULL )
X {
X Printf( "CxAllo failed\n" );
X exit( 1 );
X }
X
X /* CxReal, CxImag test */
X CxReal( bp ) = 1.0;
X CxImag( bp ) = 2.0;
X RCheck( "CxReal", CxReal( bp ), 1.0 );
X RCheck( "CxImag", CxImag( bp ), 2.0 );
X
X /* CxCons test */
X cp = CxCons( &a, -3.0, -4.0);
X CCheck( "CxCons 1", a, -3.0, -4.0 );
X CCheck( "CxCons 2", *cp, -3.0, -4.0 );
X
X /* CxNeg test */
X cp = CxNeg( &a );
X CCheck( "CxNeg 1", a, 3.0, 4.0 );
X CCheck( "CxNeg 2", *cp, 3.0, 4.0 );
X
X /* CxCopy test */
X cp = CxCopy( bp, &a );
X (void)CxCons( &a, 1.0, sqrt( 3.0 ) );
X CCheck( "CxCopy 1", *bp, 3.0, 4.0 );
X CCheck( "CxCopy 2", *cp, 3.0, 4.0 );
X
X /* CxAmpl, CxPhas test */
X RCheck( "CxAmpl 1", CxAmpl( &a ), 2.0 );
X RCheck( "CxPhas 1", CxPhas( &a ) * DEGRAD, 60.0 );
X /* try other quadrants */
X a.re = -a.re;
X RCheck( "CxAmpl 2", CxAmpl( &a ), 2.0 );
X RCheck( "CxPhas 2", CxPhas( &a ) * DEGRAD, 120.0 );
X a.im = -a.im;
X RCheck( "CxAmpl 3", CxAmpl( &a ), 2.0 );
X RCheck( "CxPhas 3", CxPhas( &a ) * DEGRAD, -120.0 );
X a.re = -a.re;
X RCheck( "CxAmpl 4", CxAmpl( &a ), 2.0 );
X RCheck( "CxPhas 4", CxPhas( &a ) * DEGRAD, -60.0 );
X /* one more for good measure */
X RCheck( "CxAmpl 5", CxAmpl( bp ), 5.0 );
X
X /* CxPhsr test */
X cp = CxPhsr( &a, 100.0, -20.0 / DEGRAD );
X RCheck( "CxPhsr 1", CxAmpl( &a ), 100.0 );
X RCheck( "CxPhsr 2", CxPhas( &a ) * DEGRAD, -20.0 );
X RCheck( "CxPhsr 3", CxAmpl( cp ), 100.0 );
X RCheck( "CxPhsr 4", CxPhas( cp ) * DEGRAD, -20.0 );
X
X /* CxConj test */
X cp = CxConj( bp );
X CCheck( "CxConj 1", *bp, 3.0, -4.0 );
X CCheck( "CxConj 2", *cp, 3.0, -4.0 );
X
X /* CxScal test */
X cp = CxScal( bp, 2.0 );
X CCheck( "CxScal 1", *bp, 6.0, -8.0 );
X CCheck( "CxScal 2", *cp, 6.0, -8.0 );
X
X /* CxAdd test */
X cp = CxAdd( CxCons( &a, -4.0, 11.0 ), bp );
X CCheck( "CxAdd 1", a, 2.0, 3.0 );
X CCheck( "CxAdd 2", *cp, 2.0, 3.0 );
X
X /* CxSub test */
X cp = CxSub( CxCons( &a, 4.0, 7.0 ), bp );
X CCheck( "CxSub 1", a, -2.0, 15.0 );
X CCheck( "CxSub 2", *cp, -2.0, 15.0 );
X
X /* CxMul test */
X cp = CxMul( CxCons( bp, -1.0, 3.0 ), CxCons( &a, 1.0, 2.0 ) );
X CCheck( "CxMul 1", *bp, -7.0, 1.0 );
X CCheck( "CxMul 2", *cp, -7.0, 1.0 );
X
X /* CxDiv test */
X cp = CxDiv( bp, &a );
X CCheck( "CxDiv 1", *bp, -1.0, 3.0 );
X CCheck( "CxDiv 2", *cp, -1.0, 3.0 );
X
X /* CxSqrt and overlapping CxMul tests */
X (void)CxCons( &a, -1.0, 2.0 );
X cp = CxSqrt( CxMul( &a, &a ) );
X CCheck( "CxSqrt 1", a, -1.0, 2.0 );
X CCheck( "CxSqrt 2", *cp, -1.0, 2.0 );
X (void)CxCons( &a, 3.0, 2.0 );
X cp = CxSqrt( CxMul( &a, &a ) );
X CCheck( "CxSqrt 3", a, 3.0, 2.0 );
X CCheck( "CxSqrt 4", *cp, 3.0, 2.0 );
X
X /* CxFree "test" */
X CxFree( bp );
X
X return errs;
X }
X
X
Xstatic void
XCCheck( s, c, r, i ) /* check complex number */
X char *s; /* message string for failure */
X complex c; /* complex to be checked */
X double r, i; /* expected real, imaginary parts */
X {
X if ( RelDif( CxReal( &c ), r ) > TOL
X || RelDif( CxImag( &c ), i ) > TOL
X ) {
X ++errs;
X Printf( "%s; s.b. (%f,%f), was (%g,%g)\n",
X s, r, i, c.re, c.im
X );
X }
X }
X
X
Xstatic void
XRCheck( s, d, r ) /* check real number */
X char *s; /* message string for failure */
X double d; /* real to be checked */
X double r; /* expected value */
X {
X if ( RelDif( d, r ) > TOL )
X {
X ++errs;
X Printf( "%s; s.b. %f, was %g\n", s, r, d );
X }
X }
X
X
Xstatic double
XRelDif( a, b ) /* returns relative difference: */
X double a, b; /* 0.0 if exactly the same,
X otherwise ratio of difference
X to the larger of the two */
X {
X double c = Abs( a );
X double d = Abs( b );
X
X d = Max( c, d );
X
X return d == 0.0 ? 0.0 : Abs( a - b ) / d;
X }
END_OF_cx_test.c
if test 4250 -ne `wc -c <cx_test.c`; then
echo shar: \"cx_test.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f cxadd.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"cxadd.c\"
else
echo shar: Extracting \"cxadd.c\" \(304 characters\)
sed "s/^X//" >cxadd.c <<'END_OF_cxadd.c'
X/*
X CxAdd -- add one complex to another
X
X last edit: 86/01/04 D A Gwyn
X
X SCCS ID: @(#)cxadd.c 1.1
X
X CxAdd( &a, &b ) adds b to a and returns &a
X*/
X
X#include <complex.h>
X
Xcomplex *
XCxAdd( ap, bp )
X register complex *ap, *bp; /* may coincide */
X {
X ap->re += bp->re;
X ap->im += bp->im;
X
X return ap;
X }
END_OF_cxadd.c
if test 304 -ne `wc -c <cxadd.c`; then
echo shar: \"cxadd.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f cxampl.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"cxampl.c\"
else
echo shar: Extracting \"cxampl.c\" \(278 characters\)
sed "s/^X//" >cxampl.c <<'END_OF_cxampl.c'
X/*
X CxAmpl -- amplitude (magnitude, modulus, norm) of a complex
X
X CxAmpl( &c ) returns |c|
X
X last edit: 86/01/04 D A Gwyn
X
X SCCS ID: @(#)cxampl.c 1.1
X*/
X
X#include <math.h>
X
X#include <complex.h>
X
Xdouble
XCxAmpl( cp )
X register complex *cp;
X {
X return hypot( cp->re, cp->im );
X }
END_OF_cxampl.c
if test 278 -ne `wc -c <cxampl.c`; then
echo shar: \"cxampl.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f cxconj.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"cxconj.c\"
else
echo shar: Extracting \"cxconj.c\" \(278 characters\)
sed "s/^X//" >cxconj.c <<'END_OF_cxconj.c'
X/*
X CxConj -- conjugate a complex
X
X CxConj( &c ) conjugates c and returns &c
X
X last edit: 86/01/04 D A Gwyn
X
X SCCS ID: @(#)cxconj.c 1.1
X*/
X
X#include <complex.h>
X
Xcomplex *
XCxConj( cp )
X register complex *cp;
X {
X /* (real part unchanged) */
X cp->im = -cp->im;
X
X return cp;
X }
END_OF_cxconj.c
if test 278 -ne `wc -c <cxconj.c`; then
echo shar: \"cxconj.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f cxcons.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"cxcons.c\"
else
echo shar: Extracting \"cxcons.c\" \(329 characters\)
sed "s/^X//" >cxcons.c <<'END_OF_cxcons.c'
X/*
X CxCons -- construct a complex from real and imaginary parts
X
X CxCons( &c, re, im ) makes c = re + i im and returns &c
X
X last edit: 86/01/04 D A Gwyn
X
X SCCS ID: @(#)cxcons.c 1.1
X*/
X
X#include <complex.h>
X
Xcomplex *
XCxCons( cp, re, im )
X register complex *cp;
X double re, im;
X {
X cp->re = re;
X cp->im = im;
X
X return cp;
X }
END_OF_cxcons.c
if test 329 -ne `wc -c <cxcons.c`; then
echo shar: \"cxcons.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f cxcopy.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"cxcopy.c\"
else
echo shar: Extracting \"cxcopy.c\" \(264 characters\)
sed "s/^X//" >cxcopy.c <<'END_OF_cxcopy.c'
X/*
X CxCopy -- copy a complex
X
X last edit: 86/01/04 D A Gwyn
X
X SCCS ID: @(#)cxcopy.c 1.1
X
X CxCopy( &a, &b ) copies b to a and returns &a
X*/
X
X#include <complex.h>
X
Xcomplex *
XCxCopy( ap, bp )
X complex *ap, *bp; /* may coincide */
X {
X *ap = *bp;
X
X return ap;
X }
END_OF_cxcopy.c
if test 264 -ne `wc -c <cxcopy.c`; then
echo shar: \"cxcopy.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f cxdiv.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"cxdiv.c\"
else
echo shar: Extracting \"cxdiv.c\" \(820 characters\)
sed "s/^X//" >cxdiv.c <<'END_OF_cxdiv.c'
X/*
X CxDiv -- divide one complex by another
X
X CxDiv( &a, &b ) divides a by b and returns &a;
X zero divisor fails
X
X last edit: 86/01/04 D A Gwyn
X
X SCCS ID: @(#)cxdiv.c 1.1 (modified for public version)
X*/
X
X#include <complex.h>
X
X#define Abs( x ) ((x) < 0 ? -(x) : (x))
X
Xcomplex *
XCxDiv( ap, bp )
X register complex *ap, *bp; /* may coincide (?) */
X {
X double r, s;
X double ap__re = ap->re;
X
X /* Note: classical formula may cause unnecessary overflow */
X r = bp->re;
X s = bp->im;
X if ( Abs( r ) >= Abs( s ) )
X {
X r = s / r; /* <= 1 */
X s = bp->re + r * s;
X ap->re = (ap->re + ap->im * r) / s;
X ap->im = (ap->im - ap__re * r) / s;
X }
X else /* Abs( s ) > Abs( r ) */
X {
X r = r / s; /* < 1 */
X s = s + r * bp->re;
X ap->re = (ap->re * r + ap->im) / s;
X ap->im = (ap->im * r - ap__re) / s;
X }
X
X return ap;
X }
END_OF_cxdiv.c
if test 820 -ne `wc -c <cxdiv.c`; then
echo shar: \"cxdiv.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f cxmul.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"cxmul.c\"
else
echo shar: Extracting \"cxmul.c\" \(424 characters\)
sed "s/^X//" >cxmul.c <<'END_OF_cxmul.c'
X/*
X CxMul -- multiply one complex by another
X
X CxMul( &a, &b ) multiplies a by b and returns &a
X
X last edit: 86/01/04 D A Gwyn
X
X SCCS ID: @(#)cxmul.c 1.1
X*/
X
X#include <complex.h>
X
Xcomplex *
XCxMul( ap, bp )
X register complex *ap, *bp; /* (may coincide) */
X {
X double ap__re = ap->re;
X double bp__re = bp->re;
X
X ap->re = ap__re * bp__re - ap->im * bp->im;
X ap->im = ap__re * bp->im + ap->im * bp__re;
X
X return ap;
X }
END_OF_cxmul.c
if test 424 -ne `wc -c <cxmul.c`; then
echo shar: \"cxmul.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f cxphas.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"cxphas.c\"
else
echo shar: Extracting \"cxphas.c\" \(406 characters\)
sed "s/^X//" >cxphas.c <<'END_OF_cxphas.c'
X/*
X CxPhas -- phase (angle, argument) of a complex
X
X CxPhas( &c ) returns arg(c) in radians (-Pi,Pi]
X
X last edit: 86/01/04 D A Gwyn
X
X SCCS ID: @(#)cxphas.c 1.1 (modified for public version)
X*/
X
X#include <math.h>
X
X#include <complex.h>
X
Xdouble
XCxPhas( cp )
X register complex *cp;
X {
X if ( cp->re == 0.0 && cp->im == 0.0 )
X return 0.0; /* can't trust atan2() */
X else
X return atan2( cp->im, cp->re );
X }
END_OF_cxphas.c
if test 406 -ne `wc -c <cxphas.c`; then
echo shar: \"cxphas.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f cxphsr.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"cxphsr.c\"
else
echo shar: Extracting \"cxphsr.c\" \(395 characters\)
sed "s/^X//" >cxphsr.c <<'END_OF_cxphsr.c'
X/*
X CxPhsr -- construct a complex "phasor" from amplitude and phase
X
X CxPhsr( &c, amp, phs ) makes c = amp exp(i phs)
X and returns &c
X
X last edit: 86/01/04 D A Gwyn
X
X SCCS ID: @(#)cxphsr.c 1.1
X*/
X
X#include <math.h>
X
X#include <complex.h>
X
Xcomplex *
XCxPhsr( cp, amp, phs )
X register complex *cp;
X double amp, phs;
X {
X cp->re = amp * cos( phs );
X cp->im = amp * sin( phs );
X
X return cp;
X }
END_OF_cxphsr.c
if test 395 -ne `wc -c <cxphsr.c`; then
echo shar: \"cxphsr.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f cxscal.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"cxscal.c\"
else
echo shar: Extracting \"cxscal.c\" \(291 characters\)
sed "s/^X//" >cxscal.c <<'END_OF_cxscal.c'
X/*
X CxScal -- multiply a complex by a scalar
X
X CxScal( &c, s ) scales c by s and returns &c
X
X last edit: 86/01/04 D A Gwyn
X
X SCCS ID: @(#)cxscal.c 1.1
X*/
X
X#include <complex.h>
X
Xcomplex *
XCxScal( cp, s )
X register complex *cp;
X double s;
X {
X cp->re *= s;
X cp->im *= s;
X
X return cp;
X }
END_OF_cxscal.c
if test 291 -ne `wc -c <cxscal.c`; then
echo shar: \"cxscal.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f cxsqrt.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"cxsqrt.c\"
else
echo shar: Extracting \"cxsqrt.c\" \(1339 characters\)
sed "s/^X//" >cxsqrt.c <<'END_OF_cxsqrt.c'
X/*
X CxSqrt -- compute square root of complex number
X
X CxSqrt( &c ) replaces c by sqrt(c) and returns &c
X
X Note: This is a double-valued function; the result of
X CxSqrt() always has nonnegative imaginary part.
X
X inspired by Jeff Hanes' version
X
X last edit: 86/01/04 D A Gwyn
X
X SCCS ID: @(#)cxsqrt.c 1.1 (modified for public version)
X*/
X
X#include <math.h>
X
X#include <complex.h>
X
X#define Sgn( x ) ((x) == 0 ? 0 : (x) > 0 ? 1 : -1)
X
Xcomplex *
XCxSqrt( cp )
X register complex *cp;
X {
X /* record signs of original real & imaginary parts */
X int re_sign = Sgn( cp->re );
X int im_sign = Sgn( cp->im );
X
X /* special cases are not necessary; they are here for speed */
X
X if ( re_sign == 0 )
X if ( im_sign == 0 )
X ; /* (0,0) already there */
X else if ( im_sign > 0 )
X cp->re = cp->im = sqrt( cp->im / 2.0 );
X else /* im_sign < 0 */
X cp->re = -(cp->im = sqrt( -cp->im / 2.0 ));
X else if ( im_sign == 0 )
X if ( re_sign > 0 )
X cp->re = sqrt( cp->re );
X/* cp->im = 0.0; /* 0 already there */
X else { /* re_sign < 0 */
X cp->im = sqrt( -cp->re );
X cp->re = 0.0;
X }
X else { /* no shortcuts */
X double ampl = CxAmpl( cp );
X
X cp->im = sqrt( (ampl - cp->re) /2.0 );
X
X if ( im_sign > 0 )
X cp->re = sqrt( (ampl + cp->re) / 2.0 );
X else /* im_sign < 0 */
X cp->re = -sqrt( (ampl + cp->re) / 2.0 );
X }
X
X return cp;
X }
END_OF_cxsqrt.c
if test 1339 -ne `wc -c <cxsqrt.c`; then
echo shar: \"cxsqrt.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f cxsub.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"cxsub.c\"
else
echo shar: Extracting \"cxsub.c\" \(320 characters\)
sed "s/^X//" >cxsub.c <<'END_OF_cxsub.c'
X/*
X CxSub -- subtract one complex from another
X
X CxSub( &a, &b ) subtracts b from a and returns &a
X
X last edit: 86/01/04 D A Gwyn
X
X SCCS ID: @(#)cxsub.c 1.1
X*/
X
X#include <complex.h>
X
Xcomplex *
XCxSub( ap, bp )
X register complex *ap, *bp; /* (may coincide) */
X {
X ap->re -= bp->re;
X ap->im -= bp->im;
X
X return ap;
X }
END_OF_cxsub.c
if test 320 -ne `wc -c <cxsub.c`; then
echo shar: \"cxsub.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
echo shar: End of shell archive.
exit 0
--
Rich $alz "Anger is an energy"
Cronus Project, BBN Labs rsalz@bbn.com
Moderator, comp.sources.unix sources@uunet.uu.net
--
Rich $alz "Anger is an energy"
Cronus Project, BBN Labs rsalz@bbn.com
Moderator, comp.sources.unix sources@uunet.uu.net