############################################################## # Scripts for local j-invariant and Grobner fan computations # ############################################################## # The following functions allow us perform several tasks needed for certifying a vertex comes from a component with good reduction, and in addition it computes potential odd primes that behave badly with respect to the local computation of leading terms (with respect to a given weight) of a given polynomial. # In particular: # 1) Generate M2 input to compute leading terms with respect to a weight vector # 2) Compute Leading terms for all primes except finitely many # 3) Compute potentially problematic primes # 4) Compute how the calculation on the leading terms must be modified for the problematic odd primes. # 5) Compute the Grobner fan of a polynomial with +1/-1 coefficients or over a field with trivial valuation. # 6) Subdivide a polyhedral cone according to a Grobner fan. ######################### # Generate the M2 input # ######################### # We create the M2 scripts to compute the leading terms of jC1Num and jC1Denom with respect to a trivial valuation on Q, we pick any weight vector in the Type V range. e.g. (36, 16, 16 9, 9, 1). # Rring = ambient ring for the polynomials # inputjCNumName = name of the variable assigned to 'inputjCNum' # inputjCDenomName = name of the variable assigned to 'inputjCDenom' # inputjCNum, inpujCDenom = name of the 2 factors we want to use, e.g. jC1Num or jC1Denom. They are assumed to be polynomials in Rring. # WeightsAs = a list of the weights of the variables in Rring (we need it to be in the Type V / Type IV canonical regions, respectively) # m2FileName = string, e.g. 'TypeV/jC1TypeV.m2' # ExplanationString = string, e.g. '-- Calculations for the left-most vertex of Type V calculations' # The output is an m2 file that computes the leading terms of both polynomials with respect to the given weight vector, assuming the residue field has characteristic 0. def writeM2Code(Rring, inputjCNumName, inputjCDenomName, inputjCNum, inputjCDenom, WeightsAs, m2FileName, ExplanationString): f = file(str(m2FileName), 'w') f.write('-- ' + str(ExplanationString) + '\n'+'\n') f.write('S=QQ['+str(R.gens())[1:-1] +', MonomialOrder=>{Weights => {'+str(WeightsAs)[1:-1]+'}},Global=>false];'+'\n'+'\n') f.write('-- We write the Numerator of the j-invariant' +'\n') f.write(str(inputjCNumName) + '=' + str(inputjCNum) +';'+'\n'+'\n') f.write('-- We write the Denominator of the j-invariant' +'\n' ) f.write(str(inputjCDenomName) + '=' + str(inputjCDenom)+';'+'\n'+'\n') f.write('-- We compute the leadingTerms of the Numerator' +'\n') f.write('toString(leadTerm(1,'+str(inputjCNumName)+'))'+'\n'+'\n') f.write('-- We compute the leadingTerms of the Denominator' +'\n') f.write('toString(leadTerm(1,'+str(inputjCDenomName)+'))'+'\n'+'\n') f.close() ################################# # Calculation of Bad Odd Primes # ################################# # For a given polynomial with rational coefficients, we compute which primes (other than 2) will lead to a different leading term when taking the p-adic valuation than the one computed with trivial valuation on the field QQ. # Rring is the ambient ring # polyToTest = polynomial in Rring to be tested # The script takes all the coefficients on 'polyToTest' and writes all the odd primes that divide some of the coefficients: def computeBadPrimesOdd(Rring, polyName): coefjC = list(set(Rring(polyName).coefficients())) len(coefjC) print 'All coefficients Factored:' print [factor(ZZ(x)) for x in coefjC] print 'All coefficients Factors with odd primes and their exponents:' print [[y for y in factor(ZZ(x)) if y[0]!=2] for x in coefjC] BadPrimesjC=[] for x in coefjC: Sprimesx=sorted(list(set([y[0] for y in factor(ZZ(x)) if y[0]!=2]))) for y in Sprimesx: if y not in BadPrimesjC: BadPrimesjC.append(y) BadPrimesjC = sorted(BadPrimesjC) return BadPrimesjC ############################################################## # Separating monomials with a fixed coefficient (up to sign) # ############################################################## # In what follows we decompose a polynomial 'polyToTest' as a sum (indexed by its coefficients) of polynomials. Each summand consists of all the monomials in 'polyToTest' with the same given coefficients. # The output is a dictionary with entries equal to a list: # s[0] = coefficient (with plus sign) # s[1] = sum of a monomials in jC1Num (resp. jC1Denom) with that coefficient # The script records the index of the coefficient list with a given value. The order of the monomials and the coefficients match, so the indices will tell us which monomials to use. We record the information in the dictionary 'indicesjC'. def separateMonomialsAndTheBadOddPrimes(Rring, polyToTest): # We record all coefficients and all monomials in polyToTest coefficientsjC = Rring(polyToTest).coefficients() absValueCoefficients = set([abs(x) for x in coefficientsjC]) monomialsjC= Rring(polyToTest).monomials() # We record all indices of the list of coefficients with a given value indicesjC ={} for j in coefficientsjC: indicesjC[j]=[k for k in range(len(coefficientsjC)) if coefficientsjC[k]==j] # We record the sum of all monomials with the given coefficient: polynomialsjCByCoefficients = {} for s in absValueCoefficients: p = [c for c in coefficientsjC if c == s] m = [d for d in coefficientsjC if d == -s] if p !=[]: pluspoly = sum([monomialsjC[k] for k in indicesjC[p[0]]]) else: pluspoly = 0 if m !=[]: negpoly = sum ([monomialsjC[k] for k in indicesjC[m[0]]]) else: negpoly = 0 polynomialsjCByCoefficients[s] = Rring(pluspoly)-Rring(negpoly) return polynomialsjCByCoefficients ############################################### # Auxiliary function to compute Leading Terms # ############################################### # We create a function to take LT of a polynomial given a weight vector. The weight vector has as many entries as the variables of the ambient polynomial ring, and they are given in the same order. The calculation assumes all the coefficients of the polynomial have valuation equal to zero. def LTFromWeights(polyP,weightVector, Rring): monPieces = polyP.monomials() Coeff = polyP.coefficients() piecesExp = [x.exponents()[0] for x in monPieces] InnerProducts=[vector(ZZ,piecesExp[k])*vector(ZZ,weightVector) for k in range(len(monPieces))] maxValueInnerProducts = max(InnerProducts) indicesForLT = [k for k in range(len(InnerProducts)) if InnerProducts[k]==maxValueInnerProducts] LTPiece = sum([Rring(monPieces[k])*Rring(Coeff[k]) for k in indicesForLT]) return LTPiece ################################################### # Computing Leading Terms and one exponent vector # ################################################### # The script takes a dictionary with polynomial values (in the ring Rring) and computes its leading term with respect to a weight vector. In case the coefficients are all 1/-1, the leading term computation is independent of the valuation of the field. def LeadingTermsFromDictionary(Rring, dictPolys, weightVector): LeadingTermsjCByCoefficients = {} for j in dictPolys.keys(): LeadingTermsjCByCoefficients[j] = LTFromWeights(Rring(dictPolys[j]),weightVector, Rring) return LeadingTermsjCByCoefficients # This script takes a dictionary of homogeneous polynomial with respect to a given weight (i.e. a leading term) and computes the exponent of its first monomial. def ExpLeadingTermFromDictionary(Rring, dictPolys): ExpLeadingTermjCByCoefficients ={} for j in dictPolys.keys(): ExpLeadingTermjCByCoefficients[j] = vector(Rring(dictPolys[j]).exponents()[0]) return ExpLeadingTermjCByCoefficients ########################### # Coefficients and primes # ########################### # We record a list of two dictionaries with the same keys (integer coefficients of a polynomial). # 1) the first dictionary takes a key and returns a list of all the bad primes for a given key # 2) the second dictionary takes a key and list of extra negative valuations that we must add for each prime in the value of the first dictionary. def ComputeBadPrimesAndExtraValuationsToAdd(keysCoefficients): CoefficientsjCBadPrimes = {} ExtraValuationToAddjC = {} for j in keysCoefficients: CoefficientsjCBadPrimes[j] = [y[0] for y in factor(ZZ(j)) if y[0]!=2] for j in keysCoefficients: if CoefficientsjCBadPrimes[j] == []: ExtraValuationToAddjC[j] = [0] else: # Since we are working with negative valuations, we must consider the negative of the exponent. ExtraValuationToAddjC[j] = [-y[1] for y in factor(ZZ(j)) if y[0]!=2] return (CoefficientsjCBadPrimes,ExtraValuationToAddjC) ######################################################################## # Comparing the values for each one of the primes and picking winners: # ######################################################################## # Computes a dictionary of dictionaries. The keys of the dictionary are odd bad primes. The value for each key is the actual p-valuation for all the BadPrimes considering a weight vector and an exponent of each leading term and the extra valuations coming from the coefficients (the keys of each entry). It is important to remark that we Always use the same weight vector def newValuationsPerLT(BadPrimesjC, CoefficientsjCBadPrimes, ExtraValuationToAddjC, ExpLeadingTermjCByCoefficients, weightVector): newValuationjC = {} for p in BadPrimesjC: newValuationjC[p] ={} for j in CoefficientsjCBadPrimes.keys(): if p in CoefficientsjCBadPrimes[j]: extraValuation = ExtraValuationToAddjC[j][CoefficientsjCBadPrimes[j].index(p)] else: extraValuation = 0 totalValuationsExponents = vector(ZZ,weightVector)*vector(ZZ,ExpLeadingTermjCByCoefficients[j]) newValuationjC[p][j] = ZZ(extraValuation) + ZZ(totalValuationsExponents) return newValuationjC # The following script picks the winners among all entries of the dicionary values obtained by the previous script. def pickWinners(dictionaryOfvaluationsPerPrimeToCompare): maxValuesjC = {} for p in dictionaryOfvaluationsPerPrimeToCompare.keys(): coefjC = dictionaryOfvaluationsPerPrimeToCompare[p].keys() maxValue = max(set(dictionaryOfvaluationsPerPrimeToCompare[p].values())) maxValuesjC[p] = [j for j in coefjC if dictionaryOfvaluationsPerPrimeToCompare[p][j] == maxValue] return maxValuesjC ################################################################################ # Reconstruct the terms on each leading term piece maximizing each p-valuation # ################################################################################ # We reconstruct the terms in LT decomposition realizing the maximal value for each prime and each coefficient where the maximal value is achieved for the given weight vector: def LTPerPrimeAndMaxCoefficients(maxValuesjC,newValuationjC,CoefficientsjCBadPrimes,ExtraValuationToAddjC, LeadingTermsjCByCoefficients): MaxTermsjC = {} for p in newValuationjC.keys(): termsToAdd = [] for j in maxValuesjC[p]: if p in CoefficientsjCBadPrimes[j]: extraValuation = ExtraValuationToAddjC[j][CoefficientsjCBadPrimes[j].index(p)] else: extraValuation = 0 # print str(extraValuation) termsToAdd.append((j,extraValuation,LeadingTermsjCByCoefficients[j])) MaxTermsjC[p] = termsToAdd return MaxTermsjC ########################### # Combining all functions # ########################### def trueValuationsPerPrimeAndWeightVector(Rring, polyToTest,WeightVector, BadPrimesjC): polynomialsjCByCoefficients = separateMonomialsAndTheBadOddPrimes(Rring, polyToTest) LeadingTermsjCByCoefficients = LeadingTermsFromDictionary(Rring, polynomialsjCByCoefficients, WeightVector) ExpLeadingTermjCByCoefficients = ExpLeadingTermFromDictionary(Rring,LeadingTermsjCByCoefficients) (CoefficientsjCBadPrimes,ExtraValuationToAddjC) = ComputeBadPrimesAndExtraValuationsToAdd(polynomialsjCByCoefficients.keys()) newValuationjC = newValuationsPerLT(BadPrimesjC, CoefficientsjCBadPrimes, ExtraValuationToAddjC, ExpLeadingTermjCByCoefficients, WeightVector) maxValuesjC = pickWinners(newValuationjC) MaxTermsjC= LTPerPrimeAndMaxCoefficients(maxValuesjC,newValuationjC,CoefficientsjCBadPrimes,ExtraValuationToAddjC, LeadingTermsjCByCoefficients) return (polynomialsjCByCoefficients, LeadingTermsjCByCoefficients, ExpLeadingTermjCByCoefficients,CoefficientsjCBadPrimes,ExtraValuationToAddjC,newValuationjC,newValuationjC,MaxTermsjC) ###################################################################### # Calculation of Grobner fans for polynomials with 1/-1 coefficients # ###################################################################### # Since we don't know if a given region lies withing a single Gr"obner cone, we need to check this is the case for the decomposition of a polynomial by their coefficients. # To simplify our calculations, we compute alternative polynomials that only contain extremal monomials def fromPolynomialsToExponentVectors(polyToTest): exponentVectorsQInZZ = [] monQ = polyToTest.monomials() coeffQ = polyToTest.coefficients() exponentVectorsQ = [x.exponents() for x in monQ] m = len(exponentVectorsQ[0][0]) for x in exponentVectorsQ: v = [ZZ(x[0][j]) for j in range(0,m)] exponentVectorsQInZZ.append(v) return exponentVectorsQInZZ # The calculations of Grobner fans only require the vertices of the Newton polytope of Q. The number of vertices of this polytope is 366 compared to the 14111 monomials in Q. The computation of the Grobner fan of the full Q does not terminate, as opposed to that of the alternative polynomial AltQ. def fromPolynomialsToExtremalMonomials(polyToTest,Rring): exponentVectorsQInZZ = fromPolynomialsToExponentVectors(polyToTest) NPQ = Polyhedron(vertices= exponentVectorsQInZZ) extremalExponents = [tuple(list(x)) for x in NPQ.vertices()] AltQMonomials = [exponentVectorsQInZZ.index(list(x)) for x in extremalExponents] monQ = polyToTest.monomials() coeffQ = polyToTest.coefficients() AltQ = Rring(0) for j in AltQMonomials: AltQ = AltQ + Rring(coeffQ[j]*monQ[j]) return AltQ def GrobnerFanComputations(polyToTest): Iq = ideal(polyToTest) GFq = Iq.groebner_fan() GFqFan = GFq.polyhedralfan() return GFqFan # We record the maximal cones with the lineality space (it is assumed to be the all ones vector) # GFan is a polyhedral fan (precomputed with the script "GrobnerFanComputations" above) def MaximalCones(GFan): RaysGFqfan = GFan.rays() # We record the maximal cones of the fan. Each cone is described by the indices of the extremal rays spanning these cones. Such indices refer to the list of 76 rays: MaximalConesPolyAlt = GFan.maximal_cones().values()[0] PolyMaximalCones = {} for cone in MaximalConesPolyAlt: allOnes = [1 for k in range(len(RaysGFqfan[0]))] negAllOnes = [-1 for k in range(len(RaysGFqfan[0]))] linSp = [allOnes,negAllOnes] # print cone raysOfCone = [RaysGFqfan[i] for i in cone] + linSp polyCone = Polyhedron(rays= raysOfCone,base_ring=QQ) PolyMaximalCones[MaximalConesPolyAlt.index(cone)] = polyCone return PolyMaximalCones # We record the maximal cones with the lineality space data. # GFan is a polyhedral fan (precomputed with the script "GrobnerFanComputations" above) # lines is a list of generators of the lineality space. def MaximalConesLines(GFan,Lines): RaysGFqfan = GFan.rays() # We record the maximal cones of the fan. Each cone is described by the indices of the extremal rays spanning these cones. Such indices refer to the list of 76 rays: MaximalConesPolyAlt = GFan.maximal_cones().values()[0] PolyMaximalCones = {} for cone in MaximalConesPolyAlt: raysOfCone = [RaysGFqfan[i] for i in cone] # + linSp polyCone = Polyhedron(rays= raysOfCone, lines = Lines, base_ring=QQ) PolyMaximalCones[MaximalConesPolyAlt.index(cone)] = polyCone return PolyMaximalCones ############################################################# # Subdivision of a polyhedral cone induced by a Grobner fan # ############################################################# # The following script computes the intersection of a polyhedral cone with all maximal cells in the Grobner fan of a polynomial jC. def IntersectionsTypeConejCGrobnerCone(TypeCone,PolyMaximalConesjC): IntersectionsjCTypeCone = {} for i in range(0,len(PolyMaximalConesjC.keys())): # print str(i) polyCone =PolyMaximalConesjC[i] # print polyCone IntersectionCone = TypeCone.intersection(polyCone) IntersectionsjCTypeCone[i] = IntersectionCone return IntersectionsjCTypeCone # The following script takes all possible intersections and computes the ones of maximal dimension. def pickIndicesMaximalIntersections(TypeCone, IntersectionsjCTypeCone): m = max(set([dim(x) for x in IntersectionsjCTypeCone.values()])) maxIntersections = [i for i in IntersectionsjCTypeCone.keys() if dim(IntersectionsjCTypeCone[i])==m] return maxIntersections # The following script combines the above scripts and returns the cone index of a maximal cones in jC intersecting the TypeCone and the actual slices. def computeSubdivisionsOfTypeConeByGrobnerFan(TypeCone,PolyMaximalConesjC): IntersectionsjCTypeCone = IntersectionsTypeConejCGrobnerCone(TypeCone,PolyMaximalConesjC) MaxIntersections = pickIndicesMaximalIntersections(TypeCone, IntersectionsjCTypeCone) # print str(MaxIntersections) SubdivisionsOfTypeConeByjC = [IntersectionsjCTypeCone[i] for i in MaxIntersections] return (MaxIntersections,SubdivisionsOfTypeConeByjC) # The following script computes the intersections of a given dimension: def pickIndicesIntersectionsGivenDimension(TypeCone, IntersectionsjCTypeCone,m): maxIntersections = [i for i in IntersectionsjCTypeCone.keys() if dim(IntersectionsjCTypeCone[i])==m] return maxIntersections # The following script takes a list of pieces subdividing a given cone, and computes the facets of the subdivision. The output is given as a tuple (the keys and the values of the intersections) def computeFacets(PiecesCones): n = dim(PiecesCones[1][0]) print n m = len(PiecesCones[1]) IntersectionsOfPieces = {} allButk= {} pickIndices = {} for k in range(0,m): print str(k) if dim(PiecesCones[1][k]) !=n: print 'The family is not equidimensional!' return False allButk[k] = {j: PiecesCones[1][k].intersection(PiecesCones[1][j]) for j in range(0,m) if j!=k} pickIndices[k] = pickIndicesIntersectionsGivenDimension(PiecesCones[1][k],allButk[k],n-1) # print pickIndicesk relevantPairs = [(k,j) for k in range(0,m-1) for j in range(k,m) if j in pickIndices[k]] for p in relevantPairs: IntersectionsOfPieces[p] = allButk[p[0]][p[1]] return (IntersectionsOfPieces.keys(), IntersectionsOfPieces.values()) ############################ # Merging identical facets # ############################ # The output of computeFacets, when applied to lower dimensional cones might lead to repetitions. If that's the case, we must re-order the data (merging indices giving the same facet). The following script performs this task. def purgeFacets(listOfIndicesAndPolyhedra): (indicesToMerge, polyhedraToMerge) = listOfIndicesAndPolyhedra allPolyhedra = list(set(polyhedraToMerge)) allIndices ={} for P in range(len(allPolyhedra)): indicesForP = [x for x in range(len(polyhedraToMerge)) if polyhedraToMerge[x] == allPolyhedra[P]] print indicesForP facetsForP = [] for x in indicesForP: for y in list(indicesToMerge[x]): if y not in facetsForP: facetsForP.append(y) allIndices[P] = facetsForP return (allIndices, allPolyhedra) #################################### # Certifying no subdivisions occur # #################################### # The following script takes a rational polyhedral cone C in R^n and the Grobner fan of a polynomial in n variables and computes the subdivision of the polyhedral cone induced by the Grobner Cone. The output is the pair of indices of maximal dimensions and a variable (True or False). The variable is 'True' if C agrees with the non-empty intersections between C and all the maximal cones in the Grobner cone. def checkNoSubdivisions(TypeCone,PolyMaximalCones): (MaxIndices,SubdivisionsOfCByMaximalCones) = computeSubdivisionsOfTypeConeByGrobnerFan(TypeCone,PolyMaximalCones) checkInclusion = all([SubdivisionsOfCByMaximalCones[i]==TypeCone for i in range(len(SubdivisionsOfCByMaximalCones))]) return (MaxIndices, checkInclusion) ################################################## # Constructing representative points from Slices # ################################################## # We record the list of rays of a slice in a Cone in M_2^trop, the polyhedron they generate, and a point in the relative interior of each cone. The latter will allow us to compute the initial form of polynomials induced by the weight given by each of these points. def findRaysPiecesSamplePoint(subdivisionsOfTypeCone): raysOfTypeConeGfan = {} TypeConeSubdividedGfan = {} pointTypeConeGfan = {} for i in range(len(subdivisionsOfTypeCone)): raysOfTypeConeGfan[i] = [list(x) for x in subdivisionsOfTypeCone[i].rays()] TypeConeSubdividedGfan[i] = Polyhedron(rays=raysOfTypeConeGfan[i], lines = subdivisionsOfTypeCone[i].lines()) pointTypeConeGfan[i]= list(2*sum([vector(r) for r in raysOfTypeConeGfan[i]])) return (raysOfTypeConeGfan,TypeConeSubdividedGfan,pointTypeConeGfan) # The following function computes a vector with non-positive coefficients from a given one by substracting a multiple of the all-ones vector and thereafter computes its primite vector def makePrimitive(rayPrelim): if list(rayPrelim) == len(rayPrelim)*[0]: print 'we have the origin' return list(rayPrelim) Q = Polyhedron(rays = [rayPrelim]) return list(Q.rays()[0].vector()) def makeNegativeAndPrimitive(intList): k = max([k for k in intList]) l =len(intList) intVector = vector(ZZ,intList) allOnes = vector(ZZ,l*[1]) rayPrelim = tuple(list(intVector - k*allOnes)) return makePrimitive(rayPrelim) ############################################## # Computation of leading terms for each cone # ############################################## def computeLTSample(weightsAs,R,Q): return LTFromWeights(R(Q),weightsAs,R) ################################################################ # Extract last n exponents of each monomial and give their sum # ################################################################ def ProjectMonomials(poly,Rxy,R): allMonomials = list(set(R(poly).monomials())) m = len(allMonomials[0].exponents()[0]) print m n = len(Rxy.gens()) print n projectedExponents = list(set([tuple(list(x.exponents()[0])[m-n:m]) for x in allMonomials])) # return projectedMonomials # Create the monomials from the exponents (with coefficient 1) and add up the results newPoly = sum([Rxy({x:1}) for x in projectedExponents]) return Rxy(newPoly) ################################################### # Re-embedding Polyhedron in larger ambient space # ################################################### # The following function re-embeds a polyhedron in R^n to a polyhedron in R^(n+2) by adding two coordinates equal to 0 and computing the list of inequalities as tuples: def inequalitiesOfReembeddedPolyhedron(PolyP): ieqsCell = [] for k in range(len(PolyP.inequalities())): oldve = list(PolyP.inequalities()[k].vector()) newve = tuple(oldve + [0,0]) ieqsCell.append(newve) return ieqsCell def equationsOfReembeddedPolyhedron(PolyP): eqsCell = [] for k in range(len(PolyP.equations())): oldve = list(PolyP.equations()[k].vector()) newve = tuple(oldve + [0,0]) eqsCell.append(newve) return eqsCell ############################################### # Projecting away from the last 2 coordinates # ############################################### # The following script takes a polyhedron and computes the polyhedron resulting from the projection away from the last 2 coordinates def projectPolyhedra(PolyhedronP): n = PolyhedronP.ambient_dim() raysP = [x.vector() for x in PolyhedronP.rays()] linesP = [x.vector() for x in PolyhedronP.lines()] negLinesP = [-x for x in linesP] allRaysP = raysP + linesP + negLinesP # project away the last 2 coordinates avoiding having the origin as a possible projection: projAllRaysP = [x[0: -2] for x in allRaysP if list(x[0:-2]) != (n-2)*[0]] return Polyhedron(rays = projAllRaysP, base_ring = QQ) ############################################# # calculation for dictionary of polynomials # ############################################# def fromDictPolyToReductedPolysWithExtremalMonomials(dictPolys,Rring): ReducedPolysByCoefficients = {} for j in dictPolys.keys(): ReducedPolysByCoefficients[j] = fromPolynomialsToExtremalMonomials(dictPolys[j],Rring) return ReducedPolysByCoefficients def fromReductedPolysToGrobnerFans(dictPolys): GrobnerFansDict = {} for j in dictPolys.keys(): print str(j) GrobnerFansDict[j] = GrobnerFanComputations(dictPolys[j]) return GrobnerFansDict