// Setting up the Polynomial ring and ideal R := PolynomialRing(RationalField(), 10, "grevlex"); MyIdeal := ideal; // Variable names, unfortunately "Names(R)" does not work for old magma versions varNames := [ "c_0011_0", "c_0011_2", "c_0011_5", "c_0011_7", "c_0101_0", "c_0101_1", "c_0101_4", "c_0101_7", "c_0101_8", "c_1001_2" ]; // Term order in which the decomposition will be reported termOrder := "lex"; // Whether to compute witnesses for 1-dim ideals computeWitnesses := true; // Whether to compute the genus for 1-dim ideals computeGenuses := true; // Data necessary to recover all Ptolemy coordinates and the triangulation print "==TRIANGULATION=BEGINS==\n% Triangulation\nL11n443\ngeometric_solution 10.14941606\noriented_manifold\nCS_unknown\n\n4 0\n torus 0.000000000000 0.000000000000\n torus 0.000000000000 0.000000000000\n torus 0.000000000000 0.000000000000\n torus 0.000000000000 0.000000000000\n\n10\n 1 1 2 3 \n 0132 1230 0132 0132\n 3 3 0 1 \n 0 0 0 0 0 0 0 0 0 -1 0 1 -1 1 0 0\n 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n 0 -1 0 1 0 0 0 0 0 0 0 0 -1 1 0 0\n 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n 0.500000000000 0.866025403784\n\n 0 4 0 5 \n 0132 0132 3012 0132\n 1 3 1 0 \n 0 0 -1 1 0 0 0 0 1 0 0 -1 0 -1 1 0\n 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n 0 0 -1 1 0 0 1 -1 1 -1 0 0 0 0 0 0\n 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n 0.500000000000 0.866025403784\n\n 6 4 7 0 \n 0132 1230 0132 0132\n 3 3 1 1 \n 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n 0.500000000000 0.866025403784\n\n 4 4 0 5 \n 0132 1302 0132 2103\n 3 3 1 0 \n 0 0 0 0 0 0 0 0 0 0 0 0 -1 2 -1 0\n 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n 0 0 -1 1 0 0 0 0 1 0 0 -1 -1 1 0 0\n 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n 0.500000000000 0.866025403784\n\n 3 1 2 3 \n 0132 0132 3012 2031\n 1 3 0 1 \n 0 0 0 0 0 0 0 0 1 1 0 -2 0 0 0 0\n 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n 0 0 0 0 0 0 0 0 1 0 0 -1 -1 1 0 0\n 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n 0.500000000000 0.866025403784\n\n 8 9 1 3 \n 0132 0132 0132 2103\n 1 3 3 1 \n 0 1 -1 0 0 0 0 0 -1 1 0 0 0 -1 1 0\n 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n 0 0 -1 1 -1 0 1 0 2 -1 0 -1 0 0 0 0\n 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n 0.500000000000 0.866025403784\n\n 2 9 8 8 \n 0132 1302 0132 0321\n 2 3 1 1 \n 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0\n 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n 0 -1 0 1 0 0 2 -2 0 0 0 0 0 0 0 0\n 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n 0.500000000000 0.866025403784\n\n 9 9 8 2 \n 0321 3201 3120 0132\n 3 3 1 2 \n 0 0 0 0 -1 0 1 0 1 -1 0 0 1 -1 0 0\n 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n 0 0 0 0 1 0 -1 0 0 0 0 0 -1 2 -1 0\n 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n 0.500000000000 0.866025403784\n\n 5 6 7 6 \n 0132 0321 3120 0132\n 2 3 1 3 \n 0 0 0 0 0 0 -1 1 0 0 0 0 1 -1 0 0\n 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n 0 -1 1 0 1 0 1 -2 0 0 0 0 -2 2 0 0\n 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n 0.500000000000 0.866025403784\n\n 7 5 7 6 \n 0321 0132 2310 2031\n 1 2 1 3 \n 0 -1 1 0 -1 0 1 0 -1 1 0 0 1 -1 0 0\n 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n 0 0 0 0 1 0 -2 1 0 0 0 0 -1 1 0 0\n 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n 0.500000000000 0.866025403784\n\n\n==TRIANGULATION=ENDS==\nPY=EVAL=SECTION=BEGINS=HERE\n{'variable_dict' : \n (lambda d: {\n 'c_0011_3' : - d['c_0011_0'],\n 'c_0011_0' : d['c_0011_0'],\n 'c_0011_1' : - d['c_0011_0'],\n 'c_1001_1' : - d['c_0011_0'],\n 'c_0011_4' : d['c_0011_0'],\n 'c_1010_4' : - d['c_0011_0'],\n 'c_0101_0' : d['c_0101_0'],\n 'c_0110_1' : d['c_0101_0'],\n 'c_0110_2' : d['c_0101_0'],\n 'c_0101_5' : d['c_0101_0'],\n 'c_0101_6' : d['c_0101_0'],\n 'c_0110_8' : d['c_0101_0'],\n 'c_1010_0' : d['c_0101_1'],\n 'c_0110_0' : d['c_0101_1'],\n 'c_0101_1' : d['c_0101_1'],\n 'c_0101_3' : d['c_0101_1'],\n 'c_1001_3' : d['c_0101_1'],\n 'c_0110_4' : d['c_0101_1'],\n 'c_0110_3' : d['c_0101_4'],\n 'c_1001_0' : d['c_0101_4'],\n 'c_1100_1' : - d['c_0101_4'],\n 'c_1010_2' : d['c_0101_4'],\n 'c_1100_5' : - d['c_0101_4'],\n 'c_0101_4' : d['c_0101_4'],\n 'c_1100_0' : - d['c_0101_8'],\n 'c_1100_2' : - d['c_0101_8'],\n 'c_1100_3' : - d['c_0101_8'],\n 'c_1100_7' : - d['c_0101_8'],\n 'c_0110_5' : d['c_0101_8'],\n 'c_0101_8' : d['c_0101_8'],\n 'c_0011_2' : d['c_0011_2'],\n 'c_0011_6' : - d['c_0011_2'],\n 'c_1010_1' : - d['c_0011_2'],\n 'c_1001_4' : - d['c_0011_2'],\n 'c_1001_5' : - d['c_0011_2'],\n 'c_1010_9' : - d['c_0011_2'],\n 'c_0101_2' : d['c_0011_5'],\n 'c_0110_6' : d['c_0011_5'],\n 'c_0110_7' : d['c_0011_5'],\n 'c_0011_5' : d['c_0011_5'],\n 'c_0011_8' : - d['c_0011_5'],\n 'c_0011_9' : - d['c_0011_5'],\n 'c_1010_3' : d['c_1001_2'],\n 'c_1001_2' : d['c_1001_2'],\n 'c_1100_4' : - d['c_1001_2'],\n 'c_1010_7' : d['c_1001_2'],\n 'c_1010_5' : - d['c_1001_2'],\n 'c_1001_9' : - d['c_1001_2'],\n 'c_0011_7' : d['c_0011_7'],\n 'c_1010_6' : - d['c_0011_7'],\n 'c_1100_9' : d['c_0011_7'],\n 'c_1001_6' : - d['c_0011_7'],\n 'c_0110_9' : - d['c_0011_7'],\n 'c_1010_8' : - d['c_0011_7'],\n 'c_1001_7' : d['c_0101_7'],\n 'c_0101_7' : d['c_0101_7'],\n 'c_0101_9' : - d['c_0101_7'],\n 'c_1100_6' : - d['c_0101_7'],\n 'c_1100_8' : - d['c_0101_7'],\n 'c_1001_8' : - d['c_0101_7'],\n 's_2_7' : d['1'],\n 's_1_7' : d['1'],\n 's_0_7' : - d['1'],\n 's_3_6' : d['1'],\n 's_2_6' : d['1'],\n 's_1_6' : d['1'],\n 's_1_5' : - d['1'],\n 's_0_5' : d['1'],\n 's_3_3' : - d['1'],\n 's_1_3' : - d['1'],\n 's_0_3' : d['1'],\n 's_2_2' : - d['1'],\n 's_1_2' : d['1'],\n 's_0_2' : d['1'],\n 's_3_1' : d['1'],\n 's_1_1' : - d['1'],\n 's_3_0' : d['1'],\n 's_2_0' : - d['1'],\n 's_1_0' : d['1'],\n 's_0_0' : - d['1'],\n 's_0_1' : - d['1'],\n 's_2_1' : d['1'],\n 's_3_2' : - d['1'],\n 's_2_3' : d['1'],\n 's_1_4' : - d['1'],\n 's_2_5' : d['1'],\n 's_0_6' : d['1'],\n 's_2_4' : d['1'],\n 's_3_7' : - d['1'],\n 's_0_4' : d['1'],\n 's_3_4' : - d['1'],\n 's_3_5' : - d['1'],\n 's_0_8' : d['1'],\n 's_1_9' : - d['1'],\n 's_3_9' : d['1'],\n 's_3_8' : d['1'],\n 's_1_8' : d['1'],\n 's_0_9' : - d['1'],\n 's_2_9' : d['1'],\n 's_2_8' : d['1']})}\nPY=EVAL=SECTION=ENDS=HERE\n"; // Various helper functions function SaturateIdeal(baseRing, theIdeal) print "Status: Computing Groebner basis..."; time Groebner(theIdeal); for i := 1 to Rank(baseRing) do saturatedIdeal := -1; print "Status: Saturating ideal (", i, "/", Rank(baseRing), ")..."; time saturatedIdeal := Saturation(theIdeal, baseRing.i); print "Status: Recomputing Groebner basis..."; time Groebner(saturatedIdeal); theIdeal := saturatedIdeal; end for; print "Status: Dimension of ideal: ", Dimension(saturatedIdeal); return saturatedIdeal; end function; function ZeroDimensionalIdealChangeOrder(theIdeal, order) if Dimension(theIdeal) le 0 then print "Status: Changing to term order ", order, "..."; time result := ChangeOrder(theIdeal, order); print "Status: Recomputing Groebner basis..."; time Groebner(result); print "Status: Confirming is prime..."; time isPrime := IsPrime(result); if not isPrime then print "Not prime!"; exit; end if; return result; else return theIdeal; end if; end function; function ZeroDimensionalIdealsChangeOrder(theIdeals, order) return [ ZeroDimensionalIdealChangeOrder(theIdeal, order) : theIdeal in theIdeals ]; end function; function FreeVariablesOfIdeal(varNames, theIdeal) D, vars := Dimension(theIdeal); return [ "\"" cat varNames[var] cat "\"" : var in vars]; end function; function FreeVariablesOfIdeals(varNames, theIdeals) return [ FreeVariablesOfIdeal(varNames, theIdeal) : theIdeal in theIdeals]; end function; function IsSuitableWitness(theIdeals, theIndex, witnessIdeal) if not Dimension(witnessIdeal) eq 0 then return false; end if; for index := 1 to #theIdeals do if not index eq theIndex then if Dimension(theIdeals[index]) gt 0 then if not Dimension(theIdeals[index] + witnessIdeal) eq -1 then return false; end if; end if; end if; end for; return true; end function; function IncrementList(theList) for index := #theList to 1 by -1 do if index eq 1 then theList[1] := theList[1] + 1; return theList; else if theList[index] lt theList[1] then theList[index] := theList[index] + 1; return theList; else theList[index] := 1; end if; end if; end for; end function; function GenerateWitnessIdeal(baseRing, theIdeal, coordinates) D, vars := Dimension(theIdeal); witnessIdeal := theIdeal; for index := 1 to #coordinates do var := vars[index]; val := coordinates[index]; witnessIdeal := witnessIdeal + ideal; end for; return SaturateIdeal(baseRing, witnessIdeal); end function; function FindWitness(baseRing, theIdeals, index) theIdeal := theIdeals[index]; D, vars := Dimension(theIdeal); if D le 0 then return ideal; end if; coordinates := [1: var in vars]; while true do witnessIdeal := GenerateWitnessIdeal(baseRing, theIdeal, coordinates); print "Status: Testing witness ", coordinates, "..."; time isSuitable := IsSuitableWitness(theIdeals, index, witnessIdeal); if isSuitable then return RadicalDecomposition(witnessIdeal)[1]; end if; coordinates := IncrementList(coordinates); end while; end function; function FindWitnesses(baseRing, theIdeals) print "Status: Finding witnesses for non-zero dimensional ideals..."; return [ FindWitness(baseRing, theIdeals, index) : index in [1..#theIdeals]]; end function; // Computations begin here saturatedIdeal := SaturateIdeal(R, MyIdeal); // Compute radical decomposition print "Status: Computing RadicalDecomposition"; time P := RadicalDecomposition(saturatedIdeal); print "Status: Number of components: ", #P; print "DECOMPOSITION=TYPE: RadicalDecomposition"; // Initialize Porder to -1 so that we can check later whether an error or user // interrupt happened. Porder := -1; // change 0-dim components to desired term order POrder := ZeroDimensionalIdealsChangeOrder(P, termOrder); print "IDEAL=DECOMPOSITION" cat "=TIME: ", Cputime(); if Type(POrder) eq RngIntElt then // An error or user interrupt occured, bail print "IDEAL=DECOMPOSITION" cat "=FAILED"; exit; end if; print "IDEAL=DECOMPOSITION" cat "=BEGINS=HERE"; POrder; print "IDEAL=DECOMPOSITION" cat "=ENDS=HERE"; print "FREE=VARIABLES=IN=COMPONENTS" cat "=BEGINS=HERE"; FreeVariablesOfIdeals(varNames, POrder); print "FREE=VARIABLES=IN=COMPONENTS" cat "=ENDS=HERE"; if computeWitnesses then witnesses := FindWitnesses(R, POrder); witnessesOrder := -1; witnessesOrder := ZeroDimensionalIdealsChangeOrder(witnesses, termOrder); if Type(witnessesOrder) eq RngIntElt then print "WITNESSES=FOR=COMPONENTS" cat "=FAILED"; exit; end if; print "==WITNESSES=FOR=COMPONENTS" cat "=BEGINS=="; for witness in witnessesOrder do print "==WITNESSES=BEGINS=="; if not IsZero(witness) then print "==WITNESS=BEGINS=="; witness; print "==WITNESS=ENDS=="; end if; print "==WITNESSES=ENDS=="; end for; print "==WITNESSES=FOR=COMPONENTS" cat "=ENDS=="; end if; if computeGenuses then print "==GENUSES=FOR=COMPONENTS" cat "=BEGINS=="; for comp in POrder do print "==GENUS=FOR=COMPONENT=BEGINS=="; D, vars := Dimension(comp); if D eq 1 then print Genus(Curve(AffineSpace(R), comp)); end if; print "==GENUS=FOR=COMPONENT=ENDS=="; end for; print "==GENUSES=FOR=COMPONENTS" cat "=ENDS=="; end if;