A Quite Different New Primitive

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

APL, especiallly extended APL, now offers a large panoply of primitives. However, one is missing for a useful extension to binary algebra. In order to avoid the design of a new symbol, it seemed interesting to present a rationale for the use of an already-existing symbol, which, for the moment, has no monadic meaning. This symbol, known by all students even outside the APL world, is the unequal ≠ symbol which, in many other programming languages, corresponds to the logical function XOR with a couple of Boolean arguments. It will be shown here why the same symbol could be adopted as the one of matrix inverse when the argument is a modulo2 integer matrix. The essential reason comes from complex arithmetics, revisited with the help of the biternion theory.

From quaternions to biternions¶

Sir Rowan Hamilton, once, invented quaternions, i.e. hypercomplex numbers with four components, one real, and three imaginary ones on 3 orthogonal imaginary axes, which are also orthogonal with the real axis.

Very useful in geometry, especially for 3D transformations, quaternions have never been implemented as a functionality of APL.

But a distinction should be made between APL as a mathematical notation - what it was, when invented by Iverson prior to its first implementation by IBM as APL\360 (without complex arithmetics) - and APL as a programming language, which has always been - exactly like all other programming languages - a slave of the hardware.

Even if we would like to use the quaternion theory with a computer, and, again, not only under APL, the present-time limitations which come from the fuzz of floating-point arithmetic would lead to unefficient implementations : although FORTRAN allows - and even obliges - the user to declare REAL constants, no FORTRAN program has ever processed any real number during the 40 years of existence of that language ! Computers only know rational numbers. So, when a complex constant is used for any classical operation such as a product or a division (and even for add and substract in the general case), the result is a complex constant containing two real numbers, each with a 50% probability of truncation on the least-significant bit. For quaternions, the probability of truncation would be the same… on each component, i.e. 4 times what it is on a single rational floating-point constant.

All microprocessors truncate information at each arithmetic operation in the general case, by rounding, always either in excess (ceiling) or in default (floor); their behaviour depends on the way bit-complements are managed in their internal cuisine (according to a recipe of bit complementing either to 1 or to 2).

Such effects, already very sensitive when a function is iterated, would become catastrophic with quaternions, although the theory is one of the most beautiful of the history of mathematics. A computer is not a mathematical processor; it is, in all cases, a bit processor first.

So, it seemed adequate to look for a computer-oriented way of representing, in general, complex constants.

Binary algebra - as exposed in the workshop presented at APL94 - may be considered as the union of two different algebras : logical algebra and modulo2 integer algebra. All operations which are necessary in modulo2 integer algebra may be performed thanks to the use of the logical Boolean primitives, because both algebras are isomorphous; only the meaning of 0 and 1 differs.

If matrix inversion a priori does not make sense in logic, it does make sense in modulo2 integer algebra; hence the need for some notation.

But binary matrices (understood as modulo2 ones) have fantastic properties, in particular the ability of representing - exactly - some complex or hypercomplex constants which, otherwise, cannot be represented in computers without fuzz.

Already, a variety of papers has shown that a simple modulo2-invertible matrix with order 2 (in APL with size 2 2) - a biternion - is able to represent - exactly - with 4 bits only, the main symmetry operator of spacetime (for me the main constant of physics), the 3-fold rotation operator around an axis, i.e. constant \(\mathbf{j}\), expressed in Euler’s notation as \(e^{2i\pi/3}\).

The reader may easily check that for any modulo2 matrix with order 2, matrix inversion is performed - without any operation - exchanging the contents of the bits which form its main diagonal; (the proof is immediate, when one knows how a regular matrix with order 2 is inverted, due to the fact that modulo 2 integer algebra is an unsigned algebra, so that the modulo2 determinant of any modulo2 integer matrix is 1 when the matrix is invertible, and 0 otherwise).

If matrix \(\begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix}\) is chosen as \(\mathbf{j}\), then, matrix \(\begin{matrix} 0 & 1 \\ 1 & 1 \end{matrix}\) is \(\div\mathbf{j}\) and also \(\mathbf{j}^2\).

The hereabove-given symmetry rule for such matrix inversion also shows that matrices :

The reader may also check that, in coherent modulo2 algebra, the algebraic sum of the three cubic roots of 1, representated by the unit matrix and the \(\mathbf{j}\) and \(\mathbf{j}^2\) matrices, then itemwise-modulo2-summed, is a null matrix, as it has to be.

The algebraic modulo2 sum is given either by +/ taken modulo 2 on a vector of matrices, or, directly, by its logical counterpart ≠/ since scalar function different (unequal ) is isomorphous with either plus modulo2 or minus modulo2.

Also, the result of ≠/ on a vector of matrices with the same size, is the last item of ≠\ applied to the same vector of matrices, independently from the order (indices) of the matrices within the vector.

Now, why can we generalise the meaning of ≠ and propose its monadic syntax as an equivalent of the one of the monadic domino for the modulo2 matrix inversion of binary matrices (i.e. rank-2 cells according to Iverson's terminology) ?

For several reasons. One comes from extended APL, the second one from quaternions.

In extended APL, monadic + means "conjugate". Already, in both versions of the unextended ISO8485 standard (e for English, and f for French with word "conjugué"), the future complex extension was foreseen.

For "real" constants, +Z is just Z; for Z a complex (a,b), then +Z is (a,−b).

So, the modulo2 counterpart of +, which, when + symbolises addition, is indeed ≠ , could also mean modulo2 conjugate for the new binary complex and hypercomplex numbers, based on biternions.

For all modulo2 invertible matrices with order 2 (the 6 matrices listed hereabove), each matrix is the complex conjugate of its modulo2 matrix inverse. j is indeed the conjugate of j2 and conversely. And self-inverse matrices are also their own conjugates. Complex conjugacy is, geometrically, a line-symmetry, the real axis being the symmetry element of the operation, seen in the complex plane.

In modulo2 matrices with order 2, the corresponding symmetry element simply becomes the 2nd diagonal of the matrices; matrix inversion and complex conjugacy have become, modulo 2, at least for all these matrices, the same operation.

Now, we have to show that the same coalescence of meanings for two (apparently) different expressions of conventional mathematics, also holds for any order i.e. for modulo2 invertible square binary matrices with any size.

Let us do it slowly - just because, at the present time, modulo2 integer matrix algebra is taught at school nowhere in the world (so that no book can help yet).

Think that you stay where you are on Earth, but that - perhaps - somebody you do not know is standing, not quite upside down - probably in a boat - not at the antipodes (this word, in Greek, does mean "upside down"), but in a point, on the Earth's surface, with the same longitude (E or W of Greenwich - or Paris), but with opposed latitude i.e. S instead of N or conversely. This guy would be your conjugate - and conversely - on Earth. (The reference to Paris meridian, in fact not much in use anymore, shows that conjugacy is independent from the location of the reference axis chosen for zero-longitude). The location of the conjugate guy, however, will depend from the choice of the Equator, or, which is equivalent, from the choice of the polar axis, orthogonal with the equatorial plane.

Now, binary matrices with order 3 (size 3 3), have the same properties as the ones of the couple formed by you and your hypothetical conjugate guy.

If their modulo2 determinant is zero, they are not invertible. The square root of the absolute value of the determinant of any regular invertible matrix does play the role of the distance to the center of the sphere (or of the hypersphere for orders larger than 3), i.e. of the module of a 3-D hypercomplex number - a ternion, of a 4-D hypercomplex number - a quaternion - or a multi-D hypercomplex number - a multi(n)ion ? with n for euphony. Mathematical properties are more general than what one may think.

So, you could not live (except if you follow Jules Verne) at the center of the Earth with a null determinant.

Longitudes and latitudes have no meaning for a hot point which lies at the center of any such system, exactly as, in the (module, argument)-representation of a regular complex number, the number is 0, for any value of the argument when the module is already 0.

This may help to show that the hyperconjugate of any hypercomplex number is its symmetric on the same hypermeridian, the symmetry element being still the unique real axis, which is the polar axis.

On Earth (if supposed spherical, as it was for the first definition of the length unit named metre) the circumference is 40,000 km. Suppose a guy is sentenced to occupy on Earth, quantised locations only, each of them being separated by one meter, on the same meridian - including the opposed part, i.e. e.g. Greenwich meridian (0°) as well as the ideal "date-change" one (180° either E or W, or π). The guy may jump on the big circle, from one location to the next one; after 40,000,000 nice entrechats - think of Nureev - he will be back on the Covent Garden stage after some time. If two guys (conjugate-sentenced) are to jump synchroneously in an endless ballet, they may meet one another - what a collision - on the poles only. Now, if the number of allowed locations were odd, they could meet only at one pole, let us say the North pole.

This is exactly what happens if you try to solve equation \(z^{999} = 1\) which has 999 different solutions, one being - obviously- 1, the unique real solution, and the 998 remaining complex solutions being pairwise-conjugated on the trigonometric circle.

The North pole of the trigonometric representation of the complex plane is number 1 (or, as a matrix, a unit matrix). The South pole is number −1 (in matrix representation any self-inverse matrix : for regular matrices, the either vertical or horizontal mirror of the unit matrix - cfr. Pauli or Dirac matrices). In complex geometry, the polar axis is taken horizontally.

In 3-D complex geometry, this is still true : Imagine, this time Charles Chaplin, in the Great Dictator, playing with the globe : spinning around the poles - as the Earth does "spontaneously", - makes every possible (discrete) meridian (i.e. drawn on the globe) occupy every day, all along the day, e.g. the location Greenwich meridian had at midnight. Successively, all discrete meridians become, for a short while, a trigonometric circle. All points on the surface of an ideal planet are defined by their module (which is 1 by definition as the radius of the planet), and by two arguments (a latitude and a longitude, measured not from 0 to 90° N or S or from 0 to 180° E or W, but from 0 to 360° or \(2\pi\) for the latitude and from 0 to 180° for the longitude.

A guy and his conjugate guy are on the same meridian, with module 1 and the same longitude. Their latitudes are indeed opposed in a plus and minus system, or, if one is Alpha, the other one is \(2\pi\)-Alpha (or \(2k\pi\) -Alpha, with \(k\) any integer).

What is true for the complex roots of 1 also holds for binary matrices.

If any binary matrix (with module 1) is invertible, it will correspond to an allowed (discrete) location on some meridian of its hypersphere with latitude \(L\).

Squaring this matrix will double the latitude (exactly as squaring a rotation matrix doubles the rotation angle in the complex plane, without making any change to the module which is 1, and with no possibility of escaping from the complex plane); similarly the matrix will stay on the same meridian. Raising the same matrix to any (integral) power does sentence the result to stay on the same meridian, so, to keep - always - the same longitude.

But, for a given order, the number of bit configurations within a matrix is finite. If one raises (modulo 2) some binary matrix \(M\) with modulo2 determinant 1 to its successive powers, one will obtain immediately a unit matrix if the matrix is self-inverse (lay at the S pole), and a different matrix, as its square, if the matrix is invertible but not self-inverse. For small matrices (there are only 6 invertible matrices with order 2), the number of different configurations is very limited. Even for huge matrices, one of these days, one of the matrix powers - let us say power \(Q\), with \(Q\) an integer - will have to reproduce an already-known configuration of bits (e.g. already found for power \(P\), another integer).

Then, matrix \(M\) raised to power \(Q+1\) will reproduce \(M\) raised to power \(P+1\), etc…

Power \(Q+K\) reproduces power \(P+K\) for any integer \(K\), negative, positive, or 0.

When \(K\) equals \(-P\), this proves :

a) that power \(Q+K\) or \(Q-P\) of \(M\) is also power 0 of \(M\) then a unit matrix,

b) that powers \(n(Q-P)\) with \(n\) any integer will be unit matrices,

c) that, if \(p\) is \(Q-P\), then power \(p+1\) of \(M\) is \(M\).

d) that power \(p-1\) of \(M\) is such that if this matrix, multiplied by \(M\) has to produce a unit matrix, power \(p-1\) of \(M\) is also the modulo2 unique matrix inverse of \(M\) (by definition of an inverse matrix).

\(p\) is the power of the cyclic groups of all the different matrices which are matrix powers of \(M\), including \(M^{-1}\) the inverse matrix of \(M\), as well as \(M^0\) the unit matrix.

Of course, \(M^{-2}\) is the square of \(M^{-1}\) and also the matrix inverse of \(M^2\) the square of \(M\); and this property will hold for any power : the matrix inverse in modulo 2 integer algebra of any invertible matrix \(M\), is, for any order, up to infinity, the conjugate of \(M\) in the geometric interpretation with hyperspperes, of hypercomplex quantities.

In other words, if \(M\) is equivalent to a rotation matrix on any meridian with angle Alpha, its inverse matrix lies on the same meridian and corresponds to angle −Alpha.

And \(M\) is some exact equivalent of \(e^{2k\pi i/p}\) with \(p\) the order of its cyclic group i.e. the power to which \(M\) raised equals a unit matrix, and with \(k\) "some" integer.

The difference with classical rotation matrices in 2D comes from the fact that \(k\) is not defined. Just word "some" is used. It is a matter of taste to choose a value for \(k\), (modulo \(p\)).

If we fix \(k\), then the value of the arguments for all the other powers of \(M\) is defined, in particular for the matrix inverse. The fact that a matrix and its inverse are conjugate in independent from the choice of \(k\).

(A similar case arises with logarithms of complex numbers which are defined modulo \(2\pi i\), and for the choice of the definition domain for inverse trigonometric functions).

Among the pair of matrices the cube of them is 1 for order 2, we may choose either of both as the complex constant \(\mathbf{j}\) (the complex cubic root of 1 with positive imaginary part) ; then the other one corresponds to \(\mathbf{j}^2\) or \(\mathbf{j}^{-1}\).

Conclusion : If monadic + produces the complex conjugate of any complex number in APL, then its modulo 2 equivalent, monadic ≠ could produce the modulo2 generalised conjugate of any hypercomplex number which may be defined with modulo 2 integer matrices; then, indeed, monadic ≠ also corresponds to a perfect monadic domino (which will, this time, always return exact results, independently from the computer and from its floating-point arithmetic; even a bugged Pentium would then compute with no error).

Algorithm¶

The explanations given so far also lead to a general algorithm for binary matrix inversion, coming both from (cyclic) group theory and from graph theory.

However, there is no immediate need for studying group theory and graph theory for a reader who just knows some trigonometry :

The path (or ballerine jump) which leads, on a meridian from the North pole to any permitted location on the meridian, e.g. exactly where matrix \(M\) (supposed invertible) lies, is exactly the same, in length and direction¹ on the meridian, than the path which leads from the conjugate of \(M\) to the North pole.

Both paths, reversed, are the same in length and direction on the meridian.

So if we take matrix \(M\) and try to transform it, using a convenient algorithm (a recipe) into a unit matrix, the same algorithm will be able to transform, automatically, a unit matrix into the conjugate of \(M\), i.e. the matrix inverse of \(M\). Simple, isn't it ?

Now, any path may be divided into two sub-paths; this holds for algorithms : instead of a single function or program, we may chose to write two successive sub-functions or subroutines which will be invoked in a cover-function or in the main program, successively. I may fly from Paris to London with a connection in Brussels. I surely will choose this path if any ticket I may buy in Brussels costs nothing since the president of the Belgian company is my brother-in-law (unfortunately a stupid and wrong hypothesis, because I have no sister).

We shall not transform matrix \(M\) directly into a unit matrix - this costs too much thinking, looping or checking - , but into some other matrix \(P\) which will be easily invertible in all cases or, at least, which may be, in turn, easily transformed into a unit matrix.

If we apply, let us say to the first row of \(M\), a transformation which will preserve the first 1 in this row, and reset the other 1s to 0, we will obtain a new row with only one 1 (of course).

Of course, again, there is always a first 1 in the first row of \(M\); otherwise, this row would contain only 0s, so that \(M\) would have a 0 determinant and would not be invertible (which does not fit the initial hypothesis).

If there is only one 1 in the first row - this will happen -, there is nothing to do (the algorithm becomes simple).

Let \(N\) be the order of the matrix.

If the number of 1s is more than one in row R1, we have to set to 0 the extra 1-bits from R1, using (as usual now) ≠. We shall apply the same recipe to the first row of a unit matrix. In fact we simply XOR the first row of the unit matrix with R1 after having reset the first 1-bit of R1 to 0 (XORing it with 1).

And we shall perform the same operations to all other rows which contain a 1 with the same index as in R1, both in \(M\) and in the unit matrix.

Then, no other row in \(M\) may contain 1 in the same column as R1 does.

If we repeat the same operation on row R2, keeping this time the first 1-bit of R2 only (which, again, must exist, otherwise \(M\) would not have been invertible), and on rows which also have a 1 in the same column as R2, the new R2 will have just one 1-bit left in a different column than R1. The same rows are XORed with the same vector also in the former unit matrix.

TheXORing elementary algorithm is now applied to row R3 in order to keep the first 1-bit, and to the rows which contain the same first 1-bit, and to the same rows, with the same XOring mask in the former unit matrix,… etc… until the last row is reached and processed.

In the end, matrix \(M\) has become a matrix in which there is only one 1 in each row, all the \(N\) 1s being located in different columns. Such matrices have a name in mathematics : they are permutation matrices.

Then, if one knows that the matrix inverse of any permutation matrix is its transposed matrix, the inverting task is almost already over (without any second sub-program).

In fact we just have to reindex the columns of the former unit matrix according to the permutation defined by the permutation matrix : this makes transposition useless (it could take time to transpose a huge matrix), so we simply reindex the columns because we have worked with the rows; of course, transposing the dimensions from the beginning would lead to the same result : keeping only the first 1 successively in all columns of \(M\), and transforming each time with the same mask the columns of a unit matrix, produces a permutation matrix which reindices the rows of the former unit-matrix, thus transformed, producing the desired modulo2 matrix-inverse or modulo hypercomplex conjugate…

As arithmetic operations, the algorithmic definition of matrix inversion or conjugation, corresponding to the monadic ≠ symbol, will have used only… XORs corresponding to classical ≠ .

It is hoped that this paper may bring some new some ideas to implementers, because APL could become an uncomparable tool also for exact hypercomplex geometry in modulo2/logic algebra with plenty of applications (even if the reader cannot imagine them still now).

MAGIC was solved using modulo2 cyclic groups - already almost 3 years ago - in half a line of APL. Many problems may be in fact solved by modulo2 matrix inversion or matrix powers, once one knows how to formulate them in this fantastic algebra (it is not still obvious in most cases…); this sentence is not just a conjecture, but a theorem, which may make, in the future, the usage of computers more fiable and accurate than ever before. The consequence of a theorem may be a new conjecture. If this conjecture is true, this will have happened thanks to APL.

Here is a possible version of the modulo2 matrix-inversion algorithm which results from the « hypersphere meridian » method :

    ∇ M←INV M;D;I;J;K;L;N;V;T;e
[1]   I←⎕io≠1 ⋄ L←⍳D←1↑S←⍴M ⋄ K←⍴I↓V←D↑1 ⋄ T←S⍴V,N←0
[2]   e←'J←V≠M[;I⍳1]⋄ J←⍴I←J/L⋄ M[I;]←M[I;]≠M[J←J/N;]⋄ T[I;]←T[I;]≠T[J;]'
[3]   'J←1∊I←,M[N←V/L;] ⋄ ⍎J/e ⋄ V←0,¯1↓V' wh 'K≠N' ⋄ V←⍋M ⋄ M←T[V;]
    ∇

(Other versions appeared e.g. in « The Power of Boolean Computing », by G. A. Langlet, SEAS94 Spring Meeting, La Hulpe, in APL94 Companion, in Les Nouvelles d’APL, in APL-CAM J.).

Note that this version (like the other ones) also « general-inverts » non-invertible square matrices with a modulo2-determinant equalling 0; (such a general inverse of a square null matrix is a conforming unit matrix, by definition). Beware of the fact that different versions may produce different general inverses; however, all versions produce the same unique inverse for a modulo2 invertible matrix.

So, one never gets DOMAIN ERROR.

For the direct solution of a modulo2 linear system, see also function SEG in the special issue for APL95 of Les Nouvelles d’APL N° 15 (May 1995), paper « In XOR we trust », by Ilya Mironov.

≠