310 IEEE TRANSACTIONS ON COMPUTERS, MARCH 1971

Short Notes_

The Relationship Between Two Fast the z = tlt2 ... tn elements of each as sequences or "strings."
Fourier Transforms When all the tj's are equal to t, we can refer to the n-dimen-
J. GOOD
sional mod t DFT. The components of a multidimensional
mod 2 DFT (of an array of real number ar) are all real be-
I.

Abstract-The purpose of this note is to show as clearly as possible cause in this case all the oj's are equal to -1. The mod 2
the mathematical relationship between the two basic fast methods
used for the calculation of discrete Fourier transforms and to general-
multidimensional DFT is also an example of a Hadamard
ize one of the methods a little further. This method applies to all those transform [12]. The matrix of a Hadamard transform is
linear transformations whose matrices are expressible as direct orthogonal and each of its elements is + 1 or - 1. It is also
products. known as an "anallagmatic pavement" [13].
Index Terms-Algorithms, circulices, direct product of matrices, A DFT, unidimensional or multidimensional, is, ofcourse,
discrete Fourier transforms, fast Fourier transforms, frequency a linear transformation of a vector (with t or tJt2 * tn com-
analysis, Hadamard transform, harmonic analysis, multidimensional
linear transformation, sparse matrices.
ponents), a multiplication of a vector by a matrix, in which
the matrix of the transformation is (Cw"rs) or (0jjrSls* * nSn)s
A fast Fourier transform (FFT) is a fast method for the In the latter matrix some convention is required for the
calculation of a discrete Fourier transform (DFT) whose ordering of the rows and columns. The simplest convention
definition we shall give below in order to make the paper is to order the t1t2 ... tn vectors r =(r1,.*. , re)'lexicographi-
adequately self-contained. For applications of FFTs and cally as if they were words in a dictionary, and similarly
DFTs see, for example, [1 ]-[4]. See also the bibliographies for the "columns" s = (s , , sj). Then the components of
[5], [6], [4]. Here we shall be concerned only with mathe- ar must be similarly ordered in order to convert it from an
matical and computational aspects. The main purpose of array into a vector.
this note is to show as clearly as possible the relationship The first FFT to be described here depends on the fact
between the two basic algorithms for FFTs [7], [8]. Qther that a one-dimensional DFT modulo T, where T = t1t2 ... tn
methods are refinements of these two. Each method has its and tl, t2, * , tn are mutually prime in pairs, can be expressed
own advantages; one method can be applied to a wider as a mod (t , t2, ,tn) multidimensional DFT. For the sake
variety of moduli, the other is more appropriate for a wide of completeness the procedure is quoted from Good [7].1
class of multidimensional transforms, not necessarily We first set up (1-1) correspondences between the scalar
Fourier. r and the vector r and also between s and s in the following
A one-dimensional mod t DFT of a vector (a sequence of manner.
t numbers ao, a,, a*,a 1 assumed to be real in this paper) is We begin by recalling the Chinese Remainder Theorem,
another vector of t complex numbers a*, a a*, de- known by'Sun-Tsu in the first century A.D. [14], [15]. It
fined by the equations states that if an integer is known modulo n integers (that is,
t- 1 its remainder is known when it is divided by each of these n
a= L acor (s=0,1,---,t- 1) integers), these integers being mutually pfime in pairs (that
S=O is, each pair has highest common factor equal to 1), then it
where o = exp (2ri/t), a tth root of unity. It has elegant can be uniquely determined modulo the product. The solu-
analogs of the properties of an ordinary Fourier transform tion can be expressed succinctly [9], [7] as follows. If
like an inversion formula, a formula for the transform of a t,,, t, are mutually prime in pairs and X = * t., and if
convolution, and even a Poisson summation formula (for the remainders are sl,,, Sn, or in the standard "congru-
example, [9]-[11]). A multidimensional mod (tl, t2, , t,) ence" notation, if
DFT of the n-dimensional array of t,t2 tn numbers ar, s sl(mod tl) *, s sn(mod tn),
where r is the vector (r,, r2, , rn) (rj=O, 1, -*, tj- 1; then
j= 1, 2,5 ..., n), is the n-dimensional array of complex num-
bers a* defined by the equations S Kr) + +
a* = E a .r.s. rnsn W
t t1Jt
1 tn Vtn
r (mod T;0 < s < ;0 < sl < t,,etc.) (3)
(rj,sj=0,l, 1 ,t -l;j= 1,5,n) (2)
where (T/tl)-', for example, means the modulo t, reciprocal
where coj =exp (27ri/tj) (j= 1, *, n). The arrays ar and a* of the integer T/t,. This reciprocal exists and is unique
can also be regarded as vectors if a rule is given for ordering modulo t1 because the integers T/t1 and t, are mutually

Manuscript received December 18, 1969; revised August 21, 1970. 'There are misprints in [7 ]. Four incorrect entries in the matrix at the
The author is with the Department of Statistics and Statistical Labora- foot of p. 363 have been found. On p. 364, line 10, 6 should be 'rV+.
tory, Virginia Polytechnic Institute and State University, Blacksburg, Va. The lower limit of the product at the foot of p. 365 should be 1. On p. 368,
24061. line 13, an asterisk was omitted.
SHORT NOTES 311

prime. Hence the first term of this expression for s is unique where we have replaced co" by K 1(r , s 1), etc., and where we
modulo T, and similarly the other terms are unique. The do not insist for the present on any constraints on t1, * t,.
reader can verify at once that the value of s does satisfy (The following discussion of (5) and (6) was given in [4].)
the n given congruences. Moreover, each of the T possible We shall describe (5) as an n-dimensional K-transform. Al-
sets of n congruences has at least one solution given by (3). though each of the functions K1, K2, ... Kn can be regarded
These solutions must be unique, since there are only - as a square matrix, the right side of (5) is, of course, not as it
residues modulo z and none of them can correspond to more stands in the form of the product of a vector by several
than one set of n congruences. Formula (3) is analogous matrices, although it can be regarded as the product of a
to Lagrange's interpolation formula. We use it to set up a vector by one large matrix whose rows are indexed by r and
(1-1) correspondence between s and s = s sn) which I columns by s. This large matrix will later be seen to be the
shall call the Sino correspondence or the Sino representation "direct product" of the smaller matrices, and also to be the
of s. For example, if t1 =3, t2 = 8, T =24, then s-16s1 + 9s2 product of large matrices which correspond in turn to K1,
(mod 24). Once we fix t1, *, tn the coefficients (here 16 and K2, -, K,. For the moment we shall point out a simple
9) are known once for all. On the other hand, the (1-1) procedure for summing (5) without matrix notation.
relationship between r and r is given by The summation can be performed by first summing over
T ~~~T
rl, then over r2, etc. The first summation is
r- r, + + rn
tI ~~~tn Eri arKI(ri, s1)
(modz0 < r < z;0 < r1 < t1,etc.) (4)
and the result is of the form L1(r2, r3, , rn, s1), that is, a
which I shall call the Ruritanian correspondence. Again, if function of r2, r3,.**, rn, s1. The next stage involves the sum
t, = 3, t2 = 8, T = 24, this gives r_ 8r1 + 3r2 (mod 24). Z L1(r2, r3,**, r,sn)K2(r2, S2)
I now assert that the equations r2

and is of the form L2(r3, * , rn, sI' 52), and so on. Thus at
T- 1
A* O=
E A rs (oj = e27ri/T)
r=0 each stage we need store only an n-dimensional array, and
and altogether we need do only (t1 + t2 + . + tn) scalar multi-
plications. This is the basis of the method of [7] for forming
tl-i tn-i
a multidimensional DFT, except that it was there expressed
a* = E. a Coijisi . .. nn
r, rn=O
in matrix notation. The method saves arithmetic as com-
pared with the multiplication of the vector by the single
O

are equivalent, where ar= ar and a* = a* by definition, and large matrix mentioned above which would require T2 scalar
where the correspondences are Sino and Ruritanian, re- multiplications. For example, with T = 5.7.8.9 = 2520, the
spectively. To see this, note that number of scalar multiplications is less than an eightieth
of what would be required (X2 = 25202) if aC were multi-
r =-(T/tv)rv (mod tv), s =_ {) (mod tQ) plied by the X x T matrix K1(rl, s1) . . . Kn(rn, sn).
tv tv tv We shall later discuss the inverse of the multidimensional
K-transform, but first we consider a generalization of the
so that rs- rVsVT/tv (mod tv) and hence, by Sun-Tsu's above simple summation procedure which leads up to the
theorem in the form given above, second method for calculating an FFT.
+
We consider the analogous multiple sum
Trlsl T (T
rs + Er s Tc/tv. tl-1 tn- I
Z E arKi(rl, s1)K2(r2, S1, S2)
Z
Moreover, r=O rn=O
* Kn(rn, S1, S21 . .
*, Sn) (6)
v

so
and we find that the above argument goes through virtually
unchanged for this more general summation. The sum-
a*
= a r, 's ....Jf,nSn = mation of arK, with respect to r1 leads to an expression of
r
the form L1(r2, .*. ,r , s1); then the summation ofL1K2 with
asclaimed. respect to r2 leads to an expression of the form L2(r3, * , rn,
A unidimensional mod t1t2 tn DFT can therefore be S1, S2), * , and finally the summation of Ln_ 1K, with respect
computed with about the same speed as a multidimensional to rn to an expression of the form Ln(Sl, S2,* *' sn). At each
mod (t1, ., tn) DFT where tl, * *, tn are mutually prime in stage only an n-dimensional array needs to be stored, but the
pairs. Accordingly we now consider the calculation of a dimensions of the array vary from stage to stage unless the
multidimensional DFT, or more generally, a linear trans- tv's are all equal.
formation from a. to a* of the form The method used by Cooley and Tukey [8 ] for computing
tl-1 tn- 1
a one-dimensional DFT is in effect to reduce it to the form
a*= arK1(rl,sl) .*. Kn(r ,sn) (5) (6). Their method applies to a modulus - = t1t2 ... t, 1, but
rj =O r,=O it was first spelled out for the case X = 2n, that is t1 = t2 =
 312 IEEE TRANSACTIONS ON COMPUTERS, MARCH 1971

= tn=2. The general case was explained in more detail the same function. Some such contexts are mentioned, for
by Bingham, Godfrey, and Tukey [16]. In this method the example, by Andrews and Caspari [17]. Also, even for the
correspondence between the one-dimensional DFT and a DFT itself, if very fast special-purpose equipment is to be
multiple sum depends on setting up (11) correspondences built, using fixed values of t1, t2, , t, the Sino-Ruritanian
between r and r and between s and s different from the representation of r and s might turn out to be more con-
Sino-Ruritanian correspondences or representations. In- venient than the mixed-radical representation.
stead r and s are each represented in mixed radix form, but In order to discuss the algorithms in matrix notation it is
each with a different kind of reversal, as follows. necessary to make use of direct products (=Kronecker
r =t1t2 tn-lr1 + t12 tn-t2r2 + ' + tlrn-1 + rn products) of matrices. (See, for example, [18]). The direct
S = t2t3 tnSn + t3t4 *.. tnSn-1 +
...
+ tns2 + Sl. (7)
product B x C of two matrices B = (bi) (i, j = 0, 1, . * ) and -

C = (cuv) (j,u v = 0, 1, ), which need not be either square nor

If we think of t1 as analogous to 10, t1t2 to 100 and so on, we of related shapes and sizes, is defined as the block matrix
see that the representation of r is mixed radical, but with the
components of r in the reverse order to the digits of r. We
-booC bo0C-
could call this the reversed-digit (mixed-radix) representa- D= bIOC blC }. (10)
tion. On the other hand, the representation of s reverses the . . . . . . . . . .

order of the radices t1, t2, * , tn but leaves the components this being of course an abbreviated notation for an ordinary
of s in the same order as the digits of s. Thus there are two larger matrix whose elements are scalars. This definition is
different kinds of reversal so the representation of r and s equivalent to saying that the elements of B x C are
together could be described as the reversed-radical or re-
actionary representation or correspondence. There are vari-
dr,s = bri,s Icr2,s2 (1 1)
ants of the method, depending on various mixed-radix where r = (r1, r2) and s = (s1, S2) are ordered lexicograph-
representations of integers, which, of course, are not mixed ically. For example, if B and C are both 2 x 2,
whentl = t2- = tn. s 00 01 10 1-1
The single summation r
Z
arcos (co = e21riI?) 00 boocoo0 boocol, bolcoo, bolcol
can be regarded as a multiple summation over the digits of 01 booclo, boocll, bolclo, bo1c11 B
C. (12)

r in its mixed-radix representation. We also substitute the

10 blocoo, b10c01, blcoO, blc1l
mixed-radix representation of s and obtain 11 F bloclo, bloc,,, blclo, blclcl J

,ar(t1t2
r
* 'tn-1r, +tlt2 ' 'tn- 2r2 + ' + tir.- I +r.)(t2 ... tnsn+
-

+tnS2+Sl) - Z arKI1(rl, sl)K2(r2, Sl, s2) .. KI(r, Sl, S2, ' ' , SO) (8)
r

where (since (0? = 1) Thus, for example, b1jcOO corresponds to r1=1, si-.,
Ki(rl,sl) =-t tlt2 tn-lrlSI
r2 = 0, S2 = 0. (The reader should now hold in mind Savage's
metamathematical principle: he should sit bolt upright in a
K2(r2, Sl, S2) _ (ollt2 ... tn- 2r2(t.S2 + SO)
hard chair with a pencil in his hand [19]). The general
definition for the direct product of several matrices is the
Kn(rn, 1.*
s ..

sn) ot)rn(t2
= tnSn + * +
..

tnS2 + SI ) (9) obvious inductive one, which turns out to be associative, so

that parentheses are not required. Moreover, formula (11)
If the elements of a. are now permuted in accordance with generalizes in the natural manner. The direct product of n
the mixed-radical representation of r, the right side of (8) matrices,
is of the form of (6) and can be computed by the simple
procedure described above. Note that it would not be strictly M(1) x M(2) X ... X M(n)
correct to write ar in place of a, in the summation because where M(V)-=(m,j) (i=0, 1, * Pv 1; j=0, 1, ... , qv-1)
we would then be using the symbol a to denote two dif- has its (r, s) element
ferent functions. N
For this Cooley-Tukey or mixed-radix method the trans- fl MrIv,sv (13)
formation from a unidimensional DFT to the multiple sum v-1

does not require that t1, t2, . , tn be mutually prime in wherein the ordering of the rows is lexicographic as is that
pairs. In this respect it has an advantage over the Sino- of the columns. This again is true even if the matrices are
Ruritanian method, being more widely applicable when not square and are of various shapes and sizes.
we are concerned with DFT's. It is perhaps used most If the ordinary matrix product M(v)N(v) can be formed,
often with the tv's all equal to 2, which is the simplest case that is, if the number of columns of M(V) is the same as the
to program for a binary computer. On the other hand, the number of rows of N(v) (v - 1, 2, * ), then
multiple sum (5) is simpler than (6) and therefore more
likely to occur (as a multidimensional K-transform) in con- (M(1) X )(N(') x N(2) X...)
texts other than the DFT, especially with K1, Kn all = (M(1)N(l)) x (M(2)N(2)) x . . . (14)
SHORT NOTES 313

In particular, if M is t x t and is nonsingular, then wise 6b=0.) Note that the nonzero elements of A(v) are,
(M x M x * )(M- 1 x M-'1 X * * ) = I, X I, X ...*= t
..
apart from replications, the same as those of M(v), but its
size depends on all the t's.
where I. is the u x u identity matrix. More briefly Thus A(') is a (P1P2 * P.) x (P2P3 ... p.qj) matrix having
- -

at most P1P2 ... Pnq1 nonzero elements. A(2) is a (P2P3 ...p* q1)
(M[n])-1 = (M- 1)[n] (15) X (p3p4. Pnqlq2) matrix having at most P2P3 ... pnqlq2
where the symbol [n] indicates the "direct nth power." nonzero elements. A(v) is a (Pv ... Pnq1 ... qv- 1) x (pv+ 1
Now the n-dimensional K-transform, as defined by (5), has . pPnql ... qv) matrix having most p,pv+1 ... Pnqlq2 ...qv
as its matrix the direct product ofthe matrices K1, K2,-- *, K,, nonzero elements. Aln) is a (Pnqlq2 ... qn- 1) x (q1q2 . qn)
where, of course, the matrix K1 has Kj(rj, s1) as its elements, matrix having at most Pnqlq2 ... qn nonzero elements.
and so on. This is true provided that we order the rows and AM ... A(n) is a (P1P2 * * * Pn) x (q,q2 ... qn) matrix, namely M.
columns lexicographically. This explains the relevance of The product a'M, where the components of a are
direct products of matrices to the previous discussion. ar(0< r, <P,, , °<rn < Pn) can thus be effected as
Now it is an interesting fact that a direct product of n
matrices can also be expressed as an ordinary product of (... ((a'A"I)A (2)) ... A (n)). (18)
sparse matrices in accordance with the following theorem. As an example, take n = 2, Pi = 1, q1 = 2, P2 = 3, q2 = 4, and
(A sparse matrix is one consisting mainly of zeros.) note that
P2q1°° 01 10 11 20 21 00 01 02 03 10 12 13 14 qlq2
P1P2 P2q,
00 [aoo ao0 0 0 0 0 boo bo1 bO2 bo3 0 0 0 00
01 0 0 aoo a01 0 0 0 0 0 0 boo bo1 bo2 bo3 01
02 0 0 0 0 aoo aoil blo bl1 b02 b13 0 0 0 0 10
0 0 0 0 blo bl1 b12 b03 11
b20 b21 b22 b23 0 0 0 20
0 0 O 0 b2O. b21 b22 b23- 21
aooboo aOObOl aolboo aolbo1 aolbO2 ao lbo3
aOObO2 aOObO3
= aOOb10 aOObll aOObl2 aOObl3 aolblo aolbl aOlbl2 aOlbl3
[aoob20 aoob21 aoob22 aOOb23 aolb20 a0lb2l aolb22 aOlb23
= A x B.

M = M() x pM(2) x x M(n) = A(1)A(2) ... A() (16) In this example, if we were to multiply a horizontal 3-vector
or a bolt-upright 8-vector by A x B it would involve 3 x 8 x 8
where = 192 scalar products, whereas multiplying first by A(') and
A(v) (aV,)) then by A , or vice versa, requires only 4 x 8+ 3 x 2 = 38
scalar products. Generally, the reduction is from
where
PlP2 ... Pnqlq2 ...qn
4V) m(V)S b6S16S2 ..bsn -I
17
to
(0 < r1< pv, 0 < r2 < pv+,,,,* * , < rn-v+1 < Png
0 < rn-v+2 < ql, ,0 < rn < qv-1 EPvPv+ 1 ... Pnqlq2 ...qv (19)
v

0 < S1 < Pv+,O

° < S2 < Pv+2, SO < Sn-v < Pn scalar multiplications.
( < Sn-v+e < ql,=1 < Sn qv)o
< The identification of this procedure with that in (5), when
PV = qv = t, (v = 1, 2,... , n), is exemplified by the following
(The "Kronecker delta" is defined by 6b = I if a = b, other- algebra for the case t1 = 2, t2 = 3.
00 01 10 11 20 21 t2t1
tlt2
[aO, aOl,aO2,alO,all,al2] coo CO, 0 000 0 0
0 0 Coo 01
Coi 0 0
0 0 0 0 Co0 Co0 02
Cll 0 0 0 0 10
0 0 CIO Ci1 0 0 11
0 0 0 0 C1o C 11 12
= [aoocoo + alOClo,aOOcOl + alOcll,aOlcOO + alclcO,aOlcOl allcll,aO2CO0 a12C10,a02CO1
+ + + a12C1 ]
314 IEEE TRANSACTIONS ON COMPUTERS, MARCH 1971

whose elements are Zriari,r2cri,sa = 1 just as after the can denote the cofactors in M(v) by MM)(i,j = 0, 1 y.. )
first summation for (5). Further, and write B(v) = b(v) where
00 01 02 10 11 12 t1t2
t2tl
Elool 1019 110. 1119 120. 12111 doo do0do2 0 0 0 00
0 0 doo do 1 dO2 01
0
d1o d11 0 0 0 10
0 0 0 d0o d d02 11
d20 d2l d22 0 20
0 0 0 d2O d2 1 d22- 21
-= [loodoo + 11odlo + 120d20,loodo, + lodl + I20d21, ,l1d02 + 111d12 + 121d22]

whose elements are Er2 Ir2,Sl dr2,s2, as in the second summa- b(v) Mpv))r Ir2 1brn-
= ...
(24)
tion for (5). The first term in this vector is (aooco0 + ajOc1O)dOO
(sl, rn=O0 19
+ (aOlcOO + aj1cjO)djO + (aO2C0O + a12C1O)d20 = Erl,r2
tv-1; S2, r, =0 1, , tv+ I-I; ; Sn-v+ 1
ari,r2Cri,O
dr2,0, and similarly the other five terms also agree with (5). rn- v = 1, , tn -1 ; ; sw, rn- I = 0 1, , tv- 1 -1).
The theorem given by (16) was proven for the case Then
Pi= =Pn=q1= =qn in [7], and the proof extends A(v)B(v) = A I
to the general case. (25)
The inverse of a multidimensional K-transform exists if
the matrices M(v), defined by M(v) = (Kv(rv, sv)), are all square This is true even if A(') is singular, in which case the right
and nonsingular. This can be easily proved by verifying side of (25) is zero. If A(V) is nonsingular, it follows from (25)
(29), which is given in the Summary and Discussion. We that (A(v))- = B(v)/A\V. A proof of (25) is
here express the proof in terms of direct products of ma- Z£ m(, ql6bq2 ... 1qn-I . M(V) bq 1 bq2 ... bqn- 1
trices. If we write a and -a* for the column vectors whose q

components are ar and a.*, these components being ar- = Em(V)qAMIV)qS26S3 ... 6Sn
qn
ranged in the lexicographic order of r and s, then (5) can be
written An ex l vo ( r2 5rn
a*' = a'(M(l) x M(2) x... x M(n)) = a'M. (20) An example of (25) is

Mo0 Mi1 0 0 Moo 0 M 0o 0 0 01

0 0 Moo Min Mo0 0 Ml, 0 0 A O O
M1o M1i1 0 0 0 Moo 0 M1o 0 0 A 0
0 0 Mi10 m1J O MO, 0 M 1lj Lo 0 0 A_

Now M(v)(M(v))- l = I,,, and we see readily from (14) that where
a' = a*{(MA(l)-fI x ..x(XM(n))-} (21) A = moomln-MO1i,MO0 = M11,M00
=-MlO,M10
= -MO1,M11 = Moo.
It may be noted that the determinant of M
Note that A(V) and its inverse are equally sparse.
IMI = Ar/t, A/t" (22)
The determinants of A(') can be neatly obtained from (22)
by supposing that A1= ..- A, 1=A,+1= . A,,=1. It
where AV = IM(v)I. This follows inductively from the case follows from this, together with the omitted argument for
n=2 which is given by MacDuffee [18] who references determining the sign, that
Hensel, Netto, and von Sterneck. In the case where M(V)
is independent of v and has order t and determinant A, we
have
IA(v)I
= Au(- 1)u(u- 1)/2 (u = T/tj) (26)
and therefore, from (25), or from a familiar property of
= IMI Antn'. (23) adjugates [20, p. 165]
If we wish to express the inverse multidimensional K- IB(v)l -. (tv 1)/rv(, -)u(u- 1)/2 (27)
transform by means of ordinary matrix multiplications, we Of course if all the M(v)'s are orthogonal (or else unitary),
can apply (16) to the inverses of the M(v)'s. Alternatively, we then their inverses are equal to their (conjugate) transposes,
SHORT NOTES 315
and the multidimensional K-transform and its inverse in- mr,s = xnl[r,sl(_l)r s (tv = 2; v = 1, 2, .n
volve nearly the same calculation. The matrix of the multi-
dimensional DFT becomes unitary if divided by z/; (or if where a is real and positive,
each M(v) is divided by t/v), just as the multidimensional n-1
Fourier transform is unitary if a suitable normalizing factor ,l= _1- 2/a, and [r,s]
s = (r,E3s,)
j=o
is used, such as 1/ 2ir per integral. When each M(v) is square,
with pv = q= tv as before, then another way of factorizing where ( denotes mod 2 addition (EXCLUSIVE OR), but the
M1) x ... x M(n) into n sparse matrices is [7] summation is ordinary summation which happens to reduce
to counting.
M)x x -M() = C(1)C(2) ... C(n) (28) 4) If, in (5), which defines a multidimensional K-trans-
form, the functions K1(rl, sl), K2(r2, s2) , when re-
where garded as elements of matrices, are square and nonsingular,
C() M( )x It2x *-x tn and if the inverse matrices have elements J1(r,, sl),
C(2) =Itl M(2) x x tn
x ... J2(r2, s2), etc., then the inverse transform can be written
c(n)-It x It x ... X M(n). ar =
I a*J,(sl, r,) .. Jn(sn, r) (29)
This follows from (14) and might sometimes be more con-
venient to use than (16) although the patterns of 1) , C(n) This can be readily proved by using Kronecker deltas
are less simple than those of A"1, -. A(n). and is another way of expressing (21).
5) When one wishes to carry out a multiplication by a
,SUMMARY AND DISCUSSION large matrix it can be worthwhile to look for a factorization
Yates [21 ] introduced a simple adding and subtracting into several sparse matrices so as to cut down on the work.
algorithm for the calculation of the interactions in 2n-fac- 6) New coding methods, generalizing the DFT, can be
torial experiments. The fact that Yates' algorithm could be invented ad lib by taking several sparse square matrices of
expressed as a multidimensional mod 2 discrete Fourier equal orders and using their product as the matrix of the
transform (DFT) was pointed out by Good [22]. Good [7] transformation or coding. In order that unique decoding
showed that the algorithm can be generalized to apply to should be possible, it is necessary that the sparse matrices
any multidimensional DFT and even more generally. This should all be nonsingular and convenient if their inverses
depends on the first three of the following facts. are sparse. Apart from the methods discussed in this paper,
1) A multidimensional mod t DFT is a linear transforma- the matrices could, for example, be taken as permutation
tion and corresponds to multiplication by a matrix (of order matrices, these being the simplest and sparsest possible
tn where n is the dimensionality) which can be expressed as nonsingular matrices. They might be too sparse for some
the nth direct power of a matrix of order t. purposes since they merely permute the elements of the
2) The nth direct power of a matrix M is equal to the vector on which they operate and so have no tendency to
nth ordinary power of a matrix A defined as flatten the variation from one component to another. Such
flattening is apt to be convenient for communication sys-
A = {M r,sn63S2 .* bSn- }
. .
tems, but the permutation by itself gives some robustness
against bursts of noise. It might also turn out to be useful
3) Although A is t' it is
by tn sparse and contains only t+1 to use functions K1, K2,.*. , K. as in (6) instead of (5), al-
nonzero elements, so that multiplication by M involves though it is presumably more difficult to find elegant
only ntn+1 ordinary multiplications instead of t2". This, and schemes of this kind.
the remarks under facts 1 and 2, extended readily to the 7) Unidimensional DFTs can be expressed as multi-
case of unequal moduli t1, t2,.*.2 tn and even more generally. dimensional ones by either the Sino-Ruritanian method or
We thus have a fast algorithm for the computation of an the "reactionary" method, the latter being more widely
n-dimensional K-transform, defined by (5), of which the applicable for this purpose. Once in the multidimensional
n-dimensional DFT is a special case. The remarks following form, the streamlining described above is applicable. The
(5) give the algorithm in nonmatrix language. paper [7] provoked the influential paper by Cooley and
When K1,l , Kn in (5) are all the same function, the Tukey in which the reactionary method was described. Their
matrices A1),...*, A () are all equal and this matrix formula- original acknowledgments to my paper were somewhat too
tion is then at its simplest. In equation (6) the K's cannot be generous since they gave the impression that the two meth-
identical unless the equation reduces to (5). ods were almost identical. Their methods more closely
A mod 2 multidimensional DFT is an example of a resembled the methods of Runge [23] and Danielson and
Hadamard transform which is discussed and applied by Lanczos [24], whose papers had been generally overlooked
Pratt, Kane, and Andrews [12]. Another example of a K- by oceanographers, X-ray crystallographers, and many
transform is the Andrews-Caspari transform [17] in which others who make frequent use of harmonic and spectral
the orthogonal matrix of the transform has elements analysis. These authors wrote before the age of electronic
316 IEEE TRANSACTIONS ON COMPUTERS, MARCH 1971

computers and were in this respect ahead of their time, elegant inverse even when m has a primitive root (whose
although a fast Fourier transform would have been of use powers take all values prime to m). But (30) can be general-
for more than 100 years. In the design of special-purpose ized in a different direction, as we now indicate.
equipment, both the Sino-Ruritanian and reactionary repre- 10) We can generalize (30) and its multidimensional form
sentations should be kept in mind. to a finite (Galois) field F. A finite field contains a prime
8) One application of the DFT is for the inversion of a power pv of elements and always has a primitive root g for
circulix, that is, a t x t matrix (cs_r) where s - r is taken which 1, g, g2,*, gp"2 are distinct elements of the field
modulo t [25]. A 1 000 000 x 1 000 000 circulix could easily [26]. If Or (r = 0, 1, , pv 2) is a vector of field elements we
be inverted since it requires only about 100 000 000 opera- can define
tions instead of some quintillion if the matrix were random. pv-2
The eigenvalues of a circulix are given by the DFT of its as =
E rgS (s =
O, l1 *- pv 2) (32)
top row, and the inverse DFT of their reciprocals gives the r=O
inverse circulix. Similarly a recursively blocked circulix, as its (unidimensional) Fourier-Galois tranform. Similarly,
defined recursively as either a circulix or a block matrix for an array caL(r = (r1,... , rj)), we can define a multidimen-
whose elements are recursively blocked circulices, or equally sional transform
defined as (a, r) where r and s are ordered lexicographically
and s - r is taken modulo (t1, t2, , t), has eigenvalues 0(* = Ea grgs* ... grnSn (33)
r
given by the multidimensional DFT of its top row, and so
can be inverted in a similar manner. (I have avoided the where g1, 92, , gn are primitive elements of F, not neces-
name "block circulix" tout court, because a block circulant sarily distinct. This can also be written in the form
is defined by Muir [20, p. 485 ] in a somewhat different sense, (X* = EOCrgclrrsl
ZJ r + + c,r,s, (34)
* * *
apart, of course, from its being a determinant and not a
matrix.) The eigenvectors of a recursively blocked circulix where c1,"., C, nare certain constants prime to pV - 1. The
are the columns of the matrix (_s1w,,,r22s2 .). inversion formula is
9) Let p be a prime number and g one its primitive roots;
that is 1, g, g22 ...,
gp 2 are all distinct modulo p. Let cxr
- = (-1)n n *gl rlsl . .. g-rnsn
,
(35)
(r = 0, 1,... , p -2) be a vector of integers and let us define a s

number-theoretic Fourier transform: and there is a convolution formula

p-2
a*= Z rgS (mod p). (30) E ( atxB gr)gsl 1 r,Sn = ac*#* (36)
= r 0
Then it is easily proved that the inverse transform is A recursively blocked circulix of group elements has, as
p-2
eigenvalues, the multidimensional Galois-Fourier trans-
Xr= - Z o*g-rs (modp) (31) form of its top row, and, as eigenvectors, the columns of the
s= matrix (g1S) x ... x (grnSn). The transform could be used for
coding "sentences" of "words" where each word, being an
and the transform of the convolution element of F, could be expressed by means of v modulo p
p-2 integers, namely the coefficients of the word when the group
E Xqfr- q is represented by means of polynomials modulo both p and
q=O
a fixed irreducible polynomial [27 ]. (Finite fields have often
is ax*f*. A (p - 1) x (p - 1) circulix of integers has eigen- been used in other ways in coding theory and in the design
values equal to the transform of its top row, and it can be of statistical experiments.)
inverted modulo p by using modulo p reciprocals of the For the multidimensional transform we can again apply
eigenvalues, almost as in fact 8. Recursively blocked cir- the simple algorithm following (5), or the equivalent matrix
culices of integers, modulo p, have their eigenvalues equal to methods, the arithmetic being in the field. I can see no ele-
multidimensional transforms of their top rows, these trans- gant transformation from unidimensional to multidimen-
forms being defined in an obvious way, modulo a fixed sional Fourier-Galois transforms.
prime p. These multidimensional transforms could be used
for a coding scheme for vectors of integers modulo p. The REFERENCES
calculation can be done by the simple method following (5), [1] IEEE Trans. Audio Electroacoust., vol. AU-15, June 1967.
[2] IEEE Trans. Audio and Electroacoust., vol. AU-15, June 1969.
or by the equivalent matrix methods, except that the arith- [3] I. J. Good, "Polynomial algebra: An application of the fast Fourier
metic must of course be modulo p. transform," Nature, vol. 222, 1969, p. 1302.
If p is replaced by a composite number m we can define [4] , "The discrete Fourier transform: Work by I. J. Good," June,
1969, pp. 1-6, mimeographed.
(r,m) = 1 [5] R. C. Singleton, "A short bibliography on the fast Fourier trans-
as* = Z crglS (mod m), form," in [2], pp. 166-169.
r<m [6] I. J. Good, "The fast Fourier transform and discrete Fourier trans-
form: Bibliography," June 1969, pp. 1-5, mimeographed.
but this seems uninteresting because it does not have an [7] , "The interaction algorithm and practical Fourier analysis,"
SHORT NOTES 317

J. Roy. Statist. Soc. Ser. B, vol. 20, 1958, pp. 361-372; vol. 22, 1960, I. INTRODUCTION
pp. 372-375. Most of the reported work on combinational cellular
[8] J. W. Cooley and J. W. Tukey, "An algorithm for the machine calcu-
lation of complex Fourier series," Math. Computation, vol. 19, 1965, arrays considers two main types of arrays, namely, the
pp. 297-301. fixed cell-function arrays [1 ] and the variable cell-function
[9] I. J. Good, "Random motion on a finite Abelian group," Proc. arrays. In fixed cell-function arrays the switching function
Cambridge Phil. Soc., vol. 47, 1951, pp. 756-762.
[10] , "The serial test and other tests for randomness," Proc. Cam- realized by each cell is fixed. Realization of a logic function
bridge Phil. Soc., vol. 49, 1953, pp. 276-284. (or a set of functions) is accomplished by modifying the
[11] , "Analogues of Poisson's summation formula," Amer. Math. interconnection structure. In variable cell-function arrays
Mon., vol. 69, 1962, pp. 259-266.
[12] W. K. Pratt, J. Kane, and H. C. Andrews, "Hadamard transform the switching function realized by each cell can be changed
image coding," Proc. IEEE, vol. 57, Jan. 1969, pp. 58-68. by means of the cutpoint technique [2] or other methods.
[13] W. W. R. Ball, Mathematical Recreations and Essays. London: Realization of a logic function or a set of functions is accom-
Macmillan, 1940. (Revised by H. M. S. Coxeter.)
[14] G. H. Hardy and E. M. Wright, Ani lttroductioni to the Theory of plished by finding an appropriate configuration of logic
Numbers. Oxford: Clarendon Press, 1938. cells in the array. In either type of array it is assumed that
[15] P. Bachmann, Niedere Zahlentheorie, pt. 1. Leipzig: Teubner, 1902, we have access to each cell in order to change the cell struc-
p. 83.
[16] C. Bingham, M. D. Godfrey, and J. W. Tukey, "Modem techniques ture, the interconnection structure, or both. This is possible
of power spectral estimation," in [1], pp. 56-66. in present circuit construction techniques where the logic
[17] H. C. Andrews and K. L. Caspari, "A generalized technique for cells (integrated circuits) are mounted on the printed-circuit
spectral analysis," IEEE Trans. Computers, vol. C-19, Jan. 1970,
pp. 16-25. board (which provides the interconnections among the
[18] C. C. MacDuffee, The Theory of Matrices. New York: Chelsea, cells). However, these techniques may not apply to future
1956, p. 82. generations of microelectronic modules.
[19] L. J. Savage, The Foundations of Statistics. New York: Wiley,
1954, p. viii. With the rapid progress that is being made in the manu-
[20] T. Muir, A Treatise on the Theory of Determinants. New York: facturing techniques for monolithic integrated circuits, it
Dover, 1960. appears that future generations of microelectronic devices
[21] F. Yates, The Design and Analysis of Factorial Experiments. Har-
penden: Imperial Bureau of Soil Science, 1937. will be the so-called large-scale integrated circuits. Since in
[22] I. J. Good, review of M. H. Quenouille's "Design and analysis of such devices each digital component is fabricated on an
experiment." in Ann. Eugen., vol. 18, 1953, p. 263. extremely small area of the substrate, it will be impractical
[23] C. Runge, "Ober die Zerlegung empirisch gegebener periodischer
Funktionen in Sinnuswellen," Z. Math. Physik, vol. 48, 1903, (i.e., very costly to implement), if not impossible, to change
pp. 443-456; vol. 52, 1905, p. 117. the cell functions or the interconnection structure of the
[24] G. C. Danielson and C. Lanczos, "Some improvements in practical array once the device is made. Thus it is deemed desirable
Fourier analysis and their application to X-ray scattering from
liquids," J. Franklin Inst., vol. 233, 1942, pp. 365-380 and 435-452. to develop a design technique such that a single type of
[25] I. J. Good, "On the inversion of circulant matrices," Biometrika, vol. array can be used to realize many different logic functions
37,1950,pp.185-186. without the necessity of changing the cell function or its
[26] G. Birkhoff and S. MacLane, A Survey of Modern Algebra. New
York: Macmillan, 1944, pp. 428-431. interconnection. One way to accomplish this is to change
[27] R. D. Carmichael, Introduction to the Theory of Groups of Finite the functional capabilities of a cellular array by appropri-
Order. Boston: Ginn, 1937, p. 255. ately setting the parameters on the edges of that array. This
concept led us to the development of the new array design
technique presented in this note. The main theoretical result
is a homogeneous cellular array that can be used to realize
A Universal Cellular Array any switching network.
JUNG-CHANG HUANG, MEMBER, IEEE
II. SYNTHESIS OF THE ARRAY
Abstract-A type of fixed-cell fixed-interconnection homo-
new
geneous cellular array that is capable of realizing any switching net-
Let us denote a sequential machine without output by
work is developed. The array is composed of identical combinational S = <Q, 1, M> as usual, where Q is the nonempty set of inter-
logic cells with three inputs and three outputs. The logic cells are nal states, I is the input alphabet, and M is the next-state
arranged in a rectangular array with uniform interconnection struc-
ture. The signal flow is unilateral in the vertical direction and bilateral
function. If III = 1, then we say S = <Q, I, M> is a one-column
in the horizontal direction. Each n-column array is capable of realizing (or autonomous) state machine (as defined in [3]).
any set of n functions of n variables. The functional capabilities of the It is well known [4] that the transition matrices of the
array can be changed by appropriately setting the parameters on the
edges of the array. Algorithms for the realization of a set of n functions machines G([), G(r), and G(r) shown in Fig. 1 constitute a
of n variables by using this type of array are presented. generator set that generates all the transition matrices of
Index Terms-Cellular arrays, combinational switching networks, all r-state, state machines. Consequently, if N1, N2, and N3
iterative switching circuits, universal arrays. are three n-input (n = [log2 rl), n-output combinational
networks that realize the next-state function of G('), G(r)
Manuscript received September 12, 1969; revised July 1, 1970, and and G(r) respectively, then the next-state function of any
August 22, 1970. This work was supported by the U. S. Army Research one-column state machine with r states' can be realized by
Office (Durham) under Contract DA-31-124-ARO(D)-97 with the
Moore School of Electrical Engineering, University of Pennsylvania,
Philadelphia, Pa. 1 If r is not a power of 2, then this state machine can be realized as a
The author is with the Department of Computer Science, University submachine of a machine having 2" states, where n= [log2rl. Thus in this
of Houston, Houston, Tex. 77004. note we only have to deal with the cases in which r is a power of 2.