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
![Pictures in GF[p]](dgfp_01.png)
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