
// ******* ENTRY DATA ****************************

  d := 7; // (AT LEAST 3) 
  all := true; 
  bound := 10^7;

// RUNNING THIS FILE WILL OUTPUT A LIST badprimes
// IF all IS SET true, THE LIST CONTAINS ALL BAD
// PRIMES FOR DEGREE d. IF IT IS SET false, THE
// LIST CONTAINS ALL BAD PRIMES FOR DEGREE d THAT
// ARE SMALLER THAN OR EQUAL TO bound.

// IF all IS false, THEN ALSO A LIST goodprimes
// IS GENERATED, CONTAINING THE PRIMES SMALLER
// THAN OR EQUAL TO BOUND THAT ARE NOT BAD.

// IF all IS SET true, THE ROUTINE IS FEASIBLE UP
// TO DEGREE 7, OR MAYBE 8.

// IF all IS SET false, THE ROUTINE IS FEASIBLE UP
// TO DEGREE 10, OR MAYBE 11 (DEPENDING ON
// THE VALUE OF bound).

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

if all then
  print "\n";
  print "*****************************************************";
  print " Computing all bad Casas-Alvero primes for degree", d;
  print "*****************************************************";
  print "";
else
  print "\n";
  print "*************************************************************************";
  print " Computing all bad Casas-Alvero primes for degree", d, "smaller than", bound;
  print "*************************************************************************";
  print "";
end if;

// ***** COMPUTING SCENARIOS *****

Write("helpfile.m", "listofscenarios := [ []: t in [1.."*IntegerToString(d-1)*"] ];": Overwrite := true);
// listofscenarios[i] will contain the scenarios of type t

Write("helpfile.m", "for i1 in [0..1] do m1 := Max([i1]) + 1;");

arraystring := "i1";
endstring := "end for; ";

for j in [2..(d-2)] do
    arraystring *:= ",i"*IntegerToString(j);
    Write("helpfile.m", "for i"*IntegerToString(j)*" in [0..m"*IntegerToString(j-1)*"] do m"*IntegerToString(j)*" := Max(["*arraystring*"]) + 1;");
    endstring *:= "end for; ";
end for;
Write("helpfile.m", "type := m"*IntegerToString(d-2)*"-1; listofscenarios[type+1] cat:= [["*arraystring*"]];");
Write("helpfile.m", endstring);

tim := Cputime();
load "helpfile.m";
print "Time needed to compute scenarios: ", Cputime(tim);

// **** COLLECTING CANDIDATE BAD PRIMES ****

tim := Cputime();
print "Collecting candidate bad primes ...";
candbadprimes := {};

if not all then
  p := 2;
  while p le bound do
    candbadprimes join:={p};
    p := NextPrime(p);
  end while;
else
print "--> Calculating Nullstellensatz expansions in characteristic 0 ...";
print "    --> Verifying type 0 ... 1 scenario(s) ... (trivial)";

for type in [1..(d-2)] do
  numofscens := #listofscenarios[type+1];
  print "    --> Verifying type", type, "...", numofscens, "scenario(s) ...";
  if not (type in [1,d-2]) then
    R := PolynomialRing(Rationals(),d-3,"elim",[type..(d-3)]);
  else
    R := PolynomialRing(Rationals(),d-3,"grevlex");
  end if;
  S<x> := PolynomialRing(R);
  // CONSTRUCT CANDIDATE COUNTEREXAMPLE f
    if type eq 1 then f := x^2*(x-1); else f := x^2*(x-1)*(&*[ x-R.j : j in [1..(type-1)]]); end if;
    if type ne d-2 then g := x^(d-2-type) + &+[ R.(j+type)*x^j : j in [0..(d-3-type)]]; f *:= g; end if;
  // print Factorization(f);
  // treat all scenarios of the given type
  counter := 0;
  for s in listofscenarios[type+1] do
//    counter +:=1;
    I := [];
    derf := Derivative(f);
    for j in [1..(d-2)] do
      derf := Derivative(derf);
      der := derf div Factorial((j+1)); // HASSE DERIVATIVES
      case s[j]:
        when 0:
          I cat:= [Evaluate(der,0)];
        when type:
          I cat:= [Evaluate(der,1)];
        else:
          I cat:= [Evaluate(der,R.s[j])];
      end case;
    end for;
    // simplify ideal with respect to leading zeroes
      index := 1;
      if s[index] eq 0 then
        I[index] := R.(type+index-1);
        index +:= 1;
      end if;
      ev := [ R.j : j in [1..(type-1)] ] cat [0 : j in [1..(index - 1)] ] cat [ R.j : j in [(type+index-1)..(d-3)] ];
      if type ge 2 then I[index] := Evaluate(I[index],ev) div R.1; end if;
    Id := IdealWithFixedBasis(I);
    co := Coordinates(Id, Parent(Basis(Id)[1]) ! 1);
    for c in co do
      for coef in Coefficients(c) do
        for fac in Factorization(Denominator(coef)) do
          candbadprimes join:={fac[1]};
        end for;
      end for;
    end for;
//    print "  ", counter, "out of", numofscens, "scenarios done ( in type", type, ")";
  end for;
end for;

end if;

print "Time needed to collect candidate bad primes:", Cputime(tim);
candbadprimes := Sort(SetToSequence(candbadprimes));
print "Number of candidate bad primes collected:", #candbadprimes;
badprimes := [];

// **** FILTERING CANDIDATE BAD PRIMES ****

tim := Cputime();
print "Filtering the list of candidate bad primes ...";

for p in candbadprimes do
  CAok := true;

  for type in [1..(d-2)] do
    numofscens := #listofscenarios[type+1];
    if not (type in [1,d-2]) then
      R := PolynomialRing(GF(p),d-3,"elim",[type..(d-3)]);
    else
      R := PolynomialRing(GF(p),d-3,"grevlex");
    end if;
    R0 := PolynomialRing(Integers(),d-3);
    S<X> := PolynomialRing(R);
    S0<x> := PolynomialRing(R0);
    // construct candidate counterexample f
      if type eq 1 then f := x^2*(x-1); else f := x^2*(x-1)*(&*[ x-R0.j : j in [1..(type-1)]]); end if;
      if type ne d-2 then g := x^(d-2-type) + &+[ R0.(j+type)*x^j : j in [0..(d-3-type)]]; f *:= g; end if;
    // print Factorization(f);
    // treat all scenarios of the given type
    counter := 0;
    for s in listofscenarios[type+1] do
      if CAok then
        I := [];
        derf := Derivative(f);
        for j in [1..(d-2)] do
          derf := Derivative(derf);
          der := S ! (derf div Factorial((j+1))); // HASSE DERIVATIVES
          case s[j]:
            when 0:
              I cat:= [Evaluate(der,0)];
            when type:
              I cat:= [Evaluate(der,1)];
            else:
              I cat:= [Evaluate(der,R.s[j])];
          end case;
        end for;
        // simplify ideal with respect to leading zeroes
          index := 1;
          if s[index] eq 0 then
            I[index] := R.(type+index-1);
            index +:= 1;
          end if;
          ev := [ R.j : j in [1..(type-1)] ] cat [0 : j in [1..(index - 1)] ] cat [ R.j : j in [(type+index-1)..(d-3)] ];
          if type ge 2 then I[index] := Evaluate(I[index],ev) div R.1; end if;
        Id := Ideal(I);
        CAok := (GroebnerBasis(Id) eq [1]);
//        counter +:= 1;
//        print "  ", counter, "out of", numofscens, "scenarios done ( in type", type, ")";
       end if;
    end for;
  end for;

if not CAok then 
  badprimes cat:= [p]; 
//  print p, "was a bad prime";
//  else print p, "was a good prime"; 
end if;
end for;

// **** FINAL OUTPUT ****

print "Time needed to filter bad primes:", Cputime(tim);
print "Number of bad primes:", #badprimes, "(stored in a variable called badprimes)";
if not all then
  goodprimes := [];
  for p in candbadprimes do
    if p in candbadprimes and not (p in badprimes) then goodprimes cat:= [p]; end if;
  end for;
  print "Number of good primes:", #goodprimes, "(stored in a variable called goodprimes)";
end if;
