// To create Table 2.1 (first block):
Ordena("1-2000.txt",8 : invers := true);
ListToTex("1-2000.ord.res" : selected := [1,2,3,4,8], factors := [1,3,8]);

// To create Table 2.1 (second block):

OrderByBiggestExponent("1-2000.txt":invers:=true);

//Clean a bit the list...
a:=Open("1-2000.max.ord.res", "r");
s:=Gets(a);
A:=[]; B:=[]; C:=[[],[],[],[],[],[],[],[],[],[]];
while not IsEof(s) do
  s:=StringToIntegerSequence(s);
  p:=Factorization(s[6]);
  s[6]:=Max([p[i][1] : i in [1..#p] | p[i][2] eq Max([p[j][2] : j in [1..#p]])]);
  if s[6] notin A then 
    Append(~A,s[6]);
    Append(~B,0);
  end if;
  if B[Index(A,s[6])] le 5 then Append(~C[s[7]],s); B[Index(A,s[6])]+:=1; end if;
  s:=Gets(a);
end while;
D:=[];
for i in [10,9,8,7,6,5,4,3,2,1] do
  D:= D cat Ordena(C[i],6:invers:=true, sortida:=[]);
end for;
delete(a);
PrintMatrixFile("1-2000.max.ord.seleccio.res",D);

ListToTex("1-2000.max.ord.seleccio.res" : selected := [1,2,3,4,5], factors := [1,3,5]);


// Compute N1eqN2.res (Table 2.3):
A:=OpenAndRead("1-2000.txt");
B:=[];
for i in [1..#A] do
  if A[i][1] eq A[i][3] and A[i][5] eq A[i][8] then
    Append(~B,A[i]);
  end if;
end for;
PrintMatrixFile("N1eqN2.res",B);  // 50275 línies
Ordena("N1eqN2.res",8 : invers := true);
ListToTex("N1eqN2.ord.res": selected := [1,2,4,5], factors := [1,5]);


// Create 2 lists: for every modular form, we compute the (probable) reducible primes, 
// and the ramification primes.

A:=Open("prob-red-primes.res","w");
B:=Open("ramif-primes.res","w");

for N in [1..1000] do
  Dec:=NewformDecomposition(NewSubspace(CuspidalSubspace(ModularSymbols(N,2,+1))));
  for i in [1..#Dec] do
    F:=Dec[i];
    f:=qEigenform(F,1000);
    k:=Parent(Coefficient(f,1));

    if Type(k) eq FldRat then 
      fprintf B,"%o %o 1\n",N, i; Flush(B);
    else   
// Depenent de si el Magma és el meu o el de l'IEM...
      P:=DefiningPolynomial(k);
//    P:=Modulus(k);
      K:=ext<Rationals()|P>;
      fprintf B,"%o %o %o\n",N, i, Discriminant(K); Flush(B);
    end if;

    p:=2;
    g:=0;
    while p lt 1000 do
      if N mod p eq 0 then p:=NextPrime(p); continue; end if;
      g:=Gcd(g,Integers()!Norm(Coefficient(f,p)^2-(p+1)^2));
      if g eq 1 then break; end if; 
      p:=NextPrime(p);
    end while;
    fprintf A,"%o %o %o\n",N, i, g; Flush(A);
  end for;
end for;
delete(A);
delete(B);





// Compute elements satisfying Theorem 3.1.3: (not optimal at all)
A:=OpenAndRead("1-2000.txt"); //204438 lines
B:=[];
reds:=OpenAndRead("prob-red-primes.res");
ramifs:=OpenAndRead("ramif-primes.res");


"Reading UpperBound2.0.res";
a:=Open("1-2000.txt","r");
aa:=Gets(a);
s:=[**];
s[1]:=[];
while IsEof(aa) ne true do
  aa:=StringToIntegerSequence(aa);
  while #s lt (aa[1]-1) do s[#s+1]:=[]; end while;
  if #s lt aa[1] and aa[2] gt 1 then s[#s+1]:=[**]; end if; //En aquest cas #s=aa[1]-1
  if #s lt aa[1] and aa[2] eq 1 then s[#s+1]:=[*[[aa[3],aa[4],aa[5]]]*]; aa:=Gets(a); continue; end if;
  while #s[aa[1]] lt (aa[2]-1) do s[aa[1]][#s[aa[1]]+1]:=[]; end while; //Aquí #s=aa[1]
  if #s[aa[1]] eq (aa[2]-1) then s[aa[1]][aa[2]]:=[[aa[3],aa[4],aa[5]]]; aa:=Gets(a); continue; end if;
  if #s[aa[1]] eq aa[2] then s[aa[1]][aa[2]]:=Append(s[aa[1]][aa[2]],[aa[3],aa[4],aa[5]]); end if;
  aa:=Gets(a);
end while;
delete(a); 

Unique:=function(N,i,M,j,l) //Suppose N>M
  h:=[];
  for w in [1..#s[N][i]] do
    if s[N][i][w][1] eq M and s[N][i][w][2] ne j and s[N][i][w][3] mod l eq 0 then 
      h:=Append(h,s[N][i][w][2]);
    end if;
  end for;

  for j2 in [1..#h] do
    if h[j2] gt j then 
      for w in [1..#s[M][h[j2]]] do
if s[M][h[j2]][w][1] eq M and s[M][h[j2]][w][2] eq j and s[M][h[j2]][w][3] mod l eq 0 then return false; end if;
      end for;
    else
      for w in [1..#s[M][j]] do
        if s[M][j][w][1] eq M and s[M][j][w][2] eq h[j2] and s[M][j][w][3] mod l eq 0 then return false; end if;
      end for;
    end if;
  end for;
  return true;
end function;

Unique2:=function(N,i,M,j,l)
  p:=Factorization(l);
  l2:=1;
  for k in p do
    if Unique(N,i,M,j,k[1]) then l2*:=k[1]^k[2]; end if;
  end for;
  return l2;
end function;

W:=AssociativeArray();
for i in [1..#ramifs] do
   if ramifs[i][1] ne reds[i][1] or ramifs[i][2] ne reds[i][2] then "ERROR"; end if;
   W[[ramifs[i][1],ramifs[i][2]]]:=Integers()!(ramifs[i][3]*reds[i][3]);
end for;


for i in [1..#A] do
if i mod 1000 eq 0 then i; end if;
  C:=[Integers()!A[i][1],Integers()!A[i][2],Integers()!A[i][3],Integers()!A[i][4],pPart(A[i][5],2),pPart(A[i][8],2)];
  p:=Factorization(Integers()!(C[1]/C[3]));
  if #p ne 1 then continue; end if;
  if C[3] mod p[1][1] eq 0 then continue; end if;  // p not divide N
  C[5] := pPart(C[5],p[1][1]); C[6] := pPart(C[6],p[1][1]);// p and l must be different
  while Gcd(C[5],C[3]) ne 1 do  // l must not divide N
    C[5]:=Integers()!(C[5]/Gcd(C[3],C[5]));
  end while;
  while Gcd(C[6],C[3]) ne 1 do  // l must not divide N
    C[6]:=Integers()!(C[6]/Gcd(C[3],C[6]));
  end while;
  if C[5] eq 1 then continue; end if;  
  w:=W[[C[3],C[4]]];
  while Gcd(w,C[5]) ne 1 do
    C[5]:=Integers()!(C[5]/Gcd(w,C[5]));
  end while;
  while Gcd(w,C[6]) ne 1 do
    C[6]:=Integers()!(C[6]/Gcd(w,C[6]));
  end while;
  if C[5] eq 1 then continue; end if;  
  if IsMinimal(C[3],C[4],C[5]) eq false then continue; end if;
  C[5]:=Unique2(C[1],C[2],C[3],C[4],C[5]);
  C[6]:=Unique2(C[1],C[2],C[3],C[4],C[6]);
  if C[5] eq 1 then continue; end if; 
  Append(~B,C);
end for;
#B;
PrintMatrixFile("Pre-Teorema",B); // 8746 lines


//Finally, we sort the results:
Ordena("Teorema.res",6 : invers := true);
ListToTex("Teorema.ord.res" : selected := [1,2,3,4,5,6], factors := [1,3,5,6]);





//Funcions auxiliars poc importants:
ChangeDirectory("c:/Documents and Settings/Xavier/Escritorio/Heckes");

Dimension:=function(N,i) // Calcula la dimensió d'f=(N,i)
  a:=Open(IntegerToString(N)*"-"*IntegerToString(i)*".txt","r");
  s:=Gets(a);
  s:=StringToIntegerSequence(s);
  delete(a);
  return(#s-1);
end function;

Dimension2:=function(N,i,j)
  return([Dimension(N,i),Dimension(N,j)]);
end function;