F := GF(2);
R<x> := PolynomialRing(F);

// g^8 + (x^3 + x)*g^4 + x

function TestIt(g)
   h := g^4;
   return IsEven(#Factorization(h^2 + (x^3 + x)*h + x));
end function;

function TestRandomOfDegree(n)
   g := R!([Random(0,1) : i in [0..n-1]] cat [1]);
   val := TestIt(g);
   if not val then
      print "Counterexample: g = ", g;
   end if;
   return val;
end function;

function TestDegree(n)
   V := VectorSpace(GF(2),n);
   for v in V do
      h := R!(Eltseq(v) cat [1])^4;
      f := h^2 + (x^3 + x)*h + x;
      F := Factorization(f);
      if not IsEven(#F) then
         print h^2 + (x^3 + x)*h + x, " is irreducible!!";
         return false;
      end if;
   end for;
   return true;
end function;

function TestDegreeRedOnly(n)
   V := VectorSpace(GF(2),n);
   for v in V do
      h := R!(Eltseq(v) cat [1])^4;
      f := h^2 + (x^3 + x)*h + x;
      if IsIrreducible(f) then
         print h^2 + (x^3 + x)*h + x, " is irreducible!!";
         return false;
      end if;
   end for;
   return true;
end function;


/////////////////////////////////////////////////////
/*


Pick an integer n = 1,2,...


and let g run over polys (not just monic)
in F_p[x] of degree \le n (say for p = 3, 5).
Count the number of times x^{p-1}g^p + 1
is irreducible, and take ratio between this
and sum_{j=1}^{n+1} (p^j)/(pj - 1).
I am interesting to know for large n
whether this ratio gets close to 1; perhaps
give ratios to three decimal places,
say for n going up to n = 20 or 15 or so
(if that is reasonable, or up to whatever
is not too time-consuming).
*/

function IrredRatio(p,n)
   k    := GF(p);
   V    := VectorSpace(k,n+1);
   R<x> := PolynomialRing(k);
   w := x^(p-1);
   cnt := &+[IsIrreducible(w*(R!Eltseq(v))^p + 1) select 1 else 0 : v in V];
   sum := &+[p^j/(p*j - 1) : j in [1..n+1]];
   return <p, n, cnt, sum, cnt/sum*1.0>;
end function;





function TestDegree5(n)
   V := VectorSpace(GF(5),n);
   for v in V do
      h := R!(Eltseq(v) cat [1])^4;
      f := x^2*h^5 + 1;
      F := Factorization(f);
      if not IsEven(#F) then
         print h^2 + (x^3 + x)*h + x, " is irreducible!!";
         return false;
      end if;
   end for;
   return true;
end function;

function TestDegree5b(n)
   p := 7;
   k := GF(p);
   R<x> := PolynomialRing(k);
   V := VectorSpace(k,n);
   minus_one_to_n := (-1)^n;
   char := [KroneckerSymbol(a,p) : a in [1..p]];
   cnt := 0;
   x_pow := x^(p-3);
   for v in V do
      for i in [1..p-1] do 
         h := R!(Eltseq(v) cat [i]);
         f := x_pow*h^p + 1;
         if char[Integers()!Discriminant(f)]*minus_one_to_n eq 1 then
            cnt := cnt + 1;
         end if;
      end for;
   end for;
   proportion := cnt/(#V*(p-1))*1.0;
   return proportion, 1-proportion;
end function;
