################################################################
#
# File author: Otto Forster
#
# This code is placed under the GNU General Public Licence
#
# Date of last change
# 2006-07-13
(**************************************************************)
(*
** Polynomial arithmetic in Fp[X]
** Fp = GF(p) = Z/pZ, p a prime number
** A polynomial F(X) = a0 + a1*X + ... + an*X**n
** is represented by an array (a0, a1, ..., an), 0 <= ak < p
** The coefficient an is supposed to be /= 0 modulo p
** The null polynomial is represented by the empty array ().
*)
(**************************************************************)
(*
** deletes leading zeroes
*)
function array_trim(F: array): array;
var
    d, i: integer;
begin
    d := length(F)-1;
    while d >= 0 and F[d] = 0 do
        dec(d);
    end;
    return F[0..d];
end;
(*------------------------------------------------------------*)
(*
** shifts array F by sh positions to the right (sh must be non-negative)
*)
function array_shift(F: array; sh: integer): array;
var
    len: integer;
    R: array;
begin
    if sh = 0 then
        return F;
    elsif sh > 0 then
        len := length(F);
        R := alloc(array,len+sh);
        R[sh..len+sh-1] := F;
        return R;
    else
        writeln("sh must be non-negative");
        return F;
    end;
end;
(*------------------------------------------------------------*)
(*
** Product of two polynomials in Fp[X]
*)
function fpX_mult(p:integer; F,G: array): array;
var
    n,m,k: integer;
    x: integer;
    R,F1: array;
begin
    n := length(F)-1;
    m := length(G)-1;
    if n<0 or m<0 then
        return ();
    end;
    R := F * G[0];
    for k := 1 to m do
        x := G[k];
        if x /= 0 then
            F1 := F*x;
            F1 := array_shift(F1,k);
            R := R + F1;
        end;
    end;
    return R mod p;
end;
(*------------------------------------------------------------*)
(*
** Calculates F(X) mod G(X) in Fp[X]
** The leading coefficient of G is supposed to be invertible
*)
function fpX_mod(p: integer; F,G: array): array;
var
    j,n,m,a,c: integer;
begin
    n := length(F)-1;
    m := length(G)-1;
    if n < m or m < 0 then
        return F;
    end;
    c := G[m];
    if c /= 1 then
        c := mod_inverse(c,p);
        if c = 0 then
            writeln("error in fpX_mod: leading coeff not invertible");
            halt(-1);
        end;
        G := (G * c) mod p;
    end;
    for j := n to m by -1 do
        a := F[j];
        if a /= 0 then
            F[j-m..j] := F[j-m..j] - G*a;
        end;
    end;
    F := F[0..m-1] mod p;
    return array_trim(F);
end;
(*------------------------------------------------------------*)
(*
** modular power F(X)**ex mod Q(X)
*)
function fpX_modpow(p: integer; F: array; ex: integer; Q: array): array;
var
    i: integer;
    R: array;
begin
    if ex < 0 then
        writeln("error");
        writeln("power: exponent must be non-negative: ",ex);
        return ();
    elsif ex = 0 then
        return {1};
    else
        R := F;
    end;
    for i := bit_length(ex)-2 to 0 by -1 do
        R := fpX_mod(p,fpX_mult(p,R,R),Q);
        if bit_test(ex,i) then
            R := fpX_mod(p,fpX_mult(p,R,F),Q);
        end;
    end;
    return R;
end;
(*------------------------------------------------------------*)
(*
** Evaluates F(x) mod p.
*)
function fpX_val(p: integer; F: array; x: integer): integer;
var
    i,z: integer;
begin
    z := 0;
    for i := length(F)-1 to 0 by -1 do
        z := (z*x + F[i]) mod p;
    end;
    return z;
end;
(**************************************************************)
(*
** Uebungen zur Vorlesung von
** Otto Forster: Endliche Koerper, SS 2006
** Aufgabe 23 c)
*)
(**************************************************************)
(***
Ein Polynom F(X) vom Grad 101 ueber dem Koerper F3 = GF(3) ist
genau dann irreduzibel, wenn es keine Nullstellen hat und F(X)
ein Teiler von X**(3**101) - X ist, d.h. X**(3**101) = X mod F(X).

Wir machen fuer das gesuchte irreduzible Polynom den Ansatz

    F(X) = X**101 + f4(X),

wobei f4(X) ein Polynom vom Grad <= 4 ist.
Da fuer x /= 0 mod 3 gilt x**101 = x mod 3, hat F(X) genau dann
keine Nullstelle in F3, wenn das konstante Glied von f4 ungleich 0 ist
und f4(1) + 1 /= 0 und f4(2) + 2 /= 0.
Die Funktion ff4(n) konstruiert ein Polynom vom Grad <= 4 aus der
3-adischen Entwicklung der ganzen Zahl n.
Die Konstruktion eines irreduziblen Polynoms vom Grad 101 ueber
dem Koerper F3 geschieht durch die Funktion aufg23c()

==> aufg23c().
,,,,,,..,,.,,.,.,,.,,.,,,.,,,,..,,,,,.,.,,,.,.,,,.,.,,..,,,,,,,.,,,.,,.,,.,.,,
.,,,.,,.,.,,.,,.,,,.,,,,,,,..,,.,.,,,.,.,,,.,.,,,,,..,,,,.,,,.,,.,,.146
-: (2, 0, 1, 2, 1)

Dies bedeutet, dass
    P1(X) = X**101 + X**4 + 2*X**3 + X**2 + 2
irreduzibel ist.
Um weitere Loesungen zu finden, ruft man aufg23c mit einem Argument auf:

==> aufg23c(147).
,.,,.,,..155
-: (2, 0, 2, 2, 1)

==> aufg23c(156).
,,,,,,,.,,.,,,.,,,,,,,..,,.,,.,.,,.,,,.,.,,,,,..,,,,,.,.,,,.,,.,.,,.,,.226
-: (1, 0, 1, 2, 2)

==> aufg23c(227).
.,,,,,,,.235
-: (1, 0, 2, 2, 2)

==> aufg23c(236).
,,,.,,.no irreducible polynomial found
-: ()

Also sind
    P2(X) = X**101 + X**4 + 2*X**3 + 2*X**2 + 2
    P3(X) = X**101 + 2*X**4 + 2*X**3 + X**2 + 1
    P4(X) = X**101 + 2*X**4 + 2*X**3 + 2*X**2 + 1
weitere irreduzible Polynome.
Es gibt keine weiteren irreduziblen Polynome ueber F3 der Gestalt
X**101 + f4(X) mit deg f4 <= 4.

***)
(*--------------------------------------------------------------*)
(*
** Erzeugt ein Polynom vom Grad <= 4 ueber F3 aus der
** 3-adischen Entwicklung der Zahl n
*)
function ff4(n: integer): array[5];
var
    i: integer;
    ff: array[5];
begin
    for i := 0 to 4 do
        ff[i] := n mod 3;
        n := n div 3;
    end;
    return ff;
end;
(*--------------------------------------------------------------*)
(*
** Sucht ein irreduzibles Polynom ueber F3 der Form
**      F(X) = X**101 + f4(X)
** mit einem Polynom f4 vom Grad <= 4.
** Die Suche wird mit dem durch die Zahl n0 dargestellten
** Polynom f4 = ff4(n0) begonnen, Default n0 = 1
** Rueckgabewert ist ein Array der Laenge 5, das das Polynom f4(X)
** darstellt.
** Falls kein irreduzibles Polynom gefunden wird, wird das
** leere Array () zurueckgegeben
*)
function aufg23c(n0 := 1): array;
var
    P, R, X, ff: array;
    i,j,k,n: integer;
begin
    P := alloc(array,102,0);
    P[101] := 1;
    (* P stellt das Polynom X**101 dar *)
    X := (0,1);
    for n := n0 to 3**5-1 do
        ff := ff4(n);
        if ff[0] = 0 or
           (fpX_val(3,ff,1) + 1) mod 3 = 0 or
           (fpX_val(3,ff,2) + 2) mod 3 = 0 then
           write(",");
           continue;
        end;
        P[0..4] := ff;
        (* nun stellt P das Polynom X**101 + ff(X) dar *)
        R := fpX_modpow(3,X,3**101,P);
        write(".");
        if R = X then
            writeln(n);
            return ff;
        end;
    end;
    writeln("no irreducible polynomial found");
    return ();
end;
(**************************************************************)