Polynomial Multiplication with Circulant Matrices: Insights Using APL
Vector 12:4, page 81
Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 81
A "Hook" to Get your Attention¶
Here is a puzzle. Show how to use nothing to make something out of nothing. Impossible? Read on and at the end of this article learn how to do it with some plain old APL.
In his book "APL Programs for the Mathematics Classroom", Norman Thomson gives the following APL one-liner to accomplish the multiplication of two polynomials,
Thomson algorithm: 1↓+/⍉(-⍳⍴L)⌽L∘.×R,L≠L
Two minor problems occur with Thomson’s code if either or both arguments are simple scalars or if ⎕io is zero. These are easily corrected however with judicious use of commas and parentheses, thus, Thomson algorithm, modified,
TAMP stands for (T)homson (A)lgorithm for (M)ultiplying (P)olynomials. (As modified — A further modification might be to replace ⍉ with ⌿.) Here in APL is a different algorithm to multiply two polynomials,
Circulant algorithm,
CAMP stands for (C)irculant (A)lgorithm for (M)ultiplying (P)olynomials. Note that to test these programs one must interpret L and R as the vectors of coefficients of two (different or the same) polynomials in ascending power order.
A two-line dyadic APL function named TAMP that incorporates the Thomson approach might have generic form and usage as follows,
Example usage: P←1 0 1 TAMP 1 2 0 3
A two-line APL function named CAMP that utilizes the Circulant algorithm would have the same generic form and usage, except only that the name of the function changes to CAMP. Both functions would return the same result in P, the list 1 2 1 5 0 3. The list is exactly the coefficients of the polynomial product in ascending power order, so that the last item in the list multiplies the highest power term of the result, etc.
Thus letting X stand for the indeterminate, the highest power term in the product polynomial is 3×X*5, the next highest is 0×X*4, the next is 5×X*3, the next is 1×X*2, the next is 2××X and the lowest power term is 1×X*0, the constant term, 1.
Circulant Matrices¶
Circulant matrices have a key property: write down the first row and you know the entire matrix. The second row is a clone of the first, except all members have been shifted one position to the right. The last element is simply "wrapped-around" into the first position. Obtain every succeeding row from its immediate predecessor by repeating the "shift and wrap" operation. The shifted clones all contain the same numbers but in different columns.
Suppose someone gives us a linear list of numbers. Then clearly one can obtain a circulant matrix from the list. The following is a useful shorthand notation for that operation; denote a circulant matrix as: circ(list) = circ(c1,c2, . . . ,cn)
Alternatively let the vector v = c1, c2, . . . , cn. Then for any vector v, circ(v;n) denotes the corresponding n by n circulant matrix. A precise definition of a circulant matrix is required at this point. In APL notation, the I, J element of the N square, possibly complex matrix C is C[I;J].
Now form: K ← N | (J-I) + 1 Consider T ← C[I;J] = C[1;K]
If boolean T is true for all I and J from 1 to N, then C is a circulant matrix. Note that the dyadic use of stile is the APL residue function, in this case the right argument modulo the left argument N. Note also that C[1;K] is just the element in the Kth position in the first row.
Time-out for a Discussion about Jargon¶
Previously X was called an indeterminate and not a variable. This choice was deliberate and the distinction is useful to algebraists. For those APLers who remember their studies in abstract algebra, the indeterminate X is fraught with non-meaning. The abstract algebraist considers it to be the "universal abstraction" — sort of a surrogate place holder that assumes whatever non-role would be useful to generalize the results of studying a particular algebraic structure. (If this sounds confusing, it is. Fortunately there are approaches to abstract algebra that dispense altogether with such a notion about X. It is no longer an indeterminate; instead it takes on concrete meaning in a formal definition of a polynomial. See for example page 127 of the text "Applied Abstract Algebra" by R. Lidl and G. Pilz.)
Undoubtedly some APLers either have never encountered abstract algebra or have been sufficiently traumatized by its study that they wish to forget the experience! Please do not be put off by references to the concepts and jargon of the discipline in what follows. If you wish, just ignore or skip over comments about "indeterminates", "rings", "ideals", "integral domains", "fields", "isomorphisms", etc. However the reader will certainly benefit from remembering several things taught in high school algebra. These include multiplication of matrices and multiplication and long division of polynomials. In any case, the listings of APL functions and their usage should be familiar and, I hope, interesting. Fortunately, given the array-like nature of the language, APLers know how to multiply two matrices with the "+.×" operator, even though they perhaps never studied matrix or linear algebra.
Back to the Discussion of Polynomial Multiplication with Functions TAMP and CAMP¶
Function TAMP has taken one polynomial with highest power 3×X*3 and has multiplied it with another polynomial of highest power X*2. Function CAMP has done the same thing. Both give the same the result, a product polynomial with highest power term 3×X*5. Further, as required, the product list contains the correct coefficients for the remaining powers. The APL expressions
give the ascending power order of the powers for each coefficient. What is going on here? Why does the Thomson program work? Why does the Circulant algorithm work?
In his book "ALGEBRA, An Algorithmic Treatment", Kenneth Iverson explains the "inner workings and hidden mechanisms" of program TAMP. He asks the reader to consider two example tables,
T1←1 0 1 ∘.× 1 2 0 3 (Use the left and right arguments from the earlier example for function TAMP.)
T2←0 1 2 ∘.+ 0 1 2 4 (Note the "plus" sign. Exponents add.).
The results are
All the counter-diagonals of table T2 are constant. Read from upper left to lower right, the constants increase from 0 to 5. This is exactly the order required for the ascending powers of the result. Scan and add along a fixed counter-diagonal in table T1. Verify that the sum gives the desired coefficient for the "power" found on the corresponding counter-diagonal in table T2. The idea behind program TAMP now stands revealed. (A caution — the preceding argument is a heuristic one; it is not a proof and may leave room for doubt. For a formal proof that function TAMP produces the polynomial product, see Iverson’s book and his convincing usage of APL notation to achieve the result.)
Can something similar be going on with the circulant approach? Notice that the primitive dyadic function "rotate" (⌽) is centrally embedded in both TAMP and CAMP. Check this out as follows. Use as inputs the previous example’s arguments (left and right). Stop both TAMP and CAMP just after the completion of dyadic ⌽ to get these results,
No clear overall pattern is apparent for table TAMP, although summing down the its columns gives the polynomial coefficients (except for an extraneous 0).
In contrast, a very useful pattern has begun to emerge for the CAMP result; it looks like the first three rows of a 7 by 7 circulant matrix. But where are the polynomial coefficients? They appear when the left argument, vector L←1 0 1, multiplies table CAMP, considered as a matrix, on the left. Carry out the vector-matrix multiplication specified for the completion of algorithm CAMP and obtain the correct result. (Ignore the extraneous 0 that appends itself to the end.)
Warning: Detour into Abstract Algebra, this Section Skippable¶
Is the apparent pattern in table CAMP a clue to the existence of a more general result? Given the title of this article, the reader will not be surprised to learn that the answer is an unqualified "yes". In the language of abstract algebra, a certain isomorphism exists. Here is the relationship in all the "terminology-as-jargon" words of abstract algebra.
"The commutative ring of n by n circulant matrices with entries in the complex field is structurally the same thing as the commutative ring of polynomials with complex coefficients, modulo the ideal generated by the indeterminate expression (X*n - 1). This latter is usually written as K[X]/(X*n - 1)"
Stated more briefly, the ring of circulants is isomorphic to the polynomial ‘residue class’ ring modulo the principal ideal (X*n - 1).
Here the notation X*n means the indeterminate raised to the nth power. Presumably the practitioners of abstract algebra know how to multiply together n times the universal abstraction embodied in X. While one might consider that such a talent is more theological than algebraic, algebraists can define away the need for a "leap of faith" that such exponentiation is possible. They do so by appealing to linear algebra and invoking the notions of a complex vector space spanned by a prescribed basis set.
The basis is taken to be the infinitely many elements X*n where n ranges over the non-negative integers. The symbol n at this point has lost its property of being an exponent. Rather it merely acts to index the basis elements. As the development proceeds, n manages to re-acquire its exponential character. This is achieved by making the product of two vectors meet the specifications of a bilinear mapping and defining the value of an ordered pair X*m, X*n to be X*(m+n). The details are a diversion, but the result is that the index n ends up acting like an exponent again. Some may find this use of definition to be egregious and the whole business rather like a ‘self-fulfilling prophecy’.
In any case the mystery can be avoided because all that really matters are the polynomial coefficients and n, the number of times X multiplies itself. However mysterious those multiplications of X might be, no details are needed, other than n, the number of factors involved. Why talk about the "ring" of circulants instead of an "integral domain" or "field"? Because in general "divisors of zero" exist for circulant matrices. This abstract algebra "jargon-phrase" means only that the cancellation law of multiplication does not hold. Later, an example is provided where the product of two non-zero circulants is the zero matrix.
Can the phrase "...modulo the ideal generated by the indeterminate expression (X*n - 1)" have some tangible implementation in the world of APL functions? To put some operational meaning to the words, first multiply two polynomials using either CAMP or TAMP. Then multiply the same two polynomials using CAMS. (CAMS is a one-line APL function to multiply two circulants. It is described in a later section.) But before using CAMS, if one list of coefficients is shorter than the other, make the matrix multiplication conformable. To do this, pad the shorter coefficient list on the right with as many zeros as it takes to get exactly as many total entries as in the longer list. Take n to be the length of either list, now both of the same length.
The first result, using CAMP or TAMP, is the product of two polynomials. The second result, using CAMS, is the product of two circulant matrices. Except possibly for zero padding, in both results the coefficient lists used as input are identical.
Next generate a coefficient list corresponding to X*n - 1. If for example, n is 10, the list will have 11 entries. All items in the list will be zero except the first and the last. The first will be -1 and the last will be 1. This list is just the coefficients of the polynomial X*10 - 1, with zeros in place of all the intermediate missing monomial terms.
Now divide the CAMP/TAMP result for the product of two polynomials by the polynomial X*10 - 1. Do this with an APL function that implements polynomial long division by using polynomial coefficient lists as inputs. This is the "modulo" operation, so focus interest on the remainder list, not the quotient list from the long division. (Note: two APL functions, one named POLYDIV and another named PDIVBY, both described later, accomplish polynomial long division. They do the same division operation, but POLYDIV does it with recursion while PDIVBY uses looping.)
Compare the remainder list of coefficients with the list from the CAMS product of the two circulants. The two lists should be identical. This result is the tangible quantification of the words and symbols describing algebraically the structural relationship between polynomials and circulants.
Why in the world should this work? In particular what "magic" is involved in identifying the ideal generator (X*n - 1) as the proper choice for the modulo calculation? In my judgement this is just one of many instances where the theorems of abstract algebra shine with considerable glory. The interested reader should consult a text on modern abstract algebra for complete details. (The matter is partially clarified in a later part of this report.)
One book on the subject that I particularly like is "Rings and Ideals" by Neal H. McCoy, 1948. It is published by the Mathematical Association of America as one of the Carus series of Mathematical Monographs. Those not familiar with the Carus Monographs might like to know that their stated purpose is "...to contribute to the dissemination of mathematical knowledge by making accessible ... a series of expository presentations of the best thoughts and the keenest researches in pure and applied mathematics."
Chapter VIII, Rings of Matrices, is most helpful, as is the discussion on page 162 of that chapter about the fundamental theorem on homomorphisms. A number of other good texts on modern abstract algebra are also available. I suggest, for further investigations, "A First Course in Abstract Algebra" by John Fraleigh and "Applied Abstract Algebra" by Rudolf Lidl and Gunter Pilz.
Back from the Land of Abstract Algebra: the Polynomial-like Structure of a Circulant Matrix¶
There is a fundamental polynomial structure to circulant matrices that is surprisingly easy to display. It supports the existence of a useful pattern in program CAMP’s table of intermediate results. It also validates the claim of structural equivalence between circ(v;n) and the residue class ring K[x]/(X*n - 1).
Consider the matrix circ(v;n) expanded into the matrix polynomial as follows,
circ(v;n) = c1×A*0 + c2×A*1 + c3×A*2 +. . . + cn×A*(n-1).
Does there exist an n by n circulant matrix A for which this tentative reformulation is true? The answer is yes. The so-called "fundamental permutation matrix" — FPM for short — handles the job. This remarkable matrix is the innocent looking circulant defined as,
FPM = A = circ(v;n) = circ(0,1,0,...0)
Note that vector v in APL is simply v←0,1,(n-2)⍴0 and that ⍴v is just n.
Here A*n is I, the n by n identity matrix of all ones along the main diagonal and all zeros elsewhere, clearly also a circulant. Further all the powers of A, from 0 to (n-1) span the space of circulants. For want of a better word the FP matrix A may be called "cyclopotent". The choice reflects the influence of the terms "idempotent" and "nilpotent" used to describe other special behaviours of the powers of matrices.
The n by n matrix of all zeros is the additive identity matrix in the ring of circulants. Note that it is also a circulant. The symbol
is taken to be the multiplicative identity matrix I in the ring of circulants. APLers beware — any matrix raised to the zeroth power in APL will replace every element, even zero elements, with ones. The result is NOT the identity matrix.
The expansion in powers of A is clearly a polynomial in A whose coefficients are the elements of the vector v (which may be complex). The left hand side of the expansion relation is just as clearly a circulant matrix generated from the coefficients of the polynomial. Its first row is identical to the vector v of coefficients.
For two "circulant polynomials" constructed as above from two different defining vectors v, the operations of addition and multiplication follow the usual rules, with one important modification.
When multiplying two such polynomials, powers of A exceeding the (n-1)st may result. In those cases, rewrite the power of A as A*n × A*m and substitute I for A*n. This leaves only the "residue" matrix A*m, multiplied by its coefficient.
Complete the multiplication by collecting all similar powers A*m together and adding up the common coefficients. Do this as necessary for every m from 0 to (n-1).
Note also that multiplication under these conditions results in powers of A that range only from m = 0 to m = n-1. Because A is cyclopotent, these n "residue classes" make up the only possibilities. In effect any power p of A is reduced modulo n, by utilizing the result that for A taken as the FPM, A*n is I, the n by n identity matrix. (Although not explicitly stated so far, one needs to use also the property that the product of two circulant matrices is also a circulant matrix.)
A clue is now manifest regarding the rather mysterious choice of X*n-1 as the ideal in the residue class ring discussed in the previous ‘Abstract Algebra’ section. Notice that A*n - I can be rewritten as A*n - I = 0.
Here A is the FPM, I is the identity matrix and the right hand side is the zero matrix. Replace A with the indeterminate X, I with the multiplicative identity unity, and the zero matrix with the additive identity 0 to get X*n - 1 = 0. Then check some theorems in the cited texts about mappings and kernels and some definitions of ideals. At this point, the perceptive reader should begin to have doubts about algorithm CAMP. Multiply two n by n circulant matrices and you get as a result another n by n circulant matrix. Assume there are no zero elements in the two original matrices. Then the coefficients in the first row of the result cannot be the polynomial product. Why? — because multiplying two polynomials with n non-zero terms will inevitably result in a polynomial with many more terms than n. Yet as just described, multiplying two circulant polynomials causes a "folding" of coefficients of powers greater than (n-1) into lower powers. No power of A greater than (n-1) occurs in the result.
Before proceeding, the reader needs to have an APL "one-liner" for multiplying two circulant matrices. The sceptic who has an APL interpreter can then check out the several calculations to come. The only tools needed are algorithms TAMP, CAMP, and CAMS — and CAMS is given next.
Circulant product algorithm, CAMS:
CAMS stands for (C)irculant (A)lgorithm for (M)ultiplying Circulant(S). Compare CAMS with CAMP to see that they are obviously closely related. CAMS does not produce the full circulant matrix, just the first row. Now here are some examples to clarify the workings of CAMS and CAMP and restore any lost faith in using function CAMP to produce polynomial products.
Use CAMP to multiply the quadratic polynomial whose coefficients are 1 1 2 by the quintic polynomial with coefficients 1 2 0 3 1 4. Get the seventh degree polynomial whose coefficients are 1 3 4 7 4 11 6 8. (Sceptics; check this with an ordinary paper and pencil multiplication of the two polynomials if you like.) Note that all coefficient lists are written in ascending order of powers and that if a power is missing its coefficient must be present as a zero in the list.
To continue, use CAMS to multiply the following 6 by 6 circulant matrices,
to get the 6 by 6 product
Take the first (CAMP) result and rewrite the single list of 8 items as two lists, folded after the sixth item and with zeros appended as necessary. This implements the necessary "powers of A modulo 6" calculation and gives table T3,
Now add down the columns of table T3 to get the list 7 11 4 7 4 11. Bingo! The summation has reproduced the first row in the circulant product, as required when the powers of a product of two circulant polynomials are reduced modulo n.
Pretty neat, huh? Yeah except that there is little likelihood that very many APLers are interested in multiplication in the commutative ring of circulant polynomials with a program like CAMS. CAMP however seems as if it might actually be useful. How does algorithm CAMP reproduce ordinary polynomial multiplication? The perceptive reader may have already noticed the trick from close examination of the above example. The list of three coefficients for the quadratic polynomial was extended by three zeros to get a circulant matrix conformable with the circulant matrix corresponding to the list of six coefficients for the quintic polynomial. One has then two 6 by 6 circulant matrices, which being now conformable, can be multiplied together with function CAMS.
What happens if both lists, the quadratic and quintic lists, are extended with enough zeros to avoid the reduction of powers of A modulo n? The answer is that one then gets the desired ordinary polynomial multiplication. Algorithm CAMP accomplishes this extension very neatly by starting its calculation with the expression R,L≠L.
This in fact is a slight over extension of zeros, one too many, but the modest inefficiency is tolerable for simplicity of the APL expression. Algorithm CAMP has taken advantage of the zero padding to simplify the expected circulant matrix times circulant matrix operation. The nature of the zero padding allows the product to be performed instead as a greatly truncated circulant matrix (a vector) times a truncated circulant matrix.
APL Functions for Long Division of Two Polynomials¶
As promised earlier, here are two APL functions for polynomial long division, or PLD for short. PLD requires only the division of monomial terms with positive integer powers, the division of scalars selected from the real or complex number fields, multiplication of two polynomials, and subtraction of two polynomials.
Recursive APL function POLYDIV implements PLD for two polynomials L and R. The left argument L is the divisor and the right argument R is the dividend. Enter as a vector only the coefficients for polynomial L, in order from lowest to highest degree. If a monomial term is absent, an entry of zero is required in proper order for the missing term. Enter the polynomial R the same way. Do NOT pad any zeros after the coefficients of terms of highest degree!
If the coefficient list in L exceeds the length of the list in R, then the degree of the polynomial divisor exceeds the degree of the polynomial dividend. In this case only the trivial result, quotient equal to 0 and remainder equal to dividend, will be returned.
POLYDIV works just as well when it is fed polynomial coefficients from the field of complex numbers, but only for those APL interpreters with complex numbers as the default variable type. The function also works with real coefficients.
POLYDIV displays results for the quotient and remainder as coefficient lists. These are in the same ascending power order as the input list; embedded zeros in these results represent absent monomials, the same usage as for input lists.
∇ L POLYDIV R;K;QV;LXQV;Z
[1] →(2=⎕nc 'Q')/NORESET ⍝ DO NOT RE-INIT GLOBAL VARS ON RECURSIVE CALLS.
[2] Q←0 ⋄ DD←¯1↑L ⋄ LL←⍴,L ⋄ RMLXQV←,R ⍝ GLOBAL VARIABLES.
[3] NORESET:
[4] →(LL>⍴,R)/EOJ ⍝ QUIT IF DIVIDEND OR REMAINDER SHORTER THAN DIVISOR.
[5] Q←(K←(¯1↑R)÷DD),Q ⍝ SCALAR DIVISION, ALWAYS WITH DD.
[6] QV←K×⌽1,LL↓(R≠R) ⍝ ZERO PADDING ON RIGHT.
[7] LXQV←¯1↓L+.×(⎕io-⍳LL)⌽(LL,⍴QV,L)⍴QV,L≠L ⍝ L×QV IS NORMAL POLY PRODUCT.
[8] RMLXQV←R-LXQV ⍝ ZERO PADS HAVE ASSURED CONFORMABLE SUBTRACTION.
[9] →(1=∧/⎕ct|RMLXQV)/EOJ ⍝ QUIT IF REMAINDER ALL ZEROS.
[10] L POLYDIV(RMLXQV←¯1↓RMLXQV) ⍝ RECURSE ON DECREASING REMAINDER.
[11] →0
[12] EOJ:
[13] 'QUOTIENT'
[14] (¯1×1≠⍴Q)↓Q ⍝ CAREFULLY DROP SUPERFLUOUS ZERO.
[15] 'REMAINDER' ⋄ RR←RMLXQV
[16] (¯1×0≠⍴RR)↓RR,(0=⍴RR)×⍴RR ⍝ CAREFULLY HANDLE RR IF A ZERO LEN VECTOR.
[17] Z←⎕ex 4 6 ⍴'Q ','DD ','LL ','RMLXQV'
∇
I wrote the preceding recursive function for this article. Because it employs global variables, the user must insure that global variable Q has been erased from the workspace before execution. Alternatively, the user can execute cover function POLY to handle all initialization and execution duties. Left and right arguments are self-documenting.
∇ DIVISOR POLY DIVIDEND
[1] ⍝ EXPUNGE GLOBAL VARIABLES BEFORE RECURSION.
[2] Z←⎕ex 4 6 ⍴'Q ','DD ','LL ','RMLXQV'
[3] DIVISOR POLYDIV DIVIDEND ⍝ POLYNOMIAL LONG DIVISION BY RECURSION.
∇
Here is a second function for PLD that does not use recursion. It is function PDIVBY, as found on page 260 of U. Grenander’s book "Mathematical Experiments on the Computer". Note: Unlike function POLYDIV, function PDIVBY uses coefficient input lists and returns result lists in DESCENDING order with respect to the powers of X. Correct answers resulted from limited tests of PDIVBY with complex input lists. It probably would produce correct results for those APLs with a complex variable type as the default, but the user should first verify that capability. Aberrant results occurred for certain special cases of no operational interest, for example
(0j1) PDIVBY (1 0j¯1) and 1 PDIVBY 1.
∇ Z←P1 PDIVBY P2;LEAD;SHIFT
[1] ⍝ RIGHT ARGUMENT: THE VECTOR OF COEFFS OF THE
[2] ⍝ DIVISOR IN DESCENDING ORDER WRT THE POWERS OF X.
[3] ⍝ LEFT ARGUMENT: VECTOR OF COEFFS OF THE DIVIDEND IN
[4] ⍝ DESCENDING ORDER WRT THE POWERS OF X.
[5] ⍝ RESULT: A 2×N ARRAY WHERE N=1+DEGREE OF DIVIDEND.
[6] ⍝ THE 1ST ROW IS THE VECTOR OF COEFFS OF THE QUOTIENT,
[7] ⍝ ORDERED THE SAME AS THE INPUTS.
[8] ⍝ THE 2ND ROW IS THE VECTOR OF COEFFS OF THE REMAINDER,
[9] ⍝ ORDERED THE SAME AS THE INPUTS.
[10] →(((⍴P1)-(P1≠0)⍳1)≥(⍴P2)-(P2≠0)⍳1)/DIVIDE
[11] Z←(2,⍴P2)⍴((⍴P2)⍴0),P2
[12] →0
[13] DIVIDE:Z←(⍴P1)⍴0
[14] LEAD←P2[DEGP2←(0≠P2←(-⍴P1)↑P2)⍳1]
[15] P2←(SHIFT←(DEGP2-(P1≠0)⍳1))⌽P2
[16] LOOP:P1←P1-P2×Z[(⍴P1)-SHIFT]←P1[DEGP2-SHIFT]÷LEAD
[17] P2←¯1⌽P2
[18] →(0≤SHIFT←SHIFT-1)/LOOP
[19] Z←(2,⍴P1)⍴Z,P1
∇
The Isomorphism in Action — It Performs as Advertised¶
Previously CAMP was used to multiply the quadratic polynomial whose coefficients are 1 1 2 by the quintic polynomial with coefficients 1 2 0 3 1 4. The result was the seventh degree polynomial whose eight coefficients are 1 3 4 7 4 11 6 8. Also previously, CAMS was used to multiply the following 6 by 6 circulant matrices,
to get the 6 by 6 product
How now to relate the seventh degree polynomial 1 3 4 7 4 11 6 8 to the circulant product 7 11 4 7 4 11? Answer — consider the sixth degree polynomial X*6 - 1 and its related coefficient list of seven items,
Use the polynomial division program POLYDIV to divide the seventh degree polynomial by the sixth degree polynomial X*6 -1, as follows
Examine the remainder, not the quotient, because a "modulo" calculation is desired. Another Bingo! The remainder is 7 11 4 7 4 11, the same as the CAMS product. This result is not accidental. It is exactly the result that one should obtain from the structural relationship — the isomorphism — between circulants and polynomials.
A Danger and an Opportunity — Divisors of Zero¶
The Chinese ideograph for danger is the same as for opportunity. One dangerous aspect of matrix multiplication is that the matrix product can be zero. In such systems, mathematicians talk about "divisors of zero". In my judgement one could equally well phrase it as the "factors of zero". Two "somethings" annihilate each other and produce "nothing". Is there opportunity here as well as danger? Can one turn things around and opportunely use "nothing" to avoid the danger of "somethings" annihilating each other? Developed next is an answer that will also answer the puzzle posed at the beginning of this article. First, as promised earlier, here are two circulant matrices that annihilate each other when multiplied together. With CAMS, check that
produces the null circulant matrix, circ(0 0 0). Hence circ(¯1 1 0) and circ(1 1 1) factor zero (are "divisors of zero"). "Something" times "something" has produced "nothing". Now pad each coefficient list with "nothing", actually a single zero, and try again. With CAMS, check that
now produces "something", the non-zero circulant circ(¯1 0 0 1). Does this construction "use nothing to make something from nothing"? Well, it certainly uses zeros ("nothing") to make a circulant matrix ("something not zero") out of the null circulant matrix (circ(0 0 0), which is "nothing") produced by the product of the circulant factors of zero.
How to Factor a Circulant Matrix¶
Polynomials have roots. Then what sense, if any, can be made out of the following procedure?
- Treat the first row of an n by n circulant as the coefficients of a polynomial of degree n-1, the coefficients being in ascending power order with respect to X.
- Assume the coefficients are from the real field. (Could equally well assume they are from the complex field. Actually almost any field will do. Finite fields will lead to some unusual happenings.)
- Factor the polynomial into its linear and possibly quadratic polynomial forms. To do this, solve for all the roots of the polynomial.
- Pad the coefficient list of two items of each linear factor with zeros on the right, until the list contains n items.
- Pad the coefficient list of three items of each quadratic factor with zeros on the right, until the list contains n items.
- With CAMS, multiply all the factors back together.
Does the product reproduce the original coefficient list for the polynomial of degree n - 1 and thus the original circulant? The answer is yes. In other words the original circulant matrix has been factored into a product of n - 1 simpler circulant matrices.
An example. Consider the 3 by 3 circulant with generating vector v = 1 0 1. This list corresponds to the quadratic X*2 + 1 with two imaginary roots i and -i
The two linear factors are (X - i) and (X + i). Padding a single zero to each list, calculate (-i 1 0) CAMS (i 1 0). The result is 1 0 1, the original circulant. The linear lists are circulant matrices whose product is the original circulant matrix. The original matrix has been factored. Note: Caution, APLers who try this using APLIWIN must write the imaginary number i as 0j1 and -i as 0j¯1.
How about a tougher problem? Consider the quintic polynomial of an earlier example. Its coefficient list is the first row of the 6 by 6 circulant generated by v=(7 11 4 7 4 11). T. C. Chen’s SCARF function in his SCARFS workspace is a robust and accurate root finder. It gives generally good results, but the iterative interaction involved requires patience and some experience in making a correct decision. The process appears to be somewhat tedious and not very practical for circulants of any significant size, say 20 by 20.
For the quintic example with coefficients 7 11 4 7 4 11, Chen’s SCARF function returned these five roots:
one real ⍝ ¯0.59948564304709
imaginary pair ⍝ 0.62914846089673 +/-i×0.83550042957063
imaginary pair ⍝ ¯0.51122382119136 +/-i×0.84205466888059
The reader may wish to check the roots with dyadic ‘base’ (⊥). Use root ⊥ 7 11 4 7 4 11, where root is one of the five listed. Results are fuzzy zeros, close to, but not quite zero. They are generally in the range of 10×⎕ct, or less, for the APLIWIN system default for comparison tolerance. How about a different approach, say finding the greatest common divisor of two circulants? Progress might be easier because there exists a known algorithm to find the GCD of two polynomials. Further, with a polynomial division function available, this Euclidean algorithm for polynomials can be programmed with little difficulty. These facts unfortunately are not enough to provide a quick answer. The quest requires, among other things, a good understanding of the meaning of division in a commutative ring.
Pressing On¶
I began a search for the (possibly non-existent) Euclidean algorithm for circulants on the Internet. The hunt turned into quite an adventure on the ‘net’. The full story is rather long and best told as a separate tale. One result is obviously the present article. Another result is the existence of a follow-on to this report. Presently only a draft, the investigation expounds upon the meaning of division in a commutative ring, affirms the existence of a Euclidean algorithm for circulants and pursues a tangible implementation thereof in APL.
The original motivation for the present article was not to find, pose and solve puzzles that seem to violate common sense. The larger purpose was to investigate ways to factor a particular kind of matrix, a circulant. One path led to the subject of Euclidean algorithms and their existence in general and in particular for circulant matrices. Another path led to finding roots of polynomials using Chen’s SCARFS workspace.
At every turn, I ran into the need to learn more about modern abstract algebra and had to engage in a crash course of self-study of the subject. The topic seems ready-made for APL. I suggest there exists an opportunity for someone to author a book, perhaps titled "APL and Modern Abstract Algebra". I believe both APLers and algebraists would benefit. (The latter to the extent that they have access to or are even aware of APL.)
Recapitulation¶
These APL programs have been displayed and discussed: long division of two polynomials, function POLYDIV and its associated cover function POLY (functions not relisted). Long division of two polynomials, function PDIVBY, (not relisted).
Thomson algorithm, modified,
TAMP stands for (T)homson (A)lgorithm for (M)ultiplying (P)olynomials.
Circulant algorithm, CAMP:
CAMP stands for (C)irculant (A)lgorithm for (M)ultiplying (P)olynomials.
Circulant product algorithm, CAMS:
CAMS stands for (C)irculant (A)lgorithm for(M)ultiplying Circulant(S).
CAMP produces the product of two polynomials as a linear list of coefficients, lowest degree term to highest degree. The internal working of the algorithm is essentially just to multiply two circulants. Because of the effective use of zero padding inside the algorithm, no reduction modulo n occurs. The result can also be interpreted as the product of two circulant matrices, but this is not generally very useful because of the extensive zero padding used to avoid reduction modulo n. CAMP does not expect that in general the input lists of polynomial coefficients will be the same length. It will produce correct results for left and right argument lists of different lengths. Because circulants and polynomials commute, CAMP produces the same answer if the left and right arguments are exchanged, even if they are of different lengths.
CAMS produces the first row of the product of two circulant matrices. Such a list is equivalent to polynomial multiplication modulo the ideal generated by (X*n - 1). The exponent n is the order of the circulant matrices. CAMS right argument list can be longer than the left argument, but not vice-versa. Unequal input list lengths are generally not a good thing. Even though the function will zero pad the left argument to match the right’s length, the user may inadvertently mix up inputs and so violate conformability requirements for matrix multiplication.
CAMS guarantees commutativity of its left and right arguments only if both lists are of equal length. Both arguments are best input with the same length, say n, which is also the exponent in the expression X*n - 1. Otherwise n is the length of the right argument.
Polynomial Long Division has produced verification in a specific example of the isomorphism between the ring of circulants and the domain of polynomials over a field, modulo the ideal generated by (X*n - 1). Judicious use of zero padding provides the key for solving the puzzle posed at the beginning of this article. If two circulant divisors of zero are so padded with the use of "nothing", they will produce "something" from "nothing".
Considered as a polynomial with real or complex coefficients, one can factor the first row of a circulant matrix. Just determine the roots of the polynomial, then judiciously pad the linear or quadratic coefficient lists. Treat these lists as circulant matrices, multiply them all together and they will reproduce the original circulant. One needs a very accurate root finder such as Chen’s SCARFS workspace and even then the procedure is not practical for larger circulants. An alternative possibility would be to develop a Euclidean algorithm for two circulant matrices and so find a greatest common divisor circulant.
Note: All APL functions are designed to run under Windows 3.1 using Iverson Software’s APLIWIN interpreter. Because almost all expressions are in upper case, the functions can be easily converted to run under IBM’s APL2. The only "case changes" required are to convert all quad functions from lower to upper case. Thus ⎕io, ⎕nc, ⎕ct and ⎕ex need to become ⎕IO, ⎕NC, ⎕CT and ⎕EX, respectively, in functions POLYDIV and POLY.
Readers interested in experimenting with the functions can type them in from "scratch" with little difficulty. The author would like to be advised of any interesting results or aberrant behaviour.
Acknowledgement¶
I thank Gerhard Niklasch of the Mathematical Institute, Technische Universität, Munich, Germany for his encouraging internet dialogues, valuable criticisms and cogent comments. As a novice, I often had to struggle with one or another of the various concepts of modern abstract algebra. I particularly appreciated his patient forbearance and helpful guidance in these instances.