Louis Baker has a Ph.D. in astronomy and has written books with code in C and Ada for McGraw-Hill. He may be reached at Mission Research Corp., 1720 Randolph SE, Albuquerque NM 87106.
FORTRAN has long been the language of choice in the scientific and engineering community. (The name FORTRAN — short for FORmula TRANslator — betrays its origins.) Over the years, FORTRAN has been extended to meet the needs of its user community. Today's FORTRANs routinely handle complex arithmetic, multi-dimensional arrays, and sophisticated graphics.
C, on the other hand, has evolved to support its user base, systems programmers who benefit from C's easy access to low-level constructs such as register variables and absolute addresses. Unfortunately, there doesn't seem to be much effort to incorporate features like complex numbers into the ANSI C standard. Most FORTRAN programmers also decry C's inability to pass a subroutine both the dimensions and values of a multi-dimensional array. Yet, there are compelling reasons for writing "engineering" code in C, especially for real-time control or data acquisition applications. C, after all, makes it easy to access ports, take advantage of DMA, and so on.
In this article I'll show you one approach, namely the use of header files, that adds complex number and matrix support to your C programs.
Roots
The header files that I'll discuss here were originally developed for a book recently published by McGraw-Hill, C Tools for Scientists and Engineers [1]. I found them quite useful for the problems addressed in that book."Why bother with C?" you may ask, evoking a mental image of Louis Baker as a Neanderthal. "Why not just use C++ and make things a lot simpler?" Although C++ provides no more support for complex numbers or multi-dimensional matrices than C does, you can at least define classes of objects such as complex numbers or matrices. In fact, new versions of C++ that support multiple inheritance make it relatively easy to define complex matrices in terms of other C++ definitions. However, C++ is less widely available than C and is even less standardized. One of C's claims to fame is its portability, but today's C++s don't insure that a port to another host or compiler will work the same way. But yes, C++ is certainly a viable option.
Complex Numbers For The Uninitiated
Complex numbers are generally represented in the form
x + iywhere i is the square root of -1. The x is usually called the "real" part of the complex number and iy the "imaginary" part. Complex numbers crop up in mathematics, for example, as the roots of polynomial equations. Every polynomial has a complex root. This means that you can plug the complex root into the polynomial for x and make the polynomial evaluate to zero. Remember factoring polynomials like x2 - 4x + 4 into (x - 2)(x - 2)? Factoring represented an easy method for solving polynomials. Factoring a higher order polynomial f(x) often involves rewriting it as (x-r)g(x) where r is a root and g(x) is a polynomial of lower degree. However, the root r may very well be a complex number — not just a "real" number such as 4.67.Complex numbers appear in engineering mostly in equations that use the complex exponential. The relationship
e(ix) = cos(x) + i sin (x)where e is 2.718281828459045... is often used to analyze AC circuits and other systems with periodic behavior. A quick glance at the equation should convince you that it enables the cosine function to be written as the real part of e(ix). So what? Well, it's often simpler to use the properties of exponents to manipulate equations than to rely on trigonometric functions.
Complex Arithmetic In C
The rules for adding, subtracting, multiplying, and dividing complex numbers can be found in most elementary algebra texts. Examine Listing 1 to see how you might write a header file to define functions that perform these common operations on complex numbers. Note that addition and subtraction of complex numbers just requires that you treat the real and complex parts of the numbers separately. Multiplication and division, as you might expect, are a bit more involved. The product
(a+bi)(c+di) = (ac-bd)+(ad+bc)ifor example, can be implemented either in a straightforward manner by CMULT() or in a more involved manner by CMLT(). The latter saves a multiplication, but the price is more additions and more storage of intermediate results. Actually, life's even more complicated than that. On some machines, particularly those that have floating-point operations performed by software emulation, CMLT() may very well be faster than CMULT(). But even if your system has a math coprocessor, where multiplication is only slightly more costly in time than addition, CMULT() should still be faster than CMLT().My method for determining the magnitude or "absolute value" of a complex number is also straightforward:
|x+iy| = sqrt(x2+y2).Some packages minimize the possibility of roundoff error by calculating this as
|x+iy| = a sqrt( 1 +(b/a)2)where
a = max(|x|,|y|)and
b = min (|x|, |y|)Listing 2 shows a C program that uses the header file definitions to calculate the Kelvin functions from approximations presented in the Handbook of Mathematical Functions [3]. The Kelvin functions are Bessel functions for which the argument is a complex number of the form sz, where z is real and s is a square root of i. The imaginary unit i has two square roots, (1+i)/sqrt(2) and (1-i)/sqrt(2). The Kelvin functions are useful, for example, in analyzing the distribution of alternating current within a wire. Listing 2 also contains two auxiliary routines, printc(), for printing complex numbers, and cexp(), for finding ez where z is a complex number.If this were C++, you could overload the "+" operator for the class of complex objects rather than use something like Listing 1's definition for complex addition, CADD(). You could use overloading in a similar fashion for most of the other operators, making complex arithmetic look much more like ordinary arithmetic. C++ would also allow you to specify that the code should be in-line, avoiding the overhead of subroutine calls.
My C Tools book uses an expanded version of this complex header file to implement root finders and an FFT (Fast Fourier Transform) package. It contains functions for the complex logarithm, conversion between the rectangular form for representing complex numbers (x + iy) and polar form
r eiqwhere r and q are real numbers, r the magnitude and q the "argument" of the complex number x+iy. I've also used the package in a number of other routines, some of which will appear in a forthcoming book called More C Tools. My pet project is to implement all the functions discussed in Abramowitz and Stegun [3] in C.
Matrices
If you use the header file given in Listing 3, you'll discover that it makes using two-dimensional arrays, or matrices, much more akin to their use in mathematics and FORTRAN. Recall that for a matrix A(m,n), the notation A(i,j) conventionally refers to the element in the i-th row and the j-th column. The matrix A declared in FORTRAN as
DIMENSION A(2,3)would have its elements stored sequentially by columns
A(1,1) A(2,1) A(1,2) A(2,2) A(1,3) A(2,3)In C, however, the declaration
double a[2][3];would give us a representation in storage of the form
a[0][0] a[0][1] a[0][2] a[1][0] a[1][1] a[1][2]The matrix in C is stored by rows, and the indexing is based at zero instead of one as the index of the lowest element. The header file supplies the INDEX() definition, among others, to allow matrix a to be referenced with a mechanism somewhat like the mathematical or FORTRAN notation, but with C conventions as to indexing starting at base zero and storage by rows. The example in Listing 4 illustrates how INDEX is used to multiply a vector by a matrix.Obviously, you can adapt INDEX to suit your taste. If you want to use indexing that starts with first rather than zeroth elements, all you have to do is replace i and j by (i-1) and (j-1). To store by columns, interchange i and j and replace the variable coln — the maximum number of columns of the matrix — by rown, the number of rows of the matrix.
The header file also supplies some constructs that make translating do loops from FORTRAN into C a little less painful.
Closing Remarks
I hope these header files and examples are of some use to you. I don't pretend that they're perfect solutions to some of the complex problems your applications need to address, but they work in most cases. Maybe they'll give you the ammunition you need to convince colleagues that C deserves a chance.[1] W. H. Press, B. P. Flannery, S. A. Teukolsky, W. T. Vetterling, Numerical Recipes In C, Cambridge University Press, 1986.
[2] L. Baker, C Tools for Scientists and Engineers, McGraw-Hill, 1989.
[3] M. Abramowitz and I. Stegun, Handbook of Mathematical Functions, U. S. Dept. of Commerce Applied Math Series Publ. 55 (also available as a paperback from Dover Publications), 1964.