ralu

<< "Algebra`PolynomialPowerMod`"
<< "FiniteFields`"

n = 10; p = 5;
r = (p^n - 1)/(p - 1);
fac = FactorInteger[r];
factor[i_] := fac[[i]][[1]];
facnum = Length[fac];
j = 0;
m[i_] := r/factor[i];


While[j < 1000, j = j + 1;
  coef = Prepend[
    Append[Table[RandomInteger[{0, p - 1}], {i, n - 1}], 1],
    RandomInteger[{1, p - 1}]];
  poly = Total[Table[coef[[i]] x^(i - 1), {i, 1, n + 1}]];
  If[MultiplicativeOrder[(-1)^n coef[[1]], p1] == p - 1, Continue[];];
   If[MemberQ[Table[Mod[poly, p] /.x -> i, {i, 1, p - 1}], 0],
   Continue[];];
  Qrows = Table[PolynomialPowerMod[x, p k, {poly, p}], {k, 0, n - 1}];
   Pd[x_] := PadRight[x, n]; Q = Pd /@ CoefficientList[Qrows, x];
  QI = Q - IdentityMatrix[n];
  ns = NullSpace[Transpose[QI], Modulus -> p];
  If[Dimensions[ns][[1]] != 1, Continue[];];
  If[PolynomialGCD[poly, \!\(
\*SubscriptBox[\(\[PartialD]\), \(x\)]poly\), Modulus -> p] != 1,
   Continue[];]; a = PolynomialPowerMod[x, r, {poly, p}];
  If[IntegerQ[a] == False, Continue[];];
  If[Mod[(-1)^n coef[[1]] - a, p] != 0, Continue[];];
  For[i = 2, i <= facnum,
   If[IntegerQ[PolynomialPowerMod[x, 80, {poly, p}]], Continue[];];
   i++]; Print[j]; Print[coef]; Print[poly];
  Break[];];