Mathematics
f(z)    Complex Numbers and Functions     Tech Note
Complex Numbers

Humans invented several number systems through the ages, including

  1. Set of integers constructed from the counting numbers, 1, 2, 3, …
  2. Rational numbers, m/n, as ratio of integers.
  3. Real numbers, x, including rational numbers (e.g., ½) and irrational numbers (e.g., ).
  4. Complex numbers, a + ib, where a is the “real part” and ib is the “imaginary part.”

 Every “new” set of invented numbers extended the previous system.

These web pages are a brief overview of complex numbers and complex math, and describe a ComplexMathLibrary unit that has a number of functions useful for numerical experiments.  Unlike many other complex math libraries, many algorithms in ComplexMathLibrary are implemented in both rectangular and polar forms, since some algorithms are easier to implement in one form than the other.

This page gives an introduction to imaginary and complex numbers both in rectangular and polar form, and an overview of the software implementation.  Other pages explain complex arithmetic, complex logarithms and powers, complex trig functions, complex hyperbolic functions, complex Bessel functions, and the complex gamma function.

For more information about complex variables, analysis and functions, see these books:

“Imaginary” Numbers

The solution to the equation, x2 = 1, is

But what is the solution to the equation, x2 = -1?  By definition, we call the “imaginary” number “i” the solution to this equation:

(Electrical engineers often use a “j” instead of an “i” to avoid confusion with the abbreviation “i” being used as “electric current” in equations such as Ohm’s Law.)

Imaginary numbers are “different” from real numbers and cannot be combined with real numbers.  So, 1 + 2 = 3, and i + 2i = 3i.  But 1 + 2i is something else, namely a “complex” number.

Complex Numbers

An ordered pair, (x, y), of a real number, x, and an imaginary number, iy, is called a “complex” number.  A complex number, z, can be written as z = (x, y), or z = x + iy, (or in other forms, which will be discussed below). 

The complex plane is a two-dimensional space, where each point can be identified in (x, y) rectangular coordinates, or (r, q) polar coordinates.

Rectangular Form

The absolute value (or sometimes called modulus or “norm”) of the complex value z is defined as follows:

Results from CAbs(z) function

Absolute value of complex number: 
CAbs = ABS(z)

     z                ABS(z)
------------  -------------------
 0.0 +  0.0i    0.000000000000000 
 0.5 +  0.5i    0.707106781186548 
-0.5 +  0.5i    0.707106781186548 
-0.5 -  0.5i    0.707106781186548 
 0.5 -  0.5i    0.707106781186548 
 1.0 +  0.0i    1.000000000000000 
 1.0 +  1.0i    1.414213562373095 
 0.0 +  1.0i    1.000000000000000 
-1.0 +  1.0i    1.414213562373095 
-1.0 +  0.0i    1.000000000000000 
-1.0 -  1.0i    1.414213562373095 
 0.0 -  1.0i    1.000000000000000 
 1.0 -  1.0i    1.414213562373095 
 5.0 +  0.0i    5.000000000000000 
 5.0 +  3.0i    5.830951894845301 
 0.0 +  3.0i    3.000000000000000 
-5.0 +  3.0i    5.830951894845301 
-5.0 +  0.0i    5.000000000000000 
-5.0 -  3.0i    5.830951894845301 
 0.0 -  3.0i    3.000000000000000 
-5.0 -  3.0i    5.830951894845301

Polar Form

Instead of (x,y) rectangular coordinates, all points can be represented in polar form with (r, q) coordinates as shown in the diagram below:

Normally the angle is limited to a symmetric range, -p < q £ p, instead of 0 < q £ 2p.

Sometimes the form, z = r[cos q + i·sin q], is abbreviated as z = r·cis q.

Software Implementation

There are several approaches that can be used in implementing a ComplexMathLibrary unit for Delphi/Kylix.  I view a complex data type much like I would an integer or real numeric data type:  a low-level primitive type.  Some may choose to use an “object oriented” approach and implement a TComplex data type as an “object.”   I prefer to use such an object-oriented approach only for a “macro” object (perhaps a matrix object), and use a more primitive “Record” to define a TComplex type, since from a mathematical standpoint a complex value “behaves” much like a primitive floating-point value.

Some complex math algorithms can be expressed more easily using polar coordinates instead of rectangular coordinates.  The TComplex data type shown below allows definition of a complex value in either rectangular or polar coordinates.  Conversions between the types are forced when necessary (like to add two complex numbers in polar form), but otherwise the programmer can control what format is appropriate for a given algorithm.  The “Case” statement in the “Record” definition allows either (x, y) rectangular coordinates, or (r, q) polar coordinates to be stored in the same data structure.

TComplex data type

TYPE
  TReal        = DOUBLE;
  TComplexForm = (cfPolar, cfRectangular);
 
  TComplex     =
    RECORD
      CASE form:  TComplexForm OF
        cfRectangular:  (x,y    :  TReal);  // z = x + i*y
        cfPolar      :  (r,theta:  TReal);  // z = r*CIS(theta)
    END;                // where CIS(theta) = COS(theta) + i*SIN(theta)

                        //       theta = -PI..PI (in canonical form)

The assignment of a TComplex variable is usually via a CSet function:

CSet function

FUNCTION CSet (CONST a,b:  TReal; CONST f:  TComplexForm):  TComplex;
BEGIN
  RESULT.form := f;

  CASE f OF
    cfRectangular:
      BEGIN
        RESULT.x := a;
        RESULT.y := b
      END;

    cfPolar:
      BEGIN
        RESULT.r := a;
        RESULT.theta := b
      END
  END
END {CSet};

The function prototype for CSet defines a default TComplexForm of cfRectangular.

The array, a, of complex values used in many of the ComplexMathVerification project routines was defined as follows:

Definition of complex array

VAR a:  ARRAY[1..21] OF TComplex;
…
a[ 1] := CSet(  0.0,  0.0);
a[ 3] := CSet( -0.5,  0.5);
a[ 4] := CSet( -0.5, -0.5);
a[ 5] := CSet(  0.5, -0.5);
a[ 6] := CSet(  1.0,  0.0);
a[ 7] := CSet(  1.0,  1.0);
a[ 8] := CSet(  0.0,  1.0);
a[ 9] := CSet( -1.0,  1.0);
a[12] := CSet(  0.0, -1.0);
a[13] := CSet(  1.0, -1.0);
a[16] := CSet(  0.0,  3.0);
a[17] := CSet( -5.0,  3.0);
a[19] := CSet( -5.0, -3.0);
a[20] := CSet(  0.0, -3.0);
a[21] := CSet( -5.0, -3.0);
 
// "Polar" definitions
a[ 2] := CSet(SQRT(0.5),       PI/4, cfPolar);   //  0.5 + 0.5i
a[10] := CSet(     1.0,          PI, cfPolar);   // -1   + 0i
a[11] := CSet( SQRT(2),      5*PI/4, cfPolar);   // -1   -  i
a[14] := CSet(     5.0,         0.0, cfPolar);   //  5   + 0i
a[15] := CSet(SQRT(34), ArcTan(3/5), cfPolar);   //  5   + 3i
a[18] := CSet(     5.0,          PI, cfPolar);   // -5   + 0i

Such complex values can be displayed in either rectangular or polar form using the CToRectStr and CToPolarStr functions (explained below).

Code to display array of complex values using CToRectStr and CToPolarStr

Memo.Lines.Add('Complex number definition and conversion:  ' +
                     'CSet, CtoRectStr, CtoPolarStr');
Memo.Lines.Add('');
Memo.Lines.Add('             rectangular          ' +
               '              polar');
Memo.Lines.Add('     ----------------------------   ' +
               '--------------------------------');
FOR k := Low(a) TO High(a) DO
BEGIN
  Memo.Lines.Add(   Format('%2d  %-s %-s',
                            [k, CToRectStr(a[k],  13,9),
                                CToPolarStr(a[k], 13,9)]))
END;

The code above results in the follow TMemo data:

Array of TComplex values

Complex number definition and conversion:  CSet, CtoRectStr, CtoPolarStr
 
             rectangular                        polar
     ----------------------------   --------------------------------
 1    0.000000000 +  0.000000000i   0.000000000 * CIS(  0.000000000)
 2    0.500000000 +  0.500000000i   0.707106781 * CIS(  0.785398163)
 3   -0.500000000 +  0.500000000i   0.707106781 * CIS(  2.356194490)
 4   -0.500000000 -  0.500000000i   0.707106781 * CIS( -2.356194490)
 5    0.500000000 -  0.500000000i   0.707106781 * CIS( -0.785398163)
 6    1.000000000 +  0.000000000i   1.000000000 * CIS(  0.000000000)
 7    1.000000000 +  1.000000000i   1.414213562 * CIS(  0.785398163)
 8    0.000000000 +  1.000000000i   1.000000000 * CIS(  1.570796327)
 9   -1.000000000 +  1.000000000i   1.414213562 * CIS(  2.356194490)
10   -1.000000000 +  0.000000000i   1.000000000 * CIS(  3.141592654)
11   -1.000000000 -  1.000000000i   1.414213562 * CIS(  3.926990817)
12    0.000000000 -  1.000000000i   1.000000000 * CIS( -1.570796327)
13    1.000000000 -  1.000000000i   1.414213562 * CIS( -0.785398163)
14    5.000000000 +  0.000000000i   5.000000000 * CIS(  0.000000000)
15    5.000000000 +  3.000000000i   5.830951895 * CIS(  0.540419500)
16    0.000000000 +  3.000000000i   3.000000000 * CIS(  1.570796327)
17   -5.000000000 +  3.000000000i   5.830951895 * CIS(  2.601173153)
18   -5.000000000 +  0.000000000i   5.000000000 * CIS(  3.141592654)
19   -5.000000000 -  3.000000000i   5.830951895 * CIS( -2.601173153)
20    0.000000000 -  3.000000000i   3.000000000 * CIS( -1.570796327)
21   -5.000000000 -  3.000000000i   5.830951895 * CIS( -2.601173153)

A dynamic array of TComplex values works just like any other Delphi (D4 or later) dynamic array:

CSet with 2D dynamic array of TComplex values

VAR  z:  ARRAY OF ARRAY OF TComplex;
. . .
BEGIN
  SetLength(z, ImageComposite.Height, ImageComposite.Width);
  FOR j := Low(z) TO High(z) DO
  BEGIN
    y := 3 - 6*(j-Low(z))/(High(z) - Low(z));
    FOR i := Low(z[j]) TO High(z[j]) DO
    BEGIN
      x := -4 + 8*(i-Low(z[j]))/(High(z[j]) - Low(z[j]));
      z[j,i] := CSet(x,y, cfRectangular);
    END
  END;
. . .

Converting rectangular coordinates to polar coordinates

Conversion from (x, y) rectangular form to polar coordinates, (r, q), is fairly simple:

To provide a range of q  to be in the interval [-p, p), the arctan2(y, x) function is used from the Delphi math library.

Converting polar coordinates to rectangular coordinates

The opposite conversion, from (r, q) polar coordinates back to the (x, y) rectangular coordinates, is also simple:

The CConvert function can be used to convert from one format to the other:

CConvert function

// complex number definition/conversion/output
FUNCTION CConvert (CONST z:  TComplex; f:  TComplexForm):  TComplex;
  VAR
    a:  TComplex;
BEGIN
  IF   IsInfinity(z.x) OR IsInfinity(z.y)
  THEN BEGIN
    CASE f OF
      cfPolar:        RESULT := CSet(Infinity, 0.0, cfPolar);
      cfRectangular:  RESULT := CSet(Infinity, Infinity)
    END
  END
  ELSE BEGIN
    IF   IsNan(z.x) OR IsNan(z.y)
    THEN RESULT := CSet(NaN, NaN)
    ELSE BEGIN
 
      IF   z.form = f
      THEN RESULT := z
      ELSE BEGIN
        CASE z.form OF
          cfPolar:            // cfPolar-to-cfRectangular conversion
            BEGIN
              a.form := cfRectangular;
              a.x := z.r * COS(z.theta);
              a.y := z.r * SIN(z.theta)
            END;
 
          cfRectangular:      // cfRectangular-to-cfPolar conversion
            BEGIN
              a.form := cfPolar;
              a.r := Cabs(z);
              a.theta := ARCTAN2(z.y, z.x);
            END;
        END;
        RESULT := a
      END
    END
  END
END {CConvert};

Using CConvert, any TComplex value can be printed either rectangular or polar format.  The default format for either form is specified in the FormatString in the function header:

CToPolarStr prototype with default string '%*.*f * CIS(%*.*f)'

FUNCTION  CToPolarStr (CONST z:  TComplex; CONST w,d:  BYTE;
                       CONST FormatString:  STRING = 
                       '%*.*f * CIS(%*.*f)'):  STRING;

CToPolar function

FUNCTION  CToPolarStr (CONST z:  TComplex;
                       CONST w,d:  BYTE;
                       CONST FormatString:  STRING):  STRING;
  VAR
    u :  TComplex;
BEGIN
  u := CConvert (z, cfPolar);
  RESULT := Format(FormatString, [w,d,u.r, w,d,u.theta]);
END {CToPolarStr};

CToRectStr prototype with default string '%*.*f +%*.*fi'

FUNCTION  CToRectStr  (CONST z:  TComplex; CONST w,d:  BYTE;
                       CONST FormatString:  STRING = 
                       '%*.*f +%*.*fi'):  STRING;

CToRectStr function

FUNCTION  CToRectStr (CONST z:  TComplex;
                      CONST w,d:  BYTE;
                      CONST FormatString:  STRING):  STRING;
  VAR
    NegativeFormat:  STRING;
    u             :  TComplex;
BEGIN
  IF   IsNan(z.x) OR IsNan(z.y)
  THEN RESULT := Format(FormatString, [w,d,NaN, w,d,NaN])
  ELSE BEGIN
    u := CConvert (z, cfRectangular);
    IF   u.y >= 0
    THEN RESULT := Format(FormatString, [w,d,u.x, w,d,u.y])
    ELSE BEGIN
      NegativeFormat := FormatString;
      IF   POS('+', NegativeFormat) > 0
      THEN NegativeFormat[ POS('+',NegativeFormat ) ] := '-';
      RESULT := Format(NegativeFormat, [w,d,u.x, w,d,ABS(u.y)])
    END
  END
END {CToRectStr};

Some extra work is made in CToRectStr so that two consecutive negative signs do not appear.

Euler’s formula and exponential form

Euler’s formula states:

One way to demonstrate the correctness of this formula is by using a series expansion of eiq, which breaks into real and imaginary series that are the cos and sin functions.

Only values on a unit circle can be described directly with Euler’s formula. Multiplication by some distance value, r, allows any complex value, z, to be described in this exponential form.  All of the following forms are equivalent:

Polar coordinates, (r, q), can be expressed in either a trigonometric or exponential form.

Download ComplexMathLibrary

Updated 18 Feb 2002


since 24 June 2001