# DEMOGFP A short demonstration of gfp-numbers, elements of GF[p]

## Initialization

This toolbox implements arithmetical and other operation in GF[p], not GF[p^n], i.e. it is basically a modulo p arithmetic. To initialize the prime number p, use

p = 7;
gfpinit(p)
Galois field prime set to 7
ans =
7

The prime number p is limited by 2*p^2<=2^53, i.e. p <= 2^26 = 67,108,864. To retrieve the current prime, use

current_prime = gfpinit
Galois field prime set to 7
current_prime =
7

## Basic operations

The four basic operations are executed as usual, for example, an addition table is as follows:

v = gfp(0:p-1);
0     1     2     3     4     5     6
1     2     3     4     5     6     0
2     3     4     5     6     0     1
3     4     5     6     0     1     2
4     5     6     0     1     2     3
5     6     0     1     2     3     4
6     0     1     2     3     4     5

Similarly, a multiplication table is obtained by

mul_table = v'*v
gfp mul_table =
0     0     0     0     0     0     0
0     1     2     3     4     5     6
0     2     4     6     1     3     5
0     3     6     2     5     1     4
0     4     1     5     2     6     3
0     5     3     1     6     4     2
0     6     5     4     3     2     1

## The reciprocal

As GF[p] being a field, every nonzero number has a unique reciprocal. One way to obtain the reciprocal is the first row of the division table:

v = gfp(1:p-1);
div_table = v'*(1./v)
gfp div_table =
1     4     5     2     3     6
2     1     3     4     6     5
3     5     1     6     2     4
4     2     6     1     5     3
5     6     4     3     1     2
6     3     2     5     4     1

or directly by

reciprocal = 1./v
gfp reciprocal =
1     4     5     2     3     6

The reciprocal can also be obtained by Fermat's small theorem x^(p-1) = 1, so that x^(p-2) must be the reciprocal of x:

v^(p-2)
gfp ans =
1     4     5     2     3     6

## Conversion of floating-point numbers

Every floating-point number is of the form m*2^e for integers m and e. Thus, they can be mapped into GF[p] by (m mod p) * 2^(e mod p):

gfp_pi = gfp(pi)
gfp_e = gfp(exp(1))
gfp gfp_pi =
0
gfp gfp_e =
5

The latter formula can be used for any odd prime; for p=2 and m=m'*2^f we use (m' mod p) * 2^(e+f). Obviously the result is zero except for odd integers.

## primitive element

Every field has at least one primitive element, i.e. a generator of the multiplicative group. Using

powers = repmat(gfp(2:p-1)',1,p-2).^repmat(1:p-2,p-2,1)
gfp powers =
2     4     1     2     4
3     2     6     4     5
4     2     1     4     2
5     4     6     2     3
6     1     6     1     6

generates all powers of 1..p-1 in the columns of powers. Each row not containing a 1 is a primitive element:

prim_elem = find(all(powers~=1,2))+1
prim_elem =
3
5

Indeed the powers are all elements 1..p-1:

gfp(prim_elem(1)) ^ (1:p-1)
gfp(prim_elem(2)) ^ (1:p-1)
gfp ans =
3     2     6     4     5     1
gfp ans =
5     4     6     2     3     1

## Matrix routines

Operations in GF[p] are always exact, neither rounding errors nor overflow or underflow may occur. Thus the term "ill-conditioned" does not apply. For arbitrarily ill-conditioned matrices it can be proved that they are regular if the determinant mod p is nonzero. Consider

gfpinit(11)
n = 5;
A = invhilb(n)
[L,U,perm] = luwpp(gfp(A));
residual = L*U - A(perm,:)
det_ = prod(diag(U))
Galois field prime set to 11
ans =
11
A =
25        -300        1050       -1400         630
-300        4800      -18900       26880      -12600
1050      -18900       79380     -117600       56700
-1400       26880     -117600      179200      -88200
630      -12600       56700      -88200       44100
gfp residual =
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
gfp det_ =
2

The routine luwpp performs an LU-decomposition of A with partial pivoting. Thus, if det_ is nonzero, the matrix is A is regular. With the small prime p=11 that approach does not work for the inverse Hilbert matrix of dimension 6 because it contains a zero row and column:

A6 = gfp(invhilb(6))
gfp A6 =
3     8     5     8     3     0
8     4     9     7     6     0
5     9     4     1     6     0
8     7     1    10     9     0
3     6     6     9     1     0
0     0     0     0     0     0

## The true Hilbert matrix in GF[p]

Alternatively one may define H to be the true Hilbert matrix in GF[p]. Consider H_5 in GF[11]:

n = 5;
gfpinit(11);
H = ( 1:n )' + ( 1:n) - 1
H = 1 ./ gfp(H)
Galois field prime set to 11
H =
1     2     3     4     5
2     3     4     5     6
3     4     5     6     7
4     5     6     7     8
5     6     7     8     9
gfp H =
1     6     4     3     9
6     4     3     9     2
4     3     9     2     8
3     9     2     8     7
9     2     8     7     5

That matrix is regular in GF[11]:

d = det(H)
gfp d =
6

Hence there is a unique inverse:

R = inv(H)
R*H
gfp R =
3     8     5     8     3
8     4     9     7     6
5     9     4     1     6
8     7     1    10     9
3     6     6     9     1
gfp ans =
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 pseudoinverse

While a matrix being regular in GF[p] has a unique inverse in GF[p], the corresponding theorem is not true for the pseudoinverse. It is not difficult to see that A=[1 3] has no pseudoinverse in GF[5]. While 0.1*[1;3] is the pseudoinverse over R, also a brute-force test demonstrates that there is no pseudoinverse in GF[5]:

gfpinit(5);
A = [1 3];
for i=0:4
for j=0:4
P = gfp([i;j]);
if A * P == 1
P
PA = P*A
disp(' ')
end
end
end
Galois field prime set to 5
gfp P =
0
2
gfp PA =
0     0
2     1

gfp P =
1
0
gfp PA =
1     3
0     0

gfp P =
2
3
gfp PA =
2     1
3     4

gfp P =
3
1
gfp PA =
3     4
1     3

gfp P =
4
4
gfp PA =
4     2
4     2

The code shows that for all matrices P with A*P=1 implies that P*A is not symmetric. The reason that no pseudoinverse exists, although A has full rank, is that A*A' is not regular:

AAt = gfp(A)*A'
gfp AAt =
0

## Regularity of extremely ill-conditioned matrices

The INTLAB routine "isregular" checks regularity on that basis for extremely ill-conditioned matrices.

Indeed, using a larger prime, very ill-conditioned matrices can be proved to be regular:

gfpinit(67108859)
n = 50;
A = randmat(n,1e200);
C = double( norm(A,inf) * norm(inv(sym(A)),inf) )
[L,U,perm] = luwpp(gfp(A));
det_ = prod(diag(U))
Galois field prime set to 67108859
ans =
67108859
C =
1.6780e+237
gfp det_ =
67108858

In that case we used the largest admissable prime for the gfp-toolbox, so that it is unlikely that the test fails by accident.

Note that the condition number C of the matrix A is computed using the symbolic toolbox of Matlab. That takes some time, but it is reliable.

## Pictures in GF[p]

For the friends of nice pictures consider

p = 101;
gfpinit(p);
ex = .5*[-1 1 1 -1];
ey = .5*[-1 -1 1 1];
close, hold on
for i=1:p
for j=1:p
patch( i+ex , j+ey , double( mod( i^2 + 3*i*j ,p))/p );
end
end
axis off
hold off
Galois field prime set to 101

## Enjoy INTLAB

INTLAB was designed and written by S.M. Rump, head of the Institute for Reliable Computing, Hamburg University of Technology. Suggestions are always welcome to rump (at) tuhh.de