// **********************************************************************
// RUNNING THIS FILE PROVES THE CASAS-ALVERO CONJECTURE IN DEGREE 12
// THE THEORETICAL JUSTIFICATION IS DESCRIBED IN THE ARTICLE
// **********************************************************************

// **********************************************************************
// PARAMETERS:

tim := Cputime();
precomp := true; // USE OUR PRECOMPUTED SCENARIOS OR NOT

p := NextPrime(10^7); // PRIME MODULO WHICH PROOF IS ATTEMPTED

types_to_be_checked := [7..7]; // SET EQUAL TO [0..8] IF YOU WANT
                               // A FULL PROOF (HOWEVER, FROM TYPE
                               // 6 ON IT IS BETTER TO SPLIT THE
                               // COMPUTATION IN PACKS)

packno := 4;                   // SET BOTH EQUAL TO 1 FOR A FULL
totalnopacks := 12;            // PROOF. THE SCENARIOS WILL BE SPLIT
                               // IN totalnopacks PACKS AND ONLY
                               // PACK packno WILL BE VERIFIED

printprogress := true;  // DO YOU WANT THE PROGRAM TO PRINT A CONFIRMATION
                        // EACH TIME A SCENARIO HAS BEEN VERIFIED?

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

// **********************************************************************
// EXTERNAL FILES:

load "precompscenarios12.m"; // THIS CAN OF COURSE BE OMITTED IF precomp
                             // IS SET false, BUT FOR SOME REASON MAGMA
                             // REFUSES TO INCLUDE THIS IN AN
                             // IF-STATEMENT

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

// 1. COMPUTE THE LIST L_res^cl
// ----------------------------------------------------------------------

print "Computing list of scenarios to be dealt with...";

// FIRST WE COMPUTE L = listofscenarios, THE LIST OF ALL SCENARIOS FOR 
// DEGREE 12

if precomp then 
  // PRECOMPSCENARIOS OMITS THE "0" WITH WHICH EVERY SCENARIO STARTS,
  // SO WE ADD IT BELOW
  for i in [1..11] do
    n := #listofscenarios[i];
    for k in [1..n] do
      listofscenarios[i][k] := [0] cat listofscenarios[i][k];
    end for;
  end for;
else
  print "precomp set false -- may take about half an hour...";
  listofscenarios := [ [] : i in [1..11] ];
  s1 := 0;
  for s2 in [0..1] do m2 := s2 + 1;
  for s3 in [0..m2] do m3 := Max([s2,s3]) + 1;
  for s4 in [0..m3] do m4 := Max([s2,s3,s4]) + 1;
  for s5 in [0..m4] do m5 := Max([s2,s3,s4,s5]) + 1;
  for s6 in [0..m5] do m6 := Max([s2,s3,s4,s5,s6]) + 1;
  for s7 in [0..m6] do m7 := Max([s2,s3,s4,s5,s6,s7]) + 1;
  for s8 in [0..m7] do m8 := Max([s2,s3,s4,s5,s6,s7,s8]) + 1;
  for s9 in [0..m8] do m9 := Max([s2,s3,s4,s5,s6,s7,s8,s9]) + 1;
  for s10 in [0..m9] do m10 := Max([s2,s3,s4,s5,s6,s7,s8,s9,s10]) + 1;
  for s11 in [0..m10] do type := Max([s2,s3,s4,s5,s6,s7,s8,s9,s10,s11]);
    listofscenarios[type+1] cat:= [[s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11]];
  end for; end for; end for; end for; end for; end for; end for; end for; end for; end for;
end if;

// NOW WE COMPUTE L_res, THE LIST OF SCENARIOS SATISFYING THE DETERMINANT
// CONDITION, MOREOVER EXCLUDING THOSE THAT MATCH WITH PROPOSITION 15

// THE FOLLOWING FUNCTION TESTS WHETHER, FOR A GIVEN SCENARIO s FOR DEGREE d = p+1,
// THE DETERMINANT CONDITION IS SATISFIED MODULO p.

function checkdet(s)
  p := #s;
  d := p + 1;
  js := {};
  for j in [2..(d-2)] do
    if s[j] eq s[d-1] then js join:= {d-j}; end if;
  end for;
  js := SetToSequence(js);
  n := #js + 1;
  matrixrowbyrow := [];
  for i in [1..(n-1)] do
    matrixrowbyrow cat:=[-1];
      for k in [1..i] do
        matrixrowbyrow cat:= [Binomial(js[i]-2,js[k]-2)*js[i]];
      end for;
      for k in [(i+1)..(n-1)] do
        matrixrowbyrow cat:= [0];
      end for;
  end for;
  matrixrowbyrow cat:= [-1];
  for i in [1..(n-1)] do
    matrixrowbyrow cat:= [(-1)^js[i]];
  end for;
  M := Matrix(Integers(),n,n,matrixrowbyrow);
  if (Determinant(M) mod p eq 0) and (#js ge 2) then return true; else return false; end if;
end function;

function prop15(s)
  // CHECK s_11 = s_1 = 0 (Proposition 15 applied to p=11, 12 = p+1)
  // OR s_4 = s_8 (Proposition 15 applied to p=2, 12 = (p+1)*p^2)
  // OR s_3 = s_9 (Proposition 15 applied to p=3, 12 = (p+1)*p)
  if s[11] eq 0 or s[4] eq s[8] or s[3] eq s[9] then return true; else return false; end if;
end function;

print "--> Filtering...";

L_res := [ [] : i in [1..11] ];
for i in [1..11] do
  for s in listofscenarios[i] do
    if checkdet(s) and (not prop15(s)) then L_res[i] cat:= [s]; end if;
  end for;
end for;

// NOW WE CLOSE L_res UNDER TAKING DESCENDANTS (THE RESULT IS STILL CALLED L_res)

// THE FOLLOWING FUNCTION COMPUTES THE SET OF DESCENDANTS OF A GIVEN SCENARIO s

function descendants(s)
  desc := {};
  d := #s + 1;
  for i in [1..Max(s)] do
    sprime := [0];
    for j in [2..(d-1)] do
      if s[j] eq i then sprime cat:= [0]; end if;
      if s[j] lt i then sprime cat:= [s[j]]; end if;
      if s[j] gt i then sprime cat:= [s[j]-1]; end if;
    end for;
    desc join:= {sprime};
  end for;
  return desc;
end function;

print "--> Closing under taking descendants...";

for i in [0..9] do
  desc := {};
  for s in L_res[11-i] do desc join:= descendants(s); end for;
  L_res[10-i] := SetToSequence(SequenceToSet(L_res[10-i]) join desc);
end for;

print "--> Done. Time needed:", Cputime(tim);

// 2. MAIN ROUTINE
// ----------------------------------------------------------------------

CAok := true;

for type in types_to_be_checked do
tim := Cputime();

if type eq 0 then 
  print "Verifying type 0 (trivial).";
else
  print "Verifying type", type;
  L := L_res[type + 1];
  print "--> Total number of scenarios to be checked:", #L;
  packsize := #L div totalnopacks;
  packrem := #L mod totalnopacks;
  if packno eq 1 then packmin := 1; packmax := packsize + packrem;
  else packmin := (packno-1)*packsize + 1 + packrem; packmax := packno*packsize + packrem; end if;
  print "--> Current pack: scenarios",packmin,"to",packmax;
  
  if type eq 1 then 
    R := PolynomialRing(GF(p),9,"grevlex");
  else
    R := PolynomialRing(GF(p),9,"elim",[type..9]);
  end if;
  R0 := PolynomialRing(Integers(),9);
  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;
    g := x^(10-type) + &+[ R0.(j+type)*x^j : j in [0..(9-type)]]; f *:= g;
  // TREAT ALL SCENARIOS IN THE PACK
  counter := 0;
  for numero in [packmin..packmax] do
    s := L[numero];
    if CAok then
      I := [];
      derf := Derivative(f);
      for j in [1..10] do
        derf := Derivative(derf);
        der := S ! (derf div Factorial((j+1))); // HASSE DERIVATIVES
        case s[j+1]:
          when 0:
            I cat:= [Evaluate(der,0)];
          when type:
            I cat:= [Evaluate(der,1)];
          else:
            I cat:= [Evaluate(der,R.s[j+1])];
        end case;
      end for;
      // SIMPLIFY IDEAL WITH RESPECT TO LEADING ZEROES
        index := 1;
        if s[index+1] 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)..9] ];
        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;
      if printprogress then print "  ", counter, "out of", packmax-packmin + 1, "scenarios done ( this was scenario", numero, ")"; end if;
     end if;
  end for;

  if totalnopacks eq 1 and packno eq 1 then 
    if CAok then
      print "--> Casas-Alvero OK for this type!";
    else 
      print "--> WARNING: CASAS-ALVERO NOT OK FOR THIS TYPE!!!"; 
    end if;
    print "--> Time needed to complete this type:", Cputime(tim);
  else
    if CAok then 
      print "Casas-Alvero OK for this pack!"; 
    else 
      print "WARNING: CASAS-ALVERO NOT OK FOR THIS PACK!!!";
    end if;
    print "--> Time needed to complete this pack:", Cputime(tim), "( This was pack no", packno, "out of", totalnopacks, ").";
  end if;

end if;
end for;

if CAok then print "Casas-Alvero conjecture succesfully verified for the given types / packs."; end if;
