//CountEmbeddings.mag //Magma code to count the number of embeddings (up to isomorphism) of the ring of integers of a quartic CM field into the endomorphism ring of a product of elliptic curves //This code was written for the paper "Igusa class polynomials, embeddings of quartic CM fields, and arithmetic intersection theory" by Helen Grundman, Jennifer Johnson-Leung, Kristin Lauter, Adriana Salerno, Bianca Viray, and Erika Wittenborn. The authors make no claims about running time of the algorithm. /********************************************************************/ //Given a list of 4 rational numbers, q1, q2, q3, q4 //return the minimal positive integer N such that //Nq1, Nq2 Nq3, Nq4 are integers function Denom(L) denominators := []; for i := 1 to 4 do //hardcoded the length Append(~denominators, Denominator(L[i])); end for; return LCM(denominators); end function; /********************************************************************/ //Check whether the element x is in the integral span of elements //in the basis B function InIdeal(x, B) A := Matrix([Coordinates(B[1]), Coordinates(B[2]), Coordinates(B[3]), Coordinates(B[4])]); v := Vector(Coordinates(x)); s := Solution(A,v); retval := true; for n := 1 to 4 do if not IsIntegral(s[n]) then retval := false; end if; end for; return retval; end function; /********************************************************************/ //Return the Rosati involution of M. //Notice this depends on the Norm of the ideal I function Rosati(M,NI) return Matrix([ [ Conjugate(M[1][1]), Conjugate(M[2][1])*NI], [ Conjugate(M[1][2])/NI, Conjugate(M[2][2])]]); end function; /********************************************************************/ //Find all elements in the integral span of the basis B with //norm = Normx and trace = Tracex function FixedNormAndTrace(Normx, Tracex, B) xvals := []; Q := Parent(B[1]); isq := -Norm(Q.1); jsq := -Norm(Q.2); A := Transpose( Matrix([ Coordinates(B[1]), Coordinates(B[2]), Coordinates(B[3]), Coordinates(B[4])])); denom1 := Denom(A[2]); denom2 := Denom(A[3]); denom3 := Denom(A[4]); x0 := Tracex/2; temp3 := Normx - x0^2; if temp3 ge 0 then for num3 := Ceiling(-Sqrt(Rationals()!(temp3*denom3^2/(isq*jsq)))) to Floor(Sqrt(temp3*denom3^2/(isq*jsq))) do temp2 := temp3 - (num3/denom3)^2*isq*jsq; for num2 := Ceiling(-Sqrt(temp2*denom2^2/(-jsq))) to Floor(Sqrt(temp2*denom2^2/(-jsq))) do temp1 := temp2 - (num2/denom2)^2*(-jsq); for num1 := Ceiling(-Sqrt(temp1*denom1^2/(-isq))) to Floor(Sqrt(temp1*denom1^2/(-isq))) do x := x0 + num1/denom1*Q.1 + num2/denom2*Q.2 + num3/denom3*Q.3; if (Norm(x) eq Normx) and (Trace(x) eq Tracex) and InIdeal(x, B) then Append(~xvals, x); end if; end for; end for; end for; end if; return xvals; end function; /********************************************************************/ function RemoveConjElts(solns, E, alpha0, alpha1, beta0, beta1, D) n := 1; while n le #solns do // Loop over automorphisms of E for r in FixedNormAndTrace(1, 0, E[1]) cat FixedNormAndTrace(1, 1, E[1]) cat FixedNormAndTrace(1, -1, E[1]) cat FixedNormAndTrace(1, 2, E[1]) cat FixedNormAndTrace(1, -2, E[1]) do // Loop over automorphisms of E' for s in FixedNormAndTrace(1, 0, E[3]) cat FixedNormAndTrace(1, 1, E[3]) cat FixedNormAndTrace(1, -1, E[3]) cat FixedNormAndTrace(1, 2, E[3]) cat FixedNormAndTrace(1, -2, E[3]) do //Calculating NormI so that we can compute Rosati involution I := LeftIdeal(QuaternionOrder(E[1]), E[2]); NormI := Norm(I); U := Matrix([[r,0],[0,s]]); Uinv := Matrix([[Conjugate(r), 0],[0, Conjugate(s)]]); //Remove embedding which is conjugate to given embedding by U // First check that U does not fix embedding if not (solns[n][1]*U eq U*solns[n][1] and solns[n][2]*U eq U*solns[n][2]) then Exclude(~solns, [U*solns[n][1]*Uinv, U*solns[n][2]*Uinv]); end if; //Remove embedding which is conjugate to the complex conjugate //of given embedding by U // First check that U does not fix complex conjugate of // embedding if not (solns[n][1]*U eq U*solns[n][1] and solns[n][2]*U eq U*Rosati(solns[n][2], NormI)) then size := #solns; Exclude(~solns, [U*solns[n][1]*Uinv, U*Rosati(solns[n][2], NormI)*Uinv]); end if; //if E and E' are isomorphic then remove embedding which are // conjugate and have the order of E, E' switched if E[1] eq E[2] then U := Matrix([[0, r], [s, 0]]); Uinv := Matrix([[0, Conjugate(s)], [Conjugate(r), 0]]); if U*Uinv ne Matrix([[1,0],[0,1]]) or Uinv*U ne Matrix([[1,0],[0,1]]) then print "Error - Uinv incorrect!!", r, s; end if; if not (solns[n][1]*U eq U*solns[n][1] and solns[n][2]*U eq U*solns[n][2]) then size := #solns; Exclude(~solns, [U*solns[n][1]*Uinv, U*solns[n][2]*Uinv]); end if; if not (solns[n][1]*U eq U*solns[n][1] and solns[n][2]*U eq U*Rosati(solns[n][2], NormI)) then size := #solns; Exclude(~solns, [U*solns[n][1]*Uinv, U*Rosati(solns[n][2], NormI)*Uinv]); end if; end if; end for; end for; n := n + 1; end while; return solns; end function; /********************************************************************/ // find Lambda1, Lambda2 satsifying embedding conditions and such that // entries have fixed norm function FindSolns(E, alpha0, alpha1, beta0, beta1, D, Norms12, Normt12, Normt22, Tracet22, Normt11, Tracet11, s11); solns := []; R1 := QuaternionOrder(E[1]); R2 := QuaternionOrder(E[3]); I := LeftIdeal(R1, E[2]); // Create list of possible values in I with Norm = Norms12 s12traces := [-Floor(Sqrt(4*Norms12))..Floor(Sqrt(4*Norms12))]; s12vals := []; for T in s12traces do s12vals := s12vals cat FixedNormAndTrace(Norms12, T, E[2]); end for; for s12 in s12vals do // loop over values in R1 with Norm = Normt11 and Trace = Tracet11 for t11 in FixedNormAndTrace(Normt11, Tracet11, E[1]) do // loop over values in R2 with Norm = Normt22 and Trace = Tracet22 for t22 in FixedNormAndTrace(Normt22, Tracet22, E[3]) do // Create list of possible values in I with Norm = Normt12 t12traces := [-Floor(Sqrt(4*Normt12))..Floor(Sqrt(4*Normt12))]; t12vals := []; for T in t12traces do t12vals := t12vals cat FixedNormAndTrace(Normt12, T, E[2]); end for; for t12 in t12vals do t21 := 1/Norm(I)*(alpha1*Conjugate(s12) - Conjugate(t12)); //define Lambda1 = omega and Lambda2 = eta omega := Matrix([[s11, s12], [Conjugate(s12)/Norm(I), D - s11]]); eta := Matrix([[t11, t12], [t21, t22]]); //check omega and eta live in correct subring // and satisfy all of the embedding conditions if InIdeal(s12, E[2]) and InIdeal(t12, E[2]) and (omega*eta eq eta*omega) and eta + Rosati(eta, Norm(I)) eq alpha0 + alpha1*omega and eta*Rosati(eta, Norm(I)) eq beta0 + beta1*omega then Append(~solns, [omega, eta]); //add to list of solutions end if; end for; end for; end for; end for; return solns; end function; /********************************************************************/ // Lists all possible pairs of matrices Lambda1, Lambda2 in R_{E, E'} // such that // Lambda1*Lambda2 = Lambda2*Lambda1 // Lambda1^2 - D*Lambda1 + 1/4*(D^2 - D) = 0 // Lambda2 + Lambda2^{\vee} = alpha0 + alpha1*Lambda1 // Lambda2*Lambda2^{\vee} = beta0 + beta1*Lambda1 function ListAllEmb(alpha0, alpha1, beta0, beta1, D, p) // Define constants associated to CM field F := QuadraticField(D); D := Discriminant(F); trw := D; //Trace(omega) nw := Integers()!(1/4*(D^2 - D)); //Norm(omega) Dtilde := 1/16*D^4*alpha1^4 + 1/2*D^3*alpha0*alpha1^3 - 1/8*D^3*alpha1^4 - D^3*alpha1^2*beta1 + 3/2*D^2*alpha0^2*alpha1^2 - 1/2*D^2*alpha0*alpha1^3 - 4*D^2*alpha0*alpha1*beta1 + 1/16*D^2*alpha1^4 - 2*D^2*alpha1^2*beta0 + D^2*alpha1^2*beta1 + 4*D^2*beta1^2 + 2*D*alpha0^3*alpha1 - 1/2*D*alpha0^2*alpha1^2 - 4*D*alpha0^2*beta1 - 8*D*alpha0*alpha1*beta0 + 4*D*alpha0*alpha1*beta1 - 2*D*alpha1^2*beta0 + 16*D*beta0*beta1 - 4*D*beta1^2 + alpha0^4 - 8*alpha0^2*beta0 + 16*beta0^2; // Give a presentation of B_{p, infinity} if p ne 2 then jsq := -p; isq := -1; while (IsSquare(FiniteField(p)!isq)) or (#RamifiedPrimes(QuaternionAlgebra) gt 1) do isq := isq - 1; end while; else jsq := -1; isq := -1; end if; QQ := Rationals(); Q := QuaternionAlgebra; ksq := -isq*jsq; // Orders = { O : O maximal order isomorphic to End(E) for E supersingular // curve over Fpbar} // Note: if E and E' are not isomorphic, but have isomorphic endomorphism // rings, Orders contains two different maximal orders, one is identified // with End(E), the other with End(E') OO := MaximalOrder(Q); LI := LeftIdealClasses(OO); Orders := []; for I in LI do Append(~Orders, RightOrder(I)); end for; // Ends = list of triples, one triple for each (ordered) pair (E, E') // Let O be a maximal order which is isomorphic to End(E), I a left ideal // of O which can be identified with Hom(E', E), and O' the right order of I // which is identified with End(E'). Then the triple associated to (E, E') // is (Basis(O), Basis(I), Basis(O')). The correspondence between // (O, I, O') and (E, E') is explained more in the paper. Ends := []; for O in Orders do for I in LeftIdealClasses(O) do Append(~Ends, [Basis(O), Basis(I), Basis(RightOrder(I))]); end for; end for; // initialize list of solutions solns := []; number := 0; //look for Lambda1, Lambda2 in R_{(E, E')} satisfying above conditions // this is equivalent to searching for s11, t11 in O, s12, t12 in I, // s22, t22 in O'. Then // Lambda1 = // [s11 s12] // [s12^{\vee} s22] // Lambda2 = // [t11 t12] // [t12^{\vee} t22] for E in Ends do NormI := Norm(LeftIdeal(QuaternionOrder(E[1]), E[2])); //possible values for s11. s22 = D - s11 s11vals := [Ceiling( (D - Sqrt(D))/2 )..Floor( (D + Sqrt(D))/2)]; Esolns := []; for s11 in s11vals do a := s11; //remnants of old notation delta := -a^2 + a*trw - nw; m := delta; Tracet11 := alpha0 + alpha1*s11; //conditions from Lambda + Lambda* Tracet22 := alpha0 + alpha1*(D - s11); for n := Ceiling(-Sqrt(Max(m^2*Dtilde - 4*D, 0))) to Floor(Sqrt(Max(m^2*Dtilde - 4*D, 0))) do //parts of old notation Normx := 1/(2*D)*(n - (1/4*a^2*D^2*alpha1^2 - 1/4*a*D^3*alpha1^2 + 1/16*D^4*alpha1^2 - 1/4*a^2*D*alpha1^2 + 1/4*a*D^2*alpha1^2 - 1/8*D^3*alpha1^2 + a^2*D*alpha1*alpha0 - a*D^2*alpha1*alpha0 + 1/4*D^3*alpha1*alpha0 + 1/16*D^2*alpha1^2 - 1/4*D^2*alpha1*alpha0 + a^2*alpha0^2 - a*D*alpha0^2 + 1/4*D^2*alpha0^2 - 2*a^2*D*beta1 + 2*a*D^2*beta1 - 1/2*D^3*beta1 - 1/4*D*alpha0^2 - 2*a*D*beta1 + 1/2*D^2*beta1 - 4*a^2*beta0 + 4*a*D*beta0 - D^2*beta0 - D*beta0)); Normu := delta*(beta0 + beta1*a) - delta*Normx; Tracexuc := beta1*delta - (D - 2*a)/delta*Normu; Normv := delta^2*Normx + delta*(D - 2*a)*Tracexuc + (D - 2*a)^2*Normu; //conditions coming from fixed n Normt11 := Normx; Normt12 := Normu*NormI/delta; Normt22 := Normv/delta^2; Norms12 := NormI*delta; //check n can give a solution which is integral if IsIntegral((m^2*Dtilde - n^2)/(4*D*p)) and IsIntegral(Normt11) and IsIntegral(Normt12) and IsIntegral(Normt22) then Esolns := Esolns cat FindSolns(E, alpha0, alpha1, beta0, beta1, D, Norms12, Normt12, Normt22, Tracet22, Normt11, Tracet11, s11); // finds all solutions with // entries of Lambda2 having // fixed norm and fixed // trace (for some entries) end if; end for; end for; // remove elements which are conjugate Esolns := RemoveConjElts(Esolns, E, alpha0, alpha1, beta0, beta1, D); if E[1] eq E[2] then number := number + #Esolns; else number := number + #Esolns/2; end if; Append(~solns, Esolns); end for; return solns, Ends, Integers()!number; end function;