// Time-stamp: "Last modified: 2021-05-19 16:23:57 (d_yasaki)" //level := 49;sqfreeL :=[7];dimbound := 2; //level := 32;sqfreeL :=[2];dimbound := 2; //level := 32;sqfreeL :=[6];dimbound := 2; //level := 32;sqfreeL :=[10];dimbound := 2; //level := 45;sqfreeL :=[5];dimbound := 8; //level := 35;sqfreeL :=[5];dimbound := 6; //level := 35;sqfreeL :=[5]; //level := 40;sqfreeL :=[5]; dimbound := 8; //level := 50;sqfreeL :=[5]; dimbound := 4; //level := 27;sqfreeL :=[3]; dimbound := 2; //level := 30;sqfreeL :=[5]; //level := 15;sqfreeL :=[5]; //level := 50;sqfreeL :=[2,3,5,6,7,10]; //level := 5; sqfreeL :=[7,10]; //level := 81; level := 169;sqfreeL :=[3]; dimbound := 24; //level := 169;sqfreeL :=[5]; dimbound := 24; //level := 169;sqfreeL :=[6]; dimbound := 24; //level := 169;sqfreeL :=[7]; dimbound := 24; //level := 169;sqfreeL :=[10]; dimbound := 24; //level := 121;sqfreeL :=[2]; dimbound := 16; //level := 121;sqfreeL :=[10]; dimbound := 16; //level := 121;sqfreeL :=[7]; dimbound := 16; //level := 121;sqfreeL :=[6]; dimbound := 16; //level := 121;sqfreeL :=[10]; dimbound := 16; //level := 73;sqfreeL :=[5]; //level := 83;sqfreeL :=[2,3,5,6,7,10]; //level := 89;sqfreeL :=[2,3,5,6,7,10]; //level := 97;sqfreeL :=[2,3,5,6,7,10]; //level := 97; fname := "data.txt"; homologybasering := RationalField(); SetSeed(3); function getho(n) fname := Sprint(n)*"-ho.txt"; L := Split(Read(fname),"\n"); levelstring := Split(L[1],":"); level := StringToInteger(levelstring[1]); dim := StringToInteger(levelstring[2]); ho1 := []; ho2 := []; prime1list := []; prime2list := []; for S in L[2..#L] do data := Split(S,":"); type := StringToInteger(data[2]); p := StringToInteger(data[1]); A := Matrix(QQ,dim,eval(data[3])); if type eq 1 then Append(~prime1list,p); Append(~ho1, A); else assert type eq 2; Append(~prime2list,p); Append(~ho2, A); end if; end for; num := Minimum({#prime1list,#prime2list}); assert prime1list[1..num] eq prime2list[1..num]; return ho1[1..num], ho2[1..num], prime1list[1..num]; end function; ho1,ho2,primes := getho(level); cellrec:= recformat< vertices : SeqEnum, stabilizer : SeqEnum, stabilizerplus : SeqEnum, facetdata : SeqEnum, mover, barycenter, type : Integers(), orientationstd : Integers(), facets : SeqEnum >; load "gl3.txt"; function q(v) // return the cone coordinates and symmetric matrix of the rank 1 // form associated to v = [a1, ..., an] ans := []; for i in [1..#v] do for j in [i..#v] do Append(~ans,v[i]*v[j]); end for; end for; return ans; end function; Qstandardorder := [ q(a) : a in standardorder ]; Hstandardorder := [Matrix(1,v)*Matrix(1,#v,v) : v in standardorder]; function standardorder_position(v) // return position of v in standardorder i := Position(Qstandardorder,q(v)); assert i ne 0; return i; end function; function matrix_stabilizer(A) G := AutomorphismGroup(LatticeWithGram(A)); G := Setseq(Set(G)); if SL eq true then G := [g: g in G | Determinant(g) eq 1]; end if; return G; end function; function matrix_is_conjugate(A,B) // return Boolean and g such that gAg^* = B, if g exists. tf,g := IsIsometric(LatticeWithGram(A),LatticeWithGram(B)); if SL and tf then if Determinant(g) eq 1 then return true, g; else flippers := [ h : h in AutomorphismGroup(LatticeWithGram(A)) | Determinant(h) eq -1]; if #flippers eq 0 then return false, _; else return true, g*flippers[1]; end if; end if; end if; if tf then return tf,g; else return tf,_; end if; end function; function bary(L) // given list of indices, return the barycenter (as a symmetric // matrix) of the cell defined by the vectors b := &+[Hstandardorder[a] : a in L]; return b; end function; function stabilizer(L) // given a list of indices, return the stabilizer of the cell // defined by those vectors return matrix_stabilizer(bary(L)); end function; function GL_type(cell, Vor) // given a list of indices and list of standard cells L, return the // type and g such that g*standard is cell // cell is list of indices, Vor is cell records for lower // dimensional cells type := 0; if not IsPositiveDefinite(bary(cell)) then print "Degenerate cell is type 0."; return 0,_; else for i in [1..#Vor] do found,g := matrix_is_conjugate(Vor[i]`barycenter,bary(cell)); if found then type := i; break i; end if; end for; assert found; return type, g; end if; end function; // projective space over ZZ/n // 7/10/12 // may be slow, but should work over ZZ for various n // Updated Dec 1, 2020 to use associative arrays function prime_power_projective_space(n,p,k) // given n, prime p, and k, return the projective space P^n // over ZZ/p^k. Should not be called on its own. // it is used in the projective_space function call. ER := [0..p^k - 1]; U := [ e : e in ER | GCD(e,p) eq 1 ]; N := [ n : n in ER | not n in U ]; Pn := []; for i in [1..n + 1] do // make leftmost invertible into 1. for v in CartesianPower(N,n + 1 - i) do for w in CartesianPower(ER, i - 1) do Append(~Pn, [a: a in v] cat [1] cat [a: a in w]); end for; end for; end for; return Pn; end function; function vector_crt(L,fact) // helper function to do the CRT pk := [ a[1]^a[2] : a in fact ]; ans := []; for i in [1..#L[1]] do num := [v[i] : v in L]; Append(~ans,CRT(num,pk)); end for; return ans; end function; function projective_space(n,level) // given n and level create the projective space P^n over ZZ/level // returns a list of representative points in ZZ^(n + 1) and a map // which takes a sequence of n+1 coprime integers and returns a // representative. fact := Factorization(level); crt := CartesianProduct([prime_power_projective_space(n,a[1],a[2]) : a in fact]); Pn := AssociativeArray(); Pn_ind := AssociativeArray(); U := [i mod level: i in [1..level] | GCD(i,level) eq 1]; stab := AssociativeArray(); for i in [k : k in [1..level - 1] | not k in U] do stab[i] := [s : s in U | s mod (level div GCD(level,i)) eq 1]; end for; stab[0] := U; function representative(v) // if v contains a unit, scale to make leftmost a 1 and // return. v := [a mod level : a in v]; for i in [ 1..#v ] do g,s,t := XGCD(v[i],level); if g eq 1 then return [ s*a mod level : a in v ]; end if; end for; // if it gets here then there are no units in v. // make first nonzero as small as possible by multiplying by unit. iszero := true; i := 0; while iszero and i lt #v do i +:= 1; if v[i] ne 0 then m,k := Minimum([s*v[i] mod level : s in U]); u := U[k]; v := [ u*a mod level : a in v ]; iszero := false; end if; end while; // now use stabilizers to make small in dictionary order currentstab := stab[v[i]]; while i lt #v and #currentstab gt 1 do i +:= 1; m,k := Minimum([s*v[i] mod level : s in currentstab]); u := currentstab[k]; v := [ u*a mod level : a in v ]; currentstab := [a: a in stab[v[i]] | a in currentstab]; end while; return v; end function; i := 1; for v in crt do w := representative(vector_crt(v,fact)); if not IsDefined(Pn_ind,w) then Pn_ind[w] := i; Pn[i] := w; i +:= 1; end if; end for; function index(v) // return the position of v in Pn ind := Pn_ind[representative(Eltseq(v))]; assert not IsZero(ind); return ind; end function; return Pn, representative, index, Pn_ind; end function; function EulerPsi(level) // the number of points in the projective line. phi := 1; for tup in Factorization(level) do NP := tup[1]; phi *:= (NP+1) * NP^(tup[2]-1); end for; return phi; end function; function opp_trans(A); // given matrix A, return backwards transpose of A (about opposite diagonal m := NumberOfRows(A); n := NumberOfColumns(A); B := ZeroMatrix(Integers(),n,m); for i in [1..n] do for j in [1..m] do B[i,j] := A[m + 1 - j, n + 1 - i]; end for; end for; return B; end function; function gammify(v) // given bottom row v, complete to a matrix in SL_n(ZZ) m := opp_trans(Matrix(#v,v)); H,T := HermiteForm(m); // T sends column v to column e_1. Ti := T^-1; // sends e_1 to v. g := opp_trans(Ti); d := Determinant(g); if d eq -1 then for i in [1..NumberOfColumns(g)] do g[1,i] := -1*g[1,i]; end for; end if; assert Determinant(g) eq 1; assert Eltseq(g[NumberOfRows(g)]) eq v; return g; end function; function number_of_points_in_projective_space(n,level) // computes the number of points in projective space n-space of Z/mZ // the number of points in the projective line. phi := 1; for tup in Factorization(level) do p := tup[1]; k := tup[2]; phi *:= &+[p^(n*k - i) : i in [0..n]]; end for; return phi; end function; /****************************************************************************** homology with level for GL_n/QQ assumes data will be loaded. ******************************************************************************/ function factorization_type(I) f := Factorization(I); return Sort([a[2] : a in f]); end function; function sum(L) if #L eq 0 then return 0; else return &+L; end if; end function; function index_action(g,i) // return j such that g acting on vertex i is vertex j return standardorder_position(Eltseq(g*Matrix(Integers(),1,standardorder[i]))); end function; // Each cell is a list (unordered) of indices. Always use Sort. function standard_simplex_indices(L) // given a list L of indices, return the indices in the // standardorder that give the standard simplex (choose from left // to right, adding vertices that increase rank. Sort(~L); ind := [L[1]]; m := Rank(Matrix([ Qstandardorder[a] : a in L ])); counter := 2; A := Matrix([ Qstandardorder[a] : a in ind ]); n := Rank(A); while n lt m do while Rank( Matrix([ Qstandardorder[a] : a in ind cat [L[counter]] ])) eq n do counter := counter + 1; end while; Append(~ind,L[counter]); A := Matrix([ Qstandardorder[a] : a in ind ]); n := Rank(A); end while; return ind; end function; function complement(L) // given indices L, return vectors that span a complement // of L I := IdentityMatrix(Rationals(),#Qstandardorder[1]); RI := [Eltseq(a) : a in Rows(I)]; A := Matrix([ Qstandardorder[i] : i in L]); n := Rank(A); ind := []; m := #Qstandardorder[1]; counter := 1; while n lt m do if Rank( Matrix([ Qstandardorder[a] : a in L] cat [ RI[i] : i in ind cat [counter] ])) gt n then Append(~ind,counter); end if; n := Rank( Matrix([ Qstandardorder[a] : a in L] cat [ RI[i] : i in ind cat [counter] ])); counter := counter +1; end while; return [ RI[i] : i in ind ]; end function; function preserves_orientation(L, g) // given list of indices, and a stabilizer group, return if g preserves // the orientation of the cell defined by L is orientable. ssi := standard_simplex_indices(L); ext := complement(ssi); A1 := Matrix([ Qstandardorder[i] : i in ssi] cat ext); newind := [standardorder_position(Eltseq(g*Matrix(Integers(),1,standardorder[i]))) : i in ssi]; A2 := Matrix([ Qstandardorder[i] : i in newind] cat ext); det := Determinant(A1)*Determinant(A2); if det lt 0 then return false; else assert det gt 0; // if it reaches here, then g preserves orientation return true; end if; end function; function ind_rank(L) return Rank(Matrix([ Qstandardorder[i] : i in L])); end function; function relative_orientation(cell, face) // given list of indices, return relative orientation. // note: the face is gotten from the Herbert method of standard // simplex. // assumes cell is std. assumes face is simplicial and ordered. std_cell := standard_simplex_indices(cell); //extra := [a : a in std_cell | not a in face]; //assert #extra eq 1; //extra := [Minimum([a : a in std_cell | not a in face])]; extra := [Minimum([a : a in std_cell | ind_rank(Append(face,a)) eq #std_cell])]; ext := complement(std_cell); A1 := Matrix([ Qstandardorder[i] : i in face cat extra] cat ext); A2 := Matrix([ Qstandardorder[i] : i in std_cell ] cat ext); d := Determinant(A1)*Determinant(A2); assert d ne 0; s := Sign(d); return s; end function; /****************************************************************************** Below here the level is fixed, and we have a projective space. ******************************************************************************/ orbitdata := recformat; printf "Creating projective space for level %o.\n", level; printf "Factorization type %o.\n",factorization_type(level); Pn,stdrep,projective_space_position, Pn_ind := projective_space(matrank-1,level); printf "Projective space has size %o.\n", #Pn; function projective_space_orbits(mats, plusmats); // compute orbits and orientation numbers N := #Keys(Pn); if #mats eq 2 then assert #plusmats eq 2; orbit_inds := [i : i in [1..N]]; orientation_numbers := [1 : i in [1..N]]; else // indices, in range [1 .. number of orbits] orbit_inds := [0 : i in [1..N]]; orientation_numbers := [0 : i in [1..N]]; k := 1; // index of the orbit being computed while exists(i){i : i in [1..N] | orbit_inds[i] eq 0} do orbit_inds[i] := k; orientation_numbers[i]:= 1; for M in mats do ii :=projective_space_position(Vector(Pn[i])*M); error if ii eq 0, "Something wrong with choice of standard reps for Pn."; orbit_inds[ii] := k; if M in plusmats then eps := 1; else eps := -1; end if; orientation_numbers[ii] := eps; end for; k +:= 1; end while; end if; return orbit_inds, orientation_numbers; end function; function orientable_orbits(c) // given a cell c, return the orientable orbits of GL-type c. // this corresponds to orientable Gamma_0-cells. G := c`stabilizer; Gp := c`stabilizerplus; orbits, orientationnum := projective_space_orbits(G, Gp); orbit_ind := Setseq(Set(orbits)); if #G eq #Gp then orientable_ind := orbit_ind; else reversers := [ g : g in G | not g in Gp ]; orientable_ind := orbit_ind; for i in orbit_ind do I := Position(orbits,i); a0 := Pn[I]; if exists(g) {g : g in reversers | projective_space_position(Vector(a0)*g) eq I} then Exclude(~orientable_ind, i); end if; end for; end if; O := rec; return O; end function; /****************************************************************************** Below here the cells are made. ******************************************************************************/ function get_facet_data(cell,Vor) // return cell with facet data attached newcell := cell; G := matrix_stabilizer(newcell`barycenter); Gp := [g : g in G |preserves_orientation(newcell`vertices,g)]; newcell`stabilizer := G; newcell`stabilizerplus := Gp; fd := []; for face in cell`facets do type, mover := GL_type(face,Vor); // use mover to get std face that face is conjugate to. if type ne 0 then simpstdface := standard_simplex_indices(Vor[type]`vertices); orientedface := [index_action(mover, k) : k in simpstdface]; o := relative_orientation(cell`vertices, orientedface); f := rec; Append(~fd, f); end if; end for; newcell`facetdata := fd; return newcell; end function; /****************************************************************************** Below here the computations start ******************************************************************************/ print "Computing orbits."; orbits := []; for i in [1..3] do print i; orbits[i] := [orientable_orbits(c): c in V[i]]; end for; //orbits := [[orientable_orbits(c): c in A] : A in V]; function boundary(dim) // gives the boundary map from n cells to n - 1 cells. // assumes everything is defined and set-up correctly. print "Computing the boundary map from",dim,"cells to",dim-1,"cells."; n := dim - matrank + 2; if n eq 1 then num_domaincells := sum([#a`orientableind : a in orbits[1]]); num_rangecells := 0; return SparseMatrix(homologybasering,num_rangecells, num_domaincells); end if; if n eq #V + 1 then num_domaincells := 0; num_rangecells := sum([#a`orientableind : a in orbits[#V]]); return SparseMatrix(homologybasering,num_rangecells, num_domaincells); end if; if n gt #V + 1 or n lt 1 then return SparseMatrix(homologybasering,0,0); end if; S := V[n]; T := V[n - 1]; S0 := orbits[n]; T0 := orbits[n - 1]; num_domaincells := sum([#a`orientableind : a in S0]); num_rangecells := sum([#a`orientableind : a in T0]); boundary_matrix := SparseMatrix(homologybasering,num_rangecells, num_domaincells); print "The matrix is",num_rangecells,"X",num_domaincells; // first fix a GL-type for sigma column := 1; for sgltype in [1..#S] do facetdata := S[sgltype]`facetdata; // record of GL-type // and mover matrix standard to there. // now fix a Gamma0-type cells of particular GL-type for sg0type in S0[sgltype]`orientableind do a0ind := Position(S0[sgltype]`orbits, sg0type); // this is the Pn index of // the a0 representative for sigma a0 := Pn[a0ind]; for L in facetdata do tgltype := L`type; assert tgltype ne 0; //c := L[2]; c := L`orientationstd; a := Vector(a0)*Matrix(matrank,matrank,L`mover); aind := projective_space_position(a); // Find a in correct boundary component of GL-type // tgltype tg0type := T0[tgltype]`orbits[aind]; b0ind := Position(T0[tgltype]`orbits, tg0type); if Position(T0[tgltype]`orientableind,tg0type) ne 0 then // Boundary face is orientable. // Find delta that sends it to b0 for that boundary // +1 if delta is orientation preserving // -1 otherwise. // use orientationnumbers eps := T0[tgltype]`orientation[aind]; // orientation number // to get the row index: row := sum([#T0[i]`orientableind : i in [1.. tgltype - 1]]) + Position( T0[tgltype]`orientableind, tg0type ); //print "Writing", c*eps,"to",[row, // column],"entry of boundary map."; boundary_matrix[row, column] := boundary_matrix[row, column] + c*eps; end if; end for; // bump to next column. if column mod 1000 eq 0 then print "Computing column",column,"of",num_domaincells; end if; column := column + 1; end for; end for; return boundary_matrix; end function; function homology(d2,d3) assert IsZero(d2*d3); r2 := Rank(d2); E3 := ElementaryDivisors(d3); e3 := [a: a in E3 | not a eq 1]; r3 := #E3; homfree := NumberOfColumns(d2) - r3 - r2; return ; end function; function cycle_basis_lifts() vecs := []; // these are bottom rows. for i in orbits[1,1]`orientableind do ind := Position(orbits[1,1]`orbits,i); assert orbits[1,1]`orientation[ind] eq 1; Append(~vecs,Pn[ind]); end for; mats := [gammify(v) : v in vecs]; return mats; end function; // Reduction of 0-sharbly [v1, v2, v3] function midrep(n,m) // return n mod m such that |n| <= m/2 //assert m gt 0; n0 := n mod m; if n0 ge m/2 then return n0 - m; else return n0; end if; end function; function reduce(Q) // deals with columns. L := []; m := Determinant(Q); if m in {1, -1} then return [Q], true; else E := EchelonForm(Q); if E[1,1] gt 1 then m := E[1,1]; t := [1,0,0]; elif E[2,2] gt 1 then //assert E[1,1] eq 1; m := E[2,2]; t := [midrep(-E[1,2],m),1,0]; else //assert E[1,1] eq 1; //assert E[2,2] eq 1; m := E[3,3]; t := [midrep(E[1,2]*E[2,3] - E[1,3],m),midrep(-E[2,3],m),1]; end if; v := [ZZ|&+[t[j]/m*Q[i,j]: j in [1..3]]:i in [1..3]]; // put v in jth column, j = 1, 2, 3 for j in [1..3] do P := Q; P[1,j] := v[1]; P[2,j] := v[2]; P[3,j] := v[3]; //if not IsZero(Determinant(P)) then if not t[j] eq 0 then Append(~L,P); end if; end for; end if; return L, false; end function; function fullreduce(Q) // reduce Q until norm determinant 1 // NOTE: Q is made of columns reducedL := []; todo := [Q]; while #todo gt 0 do L,done := reduce(todo[1]); Remove(~todo,1); if done then reducedL cat:= L; else todo cat:= L; end if; //print #todo, {Determinant(a) : a in todo}; end while; return reducedL; end function; function hecke_matrices2(p) // given prime p, return the upper triangular matrices that // act on the cusps. Taken from Paul's p. 227 of the Appendix. // assumes n = 3. // Computes for T(p,2) reps := [0..p - 1]; L := [ Matrix(3, [ p, 0, 0, 0, p, 0, 0, 0, 1 ]) ] cat [ Matrix(3,[ 1, a, b, 0, p, 0, 0, 0, p ]) : a,b in reps ] cat [ Matrix(3,[ p, 0, 0, 0, 1, a, 0, 0, p ]) : a in reps]; return L; end function; function hecke_matrices1(p) // given prime p, return the upper triangular matrices that // act on the cusps. Taken from Paul's p. 227 of the Appendix. // assumes n = 3. // Computes for T(p,1) reps := [0..p - 1]; L := [ Matrix(3,[ p, 0, 0, 0, 1, 0, 0, 0, 1 ]) ] cat [ Matrix(3,[ 1, 0, a, 0, 1, b, 0, 0, p ]) : a, b in reps ] cat [ Matrix(3,[ 1, a, 0, 0, p, 0, 0, 0, 1 ]) : a in reps]; return L; end function; function d_reps() //d runs over a set of representatives of the double cosets G\SL(3,Z)/P_3, where P_3 is the stabilizer in SL(3,Z) of (0,0,\pm 1) (acting on the right of the row vectors.) todo := [[0,1,0]] cat [v : v in Keys(Pn_ind) | IsZero(v[2])]; L := []; while #todo gt 0 do v := todo[1]; Append(~L,v); Remove(~todo,1); Exclude(~todo,stdrep([v[1],v[2],-v[3]])); end while; d := [gammify(v) : v in L]; di := [a^-1 : a in d]; return d, di; end function; time d2 := boundary(2); time d3 := boundary(3); time H0,mu := quo; d_list, di_list := d_reps(); function mat2H0(tup) // given c = coefficient, A = 3 x 3 matrix of columns, return vector in H0 A := tup[2]; v := Domain(mu)!0; Pn_ind := projective_space_position(A[3]); orbit_ind := Position(orbits[1,1]`orientableind,orbits[1,1]`orbits[Pn_ind]); if not IsZero(orbit_ind) then c := tup[1] * orbits[1,1]`orientation[Pn_ind]; v[orbit_ind] := c; end if; return mu(v); end function; function homologyaction(L) // returns a matrix which acts on ROWS! clifts := cycle_basis_lifts(); M := []; for e in Basis(H0) do v := e@@mu; terms := [ : i in Support(v),g in L]; Lv := &+[&+[mat2H0() : A in fullreduce(t[2])] : t in terms]; Append(~M,Lv); end for; return Matrix(M); end function; function makeW(m,u) // return W = [m, u\\ 0, e], e = \pm det(m), m is unital, and u in ZZ^2 W := Matrix(3,[m[1,1], m[1,2], u[1], m[2,1], m[2,2], u[2], 0, 0, Determinant(m)]); return W; end function; function next_unital(sqfree,pm : power := 0) F := QuadraticField(sqfree); ep := FundamentalUnit(F); d0L := []; while #d0L eq 0 do power +:= 1; epk := pm * ep^power; t := ZZ!Trace(epk); n := ZZ!Norm(epk); // first find roots of f mod level for d0 in Integers(level) do if IsZero(d0^2 - t*d0 + n) then Append(~d0L,ZZ!d0); end if; end for; end while; ML := []; for d0 in d0L do // now shift d0L to dL so that f(d) is small (easier to factor?) k := (t div 2 - d0) div level; left := d0 + k*level; dL := [left, left + level]; for d in dL do fd := d^2 - t*d + n; assert IsZero(fd mod level); fdN := fd div level; for c in Divisors(Abs(fdN)) do cN := c*level; M := Matrix(2,[t - d,(-fdN) div c, cN,d]); N := M; N[1,2] := -M[1,2]; N[2,1] := -M[2,1]; assert Determinant(M) eq n; assert Determinant(N) eq n; Append(~ML,M); Append(~ML,N); end for; end for; end for; bL := []; for A in ML do beta := (epk- A[2,2])/A[2,1]; Append(~bL,beta); end for; return ML,bL,power; end function; function unital(sqfree,pm,power,dbound) F := QuadraticField(sqfree); ep := FundamentalUnit(F); epk := pm * ep^power; t := ZZ!Trace(epk); n := ZZ!Norm(epk); ML := []; for d in [-dbound..dbound] do k := d^2 - t*d + n; num := AbsoluteValue(k); if num ne 1 then fac,comp := TrialDivision(num:Proof := false); if #comp ne 0 then assert #comp eq 1; fac := fac cat []; assert &*[tup[1]^tup[2]:tup in fac] eq num; end if; car := CartesianProduct([[0..a[2]] : a in fac]); else fac := [<1,1>]; car := [[1]]; end if; for tup in car do b := 1; for i in [1..#fac] do b *:= fac[i][1]^tup[i]; c := k div (-b); M := Matrix(2,[t - d,b,c,d]); N := M; N[1,2] := -M[1,2]; N[2,1] := -M[2,1]; assert Determinant(M) eq n; assert Determinant(N) eq n; Append(~ML,M); Append(~ML,N); end for; end for; end for; bL := [a[1,2]/(a[2,2] - Conjugate(epk)) : a in ML]; return ML,bL; end function; function g02g(A) // return g = A^k, where k is smallest positive integers such that g is in Gamma0(level) g := A; inGamma := (IsZero(g[3,1] mod level) and IsZero(g[3,2] mod level)); power := 1; while not inGamma do power +:= 1; g *:= A; inGamma := IsZero(g[3,1] mod level) and IsZero(g[3,2] mod level); end while; return g,power; end function; function normalize(v) // scale to make entries of v be integers d := LCM([Denominator(a) : a in v]); return [v[i] * d : i in [1..#v]]; end function; function psi_mat(m,u,i) //[f,a,gf]_G is in the image of psi, and they generated the image, at least if we work over Q. W := makeW(m,u); d := d_list[i]; g0 := d_list[i]*W*di_list[i]; g,power := g02g(g0); // a is in Q^3 with eigenvalue e for operator g e := Determinant(m); E := Nullspace(Transpose(g0) - e*IdentityMatrix(QQ,3)); assert Dimension(E) eq 1; da := normalize(Eltseq(Basis(E)[1])); // f = de_1 (where e_1 is the standard column vector (1,0,0).) f := [d[1,1], d[2,1], d[3,1]]; gf := Eltseq(g*Matrix(1,f)); A := Transpose(Matrix(ZZ,[f,gf,da])); //print "power",power; return A; end function; function psi_image(A) frA := fullreduce(A); //print Determinant(A),#frA; RA := &+[mat2H0(<1,a>) : a in frA]; return RA; end function; power := 1; //L,bL,power := next_unital(5,1: power := power); function test(sqfree,pm,power,u) V := []; L,bL := unital(sqfree,pm,power,5); for j in [1..#L] do for i in [1..#d_list] do //print [j,i]; A := psi_mat(L[j],u,i); v := psi_image(A); //if not IsZero(v) then //print i, j, bL[j],Determinant(A),v; //end if; Include(~V,v); end for; end for; return V; end function; function check(sqfree,dbound:prange := [1..10], dimbound := Dimension(H0)) vL := []; seenS := {}; P := sub; rank := 0; nextrank := 0; for power in prange do for pm in {1,-1} do for ppm in {1,-1} do L,bL := unital(sqfree,pm,ppm*power,dbound); print ppm*power, rank; for j in [1..#L] do for i in [1..#d_list] do for u in [[0,0],[1,0],[0,1],[1,1],[-1,0],[0,-1]] do //for u in [[0,0],[1,0],[0,1],[1,1],[-1,0],[0,-1],[-1,-1],[1,-1],[-1,1]] do A := psi_mat(L[j],u,i); //if #(Sprint(Determinant(A))) lt 0 then v := psi_image(A); if not v in seenS then tempvL :=[v] cat [v*T : T in ho1 cat ho2]; for w in tempvL do Include(~seenS,v); end for; nextrank := Rank(Matrix(vL cat tempvL)); end if; if nextrank gt rank then print u; print ppm*power, rank, #vL, i; P := sub; rank := nextrank; vL := Basis(P); dataL := [*level,sqfree,Dimension(H0),Dimension(P),Eltseq(EchelonForm(BasisMatrix(P)))*]; data := MyShort(dataL); Write("data3-tmp.txt",data); //if rank eq Dimension(H0) then if rank eq dimbound then return P; end if; end if; //end if; end for; end for; end for; end for; end for; end for; return P; end function; function random_check(sqfree,dbound:prange := [1..10], dimbound := Dimension(H0)) vL := []; seenS := {}; P := sub; rank := 0; nextrank := 0; while true do for pm in {1,-1} do for ppm in {1,-1} do for i in [1..#d_list] do power:=Random(prange); L,bL := unital(sqfree,pm,ppm*power,dbound); print ppm*power, rank; for j in [1..#L] do for u in [[0,0],[1,0],[0,1],[1,1],[-1,0],[0,-1]] do //for u in [[0,0],[1,0],[0,1],[1,1],[-1,0],[0,-1],[-1,-1],[1,-1],[-1,1]] do A := psi_mat(L[j],u,i); //if #(Sprint(Determinant(A))) lt 0 then v := psi_image(A); if not v in seenS then tempvL :=[v] cat [v*T : T in ho1 cat ho2]; for w in tempvL do Include(~seenS,v); end for; nextrank := Rank(Matrix(vL cat tempvL)); end if; if nextrank gt rank then print u; print ppm*power, rank, #vL, i; P := sub; rank := nextrank; vL := Basis(P); print u; L := [*level,sqfree,Dimension(H0),Dimension(P),Eltseq(EchelonForm(BasisMatrix(P)))*]; data := MyShort(L); Write("data3-tmp.txt",data); //if rank eq Dimension(H0) then if rank eq dimbound then return P; end if; end if; //end if; end for; end for; end for; end for; end for; end while; return P; end function; if Dimension(H0) ne 0 then for sqfree in sqfreeL do //P := check(sqfree,10:prange:=[3..12]); //P := check(sqfree,20:prange:=[1..15], dimbound := dimbound); //P := check(sqfree,20:prange:=[1..15]); //P := check(sqfree,20:prange:=[15..25], dimbound := dimbound); P := check(sqfree,10:prange:=[2..10], dimbound := dimbound); L := [*level,sqfree,Dimension(H0),Dimension(P),Eltseq(EchelonForm(BasisMatrix(P)))*]; data := MyShort(L); Write(fname,data); end for; end if; /* primes := [ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97 ]; fname := Sprint(level)*"-ho.txt"; Write(fname,MyShort([level,Dimension(H0)])); ho1 := []; ho2 := []; for p in primes do if GCD(p,level) eq 1 then T1 := homologyaction(hecke_matrices1(p)); Append(~ho1,T1); Write(fname,MyShort([*p,1,Eltseq(T1)*])); T2 := homologyaction(hecke_matrices2(p)); Append(~ho2,T2); Write(fname,MyShort([*p,2,Eltseq(T2)*])); end if; end for; */ quit;