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

  d := 7;
  
  // RUNNING THIS FILE PROVES THAT OUR LIST OF BAD
  // PRIMES FOR DEGREE 7 IS CORRECT.

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

// ***** 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 DENOMINATORS ****

tim := Cputime();
print "Collecting denominators ...";
denoms := {};
candbadprimes := {};

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);

    // THIS IS USED IN THE FOR-LOOP BELOW
    if type ge 3 then
      Rgrevlex := PolynomialRing(Rationals(),d-3,"grevlex");
      I := [ Rgrevlex ! f : f in I];
      Id := IdealWithFixedBasis(I);
      co1 := Coordinates(Id, Parent(Basis(Id)[1]) ! 1);
      lcm1 := LCM(&cat[ [Denominator(coef) : coef in Coefficients(c)] : c in co1 ]);
      Rglex := PolynomialRing(Rationals(),d-3,"glex");
      I := [ Rglex ! f : f in I];
      Id := IdealWithFixedBasis(I);
      co2 := Coordinates(Id, Parent(Basis(Id)[1]) ! 1);
      lcm2 := LCM(&cat[ [Denominator(coef) : coef in Coefficients(c)] : c in co2 ]);
      Rweight1 := PolynomialRing(Rationals(),d-3,"weight",[1,1,1,1,0,1,1,1,0,0,1,1,0,0,0,1]);
      I := [ Rweight1 ! f : f in I];
      Id := IdealWithFixedBasis(I);
      co3 := Coordinates(Id, Parent(Basis(Id)[1]) ! 1);
      lcm3 := LCM(&cat[ [Denominator(coef) : coef in Coefficients(c)] : c in co3 ]);
      Rweight2 := PolynomialRing(Rationals(),d-3,"weight",[1,1,1,1,0,1,1,0,0,0,1,0,0,0,0,1]);
      I := [ Rweight2 ! f : f in I];
      Id := IdealWithFixedBasis(I);
      co4 := Coordinates(Id, Parent(Basis(Id)[1]) ! 1);
      lcm4 := LCM(&cat[ [Denominator(coef) : coef in Coefficients(c)] : c in co4 ]);
      Rweight3 := PolynomialRing(Rationals(),d-3,"weight",[1,1,1,1,1,0,0,1,0,1,0,1,0,0,0,1]);
      I := [ Rweight3 ! f : f in I];
      Id := IdealWithFixedBasis(I);
      co5 := Coordinates(Id, Parent(Basis(Id)[1]) ! 1);
      lcm5 := LCM(&cat[ [Denominator(coef) : coef in Coefficients(c)] : c in co5 ]);
      Rweight4 := PolynomialRing(Rationals(),d-3,"weight",[1,1,1,1,1,2,3,4,1,4,9,16,1,8,27,64]);
      I := [ Rweight4 ! f : f in I];
      Id := IdealWithFixedBasis(I);
      co6 := Coordinates(Id, Parent(Basis(Id)[1]) ! 1);
      lcm6 := LCM(&cat[ [Denominator(coef) : coef in Coefficients(c)] : c in co6 ]);
      gcd := GCD([lcm1,lcm2,lcm3,lcm4,lcm5,lcm6]);
    end if;

    for c in co do
      for coef in Coefficients(c) do
        if type lt 3 then
          for fac in Factorization(Denominator(coef)) do
            candbadprimes join:={fac[1]};
          end for;
        else
          // IF type GE 4, THEN FACTORING THE DENOMINATORS IS TOO HARD.
          // TO GET RID OF THE FIRST FEW SUPERFLUOUS FACTORS, WE CALCULATE
          // NULLSTELLENSATZ EXPANSIONS FOR OTHER MONOMIAL ORDERINGS AND
          // TAKE GCDS
          denom := Denominator(coef);
          denomclean := GCD(denom,gcd);
          denoms join:={denomclean}; 
        end if;
      end for;
    end for;
//    print "  ", counter, "out of", numofscens, "scenarios done ( in type", type, ")";
  end for;
end for;

print "Time needed to collect denominators:", Cputime(tim);
denoms := SetToSequence(denoms);

// **** CLEANING DENOMINATORS ****

// OUR LIST OF CANDIDATE BAD PRIMES CONTAINS
// * ALL PRIMES UP TO 10^4
// * ALL PRIMES OF WHICH WE WANT TO PROVE THAT THEY ARE BAD PRIMES 
//   (TAKEN FROM badprimes7.txt)
// * A LIST OF EXTRA FACTORS THAT ARE USEFUL

load "badprimes7.txt";

candbadprimes join:= SequenceToSet(PrimesUpTo(10^4));
candbadprimes join:= SequenceToSet(badprimes7);

usefulfactors := 
{ 17159, 20089, 55441, 1845271, 2820641, 8589011, 353554063, 3685982137, 70827530059021, 143936487176023, 845538765728453, 
16110900222320863943960941, 161187396205242613484472618863, 3787912581003012853171440511037, 26944974629464206815013793743977535691, 
167613033080539302933921057826083544468012281960389771830037510186474470967, 
978468659031134217945444969222523636114093453487247216875144998215297397860991, 
95032775895156250011023088748162168418212365743244528636804387414226003768633334439106291 };

candbadprimes join:= usefulfactors;

repeat

  coulddosomething := false;

  print #denoms, "denominators"; print #candbadprimes, "candidate bad primes";
  print "filtering...";

  foundfactors := false;
  for i in [1..#denoms] do for p in candbadprimes do
    val := Valuation(denoms[i],p);
    denoms[i] := denoms[i] div p^val;
    foundfactors or:= (val gt 0);
  end for; end for;
  denoms := SetToSequence(SequenceToSet(denoms));

  // remove 1s
  repeat
    if #denoms gt 0 then
      min, pos := Min(denoms);
      if min eq 1 then denoms := Remove(denoms,pos); end if;
    end if;
  until (min gt 1) or (#denoms eq 0);

  foundsmalldenoms := false;
  for d in denoms do
    if d lt 10^90 then
      fac := Factorization(d);
      for f in fac do candbadprimes join:= {f[1]}; end for;
      foundsmalldenoms := true;
    end if;
  end for;

  foundprimes := false;
  for d in denoms do
    if IsPrime(d) then
      candbadprimes join:={d};
      foundprimes := true;
    end if;
  end for;

  foundgcds := false;
  if not (foundsmalldenoms or foundprimes or foundfactors) then
    gcds := [];
      for i in [1..#denoms] do
        for j in [(i+1)..#denoms] do
          gcd := GCD(denoms[i],denoms[j]); 
          if gcd gt 1 then
            denoms[i] := denoms[i] div gcd; denoms[j] := denoms[j] div gcd;
            gcds cat:= [gcd];
            foundgcds := true;
          end if;
        end for;
      end for;
    denoms cat:= gcds;
    denoms := SetToSequence(SequenceToSet(denoms));
  end if;

  if (foundsmalldenoms or foundprimes or foundfactors or foundgcds) then coulddosomething := true; end if;

until not coulddosomething;


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

badprimes := [];

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 ****

badprimes := Sort(badprimes);

print "Time needed to filter bad primes:", Cputime(tim);
print "Number of bad primes:", #badprimes, "(stored in a variable called badprimes)";
print "Number of untreated denominators:", #denoms, "(stored in a variable called denoms)";
