/* ************************************************************ ** ** ** Cong.mg ** ** ** ** Congruences mod ell^n ** ** ** ** Gabor Wiese & Xavier Taixés ** ** version of 17/02/09 ** ** ** ************************************************************ */ declare verbose Cong,1; SetVerbose("Cong",0); /* Computation of congruence numbers */ intrinsic CongruenceNumber (f :: RngUPolElt, g :: RngUPolElt) -> RngIntElt, RngUPolElt, RngUPolElt {The congruence number c(f,g) of the polynomials f and g is computed, as well as a linear combination c(f,g) = rf + sg. The function returns c(f,g), r, s.} /* We could assert here that f and g are coprime normed polynomials, as we have not defined the congruence number otherwise. */ // require (LeadingCoefficient(f) eq 1 and LeadingCoefficient(g) eq 1) : "The polynomials must be normed."; // require (Degree(Gcd(f,g)) lt 1) : "The polynomials must be coprime."; m := Degree(f); n := Degree(g); /* If both polynomials are normed and of degree zero, i.e. both are 1, then the congruence number is 0. */ if (n eq 0) and (m eq 0) then return 0,0,0; // treat empty sets before making Echelon forms end if; if n lt 1 then return 1,0,1; end if; if m lt 1 then return 1,1,0; end if; S := SylvesterMatrix(f,g); A,B := EchelonForm(S); v := B[m+n]; r := Polynomial([v[n-i+1] : i in [1..n]]); s := Polynomial([v[n+m-i+1] : i in [1..m]]); return A[m+n][m+n],r,s; end intrinsic; /* Computation of bounds */ // stupid function to take care of different types my_val := function (a,l) if Type(a) eq RngIntElt then out := Valuation(a,l); else out := Valuation(a); end if; return out; end function; intrinsic Bounds (f1 :: RngUPolElt, g1 :: RngUPolElt, l :: RngIntElt : factor := false, crs := <>) -> RngIntElt, RngIntElt {Computes upper and lower bounds on n for congruences of roots of f and g modulo l^n. This function returns lower bound, upper bound. If factor := true, then the polynomials f,g will be factored over pAdicRing(l) in order to improve the bounds in case they are not optimal.} if Degree(Gcd(f1,g1)) ge 1 then return []; end if; FacF:=[f[1] : f in Factorization(f1)]; FacG:=[g[1] : g in Factorization(g1)]; out_list := []; for f in FacF do for g in FacG do FX := PolynomialRing(GF(l)); /* If the polynomials are coprime mod l, then there is no congruence of zeros mod l. Hence, both the lower and the upper bound are zero. */ if Degree(Gcd(FX!f,FX!g)) eq 0 then Append(~out_list,<0,0>); continue; end if; /* If the congruence number has not yet been computed or if either f1 or g1 is reducible, then compute the congruence number. */ if (#crs eq 0) or (not IsIrreducible(f1)) or (not IsIrreducible(g1)) then c,r,s := CongruenceNumber(f,g); else c := crs[1]; r := crs[2]; s := crs[3]; end if; /* This is the upper bound. */ ub := my_val(c,l); /* If the upper bound is 1, it is trivially best possible and we can stop.*/ if ub eq 1 then vprint Cong,1: "Cor. 2.17 (b) applied."; Append(~out_list,<1,1>); continue; end if; /* If the reduction of f and the reduction of g modulo l both do not have any multiple factors, then we can apply part (c) of the corollary, telling us that the upper bound is also a lower bound. We test this condition by computing gcd mod l with the Euclidean algorithm, which is faster than computing the congruence number. */ if (Degree(Gcd(FX!f,Derivative(FX!f))) eq 0) and (Degree(Gcd(FX!g,Derivative(FX!g))) eq 0) then vprint Cong,1: "Cor. 2.17 (c)(i) applied."; Append(~out_list,); continue; else // Test conditions c(ii), d(i). if Degree(Gcd(FX!g,FX!s)) eq 0 then // gbar and sbar are coprime if Degree(Gcd(FX!g,Derivative(FX!g))) eq 0 then // gbar does not have any double factors vprint Cong,1: "Cor. 2.17 (c)(ii) applied."; b := ub; else vprint Cong,1: "Only Cor. 2.17 (d)(i) applied... hence bad result."; b := Ceiling(ub/Degree(g)); end if; else vprint Cong,1: "Only Cor. 2.17 (d)(iii) applied... trivial lower bound."; b := 1; end if; // Test conditions c(iii), d(ii). if Degree(Gcd(FX!f,FX!r)) eq 0 then // fbar and rbar are coprime if Degree(Gcd(FX!f,Derivative(FX!f))) eq 0 then // fbar does not have any double factors a := ub; vprint Cong,1: "Cor. 2.17 (c)(iii) applied."; else vprint Cong,1: "Only Cor. 2.17 (d)(ii) applied... hence bad result."; a := Ceiling(ub/Degree(f)); end if; else vprint Cong,1: "Only Cor. 2.17 (d)(iii) applied... trivial lower bound."; a := 1; end if; // the lower bound can be taken as the max of a and b. lb := Maximum(a,b); end if; if (lb ne ub) and factor then /* Improve the bounds by factoring f and g over the pAdicRing(l). */ /* As the precision we take some multiple of the max of the valuations of the congruence numbers of (f,g), (f,f'), (g,g') in order to try and make sure that we can distinguish any two zeros of either polynomial */ prec := Maximum({ub, my_val(CongruenceNumber(f,Derivative(f)),l), my_val(CongruenceNumber(g,Derivative(g)),l)})*30; ZlX := PolynomialRing(pAdicRing(l,prec)); vprint Cong,1: "Start factoring polynomials over l-adic ring..."; F := Factorisation(ZlX!f); G := Factorisation(ZlX!g); lbs := []; ubs := []; // store lower and upper bounds here // compare all pairs of factors for f1 in F do for g1 in G do here := Bounds(f1[1],g1[1],l : factor := false); if #here gt 1 then print "Stupid Error. Check!"; end if; l1 := here[1][1]; u1 := here[1][2]; Append(~lbs,l1); Append(~ubs,u1); end for; end for; // look for the maximum lower bound indices := [i : i in [1..#lbs] | lbs[i] eq Maximum(lbs)]; i := Index(ubs, Minimum([ubs[k] : k in indices])); lb := lbs[i]; ub := ubs[i]; end if; Append(~out_list, ); end for; end for; return out_list; end intrinsic; //Function to write an Integer Sequence to a file intrinsic IntegerSequenceToFile (A :: SeqEnum, S :: File) {Writes an integer sequence to a file. It puts a 0 in the non defined elements} for i in [1..#A] do if IsDefined(A,i) then fprintf S,"%o ", A[i]; else fprintf S,"0 "; end if; end for; fprintf S,"\n"; Flush(S); end intrinsic; intrinsic IntegerSequenceToFile (A :: SeqEnum, s :: MonStgElt : option:="w") {Opens the file s and writes an integer sequence to a file. It puts a 0 in the non defined elements Then it closes the file} S:=Open(s,option); IntegerSequenceToFile(A, S); delete(S); end intrinsic; //Function to put one line of a file to an integer sequence (default, the 1st line). //If the file does not exist, returns the empty sequence intrinsic FileLineToIntegerSequence ( s :: MonStgElt : n := 1) -> SeqEnum {Given a file, puts the n-th line to an integer sequence. Default: the first line. If n bigger than the number of lines of the file or the file does not exist, return the empty sequence} try a:= Open(s,"r"); for i in [1..n] do s:=Gets(a); if IsEof(s) then return []; end if; end for; delete(a); return(StringToIntegerSequence(s)); catch e return []; end try; end intrinsic; //Check if HB of f is precomputed. If not, we compute it. aux controls if we'll update "bounds.txt" or not CheckHBs := procedure(~HBs, ~f, ~aux) if not IsDefined(HBs,f`N) or HBs[f`N] eq 0 then f`M:=ModularSymbols(f`N,2,1); HBs[f`N] := HeckeBound(f`M); aux := true; end if; f`H := HBs[f`N]; end procedure; //Procedure to get the Hecke Bound and, if not previously computed, store it StoredHB := procedure(~f,~g) aux:=false; // aux controls if in one "CheckHB" a new HB is computed //We read the precomputed file. If it does not exist, we get HBs:=[]; HBs:=FileLineToIntegerSequence("bounds.txt"); //Compute HB for f and g CheckHBs(~HBs, ~f, ~aux); CheckHBs(~HBs, ~g, ~aux); if aux then //If we computed a new HB, we store it IntegerSequenceToFile(HBs,"bounds.txt"); end if; end procedure; //Function to asign f`D, if not defined (just to save some lines of code) AsignfD := procedure(~f, ~aux) //aux controls if we have computed a new D or not if not assigned f`D then if not assigned f`M then f`M:=ModularSymbols(f`N,2,1); end if; f`D := SortDecomposition(NewformDecomposition(NewSubspace(CuspidalSubspace(f`M)))); aux:=true; end if; end procedure; //Function to compute the number of newforms of level N using precomputed information NewformsOfLevelN := procedure(~f, ~NForm, ~aux) if not assigned f`d then if not assigned NForm`s then //If the file with precomputed info is not opened, we open it NForm`s:=FileLineToIntegerSequence("nforms.txt"); end if; //If the information we want is not available, we compute it: if not IsDefined(NForm`s,f`N) or NForm`s[f`N] eq 0 then AsignfD(~f, ~aux); aux:=true; // aux = true means: we will update "nforms.txt" NForm`s[f`N] := #f`D; end if; f`d := NForm`s[f`N]; end if; end procedure; //If not precomputed: computes CharpolyHecke. If precomputed: loads it. ComputeHecke := procedure(~f, ~F, ~aux1, ~file1, s1, ia1, Precomputed, p, ~Pf, ~aux) //aux1 is true means we are still looking in our precomputed file "N-i.txt" if aux1 then st1:=Gets(file1); //If "N-i.txt" has no more info, we initialize some values so that next step we can store the results if IsEof(st1) then aux1:=false; file1:=Open(s1,"a"); AsignfD(~f, ~aux); F := f`D[ia1]; // Once we initialized the variables, we can already call the same function but with aux1=false: $$(~f, ~F, ~aux1, ~file1, s1, ia1, Precomputed, p, ~Pf, ~aux); else // in case we still have not found the end of the file and we can still read the polynomials Pf:=Polynomial(StringToIntegerSequence(st1)); end if; else //if we have no more info stored, we compute the polynomial Pf := PolynomialRing(Integers())!(CharacteristicPolynomial(HeckeOperator(F,p))); if Precomputed then // if we did not have this polynomial in "N-i.txt", we save it for the next time for i in [0..Degree(Pf)] do fprintf file1, "%o ", Coefficient(Pf,i); end for; fprintf file1, "\n"; end if; end if; end procedure; HeckeInLevelNm:=function(Pf,deltap) if Type(Pf) eq SeqEnum then Pf:=Polynomial(Pf); end if; X:=Parent(Pf).1; x:=(X^2-(deltap)); P:=0; d:=Degree(Pf); for i in [0..d] do P+:=Coefficient(Pf,i)*X^(d-i)*x^i; end for; return P; end function; pPart:=function(l,p) //It computes p^Valuation(l,p) if p le 1 then return false; end if; x:=1; while l mod (p*x) eq 0 do x*:=p; end while; return x; end function; GetBiggestLB := function(l) for i in [1..#l] do if l[i][1] ne 0 then return l[i][1],l[i][2]; end if; end for; return 0,0; end function; // Application to modular forms intrinsic MFCongruences (N1 :: RngIntElt, a1 :: RngIntElt, N2 :: RngIntElt, a2 :: RngIntElt : HB := 0, factor := true, Precomputed := false) -> Any {Compute all congruences of the forms (N1,2,a1) and (N2,2,a2), i.e. the a1-th form (according to Magma's internal newspace decomposition) for level Ni and weight 2. The output is a list with entries , meaning that p^lb is a lower bound and p^ub an upper bound. Out of laziness only for Gamma_0 and weight 2 at the moment.} require (([N1,a1] ne [N2,a2]) or a1 eq 0) : "The modular forms must be different."; /* Precomputed := true : it uses/creates/modifies the following files: nforms.txt, bounds.txt and N-i.txt (where N and i are integers) */ // For the sake of simplicity, we suppose (N1,a1) > (N2,a2), and the list will be sorted: if (N1 lt N2) or (N1 eq N2 and a1 lt a2 and a1 ne 0) then aux:=N1; N1:=N2; N2:=aux; aux:=a1; a1:=a2; a2:=aux; end if; //Definition of a format to be called in procedures. FormsFormat:= recformat ; f := rec ; F:=0; // Initialize F and G with irrelevant values just to use g := rec ; G:=0; // them in some functions without being strictly defined NFormsFormat:= recformat ; // Also definition of a structure to use it with not NForm := rec ; // necessarily defined elements // Compute and sort the forms. Since it takes time, we do it just if we do not have it previously computed if not Precomputed then f`M := ModularSymbols(f`N,2,1); f`D := SortDecomposition(NewformDecomposition(NewSubspace(CuspidalSubspace(f`M)))); f`d := #f`D; if f`N ne g`N then // If N1 eq N2, we will later make g=f g`M := ModularSymbols(g`N,2,1); g`D := SortDecomposition(NewformDecomposition(NewSubspace(CuspidalSubspace(g`M)))); g`d := #g`D; end if; end if; // compute the HECKE BOUND if HB eq 0 then if Precomputed then StoredHB(~f,~g); //Compute f`H and g`H else f`H:=HeckeBound(f`M); if f`N ne g`N then // If N1 eq N2, we will later make g=f g`H:=HeckeBound(g`M); end if; end if; end if; if f`N eq g`N then // To save time, we didn't compute things twice, if they are equal g:=f; end if; HB:=Maximum(f`H,g`H); // END computation HB OUT_LIST := []; OUT_LIST_FAC := []; aux := false; //aux will tell us if we compute a new element of NForms, in order to store it later //Declare ai_ini and ai_end, deppending on the loops we want to make (i.e. if some ai eq 0) if a1 ne 0 then a1_ini:=a1; a1_end:=a1; end if; if a2 ne 0 then a2_ini:=a2; a2_end:=a2; end if; if N1 ne N2 and a1 eq 0 then //in this case, we need to know the number of newforms we will compare: NewformsOfLevelN(~f, ~NForm, ~aux); a1_ini:=1; a1_end:=f`d; end if; if N1 ne N2 and a2 eq 0 then //we need to know the number of newforms we will compare: NewformsOfLevelN(~g, ~NForm, ~aux); a2_ini:=1; a2_end:=g`d; end if; if N1 eq N2 and a1 eq 0 then //we need to know the number of newforms we will compare: NewformsOfLevelN(~f, ~NForm, ~aux); a1_end:=f`d; if a2 eq 0 then a1_ini:=2; else ai_ini:=a2+1; end if; end if; //END declare ai_ini and ai_end if aux then //if we computed new information, we store it IntegerSequenceToFile(NForm`s,"nforms.txt"); end if; //Initialize some values to call them in some functions in which they will be computed just if necessary aux:=false; Pf:=0; Pg:=0; file1:=0; file2:=0; s1:=0; s2:=0; //Start loop for all a1 we want for ia1 in [a1_ini..a1_end] do if assigned f`D then // if we computed f`D, we can assign our modular form F F := f`D[ia1]; end if; //Declare a2_ini and a2_end in a special case if N1 eq N2 and a2 eq 0 then a2_ini:=1; a2_end:=ia1-1; end if; //End declare a2_ini and a2_end //Start loop for all a2 we want for ia2 in [a2_ini..a2_end] do if assigned g`D then G := g`D[ia2]; // if we computed g`D, we can assign our modular form G end if; aux1:=Precomputed; aux2:=Precomputed; //Open "N-i.txt" files if Precomputed then s1:=IntegerToString(N1)*"-"*IntegerToString(ia1)*".txt"; s2:=IntegerToString(N2)*"-"*IntegerToString(ia2)*".txt"; file1 := Open(s1,"a"); file1 := Open(s1,"r"); file2 := Open(s2,"a"); file2 := Open(s2,"r"); end if; // entries are of the form , // p1,p2 are polys, c,r,s are the output of congruence number, p is the current prime: work_list := []; work_list_p_divide_N_1 := []; work_list_p_divide_N_2 := []; cp:=0; p_ini := 2; p_end := HB; while cp eq 0 do //this while is to avoid cases in which HB is not enough for p in [i : i in [p_ini..p_end] | IsPrime(i)] do //get all charpolys ComputeHecke(~f, ~F, ~aux1, ~file1, s1, ia1, Precomputed, p, ~Pf, ~aux); if aux and f`N eq g`N then g:=f; aux := false; end if; ComputeHecke(~g, ~G, ~aux2, ~file2, s2, ia2, Precomputed, p, ~Pg, ~aux); // compute congruence numbers and store them if Degree(Gcd(Pf,Pg)) eq 0 then // Common factors not interesting c,r,s := CongruenceNumber(Pf,Pg); // 1st step of algorithm is just for p not dividing the levels if ((N1 mod p) ne 0 and (N2 mod p) ne 0) then Append(~work_list, ); //Compute already the gcd: if the gcd is already one, we do not have to waste time computing more congruence numbers if #work_list eq 2 then // c gives no information at the prime p: cp:=Gcd(work_list[1][3]*pPart(c,work_list[1][6]),c*pPart(work_list[1][3],p)); elif #work_list gt 2 then cp:=Gcd(cp,c*pPart(cp,p)); end if; if cp eq 1 then //if cp eq 1, trivially we do not need to continue break; end if; elif p le HB then //Theorem 3.8 will use also p's dividing the levels Append(~work_list_p_divide_N_1, ); if N1 mod N2 eq 0 then rr:=Valuation(Integers()!(N1/N2),p); if rr eq 0 then Append(~work_list_p_divide_N_2, ); else delta:=0; if N2 mod p ne 0 then delta:=-p; end if; //print "Pf = ", Pf, "Pg = ", Pg, "Pgtilda = ", Parent(Pg).1^(Degree(Pg)*(rr-1))*HeckeInLevelNm(Pg,delta); Append(~work_list_p_divide_N_2, ); end if; end if; end if; end if; end for; //if HB is not enough, we continue: if cp eq 0 then p_ini:=p_end+1; p_end:=2*p_end; print "HB was not enough"; end if; end while; UB:=cp; //With cp we already have an upper bound with all possible primes of congruence: cp := [fac[1] : fac in Factorisation(cp)]; out_list := []; for p in cp do //we apply the algorithm for every different prime: //ubs are computed with step 1. lbs with steps 2 and 3: lbs1 := []; lbs2:=[]; lbs3:=[]; ubs := []; for w in work_list do number_of_tries := 0; succeeded := false; if assigned here then delete(here); end if; while (not succeeded) and (number_of_tries lt 1000) do number_of_tries := number_of_tries + 1; try here := Bounds(w[1],w[2], p : factor := factor, crs := ); succeeded := true; catch e print N1, ia1, N2, ia2; end try; end while; if not assigned here then continue ia2; end if; lb, ub := GetBiggestLB(here); if w[6] le HB then Append(~lbs1, ); // Theorem 3.8 needs all primes but just le HB Append(~lbs2, ); // Proposition 3.11 needs all primes but just le HB if lb le 1 then vprint Cong,1: "p nmid N : ia1 = ", ia1, " ia2 = ", ia2, " l= ",p," lb = ",lb, " p = ", w[6]; end if; end if; if w[6] ne p then Append(~ubs, Minimum(ub,Valuation(UB,p))); end if; end for; for w in work_list_p_divide_N_1 do // Theorem 3.8 needs also the primes dividing the levels: number_of_tries := 0; succeeded := false; if assigned here then delete(here); end if; while (not succeeded) and (number_of_tries lt 1000) do number_of_tries := number_of_tries + 1; try here := Bounds(w[1],w[2], p : factor := factor, crs := ); succeeded := true; catch e print N1, ia1, N2, ia2; end try; end while; if not assigned here then continue ia2; end if; lb, _ := GetBiggestLB(here); Append(~lbs1, ); if lb le 1 then vprint Cong,1: "Th3.8 : ia1 = ", ia1, "ia2 = ", ia2, "p= ",p,"lb = ",lb, "p = ", w[6] ; end if; end for; for w in work_list_p_divide_N_2 do // Proposition 3.11 needs also the primes dividing the levels: number_of_tries := 0; succeeded := false; if assigned here then delete(here); end if; while (not succeeded) and (number_of_tries lt 1000) do number_of_tries := number_of_tries + 1; try here := Bounds(w[1],w[2], p : factor := factor); succeeded := true; catch e print N1, ia1, N2, ia2; end try; end while; if not assigned here then continue ia2; end if; lb, _ := GetBiggestLB(here); if lb lt Infinity() then Append(~lbs2, ); if lb le 1 then vprint Cong,1: "Prop3.11 : ia1 = ", ia1, "ia2 = ", ia2, "l= ",p,"lb = ",lb, " p = ", w[3], w[3] ne p; end if; end if; end for; //lbs1,lbs2; if N1 mod N2 eq 0 then if Minimum([i[1] : i in lbs1 | i[2]]) gt Minimum([i[1] : i in lbs2 | i[2]]) then print "Th3.8 useful also for Ng|Nf in LB"; end if; if Minimum([i[1] : i in lbs1]) gt Minimum([i[1] : i in lbs2]) then print "Th3.8 useful also for Ng|Nf in LBL"; end if; lb := Maximum(Minimum([i[1] : i in lbs1 | i[2]]),Minimum([i[1] : i in lbs2 | i[2]])); lbl := Maximum(Minimum([i[1] : i in lbs1]),Minimum([i[1] : i in lbs2])); lb3 := Maximum(Minimum([i[1] : i in lbs1 | i[2] and i[4]]),Minimum([i[1] : i in lbs2 | i[2] and i[4]])); else lb := Minimum([i[1] : i in lbs1 | i[2]]); lbl := Minimum([i[1] : i in lbs1]); lb3 := Minimum([i[1] : i in lbs1 | i[2] and i[4]]); end if; //lbs1, lbs2; //if lbl eq 0 and N1 ne N2 and N1 mod N2 eq 0 then print "lbl == 0 -> Ribet Raising level could be useful?"; end if; Append(~out_list, ); end for; LBL :=1; LB:=1; LB3:=1; UB:=1; for i in [1..#out_list] do LBL*:=out_list[i][1]^out_list[i][2]; LB*:=out_list[i][1]^out_list[i][3]; LB3*:=out_list[i][1]^out_list[i][4]; UB*:=out_list[i][1]^out_list[i][5]; end for; if #out_list gt 0 then Append(~OUT_LIST_FAC,[*[N1,ia1,N2,ia2],out_list*]); // Any better way to give the results? Append(~OUT_LIST, [N1,ia1,N2,ia2,LBL,LB,LB3,UB]); end if; end for; end for; // return OUT_LIST, OUT_LIST_FAC; return OUT_LIST; end intrinsic;