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

## Contents

## 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); add_table = v'+v

gfp add_table = 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