Groups, Algorithms, Programming
www.gap-system.org
GAP website: www.gap-system.org
Original Document: http://www-circa.mcs.st-and.ac.uk/gapfinite0.php
http://www.math.colostate.edu/~hulpke
We begin by looking at ways to define a vector.
v1:= [ 1/2, 3/4, -2/3, 22/7 ];
[ 1/2, 3/4, -2/3, 22/7 ]
We can see if what we have defined is a vector
IsVector( v1 );
true
IsList( v1 );
true
We define a vector space over the rational numbers
Suppose we want a three dimensional subspace of R4
v2:= [ 1, 3, 2, 4 ];
[ 1, 3, 2, 4 ]
v3:= [ 1/2, 1/4, 1/3, 3/4 ];
[ 1/2, 1/4, 1/3, 3/4 ]
V:= VectorSpace( Rationals, [ v1, v2, v3 ] );
<vector space over Rationals, with 3 generators>
It is easy to test if vectors lie in the vector space $V$
[ 1, 1, 1, 1 ] in V;
false
[ 28, 70, 84, 45 ] in V;
true
Scalar multiplication of vectors is easy
1/2*v1;
[ 1/4, 3/8, -1/3, 11/7 ]
As is the dot product
v1*v2;
1175/84
Of course linear combinations of vectors from V lie in V. We try an example
1/2*v1 - 3/4*v2 + 5/4*v3 in V;
true
Let us construct the 3 dimensional vector space over GF(3)
V:= FullRowSpace( GF( 3 ), 3 );
( GF(3)^3 )
Let us check whether a vector is in V
[ 1, 1, 1 ] in V;
false
This gives false since 1 is not in GF(3). We need to use the identity of GF(3)
o:= One( GF( 3 ) );
Z(3)^0
[ 1, 1, 1 ]*o in V;
true
Let us look at the sum of all the vectors of V
s:= [ 0, 0, 0 ]*o;; for v in V do s:= s + v; od; Print(s, "\n");
[ 0*Z(3), 0*Z(3), 0*Z(3) ]
It is easy to see why one gets 0 since for each v in V, -v is also in V. What about a field of characteristic 2?
V:= FullRowSpace( GF( 4 ), 4 );; o:= One( GF( 4 ) );; s:= [ 0, 0, 0, 0 ]*o;; for v in V do s:= s + v; od; Print(s, "\n");
V:= FullRowSpace( GF( 4 ), 4 );; o:= One( GF( 4 ) );; s:= [ 0, 0, 0, 0 ]*o;; for v in V do s:= s + v; od; Print(s, "\n");
[ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ]
So the answer is still 0. This is not quite so obvious - what about GF(2) with dimension 1?
We set up a vector space over the rationals spanned by 4 vectors
v1:= [ 2, 2, 1, 3 ];; v2:= [ 7, 5, 5, 5 ];; v3:= [ 3, 2, 2, 1 ];; v4:= [ 2, 1, 2, 1 ];;
V:= VectorSpace( Rationals, [ v1, v2, v3, v4 ] );
<vector space over Rationals, with 4 generators>
We compute its basis
B:= Basis( V );
SemiEchelonBasis( <vector space over Rationals, with 4 generators>, ... )
What is the dimension of V ?
Length( B );
3
Let us see the basis vectors GAP has computed
BasisVectors( B );
[ [ 1, 1, 1/2, 3/2 ], [ 0, 1, -3/4, 11/4 ], [ 0, 0, 1, 3 ] ]
We can express vectors in V as a linear combination of the basis vectors B
Coefficients( B, [2,1,2,1]);
[ 2, -1, 1/4 ]
Since the coefficients are [ 2, -1, 1/4 ] we should have
[ 2, 1, 2, 1 ] = 2*B[1] - 1*B[2] + 1/4*B[3];
true
Or it can be varified by
list:=[ [ 1, 1, 1/2, 3/2 ], [ 0, 1, -3/4, 11/4 ], [ 0, 0, 1, 3 ] ];; coef:=[ 2, -1, 1/4 ];; coef*list;
[ 2, 1, 2, 1 ]
What if we try to express a vector not in the space as a linear combination?
Coefficients( B, [1, 0, 0, 0]);
fail
Not surprisingly we get "fail".
We can take linear combinations of the vectors of B with the LinearCombination command. For example
LinearCombination( B, [ 1/2, 1/3, 1/4 ] );
[ 1/2, 5/6, 1/4, 29/12 ]
produces the same result as
1/2*B[1] + 1/3*B[2] + 1/4*B[3];
[ 1/2, 5/6, 1/4, 29/12 ]
We next look at ways to define a matrix.
The following command defines a matrix as a list of lists.
m1:= [ [ 1, 2, 3 ], [ 2, 3, -1 ], [ 1, -2, 5] ];
[ [ 1, 2, 3 ], [ 2, 3, -1 ], [ 1, -2, 5 ] ]
We can display this on the screen as a more standard looking matrix.
Display( m1 );
[ [ 1, 2, 3 ], [ 2, 3, -1 ], [ 1, -2, 5 ] ]
To print the (2, 3)-entry of m1 we use
m1[2][3];
-1
To change an entry use
m1[2][2]:=300;
300
Display( m1 );
[ [ 1, 2, 3 ], [ 2, 300, -1 ], [ 1, -2, 5 ] ]
Next we define a matrix over GF(5)
We let o be the identity of GF(5) and z be the zero of GF(5)
o:= One( GF( 5 ) );
Z(5)^0
z:= Zero( GF( 5 ) );
0*Z(5)
m2:= [ [ o, z, z ], [ 2*o, 3*o, z ], [ z, o, 4*o ] ];
[ [ Z(5)^0, 0*Z(5), 0*Z(5) ], [ Z(5), Z(5)^3, 0*Z(5) ], [ 0*Z(5), Z(5)^0, Z(5)^2 ] ]
Now Display gives a more readable form
Display( m2 );
1 . . 2 3 . . 1 4
We can easily define matrices whose ( i, j ) th entry is given by a function of i and j.
m3:= List( [ 1, 2, 3 ], i -> List( [ 1, 2, 3 ], j -> i*j ) );
[ [ 1, 2, 3 ], [ 2, 4, 6 ], [ 3, 6, 9 ] ]
Display(m3);
[ [ 1, 2, 3 ], [ 2, 4, 6 ], [ 3, 6, 9 ] ]
Of course it is easy to define large matrices this way
m4:= List( [ 1 .. 12 ], i -> List( [ 1 .. 12 ], j -> i*j ) );
[ [ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 ], [ 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24 ], [ 3, 6, 9, 12, 15, 18, 21, 24, 27, 30, 33, 36 ], [ 4, 8, 12, 16, 20, 24, 28, 32, 36, 40, 44, 48 ], [ 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60 ], [ 6, 12, 18, 24, 30, 36, 42, 48, 54, 60, 66, 72 ], [ 7, 14, 21, 28, 35, 42, 49, 56, 63, 70, 77, 84 ], [ 8, 16, 24, 32, 40, 48, 56, 64, 72, 80, 88, 96 ], [ 9, 18, 27, 36, 45, 54, 63, 72, 81, 90, 99, 108 ], [ 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120 ], [ 11, 22, 33, 44, 55, 66, 77, 88, 99, 110, 121, 132 ], [ 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144 ] ]
Display( m4 );
[ [ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 ], [ 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24 ], [ 3, 6, 9, 12, 15, 18, 21, 24, 27, 30, 33, 36 ], [ 4, 8, 12, 16, 20, 24, 28, 32, 36, 40, 44, 48 ], [ 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60 ], [ 6, 12, 18, 24, 30, 36, 42, 48, 54, 60, 66, 72 ], [ 7, 14, 21, 28, 35, 42, 49, 56, 63, 70, 77, 84 ], [ 8, 16, 24, 32, 40, 48, 56, 64, 72, 80, 88, 96 ], [ 9, 18, 27, 36, 45, 54, 63, 72, 81, 90, 99, 108 ], [ 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120 ], [ 11, 22, 33, 44, 55, 66, 77, 88, 99, 110, 121, 132 ], [ 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144 ] ]
Defining the same matrix, but this time over Z(5) can be done by either
m5:= List( [ 1 .. 12 ], i -> List( [ 1 .. 12 ], j -> i*j*o ) );
[ [ Z(5)^0, Z(5), Z(5)^3, Z(5)^2, 0*Z(5), Z(5)^0, Z(5), Z(5)^3, Z(5)^2, 0*Z(5), Z(5)^0, Z(5) ], [ Z(5), Z(5)^2, Z(5)^0, Z(5)^3, 0*Z(5), Z(5), Z(5)^2, Z(5)^0, Z(5)^3, 0*Z(5), Z(5), Z(5)^2 ], [ Z(5)^3, Z(5)^0, Z(5)^2, Z(5), 0*Z(5), Z(5)^3, Z(5)^0, Z(5)^2, Z(5), 0*Z(5), Z(5)^3, Z(5)^0 ], [ Z(5)^2, Z(5)^3, Z(5), Z(5)^0, 0*Z(5), Z(5)^2, Z(5)^3, Z(5), Z(5)^0, 0*Z(5), Z(5)^2, Z(5)^3 ], [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ], [ Z(5)^0, Z(5), Z(5)^3, Z(5)^2, 0*Z(5), Z(5)^0, Z(5), Z(5)^3, Z(5)^2, 0*Z(5), Z(5)^0, Z(5) ], [ Z(5), Z(5)^2, Z(5)^0, Z(5)^3, 0*Z(5), Z(5), Z(5)^2, Z(5)^0, Z(5)^3, 0*Z(5), Z(5), Z(5)^2 ], [ Z(5)^3, Z(5)^0, Z(5)^2, Z(5), 0*Z(5), Z(5)^3, Z(5)^0, Z(5)^2, Z(5), 0*Z(5), Z(5)^3, Z(5)^0 ], [ Z(5)^2, Z(5)^3, Z(5), Z(5)^0, 0*Z(5), Z(5)^2, Z(5)^3, Z(5), Z(5)^0, 0*Z(5), Z(5)^2, Z(5)^3 ], [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ], [ Z(5)^0, Z(5), Z(5)^3, Z(5)^2, 0*Z(5), Z(5)^0, Z(5), Z(5)^3, Z(5)^2, 0*Z(5), Z(5)^0, Z(5) ], [ Z(5), Z(5)^2, Z(5)^0, Z(5)^3, 0*Z(5), Z(5), Z(5)^2, Z(5)^0, Z(5)^3, 0*Z(5), Z(5), Z(5)^2 ] ]
or
m5:= List( [1 .. 12], i -> List( [1 .. 12], j -> i*j ) )*o;
[ [ Z(5)^0, Z(5), Z(5)^3, Z(5)^2, 0*Z(5), Z(5)^0, Z(5), Z(5)^3, Z(5)^2, 0*Z(5), Z(5)^0, Z(5) ], [ Z(5), Z(5)^2, Z(5)^0, Z(5)^3, 0*Z(5), Z(5), Z(5)^2, Z(5)^0, Z(5)^3, 0*Z(5), Z(5), Z(5)^2 ], [ Z(5)^3, Z(5)^0, Z(5)^2, Z(5), 0*Z(5), Z(5)^3, Z(5)^0, Z(5)^2, Z(5), 0*Z(5), Z(5)^3, Z(5)^0 ], [ Z(5)^2, Z(5)^3, Z(5), Z(5)^0, 0*Z(5), Z(5)^2, Z(5)^3, Z(5), Z(5)^0, 0*Z(5), Z(5)^2, Z(5)^3 ], [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ], [ Z(5)^0, Z(5), Z(5)^3, Z(5)^2, 0*Z(5), Z(5)^0, Z(5), Z(5)^3, Z(5)^2, 0*Z(5), Z(5)^0, Z(5) ], [ Z(5), Z(5)^2, Z(5)^0, Z(5)^3, 0*Z(5), Z(5), Z(5)^2, Z(5)^0, Z(5)^3, 0*Z(5), Z(5), Z(5)^2 ], [ Z(5)^3, Z(5)^0, Z(5)^2, Z(5), 0*Z(5), Z(5)^3, Z(5)^0, Z(5)^2, Z(5), 0*Z(5), Z(5)^3, Z(5)^0 ], [ Z(5)^2, Z(5)^3, Z(5), Z(5)^0, 0*Z(5), Z(5)^2, Z(5)^3, Z(5), Z(5)^0, 0*Z(5), Z(5)^2, Z(5)^3 ], [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ], [ Z(5)^0, Z(5), Z(5)^3, Z(5)^2, 0*Z(5), Z(5)^0, Z(5), Z(5)^3, Z(5)^2, 0*Z(5), Z(5)^0, Z(5) ], [ Z(5), Z(5)^2, Z(5)^0, Z(5)^3, 0*Z(5), Z(5), Z(5)^2, Z(5)^0, Z(5)^3, 0*Z(5), Z(5), Z(5)^2 ] ]
Try Display with these matrices.
Display(m5);
1 2 3 4 . 1 2 3 4 . 1 2 2 4 1 3 . 2 4 1 3 . 2 4 3 1 4 2 . 3 1 4 2 . 3 1 4 3 2 1 . 4 3 2 1 . 4 3 . . . . . . . . . . . . 1 2 3 4 . 1 2 3 4 . 1 2 2 4 1 3 . 2 4 1 3 . 2 4 3 1 4 2 . 3 1 4 2 . 3 1 4 3 2 1 . 4 3 2 1 . 4 3 . . . . . . . . . . . . 1 2 3 4 . 1 2 3 4 . 1 2 2 4 1 3 . 2 4 1 3 . 2 4
The following need little explanation
The identity matrix must be square.
id:= IdentityMat( 5 );; Display( id );
[ [ 1, 0, 0, 0, 0 ], [ 0, 1, 0, 0, 0 ], [ 0, 0, 1, 0, 0 ], [ 0, 0, 0, 1, 0 ], [ 0, 0, 0, 0, 1 ] ]
The zero matrix can be rectangular
zm:= NullMat(6, 8);; Display( zm );
[ [ 0, 0, 0, 0, 0, 0, 0, 0 ], [ 0, 0, 0, 0, 0, 0, 0, 0 ], [ 0, 0, 0, 0, 0, 0, 0, 0 ], [ 0, 0, 0, 0, 0, 0, 0, 0 ], [ 0, 0, 0, 0, 0, 0, 0, 0 ], [ 0, 0, 0, 0, 0, 0, 0, 0 ] ]
zm:= NullMat(6, 7, GF(2));; Display( zm );
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A diagonal matrix is square. We give it the diagonal entries.
dm:= DiagonalMat( [ 1, 1, 2, 2, 3, 4 ] );; Display( dm );
[ [ 1, 0, 0, 0, 0, 0 ], [ 0, 1, 0, 0, 0, 0 ], [ 0, 0, 2, 0, 0, 0 ], [ 0, 0, 0, 2, 0, 0 ], [ 0, 0, 0, 0, 3, 0 ], [ 0, 0, 0, 0, 0, 4 ] ]
We can define a similar matrix over a finite field.
o:= One(GF( 3 ));; dm:= DiagonalMat( [ 1, 1, 2, 2, 3, 4 ]*o );; Display( dm );
1 . . . . . . 1 . . . . . . 2 . . . . . . 2 . . . . . . . . . . . . . 1
GAP will produce a random matrix whose entries are integers
in the range [-10..10] with a binomial distribution.
rm:= RandomMat(10, 10);; Display( rm );
[ [ 1, -4, 1, 2, -3, 2, 5, 1, 0, -2 ], [ 4, 3, 5, 3, 1, -2, 1, -1, -2, 4 ], [ -1, -3, -1, 1, 2, 0, 1, 1, -1, 1 ], [ 2, -2, -1, -4, -2, 1, -1, 1, -2, -3 ], [ -5, -4, -3, 2, 0, -4, 0, 3, 4, -1 ], [ -1, 1, 1, 2, 4, 0, -2, 2, 0, 4 ], [ 3, -1, -3, -1, -3, 7, 0, 1, -1, 2 ], [ 0, 0, 2, 0, 5, -1, 0, 1, 0, 4 ], [ -1, 0, 3, -5, 4, 0, -1, -1, 0, 1 ], [ -1, 2, 1, 3, 0, 1, 4, 0, 2, -3 ] ]
Of course, running this code again will produce a different matrix.
GAP will also produce a random matrix with rational entries -- useful for testing other routines.
Display(RandomMat(10, 10,Rationals));
[ [ 1/2, 0, -2/5, 1, -3/2, 0, 0, -1, 1, 2 ], [ -1, 1/3, 1, -3/5, -1, 2, -3, 2, 1/5, 2 ], [ -1/6, -3/2, 0, -3, 1, -1/3, -1/4, 0, -4, 0 ], [ -1/2, 1, 1, 3/2, 0, 3, 0, 1, -1, -2/3 ], [ -1, -2/3, 0, 2, -1, -1/2, -3/2, -5/4, 2, 2 ], [ 1, -3, -1/4, -1, 1, -1, 0, 2, 2/5, 0 ], [ 4, -3, -2, 2, 1/4, -1, -2, 2, -1/2, 2/5 ], [ -2, -2, -2, -1/3, -1/3, 2, -1/5, -1/2, -1, 0 ], [ -2, 1, 1/2, 0, -2/3, -1/2, 0, -1/3, 3/2, 1/3 ], [ 0, -2, 0, 2, -3, -3, 3/4, 1, 1, 3/2 ] ]
If GAP is asked to produce a random matrix with entries in a finite field, the entries will be distributed uniformly over the field elements.
rm1:= RandomMat(10, 10, GF(5));; Display( rm1 );
. 1 . 4 1 4 . 1 3 2 4 1 2 . 3 . 3 3 1 3 3 3 1 4 2 2 2 4 2 . 4 2 2 3 3 4 . 2 . 4 2 3 1 1 2 1 4 1 1 3 1 1 . 2 1 2 2 3 2 . 4 1 . 1 4 1 1 2 1 3 3 3 3 4 4 3 3 3 4 3 4 4 . 1 4 4 2 2 1 4 . 2 2 2 2 4 1 4 4 4
Addition, multiplication, etc. for matrices is straightforward. We give examples
m:= [1..8];; m[1]:= [ [1,2,3], [2,3,-1], [1,-2,5] ];; m[2]:= [ [-1,4,-3], [1,2,-1], [-1,-2,3] ];; m[3]:= m[1] + m[2];; m[4]:= m[1]*m[2];; m[5]:= 3*m[1] - 7*m[2];; m[6]:= m[1]^3;; m[7]:= m[1]^(-1);; m[8]:= m[6]*(5*m[5] + 6*m[2]*m[4])^(-1);;
for i in [1..8] do
Display(m[i]);
Print("\n");
od;
[ [ 1, 2, 3 ], [ 2, 3, -1 ], [ 1, -2, 5 ] ] [ [ -1, 4, -3 ], [ 1, 2, -1 ], [ -1, -2, 3 ] ] [ [ 0, 6, 0 ], [ 3, 5, -2 ], [ 0, -4, 8 ] ] [ [ -2, 2, 4 ], [ 2, 16, -12 ], [ -8, -10, 14 ] ] [ [ 10, -22, 30 ], [ -1, -5, 4 ], [ 10, 8, -6 ] ] [ [ 28, -10, 102 ], [ 35, 63, -4 ], [ 4, -98, 170 ] ] [ [ -13/30, 8/15, 11/30 ], [ 11/30, -1/15, -7/30 ], [ 7/30, -2/15, 1/30 ] ] [ [ 96972/194635, 288222/194635, 330503/194635 ], [ 165884/583905, 186643/194635, 495226/583905 ], [ 274034/583905, 220138/194635, 977281/583905 ] ]
To find the transpose of a matrix use TransposedMat.
m:=[ [1,2,3], [2,3,-1], [1,-2,5] ];; mdash:=TransposedMat(m);
[ [ 1, 2, 1 ], [ 2, 3, -2 ], [ 3, -1, 5 ] ]
Display(m);
[ [ 1, 2, 3 ], [ 2, 3, -1 ], [ 1, -2, 5 ] ]
Display(mdash);
[ [ 1, 2, 1 ], [ 2, 3, -2 ], [ 3, -1, 5 ] ]
We use the Determinant function in the obvious way
m1:= [ [1,2,3], [2,3,-1], [1,-2,5] ];; m2:= [ [-1,4,-3], [1,2,-1], [-1,-2,3] ];; Determinant( m1 );
-30
Determinant( m2 );
-12
Of course the determinant of a product is the product of the determinants
Determinant(m1*m2);
360
If we look at matrices over GF(2) then their determinants will be either 0 or 1. What proportion would one expect for each with $2\times 2$ matrices? Let us examine the situation.
o:= One( GF(2) );;
z:= Zero( GF(2) );;
countOne:= 0;;
countZero:= 0;;
for i in [1..1000] do
mat:= RandomMat(2, 2, GF(2));
if Determinant(mat) = o then
countOne:= countOne + 1;
else
countZero:= countZero + 1;
fi;
od;
Print("countOne=",countOne, ", countZero=", countZero, "\n");
countOne=342, countZero=658
countOne=379, countZero=621
countOne=383, countZero=617
Obviously more matrices have determinant 0 than have determinant 1. Can you see why this is so?
Show that the expected number with determinant 0 is 5/8 of the total, so 625 in this case.
What happens with $3\times 3$ matrices, $4\times 4$ matrices, $5\times 5$ matrices?
countOne:= 0;;
countZero:= 0;;
for i in [1..1000] do
mat:=RandomMat(5, 5, GF(2));
if Determinant(mat) = o then
countOne:= countOne + 1;
else
countZero:= countZero + 1;
fi;
od;
Print("countOne=",countOne, ", countZero=", countZero, "\n");
countOne=290, countZero=710
countOne=299, countZero=701
countOne=303, countZero=697
Before leaving this problem let us examine the case of matrices over GF(3), where of course there are three possible determinants
o:= One( GF(3) );;
z:= Zero( GF(3) );;
countTwo:= 0;;
countOne:= 0;;
countZero:= 0;;
for i in [1..1000] do
mat:= RandomMat(2, 2, GF(3));
if Determinant(mat) = o then
countOne:= countOne+1;
elif Determinant(mat) = z then
countZero:= countZero + 1;
else
countTwo:= countTwo + 1;
fi;
od;
Print("countTwo=",countTwo, ", countOne=",countOne, ", countZero=", countZero, "\n");
countTwo=277, countOne=313, countZero=410 countTwo=298, countOne=322, countZero=380 countTwo=298, countOne=291, countZero=411 countTwo=286, countOne=295, countZero=419 countTwo=296, countOne=304, countZero=400
To compute the eigenvalues and eigenvectors of a matrix we need to specify the field
m:= [1..3];;
m[1]:= [ [1,2,3], [2,3,-1], [1,-2,5] ];;
m[2]:= [ [-1,4,-3], [1,2,-1], [-1,-2,3] ];;
m[3]:= [ [-2,-3,-3], [-1,0,-1], [5,5,6] ];;
for i in [1,2,3] do
Print("m[",i,"]=\n");
Display(m[i]);
Print("\n");
Print("Eigenvalues of m[",i,"]: ",Eigenvalues(Rationals, m[i]), "\n");
Print("Eigenvectors of m[",i,"]: ",Eigenvectors(Rationals, m[i]), "\n");
Print("CharacPloyns of m[",i,"]: ",CharacteristicPolynomial(m[i]), "\n");
Print("\n");
od;
[ 5 ] [ [ 1, 3, -2 ] ] 30+14*x_1-9*x_1^2+x_1^3 [ ] [ ] 12-8*x_1-4*x_1^2+x_1^3 [ 2, 1 ] [ [ 1, 1, 1 ], [ 1, 0, 3/5 ], [ 0, 1, 1/5 ] ] -2+5*x_1-4*x_1^2+x_1^3
Remark: The other two eigenvalues of m[1] are $2\pm\sqrt{10}$.
m[1]= [ [ 1, 2, 3 ], [ 2, 3, -1 ], [ 1, -2, 5 ] ] Eigenvalues of m[1]: [ 5 ] Eigenvectors of m[1]: [ [ 1, 3, -2 ] ] CharacPloyns of m[1]: x_1^3-9*x_1^2+14*x_1+30 m[2]= [ [ -1, 4, -3 ], [ 1, 2, -1 ], [ -1, -2, 3 ] ] Eigenvalues of m[2]: [ ] Eigenvectors of m[2]: [ ] CharacPloyns of m[2]: x_1^3-4*x_1^2-8*x_1+12 m[3]= [ [ -2, -3, -3 ], [ -1, 0, -1 ], [ 5, 5, 6 ] ] Eigenvalues of m[3]: [ 2, 1 ] Eigenvectors of m[3]: [ [ 1, 1, 1 ], [ 1, 0, 3/5 ], [ 0, 1, 1/5 ] ] CharacPloyns of m[3]: x_1^3-4*x_1^2+5*x_1-2
Let us do the same calculation over GF(7)
o:= One( GF(7) );;
m:= [1..3];;
m[1]:= [ [1,2,3], [2,3,-1], [1,-2,5] ]*o;;
m[2]:= [ [-1,4,-3], [1,2,-1], [-1,-2,3] ]*o;;
m[3]:= [ [-2,-3,-3], [-1,0,-1], [5,5,6] ]*o;;
for i in [1, 2, 3] do
Print("m[",i,"]=\n");
Display(m[i]);
Print("\n");
Print("Eigenvalues of m[",i,"]: ",Eigenvalues(GF(7), m[i]), "\n");
Print("Eigenvectors of m[",i,"]: ",Eigenvectors(GF(7), m[i]), "\n");
Print("CharacPloyns of m[",i,"]: ",CharacteristicPolynomial(m[i]), "\n");
Print("\n");
od;
m[1]= 1 2 3 2 3 6 1 5 5 Eigenvalues of m[1]: [ Z(7)^5 ] Eigenvectors of m[1]: [ [ Z(7)^0, Z(7), Z(7)^5 ] ] CharacPloyns of m[1]: x_1^3+Z(7)^5*x_1^2+Z(7)^2 m[2]= 6 4 4 1 2 6 6 5 3 Eigenvalues of m[2]: [ Z(7) ] Eigenvectors of m[2]: [ [ Z(7)^0, Z(7)^4, 0*Z(7) ] ] CharacPloyns of m[2]: x_1^3+Z(7)*x_1^2-x_1+Z(7)^5 m[3]= 5 4 4 6 . 6 5 5 6 Eigenvalues of m[3]: [ Z(7)^0, Z(7)^2 ] Eigenvectors of m[3]: [ [ Z(7)^0, 0*Z(7), Z(7)^2 ], [ 0*Z(7), Z(7)^0, Z(7) ], [ Z(7)^0, Z(7)^0, Z(7)^0 ] ] CharacPloyns of m[3]: x_1^3+Z(7)*x_1^2+Z(7)^5*x_1+Z(7)^5
If we are given a system of linear equations $XA=B$, where $B$ is a row vector, then we can find a solution to the equations (of course only if one exists).
A:= [[1,2,1],[1,-1,2],[1,2,1]];
[ [ 1, 2, 1 ], [ 1, -1, 2 ], [ 1, 2, 1 ] ]
B:=[1,1,4/3];
[ 1, 1, 4/3 ]
C:=SolutionMat(A,B);
[ 2/3, 1/3, 0 ]
C*A=B;
true
Now this only gives one solution, even when the system of equations has infinitely many solutions. The general solution is then SolutionMat(A, B) + v for any v in NullspaceMat(A). Let us check that fact for our particular example.
NullspaceMat(A);
[ [ -1, 0, 1 ] ]
Since the null space has dimension 1 the general solution will be the particular solution SolutionMat(A, B) plus any multiple of the basis vector for NullspaceMat(A). For example
M:= SolutionMat(A,B) + 5/11*NullspaceMat(A)[1];
[ 7/33, 1/3, 5/11 ]
M*A = B;
true
Here is another example.
A:= [[1,2,1,1],[-1,1,1,4],[0,-1,1,2],[-1,-4,-4,-8],[-5,1,-6,-5]];
[ [ 1, 2, 1, 1 ], [ -1, 1, 1, 4 ], [ 0, -1, 1, 2 ], [ -1, -4, -4, -8 ], [ -5, 1, -6, -5 ] ]
B:=[-1,-1,3,8];
[ -1, -1, 3, 8 ]
If we compute
NullspaceMat(A);
[ [ 2, 1, 1, 1, 0 ], [ 3, -2, 5, 0, 1 ] ]
we see that this time it is two dimensional. The general solution will be the particular solution SolutionMat(A, B) plus any linear combination of the basis vector for NullspaceMat(A).
For example
M:= SolutionMat(A,B) + 5/11*NullspaceMat(A)[1] - 3/4* NullspaceMat(A)[2];
[ -59/44, 65/22, -57/44, 5/11, -3/4 ]
will satisfy MA = B.
M*A = B;
true
In these examples the system of equations has a solution. We could have checked this by seeing that the rank of A was equal to the rank of the augmented matrix
Rank(A);
3
To create the augmented matrix of the system we take a copy of A, then add the vector B as the final row.
augA:= ShallowCopy(A);
[ [ 1, 2, 1, 1 ], [ -1, 1, 1, 4 ], [ 0, -1, 1, 2 ], [ -1, -4, -4, -8 ], [ -5, 1, -6, -5 ] ]
Add(augA,B);
Now augA and A should have the same rank since our system had a solution.
Rank(augA) = Rank(A);
true
Of course if the rank of the augmented matrix is different from the rank of A then the system has no solution. Here is an example.
B:=[1,-1,3,8];
[ 1, -1, 3, 8 ]
augA:= ShallowCopy(A);
[ [ 1, 2, 1, 1 ], [ -1, 1, 1, 4 ], [ 0, -1, 1, 2 ], [ -1, -4, -4, -8 ], [ -5, 1, -6, -5 ] ]
Add(augA,B); Rank(augA) = Rank(A);
false
If we tried to solve XA = B, for this different row vector B, then we will find no solution to the equations and SolutionMat(A, B) will returnfail.
SolutionMat(A,B);
fail
We begin by defining R to be the ring of polynomials over the rational numbers in the indeterminate x. First we set up a list of names that GAP will use for the indeterminates; in this case it only contains one entry.
indetnames:= ["x"];
[ "x" ]
Now define R to be the polynomial ring in a single indeterminate which GAP will write as we have specified.
R:= PolynomialRing(Rationals, indetnames);
PolynomialRing(..., [ x ])
Now we want to use "x" for the name of the indeterminate so we set this up.
indets:= IndeterminatesOfPolynomialRing(R);
[ x ]
x:= indets[1];
x
We let o be the identity of R.
o:= One( R );
1
We must be careful when we define matrices over R to make sure we use the notation 2*o for the integer 2 in R. The element 2 itself is an integer, but not an element of R. Let us define a matrix over R.
mat:= [[o, x],[2*o + x, x^2]];
[ [ 1, x ], [ 2+x, x^2 ] ]
To check that mat is a matrix use the IsMatrix command.
IsMatrix( mat );
true
Note that the following will not define a matrix as the elements do not come from a ring (remember that 1 and 2 are not in R but x is in R).
mat1:= [ [1, x], [2+x, x^2] ];
[ [ 1, x ], [ 2+x, x^2 ] ]
IsMatrix(mat1);
false
We can carry out the usual matrix operations with mat.
mat^2;
[ [ 1+2*x+x^2, x+x^3 ], [ 2+x+2*x^2+x^3, 2*x+x^2+x^4 ] ]
mat^(-1);
[ [ -1/2*x, 1/2 ], [ x^-1+1/2, -1/2*x^-1 ] ]
To convert mat1 to a matrix we could multiply it by o.
mat1:= mat1*o;
[ [ 1, x ], [ 2+x, x^2 ] ]
IsMatrix(mat1);
true
mat1^(-1);
[ [ -1/2*x, 1/2 ], [ x^-1+1/2, -1/2*x^-1 ] ]
Determinant(mat1);
-2*x
We set up a polynomial ring over the rationals with 4 indeterminates. As before we first set up the names that GAP will use for the indeterminates, then set up the names that we can use for the indeterminates (which for convenience we make the same).
indetnames:= ["t","x","y","z"];; R:= PolynomialRing(Rationals, indetnames);; indets:= IndeterminatesOfPolynomialRing(R);
[ t, x, y, z ]
t:= indets[1];; x:= indets[2];; y:= indets[3];; z:= indets[4];;
Here is a $4\times 4$ matrix over R in one indeterminate
o:= One( R );; mat:= List([1..4], i -> [2^(i-1), 5^(i-1), 17^(i-1), x^(i-1)])*o;; Display(mat);
[ [ 1, 1, 1, 1 ], [ 2, 5, 17, x ], [ 4, 25, 289, x^2 ], [ 8, 125, 4913, x^3 ] ]
Factors(Determinant(mat));
[ -9180+540*x, -5+x, -2+x ]
GAP will factorise polynomials in one indeterminate but not polynomials with more than one indeterminate
o:= One(R);; mat:= List([1..4], i -> [t^(i-1), x^(i-1), y^(i-1), z^(i-1)])*o;; Display(mat);
[ [ 1, 1, 1, 1 ], [ t, x, y, z ], [ t^2, x^2, y^2, z^2 ], [ t^3, x^3, y^3, z^3 ] ]
Determinant(mat);
t^3*x^2*y-t^3*x^2*z-t^3*x*y^2+t^3*x*z^2+t^3*y^2*z-t^3*y*z^2-t^2*x^3*y+t^2*x^3* z+t^2*x*y^3-t^2*x*z^3-t^2*y^3*z+t^2*\ y*z^3+t*x^3*y^2-t*x^3*z^2-t*x^2*y^3+t*x^2*z^3+t*y^3*z^2-t*y^2*z^3-x^3*y^2*z+x^3* y*z^2+x^2*y^3*z-x^2*y*z^3-x*y^3*z^\ 2+x*y^2*z^3
We define a matrix using a function
f:=function(i, j)
if j < i then return x;
else
return -2;
fi;
end;<
function( i, j ) ... end
n:= 6;; mat:= List([1 .. n], i -> List([1 .. n], j -> f(i,j)))*o;; Display(mat);
[ [ -2, -2, -2, -2, -2, -2 ], [ x, -2, -2, -2, -2, -2 ], [ x, x, -2, -2, -2, -2 ], [ x, x, x, -2, -2, -2 ], [ x, x, x, x, -2, -2 ], [ x, x, x, x, x, -2 ] ]
Determinant(mat);64+160*x+160*x^2+80*x^3+20*x^4+2*x^5
Factors(Determinant(mat));[ 4+2*x, 2+x, 2+x, 2+x, 2+x ]Clearly there is a theorem here!
for n in [2..12] do mat:= List([1 .. n], i -> List([1 .. n], j -> f(i,j)))*o;; Display(mat); Print(Determinant(mat), "\n"); Print(Factors(Determinant(mat)), "\n"); od;[ [ -2, -2 ], [ x, -2 ] ] 4+2*x [ 4+2*x ] [ [ -2, -2, -2 ], [ x, -2, -2 ], [ x, x, -2 ] ] -8-8*x-2*x^2 [ -4-2*x, 2+x ] [ [ -2, -2, -2, -2 ], [ x, -2, -2, -2 ], [ x, x, -2, -2 ], [ x, x, x, -2 ] ] 16+24*x+12*x^2+2*x^3 [ 4+2*x, 2+x, 2+x ] [ [ -2, -2, -2, -2, -2 ], [ x, -2, -2, -2, -2 ], [ x, x, -2, -2, -2 ], [ x, x, x, -2, -2 ], [ x, x, x, x, -2 ] ] -32-64*x-48*x^2-16*x^3-2*x^4 [ -4-2*x, 2+x, 2+x, 2+x ] [ [ -2, -2, -2, -2, -2, -2 ], [ x, -2, -2, -2, -2, -2 ], [ x, x, -2, -2, -2, -2 ], [ x, x, x, -2, -2, -2 ], [ x, x, x, x, -2, -2 ], [ x, x, x, x, x, -2 ] ] 64+160*x+160*x^2+80*x^3+20*x^4+2*x^5 [ 4+2*x, 2+x, 2+x, 2+x, 2+x ] [ [ -2, -2, -2, -2, -2, -2, -2 ], [ x, -2, -2, -2, -2, -2, -2 ], [ x, x, -2, -2, -2, -2, -2 ], [ x, x, x, -2, -2, -2, -2 ], [ x, x, x, x, -2, -2, -2 ], [ x, x, x, x, x, -2, -2 ], [ x, x, x, x, x, x, -2 ] ] -128-384*x-480*x^2-320*x^3-120*x^4-24*x^5-2*x^6 [ -4-2*x, 2+x, 2+x, 2+x, 2+x, 2+x ] [ [ -2, -2, -2, -2, -2, -2, -2, -2 ], [ x, -2, -2, -2, -2, -2, -2, -2 ], [ x, x, -2, -2, -2, -2, -2, -2 ], [ x, x, x, -2, -2, -2, -2, -2 ], [ x, x, x, x, -2, -2, -2, -2 ], [ x, x, x, x, x, -2, -2, -2 ], [ x, x, x, x, x, x, -2, -2 ], [ x, x, x, x, x, x, x, -2 ] ] 256+896*x+1344*x^2+1120*x^3+560*x^4+168*x^5+28*x^6+2*x^7 [ 4+2*x, 2+x, 2+x, 2+x, 2+x, 2+x, 2+x ] [ [ -2, -2, -2, -2, -2, -2, -2, -2, -2 ], [ x, -2, -2, -2, -2, -2, -2, -2, -2 ], [ x, x, -2, -2, -2, -2, -2, -2, -2 ], [ x, x, x, -2, -2, -2, -2, -2, -2 ], [ x, x, x, x, -2, -2, -2, -2, -2 ], [ x, x, x, x, x, -2, -2, -2, -2 ], [ x, x, x, x, x, x, -2, -2, -2 ], [ x, x, x, x, x, x, x, -2, -2 ], [ x, x, x, x, x, x, x, x, -2 ] ] -512-2048*x-3584*x^2-3584*x^3-2240*x^4-896*x^5-224*x^6-32*x^7-2*x^8 [ -4-2*x, 2+x, 2+x, 2+x, 2+x, 2+x, 2+x, 2+x ] [ [ -2, -2, -2, -2, -2, -2, -2, -2, -2, -2 ], [ x, -2, -2, -2, -2, -2, -2, -2, -2, -2 ], [ x, x, -2, -2, -2, -2, -2, -2, -2, -2 ], [ x, x, x, -2, -2, -2, -2, -2, -2, -2 ], [ x, x, x, x, -2, -2, -2, -2, -2, -2 ], [ x, x, x, x, x, -2, -2, -2, -2, -2 ], [ x, x, x, x, x, x, -2, -2, -2, -2 ], [ x, x, x, x, x, x, x, -2, -2, -2 ], [ x, x, x, x, x, x, x, x, -2, -2 ], [ x, x, x, x, x, x, x, x, x, -2 ] ] 1024+4608*x+9216*x^2+10752*x^3+8064*x^4+4032*x^5+1344*x^6+288*x^7+36*x^8+2*x^9 [ 4+2*x, 2+x, 2+x, 2+x, 2+x, 2+x, 2+x, 2+x, 2+x ] [ [ -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2 ], [ x, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2 ], [ x, x, -2, -2, -2, -2, -2, -2, -2, -2, -2 ], [ x, x, x, -2, -2, -2, -2, -2, -2, -2, -2 ], [ x, x, x, x, -2, -2, -2, -2, -2, -2, -2 ], [ x, x, x, x, x, -2, -2, -2, -2, -2, -2 ], [ x, x, x, x, x, x, -2, -2, -2, -2, -2 ], [ x, x, x, x, x, x, x, -2, -2, -2, -2 ], [ x, x, x, x, x, x, x, x, -2, -2, -2 ], [ x, x, x, x, x, x, x, x, x, -2, -2 ], [ x, x, x, x, x, x, x, x, x, x, -2 ] ] -2048-10240*x-23040*x^2-30720*x^3-26880*x^4-16128*x^5-6720*x^6-1920* x^7-360*x^8-40*x^9-2*x^10 [ -4-2*x, 2+x, 2+x, 2+x, 2+x, 2+x, 2+x, 2+x, 2+x, 2+x ] [ [ -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2 ], [ x, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2 ], [ x, x, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2 ], [ x, x, x, -2, -2, -2, -2, -2, -2, -2, -2, -2 ], [ x, x, x, x, -2, -2, -2, -2, -2, -2, -2, -2 ], [ x, x, x, x, x, -2, -2, -2, -2, -2, -2, -2 ], [ x, x, x, x, x, x, -2, -2, -2, -2, -2, -2 ], [ x, x, x, x, x, x, x, -2, -2, -2, -2, -2 ], [ x, x, x, x, x, x, x, x, -2, -2, -2, -2 ], [ x, x, x, x, x, x, x, x, x, -2, -2, -2 ], [ x, x, x, x, x, x, x, x, x, x, -2, -2 ], [ x, x, x, x, x, x, x, x, x, x, x, -2 ] ] 4096+22528*x+56320*x^2+84480*x^3+84480*x^4+59136*x^5+29568*x^6+10560* x^7+2640*x^8+440*x^9+44*x^10+2*x^11 [ 4+2*x, 2+x, 2+x, 2+x, 2+x, 2+x, 2+x, 2+x, 2+x, 2+x, 2+x ]