#### Arithmetic in the Galois Field GF(p^m)

These procedures require the Modular arithmetic and Polynomial Operations Modulo a Prime packages.

```/* gf(P, M, AddTable, MultTable) is true if AddTable and MultTable are       */
/*   addition and multiplication tables in GF(P^M) for prime P and M>=2.     */
/* e.g. gf(2, 2, X, Y) gives                                                 */
/*  X=[0,1,2,3,1,0,3,2,2,3,0,1,3,2,1,0], Y=[0,0,0,0,0,1,2,3,0,2,3,1,0,3,1,2] */
/* i.e. X is 0 1 2 3  Y is 0 0 0 0                                           */
/*           1 0 3 2       0 1 2 3                                           */
/*           2 3 0 1       0 2 3 1                                           */
/*           3 2 1 0       0 3 1 2                                           */
gf(P, M, AddTable, MultTable):-
M > 0,
irreducible_polynomial(M, P, Qs), !,
P_1 is P - 1,
findall(I, for(0, P_1, I), Is),
findall(Xs, with_replacement(M, Is, Xs), Xss),
% For P=2, M=3:
% Qs=[1,0,1,1], i.e. x^3+x+1
% Is=[0,1]
% Xss=[[0,0,0],[0,0,1],[0,1,0],[0,1,1],[1,0,0],[1,0,1],[1,1,0],[1,1,1]]
%           0,      1,      x,    x+1,    x^2,  x^2+1,  x^2+x, x^2+x+1
findall(Z, add_table(Xss, P, Z), AddTable),
findall(Z, mult_table(Xss, P, Qs, Z), MultTable).

add_table(Xss, P, Z):-
member(Xs, Xss),
member(Ys, Xss),
p_addmod(Xs, Ys, P, Z0),
evaluate(Z0, P, Z).

mult_table(Xss, P, Qs, Z):-
member(Xs, Xss),
member(Ys, Xss),
p_multmod(Xs, Ys, P, Ws),
p_divmod(Ws, Qs, P, _, Z0),
evaluate(Z0, P, Z).

/* evaluate(Polynomial, X, Y) is true if Y is the value of the Polynomial    */
/*   evaluated for the value X.                                              */
/* e.g. evaluate([2,3,4,5], 2, 41).                                          */
evaluate([], _, 0).
evaluate([A|As], X, Y):-
evaluate_1(As, X, A, Y).

evaluate_1([], _, Y, Y).
evaluate_1([A|As], X, Y0, Y):-
Y1 is Y0 * X + A,
evaluate_1(As, X, Y1, Y).

/* is_irreducible(+Qs, +P) is true if Qs is an irreducible polynomial        */
/*   modulo P.                                                               */
/* A polynomial q(x) of degree d is irreducible modulo p if and only if      */
/*   x^(p^d)-x is divisible by q(x), and q(x) and x^(p^(d/r))-x are          */
/*   mutually prime for all primes r that divide d.                          */
is_irreducible(Qs, P):-
length_1(Qs, -1, D),
power(P, D, PD),
power_of_x(PD, XPD),
p_subtmod(XPD, [1,0], P, XPDX),
p_divmod(XPDX, Qs, P, _, [0]),
is_irreducible_1(2, D, Qs, P).

is_irreducible_1(R, D, _, _):-
R > D, !.
is_irreducible_1(R, D, Qs, P):-
is_prime(R),
0 =:= D mod R, !,
DR is D // R,
power(P, DR, PDR),
power_of_x(PDR, XPDR),
p_subtmod(XPDR, [1,0], P, XPDRX),
p_gcdmod(XPDRX, Qs, P, [1]),
R1 is R + 1,
is_irreducible_1(R1, D, Qs, P).
is_irreducible_1(R, D, Qs, P):-
R1 is R + 1,
is_irreducible_1(R1, D, Qs, P).

/* power_of_x(P, Rs) is true if Rs is the polynomial x^P.                    */
/* e.g. power_of_x(3, [1,0,0,0]).                                            */
power_of_x(P, [1|Rs]):-power_of_x_1(P, Rs).

power_of_x_1(0, []):-!.
power_of_x_1(I, [0|Rs]):-
I1 is I - 1,
power_of_x_1(I1, Rs).

/* irreducible_polynomial(+D, +P, -Qs) is true if Qs is an irreducible       */
/*   monic polynomial of degree D, modulo a prime P.                         */
/* The number of irreducible polynomials (modulo 2) for D=1,2,... are        */
/*   2,1,2,3,6,9,18,30,56,99,186,335,630,1161,2182,...                       */
/* e.g. irreducible_polynomial(4, 3, [1,0,0,1,2]). and 17 other solutions    */
irreducible_polynomial(D, P, [1|Qs]):-
is_prime(P),
P_1 is P - 1,
findall(I, for(0, P_1, I), Is),
with_replacement(D, Is, Qs),
is_irreducible([1|Qs], P).

/* for(I, J, K) is true if K is an integer between I and J inclusive.        */
for(I, J, I):-I =< J.
for(I, J, K):-I < J, I1 is I + 1, for(I1, J, K).

/* is_prime(P) is true if P is a prime number.                               */
is_prime(2):-!.
is_prime(P):-
P > 2,
P mod 2 =:= 1,
SqrtP is sqrt(P),
is_prime_1(3, P, SqrtP).

is_prime_1(I, _, SqrtP):-
I > SqrtP, !.
is_prime_1(I, P, SqrtP):-
P mod I > 0,
I1 is I + 2,
is_prime_1(I1, P, SqrtP).

/* power(X, N, Y) is true if Y is the result of raising X to the power N.    */
power(X, N, Y):-N >= 0, power_1(N, X, 1, Y).

power_1(0, _, Y, Y):-!.
power_1(1, X, Y0, Y):-!,
Y is Y0 * X.
power_1(N, X, Y0, Y):-
1 =:= N mod 2, !,
N1 is N // 2,
X1 is X * X,
Y1 is Y0 * X,
power_1(N1, X1, Y1, Y).
power_1(N, X, Y0, Y):-
N1 is N // 2,
X1 is X * X,
power_1(N1, X1, Y0, Y).

/* with_replacement(N, Xs, Ys) is true if Ys is the result of picking N      */
/*   things out of Xs with replacement.                                      */
with_replacement(0, _, []):-!.
with_replacement(I, Xs, [X|Ys]):-
member(X, Xs),
I1 is I - 1,
with_replacement(I1, Xs, Ys).

/* length(Xs, L) is true if L is the number of elements in the list Xs.      */
%length(Xs, L):-length_1(Xs, 0, L).

length_1([], L, L).
length_1([_|Xs], L0, L):-L1 is L0 + 1, length_1(Xs, L1, L).

/* member(X, Xs) is true if the element X is contained in the list Xs.       */
%member(X, [X|_]).
%member(X, [_|Xs]):-member(X, Xs).
```