CMAT
CMAT is a matrix calculator program, written in C.
Calculations can be performed on matrices with complex rational coefficients
using exact arithmetic routines, as well as on matrices with elements mod p.
There are Unix, DOS and Windows XP versions.
To run under windows2000+ from a dos window, the user should first download a dos emulator from http://dosbox.sourceforge.net/.
The Unix, DOS and Windows XP versions, together with the C source, can be downloaded.
People using the UNIX version have to create a .cmatrc file in the
working directory, consisting of the two lines
- datpath=.
- cmat=.
(Also see cmat_bugfix.html for bug reports.)
There are three calculator programs within CMAT: CMATR, CMATCR
and CMATM.
CMATR works over the rationals, CMATCR works the field of complex
rationals and CMATM works over the field of p elements, where p is any prime
less than 216 = 65536.
The programs use multiple precision arithmetic routines based on those in
[Fla][342-357,175-185]. (See documentation.)
- Up to M0 (=30) objects of each type can be created and stored for use
in future sessions.
- (rational) numbers: R[0],...,R[M0-1];
- (rational) matrices: RM[0],...,RM[M0-1];
- polynomials (rational): PR[0],...,PR[M0-1];
- polynomial (rational) matrices: PRM[0],...,PRM[M0-1];
- (complex rational) numbers CR[0],...,CR[M0-1];
- (complex rational) matrices CRM[0],...,CRM[M0-1];
- polynomials (complex rational): PCR[0],...,PCR[M0-1];
- polynomial (complex rational) matrices: PCRM[0],...,PCRM[M0-1];
- (mod p) matrices: mM[0],...,mM[M0-1];
- polynomials (mod p): Pm[0],...,Pm[M0-1];
- polynomial (mod p) matrices: PmM[0],...,PmM[M0-1].
The program can handle integers of arbitrary size. Matrix size is
restricted only by the time taken to create the matrix in question and to
compute with it. Execution times are much less when calculating with matrices
mod p.
RATIONAL/COMPLEX/MOD P NUMERICAL ROUTINES AVAILABLE.
Conjugation, addition, subtraction, multiplication, exponentiation,
ratio and inverse. Also the m-th root of rational can be computed to a given
number of decimal places.
POLYNOMIAL ROUTINES AVAILABLE.
Addition, subtraction, multiplication, scalar-multiplication,
exponentiation, derivative, evaluation, quotient, remainder, greatest common
divisor, expressing the gcd as a linear combination of the polynomials involved,
square-free decomposition. The resultant of two polynomials over ℤ[x] and the
discriminant of a polynomial over ℤ[x] can be found.
There is modular exponentiation for polynomials mod p.
The Berlekamp matrix of a squarefree polynomial mod p. (Our Berlekamp
matrix is the transpose of the usual Berlekamp matrix.) Factorization of
polynomials with rational, complex rational or mod p coefficients into
irreducibles. Finding an irreducible polynomial of given degree mod p.
Finding the nth cyclotomic polynomial. (mod p)
MATRIX ROUTINES AVAILABLE.
- Addition, negation, complex-conjugate transpose, multiplication by a
scalar, matrix multiplication, exponentiation, Kronecker product, evaluating
matrix polynomials at a square matrix.
- The standard elementary row and column operations can be performed on all
types of matrices.
- The entries of a stored matrix can be changed, either individually, or by
row or column. A single row or column can be deleted.
- A single row or column vector can be selected from a stored matrix.
More generally, a submatrix of a stored matrix can be formed by selecting a
subfamily of rows or columns, respectively.
- Two matrices can be joined by rows or columns. The partitioned matrix
[A|I] can be formed from A. A matrix can be unwrapped into a column vector
(a process useful for finding the minimum polynomial of a matrix).
- Forming the matrix A-tI, for use in eigenvector calculations.
- Creating special matrices such as the identity matrix, elementary Jordan
matrices and companion matrices as well as band matrices.
- Creating matrices of scalars whose (i-j)th element is given by simple
arithmetical expressions such as i + j - i^2 + 3*j for the numerators
and denominators, may be created wholly within CMAT.
- Finding the inverse, adjoint, determinant, characteristic polynomial
and minimum polynomial of a matrix of scalars, P-1AP. The Cholesky
decomposition of a positive definite matrix is found: A=L*DL, where L is
upper unit triangular and D is a diagonal matrix. The output has the diagonal
elements of D on its diagonal and its above diagonal elements are those of L.
The inverse of a matrix M of integers mod m can be found, when m > 1 is
a positive integer which is coprime to the determinant of M.
- Finding the reduced row-echelon form, solving linear equations,
finding bases for the null-space N(A) and column-space C(A) of a matrix A.
These procedures are very useful in algorithms for finding the Jordan canonical
form, or more generally, the rational canonical form of a square matrix.
(See [paper].)
- Calculating dot product, length and using the Gram-Schmidt process
to find an orthogonal basis for the column-space of a matrix.
Finding the cross-product of the row vectors of an (n - 1) x n matrix.
- Finding the determinant and adjoint of a matrix with polynomial entries.
- Finding the Smith canonical form of a matrix B with polynomial entries.
(e.g. when applied to the matrix B = xI - A, where A is a matrix with
rational or complex rational entries, this gives the similarity invariants
of A, in particular we get the minimum polynomial of A as the invariant of
highest degree.)
We also get polynomial matrices P and Q, whose determinants are constant
polynomials, with the property that PBQ is in Smith canonical form.
- Elementary row/column operations for matrices with polynomial entries.
- Printing matrices of rationals to a file or the screen, accurate to a
given number of decimals.
- Sending data to the line-printer.
- A file of r (≤ M0-j) prepared matrices of the same type
(rational, floating point, complex rational or modular) can be entered into
array positions j,...j+r-1.
RUNNING CMAT.
After entering cmat, the user can now descend through various levels
by selecting from the displayed screen options. To ascend through the various
levels, type q when this is an option. When leaving the program, the user
is given the option of saving data for reuse in future sessions.
CALCULATING WITH CMAT.
To enter a rational number and store the number as R[0] for example,
select the appropriate menu and type n 0.
NOTE: Use CONTROL-H, not the backspace key, to backspace over characters.
Integers are entered as usual, while non-integer rationals are entered
with a slash as integer/positive integer: e.g. 1/2, while complex numbers
with non-zero imaginary part must end with an i: e.g. 1/2i represents (1/2)i.
No brackets are to be used when entering numbers. Rational numbers are
outputted in 'lowest terms' with positive denominator.
Some latitude is allowed as regards spacing. For example,
the expressions 1-i, 1- i, 1 - i represent the same complex rational, while
1 -i represents the two numbers 1 and -i.
Matrices are entered row by row, the end of a row being marked by a
carriage return. Spaces separate entries of a row. Having entered matrices
RM[0] and RM[1], the product RM[0] * RM[1] is sent to RM[2] (say) or RM[0],
RM[1]) by selecting the multiplication option t from the menu below and
typing t 0 1 2.
Rational Matrix Arithmetic
------------------------------------
a j k l : RM[j] + RM[k] -> RM[l]
t j k l : RM[j] * RM[k] -> RM[l]
s j k l : RM[j] - RM[k] -> RM[l]
m j k : -RM[j] -> RM[k]
f j k l : R[j](RM[k]) -> RM[l] (scalar multiple)
* j k : (RM[j])* -> RM[k] (transpose)
e n j k : RM[j] ^ n -> RM[k] (exponentiation)
z j k l : RM[j] - R[k]I -> RM[l]
p : print numbers and matrices
x : to enter data
q : to exit
------------------------------------
To discontinue entering a polynomial or matrix from the keyboard,
type q or any non-alphanumeric character when entering a coefficient. A
default polynomial or matrix is then stored.
To calculate the 10th power of RM[0], select the exponentiation option e
from the above menu, typing e 10 0 3 to send the output to RM[3].
The situation is somewhat simpler with modular arithmetic calculations. Here
no numbers are stored and non-negative numbers less than the relevant
modulus are entered directly, instead of being stored. For example to multiply
the matrix mM[0] by 5 (mod 7) and to store the result as mM[1], select the
appropriate menu with the menu line
f,j k l p :j * (mM[k]) -> mM[l](mod p)(scalar multiple)
and type f 5 0 1 7.
It is important to remember that when using the modular arithmetic section of
CMAT, operations should be carried out only on stored objects having the same
modulus.
To terminate input prematurely after entering the first letter of a menu
option, type q followed by RETURN.
Finally, to input a file of r (≤ M0-j) prepared matrices of the same type
(rational, floating point, complex rational or modular) into array positions
j,...j+r-1, the first line of the file should contain the number of matrices;
the next line consists of the row and column number of the first matrix,
followed by the respective rows of the matrix. For example, the following file
represents two matrices of rationals, a 3 x 3 matrix followed by a 2 x 2 matrix:
- /* rational matrix file */
- 2
- 3 3
- 2/3 5 -7/8
- 1/2 12 -5
- 7 6 4
- 2 2
- 1 0
- 0 1
Comments on algorithms used in CMAT
- The tricks of P. Henrici mentioned in [Buc][200-201]
are used to speed up addition and multiplication of rationals.
-
The Gauss-Jordan method is used to find the reduced row-echelon form and to
solve system of linear equations over all three fields.
- The adjoint, inverse, determinant and characteristic polynomial of a
matrix of rationals or complex rationals are found by the
Faddeeva-Leverrier method. (See [Fad][177-181].)
For matrices over ℤp, we use a modification of an algorithm due to Frobenius
(see [McW][168-175]) to find the characteristic
polynomial.
Over ℤp, the Gauss-Jordan algorithm is used to find the inverse and
determinant of a matrix. The adjoint is found using the fact that adj(A)=0 if
rank(A) ≤ n-2, rank(adj(A))=1 if rank(A)=n-1 and adj(A)=(det(A))A-1
otherwise.
-
The inverse of an integer matrix mod m is found by using the adjoint matrix
and multiplying mod m by the inverse of determinant mod m.
- The minimum polynomial mA(x) of an nxn matrix A is found by
finding the least positive integer m such that Am is expressible
as a linear combination of the matrices In,...,Am-1.
This is done by unwrapping the matrices In,...,Ar into column vectors for 1 ≤ r ≤ m and solving the equation
a[0]I+···+a[r-1]Ar-1 = Ar
as a system of n2 equations in r unknowns. The system will be inconsistent
for 1 ≤ r ≤ m, but will have a unique solution when r = m, giving the minimum
polynomial mA(x) = xm-a[m-1]xm-1-···-a[0].
- The Smith canonical form is found by the algorithm in [Har][111-113].
Faster algorithms exist (see [Lu2]).
-
Factorization of polynomials in ℤp[x] is accomplished using a method of
finding a non-trivial factor due to P. Camion (see [Cam]).
This is used in conjunction with the Berlekamp matrix of a polynomial. For
information on this matrix see [Knu][420-429].
(Our Berlekamp matrix is the transpose of Knuth's.)
- Factorization of a polynomial in ℤ[x] is accomplished using an algorithm
outlined in [Mus]. The algorithm uses the degree-set testing procedure which
sometimes uncovers irreducibility before more complicated tests are employed.
The Hensel-Zassenhaus quadratic and linear liftings are used to obtain a
factorization of f(x) modulo a high power of a prime.
The degree-set test also has the effect of reducing the number of test
divisions needed to determine possible factors over ℤ[x] after the liftings
are completed.
- Factorization of a polynomial in ℚ(i)[x] is accomplished using an
algorithm outlined in [Tra].
- A squarefree decomposition algorithm for polynomials due to D.Y.Y. Yun
and described in [Knu][631-632] is used.
- Irreducible polynomials of given degree (mod p) are constructed using
a probabilistic algorithm from [Lu2][145-149].
- The algorithm used for calculating the resultant of two polynomials
with integer coefficients is outlined in the article by R. Loos in
[Buc][130].
- The cyclotomic polynomial Φn(x) is calculated using an algorithm
from [Lu2][354-356], which in turn is based on an account
in [Lu1][58-63].
Factoring Φpn-1(x) over ℤp[x] gives the {φ(pn-1)}/n monic primitive
polynomials of degree n over ℤp.
- The m-th root of a rational number is calculated using the algorithm
from [Ma1].
REFERENCES
- [Buc] B. Buchberger, G.E. Collins and R. Loos (Editors)
Computer Algebra, Symbolic and Algebraic Computation.
1972, Springer-Verlag Wien, New York.
- [Cam] P. Camion. A Deterministic Algorithm for Factorizing Polynomials of Fq[x].
Annals of Discrete Mathematics 17 (1983) 149-157.
- [Fad] V.N. Faddeeva. Computational Methods in Linear Algebra.
1959, Dover, New York.
- [Fla] H. Flanders. Scientific Pascal.
1984, Reston Publishing Company, Reston, Virginia. (Second Edition Birkhäuser, 1995).
- [Har] B. Hartley and T.O. Hawkes. Rings, Modules and Linear Algebra.
1970, Chapman and Hall, London.
- [Knu] D.E. Knuth. The Art of Computer Programming, Volume 2.
Second Edition 1981, Addison-Wesley Publishing Company, Reading, Massachusetts.
- [Lu1] H. Lüneburg. Galoisfelder, Kreisteilungskörper und Schieberegisterfolgen.
1979, B.I. Wissenschaftsverlag, Mannheim/Wien/Zürich.
- [Lu2] H. Lüneburg. On the Rational Normal Form of Endomorphisms.
1987, B.I. Wissenschaftsverlag, Mannheim/Wien/Zürich.
- [Ma1] K.R. Matthews. Computing mth roots.
College Mathematics Journal 19 (1988) 174-176.
- [Ma2] K.R. Matthews. A Rational Canonical Form Algorithm.
Mathematica Bohemica 117 (1992) 315-324 (pdf file)
- [McW] W.A. McWorter, Jr. An Algorithm for the Characteristic Polynomial.
Mathematics Magazine 56 (1983) 168-175.
- [Mus] D.R. Musser. Multivariate Polynomial Factorization.
J.A.C.M. 22 (1975) 291-308.
- [Tra] B.M. Trager. Algebraic Factoring and Rational Function Integration.
Proceedings of the 1976 ACM Symposium on Symbolic and Algebraic Computation,
219-226.
ACKNOWLEDGMENT
I am indebted over the years for help in fixing bugs and improving program
design generously provided by Michael Forbes (in the early years), Peter Adams
and Sean Vickery.
Email
http://www.numbertheory.org/keith.html
Last modified 28th June 2014