Complex numbers in APL

Alan G. Hawkes, Department of Management Science and Statistics, University College of Swansea, U.K.

mathematics and statistics programming techniques

Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 107

This paper describes a systematic method of representing complex numbers in standard APL in such a way that basic complex arithmetic and complex matrix algebra may be simply done. A listing of functions is given in section five.

The VECTOR typesetters were persuaded to resist the challenge of setting this article in print, and consequently it is reproduced directly from Alan’s originals.

1. Introduction: Representation of complex numbers¶

One of the strengths of FORTRAN as a language for scientific and engineering computing is that it makes provision for declaring COMPLEX variables which may then be manipulated according to the rules of complex arithmetic by the standard algebraic symbols +, −, *, /. In addition, special functions or subroutines such as CEXP are available: in particular complex matrix algebra can be done by calling the appropriate subroutines.

Standard APL makes no provision for complex variables. After several proposals and discussions, Penfield (1979) produced proposals to extend the APL language into the complex field in a consistent manner. McDonnell (1981) describes such an extension which has been implemented in I.P. Sharp’s APL. I would urge all writers of APL interpreters to follow suit as soon as possible; the scientists would then be more readily moved to change to APL.

In the meantime, it is quite easy to develop a simple means of handling complex numbers in standard APL. It could no doubt be done more elegantly in versions supporting nested arrays. I have not attempted to be exhaustive, but have restricted myself to the basis found necessary in my own applications. The functions mentioned in the text are listed in section five. They have been knocked together quite rapidly to solve a particular problem. The author would welcome amendments which would render them move elegant, efficient or accurate.

A single complex number \(z = x + iy\) has real part \(x\) and imaginary part \(y\) (\(x\) and \(y\) both real), so that it may be thought of as an ordered pair of real numbers. We adopt the convention that a complex array of shape S will be represented by a real array of shape 2,S such that Z← X, [0.5]Y where \(X\) is the real part and \(Y\) the imaginary part (each of shape S).

Thus, for example A←7 8 is a complex scalar \(7+i8\) with

      REAL A
7
      IMAG A
8

are real scalars corresponding to the real and imaginary parts of A (here S is the empty vector). In contrast

      ⎕←B←2 1⍴7 8
7
8

is a complex vector with one element. Superficially

      REAL B
7
      IMAG B
8

appear to be the same as the previous case but this time 7 and 8 are one-element vectors (⍴REAL B and ⍴IMAG B are both 1).

More generally, the complex numbers \(1+i5\), \(2+i6\), \(3+i7\), \(4+i8\) may be assigned to an array of shape 2 4; so that

      C
1 2 3 4
5 6 7 8

is a four-element complex vector with

      REAL C
1 2 3 4
      IMAG C
5 6 7 8

An example of a complex 3 by 4 matrix, whose representation is an array of shape 2 3 4 is

      D
 1  2  3  4
 5  6  7  8
 9 10 11 12

13 14 15 16
17 18 19 20
21 22 23 24

with real and imaginary parts

      REAL D
 1  2  3  4
 5  6  7  8
 9 10 11 12

      IMAG D
13 14 15 16
17 18 19 20
21 22 23 24

Higher-dimensional arrays could be defined, but do not seem likely to be much used.

An alternative polar representation of a complex number

\[\begin{gather*} x + iy = re^{i\theta} = r\cos\theta + i\sin\theta \\ r = (x^2 + y^2)^{\frac{1}{2}} \end{gather*}\]

is often useful. The modulus \(r\) and argument \(\theta\) are obtained by functions MOD Z and ARG Z, while R CART THETA does the inverse transformation from polar to Cartesian co-ordinates. In this connection the functions RAD and DEG are useful to convert between angles measured in degrees and radians. CONJ Z yields the complete conjugate \(X - iY\), while COMPLEX X converts a real array into a complex array with zero imaginary part, so that for example

      COMPLEX 3
3 0

is a complex scalar and

      COMPLEX 7 5
7 5
0 0

is a complex vector of length two.

2. Complex arithmetic¶

Two complex arrays of the same shape may be added or subtracted by Z + W or Z − W in the usual say. However, this does not satisfy the rule that if one of the arguments is a (complex) scalar it is repeated to become the same shape as the other argument and thus added to, or subtracted from, each element of the other argument. This property is a feature of the functions CADD and CSUB so, for example, using C from the previous section we have

      C CSUB 1 2
0 1 2 3
3 4 5 6

and we see that the complex scalar \(1 + i2\) is subtracted from each element of the complex vector C. The same is true for elementwise multiplication and division in the functions Z CTIMES W and Z CDIV W respectively. In each case Z and W must be the same shape or be complex scalars. As in the real case, (complex) vectors of length 1 and 1 by 1 matrices may act like scalars in this respect.

Simple monadic functions are CRECIP Z, the complex version of monadic ÷, CEXP Z and CLOG Z.

In principle it would be possible to write functions to manipulate complex arrays, such as CTAKE or CDROP etc., but one can always achieve the same ends by appropriate real manipulations.

3. Complex matrix algebra¶

The addition and subtraction of matrices whose elements are complex numbers is just a special case of the addition and subtraction of general arrays, as discussed in the previous section. Matrix multiplication is done by Z CMM W, which is the complex equivalent of Z+.×W, having equivalent rules about conformability and (complex) shape of the result.

This leaves three rather more difficult problems. The most important is that of finding the inverse of a complex matrix. This is an interesting problem which I discuss in this section. In the following section I deal with the computation of the exponential of a matrix. The third problem is that of finding eigenvalues and eigenvectors of a complex matrix. I have not yet attempted this because it is the most difficult and because I have as yet not had any need to do so. However, it should be done and I would be pleased to hear from anyone who may already have done it.

The inverse of a complex matrix¶

In principle one could write a function which carries out some form of Gaussian elimination to find the inverse of a complex matrix. However, since APL is very good at inverting real matrices with ⌹ it seems a good idea to make use of that facility.

Denote the inverse of the complex square matrix \(Z = X + iY\) by \(U + iV\). Then, since their product must be the identity \(I\), we have

\[\begin{align*} XU - YV &= I \tag{3.1} \\ XV + YU &= 0 \tag{3.2} \end{align*}\]

Suppose there exist real \(\lambda_1, \lambda_2\) such that

\[\begin{alignat*}{2} B &= \lambda_1X - \lambda_2Y &&\quad\text{is invertible. Let} \\ C &= \lambda_2X + \lambda_1Y\;. && \end{alignat*}\]

Then \(\lambda_1 \times (3.1) - \lambda_2 \times (3.2)\) and \(\lambda_2 \times (3.1) + \lambda_1 \times (3.2)\) yield

\[\begin{align*} BU - CV &= \lambda_1I \\ CU + BV &= \lambda_2I, \end{align*}\]

which can be solved to yield

\[\begin{align*} U &= (B + CB^{-1}C)^{-1}\,(\lambda_1I + \lambda_2CB^{-1}) \tag{3.3} \\ V &= (B + CB^{-1}C)^{-1}\,(\lambda_2I - \lambda_1CB^{-1}) \tag{3.4} \end{align*}\]

If Z is non-singular there are at most 1↑⍴X values of the ratio \(\lambda_1/\lambda_2\) for which \(B\) is singular. Our general routine CINV Z therefore chooses \(\lambda_1,\lambda_2\) at random. There will be a very small probability of hitting one of these ‘eigenratios’ . However, to be on the safe side, an error trap is set in line [7] which returns to choose another pair of values for \(\lambda_1,\lambda_2\) if \(B\) should be singular. Of course, if Z is singular, \(B\) will always be singular and we would get into an infinite loop. Thus we set a counter and switch off the error trap in line [10] after five attempts to find a non-singular \(B\). Once a suitable \(B\) is found, the inverse \(U + iV\) is constructed from (3.3) and (3.4) using only real matrix algebra. Error-trapping routines vary from one interpreter to another, so anyone using this function may need to make appropriate minor adjustments.

For particular problems one sometimes knows that the real part \(X\) of a matrix must be invertible. In that case one can take \(\lambda_1 = 1\), \(\lambda_2 = 0\) which imply that

\[\begin{alignat*}{2} &B = X\;,\ C = Y \text{ and thus} && \\ &U = (X + YX^{-1}Y)^{-1} && \tag{3.5} \\ &V = -(X + YX^{-1}Y)^{-1}\,YX^{-1} && \tag{3.6} \end{alignat*}\]

This solution is implemented in the function CINVR Z. Similarly, if it is known that \(Y\) is non-singular, one may take \(\lambda_1 = 0\), \(\lambda_2 = -1\) so that \(B = Y\), \(C = -X\) which imply that

\[\begin{alignat*}{2} &U = (Y + XY^{-1}X)^{-1}\,XY^{-1} && \tag{3.7} \\ &V = -(Y + XY^{-1}X)^{-1} && \tag{3.8} \end{alignat*}\]

This solution is implemented in the function CINVI Z. A simple example in which neither the real nor imaginary parts are invertible is

\[\begin{equation*} \begin{pmatrix} 1 & 1 \\ 0 & i \end{pmatrix}^{-1} = \begin{pmatrix} 1 & i \\ 0 & -i \end{pmatrix} \;. \end{equation*}\]

Any values of \(\lambda_1,\lambda_2 \ne 0\) would suffice to invert this one. The general routine CINV does the job.

      CINV   2 2 2 ⍴ 1 1 0 0 0 0 0 1
1  0
0  0

0  1
0 ¯1

The matrix \(\begin{pmatrix} 1 + i & 2 \\ 1 & 1-i \end{pmatrix}\) ,however, is singular so CINV 2 2 2 ⍴ 1 2 1 1 1 0 0 ¯1 yields the dreaded DOMAIN ERROR.

Note that any of these three inverse functions are generalisations of the real monadic ⌹ , except that we allow them to apply only to complex square matrices, vectors of length 1, or scalars. In the real case, of course, one gets a least squares fit for the non-square case. I have not attempted any generalisation of this feature.

4. The exponential of a matrix¶

In applied probability, and other areas of applied mathematics, one often needs to evaluate the exponential of a square matrix, defined by the infinite series

\[\begin{equation*} e^A = I + \frac{A^2}{2!} + \frac{A^3}{3!} + \ldots\;. \end{equation*}\]

The accurate and efficient computation of this is a rather delicate matter, see Moler and Van Loan (1978). I have adopted the simplest, and probably most effective, method of scaling and squaring. Consider the case of a real matrix \(A\), which is implemented in the function MEXP A.

The most obvious thing to do is to truncate the series

\[\begin{equation*} e^A \simeq T_k(A) = \sum_{r=0}^{k} A^r \div r! \quad. \tag{4.1} \end{equation*}\]

Bounds for the error of truncation are given by

\[\begin{equation*} \lVert T_k(A) - e^A \rVert \le \frac{\lVert A\rVert^{k+1}}{k+1!}\,\bigl[1 - \lVert A\rVert/(k+2)\bigr]^{-1} \;, \tag{4.2} \end{equation*}\]

where \(\lVert A\rVert\) is the norm of the matrix. A number of different norms may be used. For simplicity of calculation use the 1-norm

\[\begin{equation*} \lVert A\rVert = \operatorname*{Max}_i \sum_j \lvert a_{ij}\rvert \quad. \tag{4.3} \end{equation*}\]

The scaling and squaring method, as implemented here, proceeds in three stages:-

(i): find the smallest non-negative integer \(M\) such that

\[\begin{equation*} \lVert A \div 2^M\rVert = \lVert A\rVert \div 2^M < \tfrac{1}{2}\;, \end{equation*}\]
(ii): Approximate \(e^{(A \div 2^M)}\) by the finite series (4.1), replacing \(A\) by \(A \div 2^M\) and with \(k = 13\) (truncation error less than \(10^{-15}\) ),
(iii): find \(e^A = \bigl(e^{A \div 2^M}\bigr)^{2^M}\) by squaring \(M\) times.

For a complex matrix one simply needs to interpret \(\lvert a_{ij}\rvert\) in (4.3) by the complex modulus, MOD, and replace +.× by CMM when performing matrix multiplication. This is done in the function CMEXP. For MEXP and CMEXP the inputs must be real, or complex as appropriate, square matrices, vectors of length 1 or scalars. In each case, the output has the same shape as the input.

Sometimes there are short cuts. For example, if \(I\) is a real identity matrix of the same shape as a real matrix \(Q\) and \(s\) a complex scalar, then

\[\begin{equation*} e^{-sI + Q} = e^{-s}\,e^Q \end{equation*}\]

so that the exponential of this particular complex matrix can be written as the product of the exponential of a complex scalar and the exponential of a real matrix. In terms of the functions described in this paper this is equivalent to the rather less elegant form CMEXP (COMPLEX Q) CSUB S CTIMES COMPLEX I is equal to (CEXP −S) CTIMES COMPLEX MEXP Q.

Oh well, at least it’s better than doing it in FORTRAN!

Finally, a word of warning. The obvious result

\[\begin{equation*} e^{A + B} = e^A\,e^B \end{equation*}\]

is only true for matrices if \(A\) and \(B\) commute.

References¶

FORBES, D. (1981), Complex floor revisited, APL81, APLQQ vol. 12 No. 1 September 1981, 107-111.

McDONNELL (1981), An implementation of complex APL, APLQQ vol. 11 No. 3, March, 1981, 19-22.

MOLER, C. and VAN LOAN, C. (1978), Ninteen dubious ways to compute the exponential of a matrix, SIAM Review, vol. 20 No. 4, 801-836.

PENFIELD, P. (1979), Proposal for a complex APL, APL79, APLQQ vol. 9 No. 4, June 1979, 47-53.

PENFIELD, P. (1981), Principal values and branch cuts in complex APL, APL81, APLQQ, vol. 12 No. 1, September 1981, 248-256.

Appendix: Function Listings¶

      ∇R←ARG X;I
[1]   ⍝ ARGUMENT OF COMPLEX X, -PI<THETA≤PI
[2]    R←REAL X
[3]    I←IMAG X
[4]    R←(2×((I≥0)-0.5))×¯2○¯1⌈1⌊R÷((R*2)+(I*2))*0.5
      ∇

      ∇R←X CADD Y
[1]   ⍝ ADDITION OF COMPLEX X AND Y
[2]    R←((REAL X)+REAL Y),[0.5](IMAG X)+IMAG Y
      ∇

      ∇Z←R CART THETA
[1]   ⍝ POLAR TO CARTESIAN(COMLEX); THETA IN RADIANS
[2]    Z←(R×2○THETA),[0.5]R×1○THETA
      ∇

      ∇R←X CDIV Y
[1]   ⍝ DIVISION OF COMPLEX X BY Y
[2]    R←X CTIMES CRECIP Y
      ∇

      ∇R←X CTIMES Y;RX;IX;RY;IY
[1]   ⍝ MULTIPLIES ELEMENTWISE; X AND Y MUST BE SAME SHAPE OR COMPLEX SCALAR
[2]    RX←REAL X
[3]    IX←IMAG X
[4]    RY←REAL Y
[5]    IY←IMAG Y
[6]    R←((RX×RY)-(IX×IY)),[0.5](RX×IY)+(IX×RY)
      ∇

      ∇R←CRECIP X;RX;IX;MODSQ
[1]   ⍝ COMPLEX RECIPROCAL OF ARRAY X
[2]    RX←REAL X
[3]    IX←IMAG X
[4]    R←(RX÷MODSQ),[0.5]-IX÷MODSQ←(RX*2)+(IX*2)
      ∇

      ∇R←CEXP Z
[1]   ⍝ EXPONENTIAL OF COMPLEX Z
[2]    R←(*REAL Z)CART IMAG Z
      ∇

      ∇R←CINV Z;X;Y;B;C;I;U;LAMDA;D;S;J;S∆CINV
[1]   ⍝ INVERSE OF COMPLEX MATRIX WITH RANDOM TRANSFORM: AGH 1984
[2]    →(2≠1↑S←⍴Z)/MESS
[3]    →((1=⍴S)∨(2=⍴S)∧1=¯1↑S)/SCALAR
[4]    →((3=⍴S)∧0=-/1↓S)/CONTINUE
[5]   MESS:→0,⍴⎕←'ARGUMENT MUST BE COMPLEX SQUARE MATRIX'
[6]   SCALAR:Z←2 1 1⍴Z
[7]   CONTINUE:S∆CINV←1 2⍴COUNT,1
[8]    J←1
[9]   COUNT:→(5≥J←J+1)/RAND
[10]   S∆CINV←0 2⍴1
[11]  RAND:LAMDA←LAMDA÷+/LAMDA←?2⍴1000000000
[12]   B←(LAMDA[1]×X←REAL Z)-LAMDA[2]×Y←IMAG Z
[13]   C←(LAMDA[1]×Y)+LAMDA[2]×X
[14]   I←I∘.=I←⍳1↑⍴X
[15]   D←⌹B+(R←C+.×⌹B)+.×C
[16]   U←D+.×(LAMDA[1]×I)+LAMDA[2]×R
[17]   R←S⍴U,[0.5]D+.×(LAMDA[2]×I)-LAMDA[1]×R
      ∇

      ∇R←CINVI Z;Y;D;X;S
[1]   ⍝ INVERSE OF COMPLEX Z WHEN IMAG PART INVERTIBLE
[2]    →(2≠1↑S←⍴Z)/MESS
[3]    →((1=⍴S)∨(2=⍴S)∧1=¯1↑S)/SCALAR
[4]    →((3=⍴S)∧0=-/1↓S)/CONTINUE
[5]   MESS:→0,⎕←'ARGUMENT MUST BE COMPLEX SQUARE MATRIX'
[6]   SCALAR:Z←2 1 1⍴Z
[7]   CONTINUE:D←⌹Y+(R←X+.×⌹Y←IMAG Z)+.×X←REAL Z
[8]    R←S⍴(D+.×R),[0.5]-D
      ∇

      ∇R←CINVR Z;Y;D;X;S
[1]   ⍝ INVERSE OF COMPLEX Z WHEN REAL PART INVERTIBLE
[2]    →(2≠1↑S←⍴Z)/MESS
[3]    →((1=⍴S)∨(2=⍴S)∧1=¯1↑S)/SCALAR
[4]    →((3=⍴S)∧0=-/1↓S)/CONTINUE
[5]   MESS:→0,⎕←'ARGUMENT MUST BE COMPLEX SQUARE MATRIX'
[6]   SCALAR:Z←2 1 1⍴Z
[7]   CONTINUE:D←⌹X+(R←Y+.×⌹X←REAL Z)+.×Y←IMAG Z
[8]    R←S⍴D,[0.5]D+.×-R
      ∇

      ∇R←CLOG Z
[1]   ⍝ LOGARITHM OF COMPLEX Z
[2]    R←(⍟MOD Z),[0.5]ARG Z
      ∇

      ∇R←CMEXP A;J;T;M;S
[1]   ⍝ EXPONENTIAL OF COMPLEX MATRIX: AGH84
[2]    →(2≠1↑S←⍴A)/MESS
[3]    →((1=⍴S)∨(2=⍴S)∧1=¯1↑S)/SCALAR
[4]    →((3=⍴⍴A)∧0=-/1↓⍴A)/CONTINUE
[5]   MESS:→0,⍴⎕←'ARGUMENT MUST BE COMPLEX SQUARE MATRIX'
[6]   SCALAR:A←2 1 1⍴A
[7]   CONTINUE:A←A÷2*M←0⌈⌈1+2⍟1E¯10+⌈/+/MOD A
[8]    R←T←COMPLEX T∘.=T←⍳(J←1)↑1↓⍴A
[9]   SERIES:R←R+T←T CMM A÷J
[10]   →(13≥J←J+1)/SERIES
[11]   J←0
[12]  SQUARE:→(M<J←J+1)/EXIT
[13]   R←R CMM R
[14]   →SQUARE
[15]  EXIT:R←S⍴R
      ∇

      ∇R←COMPLEX X
[1]   ⍝ CONVERTS REAL X TO COMPLEX
[2]    R←X,[0.5](⍴X)⍴0
      ∇

      ∇R←X CMM Y;RX;RY;IX;IY
[1]   ⍝ COMPLEX MATRIX MULTIPLICATION
[2]    RX←REAL X
[3]    IX←IMAG X
[4]    RY←REAL Y
[5]    IY←IMAG Y
[6]    R←((RX+.×RY)-(IX+.×IY)),[0.5](RX+.×IY)+(IX+.×RY)
      ∇

      ∇R←CONJ X
[1]   ⍝ CONJUGATE OF COMPLEX X
[2]    R←(REAL X),[0.5]-IMAG X
      ∇

      ∇R←X CSUB Y
[1]   ⍝ SUBTRACTION OF COMPLEX Y FROM X
[2]    R←((REAL X)-REAL Y),[0.5](IMAG X)-IMAG Y
      ∇

      ∇R←DEG X
[1]   ⍝ CONVERTS X IN RADIANS TO DEGREES
[2]    R←X×180÷○1
      ∇

      ∇R←IMAG X;S
[1]   ⍝ EXTRACTS IMAG PART OF COMPLEX X
[2]    →(2=1↑⍴X)/CONTINUE
[3]    →0,⍴⎕←'COMPLEX X MUST HAVE 2=1↑⍴X'
[4]   CONTINUE:R←S⍴(¯1,S←1↓⍴X)↑X
      ∇

      ∇R←MOD X
[1]   ⍝ MODULUS OF COMPLEX X
[2]    R←(((REAL X)*2)+((IMAG X)*2))*0.5
      ∇

      ∇R←MEXP A;J;T;M;S
[1]   ⍝ EXPONENTIAL OF REAL MATRIX: AGH84
[2]    S←⍴A
[3]    →(1=⍴,A)/SCALAR
[4]    →((2=⍴⍴A)∧(0=-/⍴A))/CONTINUE
[5]    →0,⍴⎕←'ARGUMENT MUST BE REAL SQUARE MATRIX'
[6]   SCALAR:A←,A
[7]   CONTINUE:A←A÷2*M←0⌈⌈1+2⍟1E¯10+⌈/+/|A
[8]    R←T←T∘.=T←⍳(J←1)↑⍴A
[9]   SERIES:R←R+T←T+.×A÷J
[10]   →(13≥J←J+1)/SERIES
[11]   J←0
[12]  SQUARE:→(M<J←J+1)/EXIT
[13]   R←R+.×R
[14]   →SQUARE
[15]  EXIT:R←S⍴R
      ∇

      ∇R←REAL X;S
[1]   ⍝ EXTRACTS REAL PART OF COMPLEX X
[2]    →(2=1↑⍴X)/CONTINUE
[3]    →0,⍴⎕←'COMPLEX X MUST HAVE 2=1↑⍴X'
[4]   CONTINUE:R←S⍴(1,S←1↓⍴X)↑X
      ∇

      ∇R←RAD X
[1]   ⍝ CONVERTS X IN DEGREES TO RADIANS
[2]    R←○X÷180
      ∇