//////////////////////////////////////////////////////////////////////////////// version="$Id: compcox.lib 8 2015-02-25 21:27:46Z ami $"; category="Algebraic Geometry"; info=" LIBRARY: compcox.lib PURPOSE: Modifying canonically embedded Mori Dream Spaces (CEMDS). AUTHORS: Ulrich Derenthal, Juergen Hausen, Armand Heim, Simon Keicher, Antonio Laface. OVERVIEW: This library provides a framework for modifying canonically embedded Mori Dream Spaces (CEMDS). NOTE: Sometimes variables need to be exported. They carry names of the structure g*Result, where * may be any combination of prefixes and types defined below. PROCEDURES: //CEMDS type and basic methods definition assignCEMDS(list #); Overrides the standard singular assignment for more comfortable CEMDS assignments. printCEMDS(CEMDS cemds); Prints the CEMDS's data to stdout. createCEMDS(intmat rvcvzP, fan scnSigma, ideal spG); Composes a CEMDS from a matrix of fan's rays, the fan itself and an ideal embedding the MDS in an ambient toric variety. encodeCEMDS(CEMDS cemds); Prints SINGULAR code defining the passed CEMDS on execution. fetchCEMDS(CEMDS cemdsRight); Associates a CEMDS associated with an arbitrary ring with the current basering, especially fetching the ideal \"spG\" embedding the MDS in its ambient toric variety. //Matrix and vector operations gale(matz); Computes a matrix whose rows generate the integer kernel of the passed matrix. galeExtension(matzQ, matzPOld); Computes a matrix whose rows generate the integer kernel of the first passed matrix, assuming this matrix is some gale dual of the second passed matrix extended by some columns. intmatAppendCol(rvcvzInput, intvec cvz); Appends a column to a bigintmat/intmat/intvec. intmatAppendCols(rvcvzA, rvcvzB); Appends the second matrix to the first matrix by generating new columns. intmatAppendRow(cvrvzInput, intvec rvz); Appends a row to a bigintmat/intmat/intvec. intmatAppendRows(cvrvzA, cvrvzB); Appends the second matrix to the first matrix by generating new rows. intmatContainsRow(cvrvz, intvec rvz); Checks whether a bigintmat/intmat contains a specific row. intmatDeleteRow(cvrvz, int irvzDel); Deletes a specified row from a bigintmat/intmat. intmatNonNegativeCols(rvcvz); Adds an individual constant vector to each of the bigintmat's/intmat's columns so that no entry is negative and at least one entry equals zero. intmatReplaceCol(rvcvz, int icvz, intvec cvzReplace); Replaces a bigintmat's/intmat's column by another one. intmatTakeCol(rvcvz, int icvzTake); Extracts a specified column from a bigintmat/intmat. intmatTakeCols(rvcvz, intvec rvfTakeCol); Generates a submatrix of a bigintmat/intmat by only taking the columns specified. intmatTakeRow(cvrvz, int irvzTake); Extracts a specified row from a bigintmat/intmat. intmatTakeRows(cvrvz, intvec cvfTakeRow); Generates a submatrix of a bigintmat/intmat by only taking the rows specified. intvecComponentAnd(intvec vfV, intvec vfW); Performs the logic AND operation (&&) on two vectors component by component. intvecComponentNot(intvec vfV); Performs the logic NOT operation (!) on a vector component by component. intvecComponentOr(intvec vfV, intvec vfW); Performs the logic OR operation (||) on two vectors component by component. intvecMin(intvec vz); Gets the minimal entry of an intvec. intvecSumFlags(intvec vfV); Determines the number of nonzero entries in an intvec. isRowLinearlyDependent(cvrvz, intvec rvz); Checks whether the collection consisting of the row vectors of a bigintmat/intmat and another intvec is linearly dependent. lusolveInvertible(matrix P, matrix L, matrix U, matrix cvnumB); Solves the system of linear equations Ax = cvnumB <=> matnumL*matnumU*x = matnumP*cvnumB for invertible A. lusolveUnderdetFullRank(matrix P, matrix L, matrix U, matrix cvnumB);Solves the underdetermined system of linear equations Ax = cvnumB <=> L*U*x = P*cvnumB. unitMatrix(int n); Builds an n x n unit matrix with integer entries. vectorComponentZero(vnum, int nnum); Marks all those entries of a vector positively which equal zero. vectorFromIntvec(intvec vz); Converts an intvec to a vector. vectorTake(vector vz, intvec vf); Generates a vector by taking specified entries of an original vector. //Polynomial handling containedVariables(poly p); Collects the indices of all variables that are contained within the passed polynomial with a positive exponent. evaluatePoly(poly p, vz); Calculates p(vz). exponentMatrixFromPoly(poly p); Lists the exponents of all polynomial's terms as rows of a matrix. getDegree(intmat rvcvzQ, poly hpQ); Calculates the degree of the passed polynomial homogeneous w.r.t. to the grading given by the passed matrix. intersectContainedVariables(list #); Collects the indices of all variables that are contained within all passed polynomials with a positive exponent. isOfTotalDegree(int nDeg, poly p); Checks whether a polynomial is homogeneous and of a specific degree w.r.t. standard grading. isTotalDegreeHomogeneous(poly p); Determines whether a polynomial is homogeneous w.r.t. standard grading. veronese(int nDeg, vector vnum); Applies the veronese embedding of a certain degree on a vector. //Geometric functions isPointInVariety(ideal spI, vz); Checks whether a point is contained in the variety V(K^n; spI) of the passed ideal spI. hyperplanes(list svnumPts); Computes all hyperplanes through the passed points. quadrics(list svnumPts); Computes all quadrics through the passed points. //Cone and fan operations containsRay(fan scn, intvec vzRay); Checks whether a ray is contained in a fan. contractRay(fan scn, intvec vzRay); Computes the fan obtained by contracting one of the rays of the original fan. generateOrthantFace(intvec vfIndices); Calculates the orthant face generated by the canonical basis vectors flagged true in the passed vector. irrelevantIdealFromFan(fan sof); Calculates the irrelevant ideal that is associated with a fan consisting of faces of the positive orthant. isRayContractible(rvcvzP, int icvzRay); Checks whether a ray from a list of rays is contractible, i.e. whether its corresponding ray in the gale dual version is extremal. mapFan(intmat rvcvzP, fan scn); Maps a fan under a matrix by mapping each cone under that matrix. movingConeFromWeights(intmat rvcvzQ); Calculates the moving cone of the cone defined by the row vectors of the passed intmat. orthantFanFromWeight(intmat rvcvzQ, intvec cvzWeight); Calculates the fan corresponding to a grading matrix and a weight vector. stellarSubdivision(fan scn, intvec vzNewRay); Calculates the stellar subdivision of a fan w.r.t. a new ray to be inserted. //Ideal and poly generation binomialIdeal(intmat rvcvzP); Calculates the lattice ideal generated by the rows of a matrix. minimalizeIdeal(ideal sp); Computes a minimal representation of an ideal. toricMorphismPullback(intmat matz, ideal sp); Pulls back an ideal under a morphism of standard tori determined by a martix. varproductOrbitClosureIdeal(intmat rvcvzP, intvec rvzPointInOrbit); Calculates the ideal of a point's orbit closure under a torus action implied by the rays of the considered toric variety's fan. varproductIdealSaturation(ideal spI); Calculates the saturation of an ideal w.r.t. the monomial var(1)*...*var(nvars). xIdealSaturation(ideal spI, list simonLexVars, intvec vfSatVars); Calculates the saturation of an ideal w.r.t. the monomial defined by vfSatVars. In contrast to working with the standard elimination order, some variables can be ordered lexicographically. veronesePoly(int nDeg, vector vnum); Expands a vector of coefficients for the monomials of a certain degree to the sum of those monomials multiplied by their corresponding coefficient. //Algorithms on CEMDSs blowupCEMDS(CEMDS cemds, list spC, intvec spD, list #); Blows up a CEMDS in a subvariety contained in the smooth locus. blowupCEMDSpoints(CEMDS cemds, list #); Blows up a CEMDS in some points contained therein. compressCEMDS(CEMDS cemds, list #); Embeds a CEMDS in a new ambient toric variety such that abundant relations in the ideal defining the CEMDS (i.e. of the type var(i) - p) are erased. Optionally checks whether the result really is a CEMDS and transfers points from the uncompressed CEMDS to the compressed one. contractCEMDS(CEMDS cemds, intvec vfKeepRay, list #); Blows down the ambient toric variety of a CEMDS; all rays not to be kept are deleted. linearBlowup(list svnumPts, list #); Blows up the projective space of dimension (n-1) in some of its points using linear relations for the stretching part only. linearQuadricBlowupP2(list svnumPts, list #); Blows up the projective space of dimension 2 in some of its points using relations of degree 1 and 2 for the stretching part only. modifyCEMDS(CEMDS cemds, list sRaysFan, list #); Blows up the ambient toric variety of a CEMDS by a refining of its fan. Optionally checks whether the result really is a CEMDS and transfers points from the unmodified CEMDS to the modified one. modifyCEMDS2(CEMDS cemds, intvec rvnVanishingDegrees, list #); Blows up the ambient toric variety of a CEMDS by adding some rays to its fan and performing stellar subdivision. Optionally checks whether the result really is a CEMDS and transfers points from the unmodified CEMDS to the modified one. stretchCEMDS(list cemds, list shpQPrimes, list #); Embeds a CEMDS in a new ambient toric variety such that each passed polynomial is set in relation to a new variable. Optionally transfers points from the unstretched CEMDS to the stretched one. verifyAlmostFree(intmat rvcvzP); Checks whether the grading defined by a gale dual matrix to the passed matrix is an almost free grading. verifyDimension(ideal spG); Checks whether dim(spG) - dim(spG + + ) >= 2 for all i != j. verifyVarPrimality(ideal spG); Checks whether each ring variable's class in the basering modulo spG defines a prime ideal. KEYWORDS: "; /* Hungarian notation: * * Prefixes: * cv: column vector * g: global variable, used for exporting * i: index in a list/vector/... * mat: matrix with row/column priority unspecified * n: count/number of entries in a list/vector/... * rv: row vector * s: set, e.g. a list (unsorted) * st: stack * v: unspecified row or column vector * * Types: * bin: binomial * cemds: list containing a ray matrix, a fan and embedding equations defining a CEMDS * cn: general cone * f: flag/boolean * hp*: homogeneous polynomial (w.r.t. to grading matrix *) * map: map * mon: monomial * n: unsigned integer (0, 1, ...) * num: number in the current basering * of: A cone being a face of the positive orthant * p: polynomial * str: string * z: signed integer (0, +-1, ...) */ /* Printlevel settings: * 0: Silent mode * 1: Show status information * 2: Debugging information */ //Included libraries LIB "matrix.lib"; LIB "multigrading.lib"; LIB "primdec.lib"; /*********************************************** * CEMDS type and basic methods definition * ***********************************************/ //CEMDS struct definition on loading the library proc mod_init() { newstruct("CEMDS", "def R, intmat rvcvzP, fan scnSigma, ideal spG"); system("install", "CEMDS", "=", assignCEMDS, 1); system("install", "CEMDS", "print", printCEMDS, 1); system("--ticks-per-sec", 1000); // set timer resolution to ms } proc assignCEMDS(list #) "USAGE: cemdsLeft = list(rvcvzP, scnSigma, spG) with rvcvzP: intmat, scnSigma: fan, spG: ideal. ASSUME: A basering is defined. PURPOSE: Overrides the standard singular assignment for more comfortable CEMDS assignments. RETURN: The CEMDS fetched or composed from the passed parameter(s). NOTE: cemdsLeft = cemdsRight; with cemdsRight: CEMDS may also work as soon as Singular allows overriding the standard assignment, too. EXAMPLE: None." { //Empty assignment -> empty CEMDS if (size(#) == 0) { CEMDS cemdsLeft; cemds.R = basering; return(cemdsLeft); } //One argument -> has to be a CEMDS already. Fetch it to the current basering. //NOTE: This case treats CEMDS assignments under the same preconditions that the = standard assignment operator assumes. This standard assignment operator under standard preconditions cannot be overridden at the moment, so the code below may only be interpreted in future releases of Singular that allow overriding the standard assignment operator under standard preconditions. if (size(#) == 1 && typeof(#) == "CEMDS") { return(fetchCEMDS(#)); } //Three arguments -> individual assignment of matrix, fan and generators if (size(#) == 3 && typeof(#[1]) == "intmat" && typeof(#[2]) == "fan" && typeof(#[3]) == "ideal") { return(createCEMDS(#[1], #[2], #[3])); } ERROR("A CEMDS can only be assigned another CEMDS or a list of its matrix, fan and embedding ideal."); } proc printCEMDS(CEMDS cemds) "USAGE: printCEMDS(cemds); OR print(cemds); OR cemds; with cemds: CEMDS. ASSUME: The CEMDS is associated with some ring. PURPOSE: Prints the CEMDS's data to stdout. RETURN: Nothing. EXAMPLE: None." { //The CEMDS to print must be associated with a ring. if (!defined(cemds.R)) { print("The CEMDS is not associated with any ring (Member \"R\" not defined)."); } if (typeof(cemds.R) != "ring") { print("The CEMDS is not associated with any ring (Member \"R\" not of type \"ring\")."); } //Switch to the CEMDS's ring def RCemds = cemds.R; setring RCemds; //Print the CEMDS //Its ring print(""); print("The CEMDS's ring:"); print(cemds.R); //Its rays print(""); print("The column matrix P of the CEMDS's fan's rays:"); print(cemds.rvcvzP); //Its fan print(""); print("The CEMDS's fan via its maximal cones, each one denoted by a column matrix of its rays:"); int iDim; if (size(fVector(cemds.scnSigma)) == 0) { int nDim = 0; print("Empty fan of ambient dimension " + string(ambientDimension(cemds.scnSigma))); } else { int nDim = dimension(cemds.scnSigma); } int icnMax; int ncnMax; cone cnMax; bigintmat rvcvzRays; //Find all cones and export them first //Iterate over the cone dimension for (iDim = 1; iDim <= nDim; iDim++) { //Iterate over the maximal cones of a dimension ncnMax = numberOfConesOfDimension(cemds.scnSigma, iDim, 0, 1); if (ncnMax > 0) { print("Dimension " + string(iDim) + ":"); } for (icnMax = 1; icnMax <= ncnMax; icnMax++) { cnMax = getCone(cemds.scnSigma, iDim, icnMax, 0, 1); rvcvzRays = transpose(rays(cnMax)); print(ordStrFromInt(icnMax) + " maximal cone:"); print(rvcvzRays); print(""); } } //Its ideal print(""); if (defined(cemds.spG)) { print("The equations' ideal G embedding the MDS into its ambient toric variety:"); print(cemds.spG); } else { print("The equations' ideal G embedding the MDS into its ambient toric variety is not defined (Member \"spG\" not defined)."); } } proc createCEMDS(intmat rvcvzP, fan scnSigma, ideal spG) "USAGE: createCEMDS(rvcvzP, scnSigma, spG); with rvcvzP: intmat, scnSigma: fan, spG: ideal. ASSUME: A basering with spG local to it is defined. The columns of rvcvzP create the rays of scnSigma. PURPOSE: Composes a CEMDS from a matrix of fan's rays, the fan itself and an ideal embedding the MDS in an ambient toric variety. RETURN: A CEMDS that is CEMDS defined by the passed parameters. EXAMPLE: example createCEMDS; shows an example." { //Make sure a basering is defined if (!defined(basering)) { ERROR("A basering needs to be defined for creating CEMDSs."); } //Compose and return the CEMDS CEMDS cemdsLeft; cemdsLeft.R = basering; cemdsLeft.rvcvzP = rvcvzP; cemdsLeft.scnSigma = scnSigma; cemdsLeft.spG = spG; return(cemdsLeft); } example { "EXAMPLE:"; echo=2; //A CEMDS representing the projective plane P2 //Define input parameters first ring R = 0,T(1..3),dp; intmat rvcvzP[2][3] = 1,0,-1, 0,1,-1; intmat cvrvz1[2][2] = 1,0, 0,1; intmat cvrvz2[2][2] = 0,1, -1,-1; intmat cvrvz3[2][2] = -1,-1, 1,0; cone cn1 = coneViaPoints(cvrvz1); cone cn2 = coneViaPoints(cvrvz2); cone cn3 = coneViaPoints(cvrvz3); fan scnSigma = fanViaCones(cn1, cn2, cn3); ideal spG = 0; //Then create a new CEMDS. CEMDS cemds = createCEMDS(rvcvzP, scnSigma, spG); print(cemds); //Should yield the CEMDS arising from the parameters passed associated with the ring defined above. } proc encodeCEMDS(CEMDS cemds) "USAGE: encodeCEMDS(cemds); with cemds: CEMDS. ASSUME: The CEMDS is associated with some ring. PURPOSE: Prints SINGULAR code defining the passed CEMDS on execution. RETURN: Nothing. EXAMPLE: example encodeCEMDS; shows an example." { //The CEMDS to print must be associated with a ring. if (!defined(cemds.R)) { print("The CEMDS is not associated with any ring (Member \"R\" not defined)."); } if (typeof(cemds.R) != "ring") { print("The CEMDS is not associated with any ring (Member \"R\" not of type \"ring\")."); } //Switch to the CEMDS's ring def RCemds = cemds.R; setring RCemds; print("//---------------START CEMDS ENCODING---------------"); //Export generation code for ring print("ring R = (" + charstr(cemds.R) + "),(" + varstr(cemds.R) + "),(" + ordstr(cemds.R) + ");"); //Export generation code for rays' matrix print("intmat rvcvzP = intmat(intvec("); cemds.rvcvzP; print("), " + string(nrows(cemds.rvcvzP)) + ", " + string(ncols(cemds.rvcvzP)) + ");"); //Export generation code for fan print("intmat cvrvz;"); //Make sure the passed CEMDS's fan is defined if (size(fVector(cemds.scnSigma)) == 0) { print("//WARNING: The passed CEMDS's fan is not defined."); print("fan scnSigma;"); } else { print("list scn = list();"); print("list scnList = list();"); //Helper variables int iDim; int nDim = dimension(cemds.scnSigma); int icnMax; int ncnMax; cone cnMax; bigintmat cvrvzRays; //Find all cones and export them first //Iterate over the cone dimension for (iDim = 1; iDim <= nDim; iDim++) { //Iterate over the maximal cones of a dimension ncnMax = numberOfConesOfDimension(cemds.scnSigma, iDim, 0, 1); for (icnMax = 1; icnMax <= ncnMax; icnMax++) { cnMax = getCone(cemds.scnSigma, iDim, icnMax, 0, 1); cvrvzRays = rays(cnMax); print("cvrvz = intmat(intvec("); print(cvrvzRays); print("), " + string(nrows(cvrvzRays)) + ", " + string(ncols(cvrvzRays)) + ");"); print("scn = coneViaPoints(cvrvz);"); print("scnList = scnList + scn;"); } } print("fan scnSigma = fanViaCones(scnList);"); } //Export generation code for ideal print("ideal spG = "); print(cemds.spG); print(";"); //Export generation code for creationg CEMDS print("CEMDS cemds = createCEMDS(rvcvzP, scnSigma, spG);"); print("print(\"\");"); print("print(\"The CEMDS was successfully imported to the variable \\\"cemds\\\"!\");"); print("//--------------- END CEMDS ENCODING ---------------"); return(); } example { "EXAMPLE:"; echo=2; //A CEMDS representing the projective plane P2 //Define input parameters first ring R = (0,a),T(1..3),dp; intmat rvcvzP[2][3] = 1,0,-1, 0,1,-1; intmat cvrvz1[2][2] = 1,0, 0,1; intmat cvrvz2[2][2] = 0,1, -1,-1; intmat cvrvz3[2][2] = -1,-1, 1,0; cone cn1 = coneViaPoints(cvrvz1); cone cn2 = coneViaPoints(cvrvz2); cone cn3 = coneViaPoints(cvrvz3); fan scnSigma = fanViaCones(cn1, cn2, cn3); ideal spG = 0; //Then create a new CEMDS. CEMDS cemds = createCEMDS(rvcvzP, scnSigma, spG); //Encode the CEMDS. This should yield SINGULAR code defining the CEMDS and thus similar to the code used above. encodeCEMDS(cemds); } proc fetchCEMDS(CEMDS cemdsRight) "USAGE: fetchCEMDS(cemdsRight); with cemdsRight: CEMDS. ASSUME: A basering is defined. The passed CEMDS is associated with a ring. PURPOSE: Associates a CEMDS associated with an arbitrary ring with the current basering, especially fetching the ideal \"spG\" embedding the MDS in its ambient toric variety. RETURN: A CEMDS that is the passed CEMDS associated with the current basering with the ideal \"spG\" embedding the MDS in its ambient toric variety fetched. EXAMPLE: example fetchCEMDS; shows an example." { //Make sure a basering is defined if (!defined(basering)) { ERROR("Cannot fetch the CEMDS as there is no basering defined."); } //Make sure the CEMDS to fetch is associated with a ring if (!defined(cemdsRight.R)) { ERROR("Cannot fetch the CEMDS because it is not associated with any ring."); } //Name the participating rings def RLeft = basering; def RRight = cemdsRight.R; //Compose a new CEMDS by copying the old CEMDS's matrix and fan and fetching its ideal to the current basering CEMDS cemdsLeft; cemdsLeft.R = basering; cemdsLeft.rvcvzP = cemdsRight.rvcvzP; cemdsLeft.scnSigma = cemdsRight.scnSigma; setring RRight; ideal spGRightClone = cemdsRight.spG; setring RLeft; ideal spGLeftClone = fetch(cemdsRight.R, spGRightClone); cemdsLeft.spG = spGLeftClone; return(cemdsLeft); } example { "EXAMPLE:"; echo=2; //Create a CEMDS associated with a certain basering ring R = 0,T(1..5),dp; intmat rvcvzP[4][5] = 1,0,-1,0,0, 0,1,-1,0,0, 0,0,-1,1,0, 0,0,-1,0,1; intmat cvrvz1[4][4] = -1,-1,-1,-1, 0,0,0,1, 0,0,1,0, 0,1,0,0; intmat cvrvz2[4][4] = -1,-1,-1,-1, 0,0,0,1, 0,0,1,0, 1,0,0,0; intmat cvrvz3[4][4] = -1,-1,-1,-1, 0,0,0,1, 0,1,0,0, 1,0,0,0; intmat cvrvz4[4][4] = -1,-1,-1,-1, 0,0,1,0, 0,1,0,0, 1,0,0,0; intmat cvrvz5[4][4] = 0,0,0,1, 0,0,1,0, 0,1,0,0, 1,0,0,0; cone cn1 = coneViaPoints(cvrvz1); cone cn2 = coneViaPoints(cvrvz2); cone cn3 = coneViaPoints(cvrvz3); cone cn4 = coneViaPoints(cvrvz4); cone cn5 = coneViaPoints(cvrvz5); fan scnSigma = fanViaCones(cn1, cn2, cn3, cn4, cn5); poly p1 = T(1)-T(2)+T(4); poly p2 = T(2)-T(3)+T(5); ideal spG = p1, p2; CEMDS cemds = createCEMDS(rvcvzP, scnSigma, spG); //Then switch to another ring and fetch the CEMDS. ring S = 0,T(1..5),lp; CEMDS cemdsFetched = fetchCEMDS(cemds); print(cemdsFetched); //Should yield the original CEMDS defined first but associated with the new ring. } /*********************************************** * Helpers * ***********************************************/ static proc ordStrFromInt(int n) "USAGE: ordStrFromInt(n); with n: int. PURPOSE: Converts an int to its ordinal string (1 -> \"1st\", 2 -> \"2nd\", ...). RETURN: The ordinal string of the passed int. EXAMPLE: None." { int nMod100 = n % 100; int nMod10 = nMod100 % 10; string strResult; if (nMod10 == 1 && nMod100 != 11) { strResult = string(n) + "st"; } else { if (nMod10 == 2 && nMod100 != 12) { strResult = string(n) + "nd"; } else { if (nMod10 == 3 && nMod100 != 13) { strResult = string(n) + "rd"; } else { strResult = string(n) + "th"; } } } return(strResult); } /*********************************************** * Matrix and vector operations * ***********************************************/ proc gale(matz) "USAGE: gale(matz); with matz: bigintmat/intmat. ASSUME: intmat limitations are not exceeded by bigintmats. PURPOSE: Computes a matrix whose rows generate the integer kernel of the passed matrix. RETURN: A bigintmat/intmat whose rows generate the integer kernel of the passed matrix. EXAMPLE: example gale; shows an example." { //Input type validation if (typeof(matz) != "intmat" && typeof(matz) != "bigintmat") { ERROR("The parameter passed must be of type \"bigintmat\" or \"intmat\"."); } return(transpose(kernelLattice(matz))); } example { "EXAMPLE:"; echo=2; //bigintmat example bigintmat rvcvzP1[2][4] = 1,0,-1,-1, 0,1,-1,0; bigintmat rvcvzQ1 = gale(rvcvzP1); print(rvcvzQ1); //Should yield // (1 1 1 0) // (1 0 0 1) //or something in the linear span of these two rows. //intmat example intmat rvcvzP2[2][4] = 1,0,-1,-1, 0,1,-1,0; intmat rvcvzQ2 = gale(rvcvzP2); print(rvcvzQ2); //Same as above, but as intmat. } proc galeExtension(matzQ, matzPOld) "USAGE: galeExtension(matzQ, matzPOld); with matzQ: bigintmat/intmat, matzPOld: bigintmat/intmat. ASSUME: intmat limitations are not exceeded by bigintmats. PURPOSE: Computes a matrix whose rows generate the integer kernel of the first passed matrix, assuming this matrix is some gale dual of the second passed matrix extended by some columns. RETURN: An intmat whose rows generate the integer kernel of the first passed matrix such that the second passed matrix is a submatrix of the resulting matrix located in the upper left corner. EXAMPLE: example galeExtension; shows an example." { //Input type check if (typeof(matzQ) != "intmat" && typeof(matzQ) != "bigintmat") { ERROR("The first parameter passed must be of type \"bigintmat\" or \"intmat\"."); } if (typeof(matzPOld) != "intmat" && typeof(matzPOld) != "bigintmat") { ERROR("The second parameter passed must be of type \"bigintmat\" or \"intmat\"."); } //First, compute a gale dual for matzQ. if (typeof(matzQ) == "bigintmat") { bigintmat matzQGale = gale(matzQ); } else { intmat matzQGale = gale(matzQ); } int icvzQGale; int ncvzQGale = ncols(matzQGale); int irvzQGale; int nrvzQGale = nrows(matzQGale); int ncvzPOld = ncols(matzPOld); int nrvzPOld = nrows(matzPOld); //Check whether the old gale dual matrix can be placed within the new one. if (ncvzPOld > ncvzQGale || nrvzPOld > nrvzQGale) { ERROR("The old gale matrix must not be larger than the current one."); } //Copy the old gale dual matrix to the resulting one //First, compute a gale dual for matzQ. intvec cvzZeros = 0:nrvzPOld; for (icvzQGale = ncvzPOld + 1; icvzQGale <= ncvzQGale; icvzQGale++) { matzPOld = intmatAppendCol(matzPOld, cvzZeros); } //Check whether the old matrix was indeed gale dual to the unextended matrix intmat matzZero[nrows(matzQ)][nrvzPOld]; if (matzQ * transpose(matzPOld) != matzZero) { ERROR("The matrix passed as old gale dual matrix was not gale dual to the original matrix."); } //Add new rows from the calculated gale matrix that do not linearly depend from the already added rows. intvec rvzNewRow; for (irvzQGale = 1; irvzQGale <= nrvzQGale; irvzQGale++) { rvzNewRow = intmatTakeRow(matzQGale, irvzQGale); if (!isRowLinearlyDependent(matzPOld, rvzNewRow)) { matzPOld = intmatAppendRow(matzPOld, rvzNewRow); } } return(matzPOld); } example { "EXAMPLE:"; echo=2; //bigintmat example bigintmat rvcvzP1[2][6] = 1,0,-1,1,0,-1, 0,1,-1,1,-1,0; bigintmat rvcvzQ1 = gale(rvcvzP1); intmat rvcvzNewColCoefs[6][3] = 0,1,1, 1,0,0, 0,0,0, 1,1,0, 0,0,1, 0,0,0; bigintmat rvcvzQ1Dash = intmatAppendCols(rvcvzQ1, rvcvzQ1 * rvcvzNewColCoefs); bigintmat rvcvzP1Dash = galeExtension(rvcvzQ1Dash, rvcvzP1); print(rvcvzP1Dash); //This matrix gale dual to rvcvzQ1Dash should have the original matrix as a submatrix in its top left corner. //intmat example intmat rvcvzP2[2][6] = 1,0,-1,1,0,-1, 0,1,-1,1,-1,0; intmat rvcvzQ2 = gale(rvcvzP2); intmat rvcvzQ2Dash = intmatAppendCols(rvcvzQ2, rvcvzQ2 * rvcvzNewColCoefs); intmat rvcvzP2Dash = galeExtension(rvcvzQ2Dash, rvcvzP2); print(rvcvzP2Dash); //This matrix gale dual to rvcvzQ2Dash should have the original matrix as a submatrix in its top left corner. intmat rvcvzP2DashDash = gale(rvcvzQ2Dash); print(rvcvzP2DashDash); //Cf. the "standard" matrix gale dual to rvcvzQ2Dash, which need not have the original matrix as a submatrix. } proc intmatAppendCol(rvcvzInput, intvec cvz) "USAGE: intmatAppendCol(rvcvzInput, cvz); with rvcvzInput: bigintmat/intmat/intvec, cvz: intvec. ASSUME: intmat limitations are not exceeded by bigintmats. PURPOSE: Appends a column to a bigintmat/intmat/intvec. RETURN: An bigintmat/intmat which is the passed bigintmat/intmat/intvec extended by the passed column. EXAMPLE: example intmatAppendCol; shows an example." { //Input type validation if (typeof(rvcvzInput) == "intvec") { intmat rvcvz = intmat(rvcvzInput, size(rvcvzInput), 1); } else { if (typeof(rvcvzInput) == "intmat") { intmat rvcvz = rvcvzInput; } else { if (typeof(rvcvzInput) == "bigintmat") { bigintmat rvcvz = rvcvzInput; } else { ERROR("The first parameter passed must be a of type \"bigintmat\", \"intmat\" or \"intvec\"."); } } } //Input dimension validation int icvz; int ncvz = ncols(rvcvz); int iz; int nz = nrows(rvcvz); if (size(cvz) != nz) { ERROR("The passed matrix and vector are of different size."); } //Set the output type if (typeof(rvcvz) == "bigintmat") { bigintmat rvcvzResult[nz][ncvz + 1]; } else { intmat rvcvzResult[nz][ncvz + 1]; } //Actually append the passed vector to the original matrix/vector for (iz = 1; iz <= nz; iz++) { for (icvz = 1; icvz <= ncvz; icvz++) { rvcvzResult[iz, icvz] = rvcvz[iz, icvz]; } rvcvzResult[iz, ncvz + 1] = cvz[iz]; } return(rvcvzResult); } example { "EXAMPLE:"; echo=2; //bigintmat example bigintmat rvcvz1[2][3] = 1,2,3, 4,5,6; intvec vz = -1,9; bigintmat rvcvzResult1 = intmatAppendCol(rvcvz1, vz); print(rvcvzResult1); //Should yield the intmat // (1 2 3 -1) // (4 5 6 9) //intmat example intmat rvcvz2[2][3] = 1,2,3, 4,5,6; intmat rvcvzResult2 = intmatAppendCol(rvcvz2, vz); print(rvcvzResult2); //Should yield the intmat // (1 2 3 -1) // (4 5 6 9) //intvec example intvec cvz1 = 1,2,3; intvec cvz2 = 7,8,9; intmat rvcvzResult3 = intmatAppendCol(cvz1, cvz2); print(rvcvzResult3); //Should yield // (1 7) // (2 8) // (3 9) } proc intmatAppendCols(rvcvzA, rvcvzB) "USAGE: intmatAppendCols(rvcvzA, rvcvzB); with rvcvzA: bigintmat/intmat, rvcvzB: bigintmat/intmat. ASSUME: intmat limitations are not exceeded by bigintmats. PURPOSE: Appends the second matrix to the first matrix by generating new columns. RETURN: An bigintmat/intmat which is the passed first matrix extended by appending the second passed matrix as new columns. EXAMPLE: example intmatAppendCols; shows an example." { //Input type validation if (typeof(rvcvzA) != "intmat" && typeof(rvcvzA) != "bigintmat") { ERROR("The first parameter passed must be of type \"bigintmat\" or \"intmat\"."); } if (typeof(rvcvzB) != "intmat" && typeof(rvcvzB) != "bigintmat") { ERROR("The second parameter passed must be of type \"bigintmat\" or \"intmat\"."); } //Input dimension validation int icvzA; int ncvzA = ncols(rvcvzA); int icvzB; int ncvzB = ncols(rvcvzB); int izA; int nzA = nrows(rvcvzA); if (nzA != nrows(rvcvzB)) { ERROR("The passed matrices have different numbers of rows."); } //Set the output type if (typeof(rvcvzA) == "bigintmat" || typeof(rvcvzB) == "bigintmat") { bigintmat rvcvzResult[nzA][ncvzA + ncvzB]; } else { intmat rvcvzResult[nzA][ncvzA + ncvzB]; } //Actually append the second to the first matrix for (izA = 1; izA <= nzA; izA++) { for (icvzA = 1; icvzA <= ncvzA; icvzA++) { rvcvzResult[izA, icvzA] = rvcvzA[izA, icvzA]; } for (icvzB = 1; icvzB <= ncvzB; icvzB++) { rvcvzResult[izA, ncvzA + icvzB] = rvcvzB[izA, icvzB]; } } return(rvcvzResult); } example { "EXAMPLE:"; echo=2; //bigintmat example bigintmat rvcvzA1[2][3] = 1,2,3, 4,5,6; intmat rvcvzB[2][2] = 7,8, 9,0; bigintmat rvcvzResult1 = intmatAppendCols(rvcvzA1, rvcvzB); print(rvcvzResult1); //Should yield // (1 2 3 7 8) // (4 5 6 9 0) //intmat example intmat rvcvzA2[2][3] = 1,2,3, 4,5,6; intmat rvcvzResult2 = intmatAppendCols(rvcvzA2, rvcvzB); print(rvcvzResult2); //Should yield // (1 2 3 7 8) // (4 5 6 9 0) } proc intmatAppendRow(cvrvzInput, intvec rvz) "USAGE: intmatAppendRow(cvrvzInput, rvz); with cvrvzInput: bigintmat/intmat/intvec, rvz: intvec. ASSUME: intmat limitations are not exceeded by bigintmats. PURPOSE: Appends a row to a bigintmat/intmat/intvec. RETURN: An bigintmat/intmat which is the passed bigintmat/intmat/intvec extended by the passed row. EXAMPLE: example intmatAppendRow; shows an example." { //Input type validation if (typeof(cvrvzInput) == "intvec") { intmat cvrvz = intmat(cvrvzInput, 1, size(cvrvzInput)); } else { if (typeof(cvrvzInput) == "intmat") { intmat cvrvz = cvrvzInput; } else { if (typeof(cvrvzInput) == "bigintmat") { bigintmat cvrvz = cvrvzInput; } else { ERROR("The first parameter passed must be a of type \"bigintmat\", \"intmat\" or \"intvec\"."); } } } //Input dimension validation int irvz; int nrvz = nrows(cvrvz); int iz; int nz = ncols(cvrvz); if (size(rvz) != nz) { ERROR("The passed matrix and vector are of different size."); } //Set the output type if (typeof(cvrvzInput) == "bigintmat") { bigintmat cvrvzResult[nrvz + 1][nz]; } else { intmat cvrvzResult[nrvz + 1][nz]; } //Actually append the passed vector to the original matrix/vector for (iz = 1; iz <= nz; iz++) { for (irvz = 1; irvz <= nrvz; irvz++) { cvrvzResult[irvz, iz] = cvrvz[irvz, iz]; } cvrvzResult[nrvz + 1, iz] = rvz[iz]; } return(cvrvzResult); } example { "EXAMPLE:"; echo=2; //bigintmat example bigintmat cvrvz1[2][3] = ( 1,2,3, 4,5,6); intvec rvz = 7,8,9; bigintmat cvrvzResult1 = intmatAppendRow(cvrvz1, rvz); print(cvrvzResult1); //Should yield // (1 2 3) // (4 5 6) // (7 8 9) //intmat example intmat cvrvz2[2][3] = ( 1,2,3, 4,5,6); intmat cvrvzResult2 = intmatAppendRow(cvrvz2, rvz); print(cvrvzResult2); //Should yield // (1 2 3) // (4 5 6) // (7 8 9) //intvec example intvec rvz1 = 1,2,3; intvec rvz2 = 7,8,9; intmat cvrvzResult3 = intmatAppendRow(rvz1, rvz2); print(cvrvzResult3); //Should yield // (1 2 3) // (7 8 9) } proc intmatAppendRows(cvrvzA, cvrvzB) "USAGE: intmatAppendRows(cvrvzA, cvrvzB); with cvrvzA: bigintmat/intmat, cvrvzB: bigintmat/intmat. ASSUME: intmat limitations are not exceeded by bigintmats. PURPOSE: Appends the second matrix to the first matrix by generating new rows. RETURN: An bigintmat/intmat which is the passed first matrix extended by appending the second passed matrix as new rows. EXAMPLE: example intmatAppendRows; shows an example." { //Input type validation if (typeof(cvrvzA) != "intmat" && typeof(cvrvzA) != "bigintmat") { ERROR("The first parameter passed must be of type \"bigintmat\" or \"intmat\"."); } if (typeof(cvrvzB) != "intmat" && typeof(cvrvzB) != "bigintmat") { ERROR("The second parameter passed must be of type \"bigintmat\" or \"intmat\"."); } //Input dimension validation int irvzA; int nrvzA = nrows(cvrvzA); int irvzB; int nrvzB = nrows(cvrvzB); int izA; int nzA = ncols(cvrvzA); if (nzA != ncols(cvrvzB)) { ERROR("The passed matrices have different numbers of columns."); } //Set the output type if (typeof(cvrvzA) == "bigintmat" || typeof(cvrvzB) == "bigintmat") { bigintmat cvrvzResult[nrvzA + nrvzB][nzA]; } else { intmat cvrvzResult[nrvzA + nrvzB][nzA]; } //Actually append the second to the first matrix for (izA = 1; izA <= nzA; izA++) { for (irvzA = 1; irvzA <= nrvzA; irvzA++) { cvrvzResult[irvzA, izA] = cvrvzA[irvzA, izA]; } for (irvzB = 1; irvzB <= nrvzB; irvzB++) { cvrvzResult[nrvzA + irvzB, izA] = cvrvzB[irvzB, izA]; } } return(cvrvzResult); } example { "EXAMPLE:"; echo=2; //bigintmat example bigintmat cvrvzA1[3][2] = 1,2, 3,4, 5,6; intmat cvrvzB[2][2] = 7,8, 9,0; bigintmat cvrvzResult1 = intmatAppendRows(cvrvzA1, cvrvzB); print(cvrvzResult1); //Should yield // (1 2) // (3 4) // (5 6) // (7 8) // (9 0) //intmat example intmat cvrvzA2[3][2] = 1,2, 3,4, 5,6; intmat cvrvzResult2 = intmatAppendRows(cvrvzA2, cvrvzB); print(cvrvzResult2); //Should yield // (1 2) // (3 4) // (5 6) // (7 8) // (9 0) } proc intmatContainsRow(cvrvz, intvec rvz) "USAGE: intmatContainsRow(cvrvz, rvz); with cvrvz: bigintmat/intmat, rvz: intvec. ASSUME: intmat limitations are not exceeded by bigintmats. PURPOSE: Checks whether a bigintmat/intmat contains a specific row. RETURN: The int 1 if the passed bigintmat/intmat contains the passed row, 0 else. EXAMPLE: example intmatContainsRow; shows an example." { //Input type validation if (typeof(cvrvz) != "intmat" && typeof(cvrvz) != "bigintmat") { ERROR("The first parameter passed must be of type \"bigintmat\" or \"intmat\"."); } //Input dimension validation int irvz; int nrvz = nrows(cvrvz); int iz; int nz = ncols(cvrvz); if (nz != size(rvz)) { ERROR("The passed matrix's column count and the passed vectors size differ."); } //Search for the passed vector row by row for (irvz = 1; irvz <= nrvz; irvz++) { //Check entries of a row for (iz = 1; iz <= nz; iz++) { //Check for a mismatch in the current row; if one was found, continue with the next row. if (cvrvz[irvz, iz] != rvz[iz]) { break; } //If no mismatch was found in a row, the passed vector is contained in the passed matrix if (iz == nz) { return(1); } } } return(0); } example { "EXAMPLE:"; echo=2; //bigintmat example bigintmat cvrvz1[2][3] = 1,2,3, 4,5,6; intvec rvz1 = 4,5,6; int fResult1 = intmatContainsRow(cvrvz1, rvz1); print(fResult1); //Should be true intvec rvz2 = 7,8,9; int fResult2 = intmatContainsRow(cvrvz1, rvz2); print(fResult2); //Should be false //intmat example intmat cvrvz2[2][3] = 1,2,3, 4,5,6; int fResult3 = intmatContainsRow(cvrvz2, rvz1); print(fResult3); //Should be true int fResult4 = intmatContainsRow(cvrvz2, rvz2); print(fResult4); //Should be false } proc intmatDeleteRow(cvrvz, int irvzDel) "USAGE: intmatDeleteRow(cvrvz, irvzDel); with cvrvz: bigintmat/intmat, irvzDel: int. ASSUME: intmat limitations are not exceeded by bigintmats. PURPOSE: Deletes a specified row from a bigintmat/intmat. RETURN: An bigintmat/intmat which is the passed bigintmat/intmat missing the passed index's row. EXAMPLE: example intmatDeleteRow; shows an example." { //Input type validation if (typeof(cvrvz) != "intmat" && typeof(cvrvz) != "bigintmat") { ERROR("The first parameter passed must be of type \"bigintmat\" or \"intmat\"."); } //Index in bounds validation int irvz; int nrvz = nrows(cvrvz); int iz; int nz = ncols(cvrvz); if (irvzDel <= 0 || irvzDel > nrvz) { ERROR("Index out of bounds."); } //Set the output type if (typeof(cvrvz) == "bigintmat") { bigintmat cvrvzResult[nrvz - 1][nz]; } else { intmat cvrvzResult[nrvz - 1][nz]; } //Simply copy all rows not to be deleted to the result matrix. for (irvz = 1; irvz <= nrvz; irvz++) { if (irvz != irvzDel) { for (iz = 1; iz <= nz; iz++) { if (irvz < irvzDel) { cvrvzResult[irvz, iz] = cvrvz[irvz, iz]; } else { cvrvzResult[irvz-1, iz] = cvrvz[irvz, iz]; } } } } return(cvrvzResult); } example { "EXAMPLE:"; echo=2; //bigintmat example bigintmat cvrvz1[3][3] = 1,2,3, 4,5,6, 7,8,9; bigintmat cvrvzResult1 = intmatDeleteRow(cvrvz1, 2); print(cvrvzResult1); //Should yield // (1 2 3) // (7 8 9) //intmat example intmat cvrvz2[3][3] = 1,2,3, 4,5,6, 7,8,9; intmat cvrvzResult2 = intmatDeleteRow(cvrvz2, 2); print(cvrvzResult2); //Should yield // (1 2 3) // (7 8 9) } proc intmatNonNegativeCols(rvcvz) "USAGE: intmatNonNegativeCols(rvcvz); with rvcvz: bigintmat/intmat. ASSUME: intmat limitations are not exceeded by bigintmats. PURPOSE: Adds an individual constant vector to each of the bigintmat's/intmat's columns so that no entry is negative and at least one entry equals zero. RETURN: The bigintmat/intmat with no negative entries and at least one entry being zero in each column resulting from adding constant vectors to the columns of the passed matrix. EXAMPLE: example intmatNonNegativeCols; shows an example." { //Input type validation if (typeof(rvcvz) != "bigintmat" && typeof(rvcvz) != "intmat") { ERROR("The parameter passed must be of type \"bigintmat\" or \"intmat\"."); } int icvz; int ncvz = ncols(rvcvz); int nz = nrows(rvcvz); //Iterate through the matrix's columns; add a constant vector to them so the column's minimal entry plus the constant vector will equal 0. intvec cvzCurrent; intvec vzColMin; for (icvz = 1; icvz <= ncvz; icvz++) { cvzCurrent = intmatTakeCol(rvcvz, icvz); vzColMin = intvecMin(cvzCurrent):nz; rvcvz = intmatReplaceCol(rvcvz, icvz, cvzCurrent - vzColMin); } return(rvcvz); } example { "EXAMPLE:"; echo=2; //bigintmat example bigintmat rvcvz1[2][8] = 2,1,0,-1,1,0,-1,-2, 1,0,-1,-2,2,1,0,-1; bigintmat rvcvzResult1 = intmatNonNegativeCols(rvcvz1); print(rvcvzResult1); //Should yield the bigintmat // (1 1 1 1 0 0 0 0) // (0 0 0 0 1 1 1 1) //intmat example intmat rvcvz2[2][8] = 2,1,0,-1,1,0,-1,-2, 1,0,-1,-2,2,1,0,-1; intmat rvcvzResult2 = intmatNonNegativeCols(rvcvz2); print(rvcvzResult2); //Should yield the intmat // (1 1 1 1 0 0 0 0) // (0 0 0 0 1 1 1 1) } proc intmatReplaceCol(rvcvz, int icvz, intvec cvzReplace) "USAGE: intmatReplaceCol(rvcvz, icvz, cvzReplace); with rvcvz: bigintmat/intmat, icvz: int, cvzReplace: intvec. ASSUME: intmat limitations are not exceeded by bigintmats. PURPOSE: Replaces a bigintmat's/intmat's column by another one. RETURN: The original bigintmat/intmat with the column corresponding to the passed index replaced by the passed intvec. EXAMPLE: example intmatReplaceCol; shows an example." { //Input type validation if (typeof(rvcvz) != "intmat" && typeof(rvcvz) != "bigintmat") { ERROR("The first parameter passed must be of type \"bigintmat\" or \"intmat\"."); } //Index in bounds validation if (icvz <= 0 || icvz > ncols(rvcvz)) { ERROR("The column index passed was out of bounds."); } //Passed vector's size validation int iz; int nz = nrows(rvcvz); if (size(cvzReplace) != nz) { ERROR("The passed intvec to replace the old column must have exactly as many entries as the passed matrix has rows."); } //Actually replace the matrix's column with index icvz by the other column verctor passed. for (iz = 1; iz <= nz; iz++) { rvcvz[iz, icvz] = cvzReplace[iz]; } return(rvcvz); } example { "EXAMPLE:"; echo=2; //bigintmat example bigintmat rvcvz1[2][3] = 1,2,0, 4,5,5; intvec cvzReplace = 3,6; bigintmat rvcvzResult1 = intmatReplaceCol(rvcvz1, 3, cvzReplace); print(rvcvzResult1); //Should yield the bigintmat // (1 2 3) // (4 5 6) //intmat example intmat rvcvz2[2][3] = 1,2,0, 4,5,5; intmat rvcvzResult2 = intmatReplaceCol(rvcvz2, 3, cvzReplace); print(rvcvzResult2); //Should yield the intmat // (1 2 3) // (4 5 6) } proc intmatTakeCol(rvcvz, int icvzTake) "USAGE: intmatTakeCol(rvcvz, icvzTake); with rvcvz: bigintmat/intmat, icvzTake: int. ASSUME: intmat limitations are not exceeded by bigintmats. PURPOSE: Extracts a specified column from a bigintmat/intmat. RETURN: An intvec that is the icvzTake-th column of the passed bigintmat/intmat. EXAMPLE: example intmatTakeCol; shows an example." { //Input type validation if (typeof(rvcvz) != "intmat" && typeof(rvcvz) != "bigintmat") { ERROR("The first parameter passed must be of type \"bigintmat\" or \"intmat\"."); } //Index in bounds validation int ncvz = ncols(rvcvz); if (icvzTake > ncvz || icvzTake < 1) { ERROR("Index out of bounds."); } //Take the column specified via a matrix multiplication. intvec cvTake = 0:ncvz; cvTake[icvzTake] = 1; return(rvcvz * cvTake); } example { "EXAMPLE:"; echo=2; //bigintmat example bigintmat rvcvz1[2][3] = 1,2,3, 4,5,6; intvec cvzResult1 = intmatTakeCol(rvcvz1, 2); print(cvzResult1); //Should yield the vector (2,5) //intmat example intmat rvcvz2[2][3] = 1,2,3, 4,5,6; intvec cvzResult2 = intmatTakeCol(rvcvz2, 2); print(cvzResult2); //Should yield the vector (2,5) } proc intmatTakeCols(rvcvz, intvec rvfTakeCol) "USAGE: intmatTakeCols(rvcvz, rvfTakeCol); with rvcvz: bigintmat/intmat, rvfTakeCol: intvec of flags. ASSUME: intmat limitations are not exceeded by bigintmats. PURPOSE: Generates a submatrix of a bigintmat/intmat by only taking the columns specified. RETURN: An bigintmat/intmat that is the submatrix of the passed bigintmat/intmat consisting of the columns flagged true in the passed vector. EXAMPLE: example intmatTakeCols; shows an example." { //Input type validation if (typeof(rvcvz) != "intmat" && typeof(rvcvz) != "bigintmat") { ERROR("The first parameter passed must be of type \"bigintmat\" or \"intmat\"."); } //Vector size validation int ifTakeCol; int nfTakeCol = size(rvfTakeCol); if (nfTakeCol != ncols(rvcvz)) { ERROR("The passed vector determining which columns to take must have as many entries as the passed matrix has columns."); } //Check whether there is at least one column to take specified. if (rvfTakeCol == 0:nfTakeCol) { ERROR("No column to take was specified."); } //Set the output type if (typeof(rvcvz) == "bigintmat") { bigintmat rvcvzTake[nfTakeCol][intvecSumFlags(rvfTakeCol)]; } else { intmat rvcvzTake[nfTakeCol][intvecSumFlags(rvfTakeCol)]; } //Generate a column selection matrix to multiply the input matrix with int icvzTake = 1; for (ifTakeCol = 1; ifTakeCol <= nfTakeCol; ifTakeCol++) { if (rvfTakeCol[ifTakeCol]) { rvcvzTake[ifTakeCol, icvzTake] = 1; icvzTake++; } } //Take the columns requested by a matrix multiplication. return(rvcvz * rvcvzTake); } example { "EXAMPLE:"; echo=2; //bigintmat example bigintmat rvcvz1[2][4] = 1,2,3,4, 5,6,7,8; intvec rvfTakeCol = 1,0,1,1; bigintmat rvcvzResult1 = intmatTakeCols(rvcvz1, rvfTakeCol); print(rvcvzResult1); //Should yield // (1 3 4) // (5 7 8) //intmat example intmat rvcvz2[2][4] = 1,2,3,4, 5,6,7,8; intmat rvcvzResult2 = intmatTakeCols(rvcvz2, rvfTakeCol); print(rvcvzResult2); //Should yield // (1 3 4) // (5 7 8) } proc intmatTakeRow(cvrvz, int irvzTake) "USAGE: intmatTakeRow(cvrvz, irvzTake); with cvrvz: bigintmat/intmat, irvzTake: int. ASSUME: intmat limitations are not exceeded by bigintmats. PURPOSE: Extracts a specified row from a bigintmat/intmat. RETURN: An intvec that is the irvzTake-th row of the passed bigintmat/intmat. EXAMPLE: example intmatTakeRow; shows an example." { //Input type validation if (typeof(cvrvz) != "intmat" && typeof(cvrvz) != "bigintmat") { ERROR("The first parameter passed must be of type \"bigintmat\" or \"intmat\"."); } //Index in bounds validation int nrvz = nrows(cvrvz); if (irvzTake > nrvz) { ERROR("Index out of bounds."); } //Take the row specified via a matrix multiplication. intvec rvzTake = 0:nrvz; rvzTake[irvzTake] = 1; return(intvec(transpose(rvzTake) * cvrvz)); } example { "EXAMPLE:"; echo=2; //bigintmat example bigintmat cvrvz1[3][3] = 1,2,3, 4,5,6, 7,8,9; intvec rvzResult1 = intmatTakeRow(cvrvz1, 2); print(rvzResult1); //Should yield the vector (4,5,6) //intmat example intmat cvrvz2[3][3] = 1,2,3, 4,5,6, 7,8,9; intvec rvzResult2 = intmatTakeRow(cvrvz2, 2); print(rvzResult2); //Should yield the vector (4,5,6) } proc intmatTakeRows(cvrvz, intvec cvfTakeRow) "USAGE: intmatTakeRows(cvrvz, cvfTakeRow); with cvrvz: bigintmat/intmat, cvfTakeRow: intvec of flags. ASSUME: intmat limitations are not exceeded by bigintmats. PURPOSE: Generates a submatrix of a bigintmat/intmat by only taking the rows specified. RETURN: An bigintmat/intmat which is the submatrix of the passed bigintmat/intmat consisting of the rows flagged true in the passed vector. EXAMPLE: example intmatTakeRows; shows an example." { //Input type validation if (typeof(cvrvz) != "intmat" && typeof(cvrvz) != "bigintmat") { ERROR("The first parameter passed must be of type \"bigintmat\" or \"intmat\"."); } //Vector size validation int ifTakeRow; int nfTakeRow = size(cvfTakeRow); if (nfTakeRow != nrows(cvrvz)) { ERROR("The passed vector determining which rows to take must have as many entries as the passed matrix has columns."); } //Check whether there is at least one row to take specified. if (cvfTakeRow == 0:nfTakeRow) { ERROR("No row to take was specified."); } //Generate a row selection matrix to multiply the input matrix with int irvzTake = 1; intmat cvrvzTake[intvecSumFlags(cvfTakeRow)][nfTakeRow]; for (ifTakeRow = 1; ifTakeRow <= nfTakeRow; ifTakeRow++) { if (cvfTakeRow[ifTakeRow]) { cvrvzTake[irvzTake, ifTakeRow] = 1; irvzTake++; } } //Take the rows requested by a matrix multiplication. return(cvrvzTake * cvrvz); } example { "EXAMPLE:"; echo=2; //bigintmat example bigintmat cvrvz1[4][2] = 1,5, 2,6, 3,7, 4,8; intvec cvfTakeRow = 1,0,1,1; bigintmat cvrvzResult1 = intmatTakeRows(cvrvz1, cvfTakeRow); print(cvrvzResult1); //Should yield // (1 5) // (3 7) // (4 8) //intmat example intmat cvrvz2[4][2] = 1,5, 2,6, 3,7, 4,8; intmat cvrvzResult2 = intmatTakeRows(cvrvz2, cvfTakeRow); print(cvrvzResult2); //Should yield // (1 5) // (3 7) // (4 8) } proc intvecComponentAnd(intvec vfV, intvec vfW) "USAGE: intvecComponentAnd(vfV, vfW); with vfV: intvec of flags, vfW: intvec of flags. PURPOSE: Performs the logic AND operation (&&) on two vectors component by component. RETURN: An intvec whose components are: result[i] = vfV[i] && vfW[i]. EXAMPLE: example intvecComponentAnd; shows an example." { //Check whether the sizes of the two passed vectors match. int ifV; int nfV = size(vfV); if (nfV != size(vfW)) { ERROR("The two passed vectors are not of equal size."); } //Perform the logic AND operation component by component. intvec vfResult = 0:nfV; for (ifV = 1; ifV <= nfV; ifV++) { vfResult[ifV] = vfV[ifV] && vfW[ifV]; } return(vfResult); } example { "EXAMPLE:"; echo=2; intvec rvfV = 1,1,0,0; intvec rvfW = 1,0,1,0; intvec rvfResult = intvecComponentAnd(rvfV, rvfW); print(rvfResult); //Should yield the vector (1,0,0,0). } proc intvecComponentNot(intvec vfV) "USAGE: intvecComponentNot(vfV); with vfV: intvec of flags. PURPOSE: Performs the logic NOT operation (!) on a vector component by component. RETURN: An intvec whose components are: result[i] = !vfV[i]. EXAMPLE: example intvecComponentNot; shows an example." { int ifV; int nfV = size(vfV); //Perform the logic NOT operation component by component. intvec vfResult = 0:nfV; for (ifV = 1; ifV <= nfV; ifV++) { vfResult[ifV] = !vfV[ifV]; } return(vfResult); } example { "EXAMPLE:"; echo=2; intvec vf = 1,0; intvec vfResult = intvecComponentNot(vf); print(vfResult); //Should yield the vector (0,1). } proc intvecComponentOr(intvec vfV, intvec vfW) "USAGE: intvecComponentOr(vfV, vfW); with vfV: intvec of flags, vfW: intvec of flags. PURPOSE: Performs the logic OR operation (||) on two vectors component by component. RETURN: An intvec whose components are: result[i] = vfV[i] || vfW[i]. EXAMPLE: example intvecComponentOr; shows an example." { //Check whether the sizes of the two passed vectors match. int ifV; int nfV = size(vfV); if (nfV != size(vfW)) { ERROR("The two passed vectors are not of equal size."); } //Perform the logic OR operation component by component. intvec vfResult = 0:nfV; for (ifV = 1; ifV <= nfV; ifV++) { vfResult[ifV] = vfV[ifV] || vfW[ifV]; } return(vfResult); } example { "EXAMPLE:"; echo=2; intvec vfV = 1,1,0,0; intvec vfW = 1,0,1,0; intvec vfResult = intvecComponentOr(vfV, vfW); print(vfResult); //Should yield the vector (1,1,1,0). } proc intvecMin(intvec vz) "USAGE: intvecMin(vz); with vz: intvec. PURPOSE: Gets the minimal entry of an intvec. RETURN: An int that is the minimal entry of the passed intvec. EXAMPLE: example intvecMin; shows an example." { int iz; int nz = size(vz); //Iterate through the invec to find its minimal entry. int zMin = vz[1]; for (iz = 2; iz <= nz; iz++) { if (vz[iz] < zMin) { zMin = vz[iz]; } } return(zMin); } example { "EXAMPLE:"; echo=2; intvec vz = 1,8,-5,0; int zMin = intvecMin(vz); zMin; //Should be -5. } proc intvecSumFlags(intvec vfV) "USAGE: intvecSumFlags(vfV); with vfV: intvec of flags. PURPOSE: Determines the number of nonzero entries in an intvec. RETURN: An int that is the count of all nonzero entries of the passed vector. EXAMPLE: example intvecSumFlags; shows an example." { int ifV; int nfV = size(vfV); int nResult = 0; for (ifV = 1; ifV <= nfV; ifV++) { if (vfV[ifV]) { nResult++; } } return(nResult); } example { "EXAMPLE:"; echo=2; intvec vf = 0,1,2,3; int nFlags = intvecSumFlags(vf); print(nFlags); //Should return 0 + 1 + 1 + 1 = 3. } proc isRowLinearlyDependent(cvrvz, intvec rvz) "USAGE: isRowLinearlyDependent(cvrvz, intvec rvz); with cvrvz: bigintmat/intmat, rvz: intvec. ASSUME: intmat limitations are not exceeded by bigintmats. PURPOSE: Checks whether the collection consisting of the row vectors of a bigintmat/intmat and another intvec is linearly dependent. RETURN: An int pointing out whether the collection consisting of the row vectors of the passed bigintmat/intmat and the passed intvec is linearly dependent (1) or not (0). EXAMPLE: example isRowLinearlyDependent; shows an example." { //Input type validation if (typeof(cvrvz) != "intmat" && typeof(cvrvz) != "bigintmat") { ERROR("The first parameter passed must be of type \"bigintmat\" or \"intmat\"."); } cvrvz = intmatAppendRow(cvrvz, rvz); intvec vzKernelEntries = intvec(kernelLattice(transpose(cvrvz))); return(vzKernelEntries != 0:size(vzKernelEntries)); } example { "EXAMPLE:"; echo=2; //bigintmat example bigintmat cvrvz1[2][3] = 1,0,0, 0,1,0; intvec rvz1 = 0,0,1; int f1 = isRowLinearlyDependent(cvrvz1, rvz1); print(f1); //Should be false intvec rvz2 = 1,1,0; int f2 = isRowLinearlyDependent(cvrvz1, rvz2); print(f2); //Should be true //intmat example intmat cvrvz2[2][3] = 1,0,0, 0,1,0; int f3 = isRowLinearlyDependent(cvrvz2, rvz1); print(f3); //Should be false int f4 = isRowLinearlyDependent(cvrvz2, rvz2); print(f4); //Should be true } proc lusolveInvertible(matrix P, matrix L, matrix U, matrix cvnumB) "USAGE: lusolveInvertible(P, L, U, cvnumB); with P: matrix of numbers, L: matrix of numbers, U: matrix of numbers, cvnumB: matrix of numbers. ASSUME: P, L, U arise from the LU-decomposition of an invertible matrix A (with P*A = L*U). PURPOSE: Solves the system of linear equations Ax = cvnumB <=> L*U*x = P*cvnumB for invertible A. RETURN: A matrix of numbers containing exactly one column with the solution to the system of linear equations Ax = cvnumB. EXAMPLE: example lusolveInvertible; shows an example." { //Check for correct input dimensions: P, L, U are square with identical count of rows and columns, B has as many entries as another matrix has rows/columns int ncvnumP = ncols(P); int nrvnumP = nrows(P); int ncvnumL = ncols(L); int nrvnumL = nrows(L); int ncvnumU = ncols(U); int nrvnumU = nrows(U); int inumB; int inumB2; int nnumB = nrows(cvnumB); if (ncvnumP != nrvnumP || ncvnumP != ncvnumL || ncvnumP != nrvnumL || ncvnumP != ncvnumU || ncvnumP != nrvnumU || ncvnumP != nnumB) { ERROR("lusolveInvertible: Input matrix dimension mismatch.") } //Solve Ly = Pb matrix cvnumY[nnumB][1]; cvnumB = P * cvnumB; for (inumB = 1; inumB <= nnumB; inumB++) { cvnumY[inumB, 1] = cvnumB[inumB, 1] / L[inumB, inumB]; for (inumB2 = 1; inumB2 < inumB; inumB2++) { cvnumY[inumB, 1] = cvnumY[inumB, 1] - L[inumB, inumB2] * cvnumY[inumB2, 1]; } } //Solve Ux = y matrix cvnumX[nnumB][1]; for (inumB = nnumB; inumB >= 1; inumB--) { cvnumX[inumB, 1] = cvnumY[inumB, 1]; for (inumB2 = nnumB; inumB2 > inumB; inumB2--) { cvnumX[inumB, 1] = cvnumX[inumB, 1] - U[inumB, inumB2] * cvnumX[inumB2, 1]; } cvnumX[inumB, 1] = cvnumX[inumB, 1] / U[inumB, inumB]; } return(cvnumX); } example { "EXAMPLE:"; echo=2; ring R = 0,x,dp; matrix matnumA[2][2] = 4,3,6,3; matrix cvnumB[2][1] = 7,9; list smatnumPLU = ludecomp(matnumA); matrix cvnumX = lusolveInvertible(smatnumPLU[1], smatnumPLU[2], smatnumPLU[3], cvnumB); print(cvnumX); //Should yield the 2x1-matrix [1,1]. } proc lusolveUnderdetFullRank(matrix P, matrix L, matrix U, matrix cvnumB) "USAGE: lusolveUnderdetFullRank(P, L, U, cvnumB); with P: matrix of numbers, L: matrix of numbers, U: matrix of numbers, cvnumB: matrix of numbers. ASSUME: P, L, U arise from the LU-decomposition of a matrix A of full rank (with P*A = L*U). Ax = cvnumB is an underdetermined system of linear equations. PURPOSE: Solves the underdetermined system of linear equations Ax = cvnumB <=> L*U*x = P*cvnumB. RETURN: A list containing (i) a matrix of numbers containing exactly one column with a specific solution to the system of linear equations Ax = cvnumB. (ii) a matrix of numbers whose columns span the affine linear space of all solutions to the system of linear equations Ax = cvnumB. EXAMPLE: example lusolveUnderdetFullRank; shows an example." { //Check for correct input dimensions: P, L, U are square with identical count of rows and columns, B has as many entries as another matrix has rows/columns int ncvnumP = ncols(P); int nrvnumP = nrows(P); int ncvnumL = ncols(L); int nrvnumL = nrows(L); int icvnumU; int ncvnumU = ncols(U); int irvnumU; int nrvnumU = nrows(U); int nnumB = nrows(cvnumB); if (ncvnumP != nrvnumP || ncvnumP != ncvnumL || ncvnumP != nrvnumL || ncvnumP >= ncvnumU || ncvnumP != nrvnumU || ncvnumP != nnumB) { ERROR("lusolveUnderdetFullRank: Input matrix dimension mismatch.") } //Find appropriate square upper triangular submatrix U' of U for solving LU'x = Pb first. irvnumU = 1; for (icvnumU = 1; icvnumU <= ncvnumU; icvnumU++) { if (irvnumU <= nrvnumU) { if (U[irvnumU, icvnumU] != 0) { if (!defined(rvicvnumUDash)) { intvec rvicvnumUDash = icvnumU; } else { rvicvnumUDash = rvicvnumUDash, icvnumU; } irvnumU++; } else { if (!defined(rvicvnumULinear)) { intvec rvicvnumULinear = icvnumU; } else { rvicvnumULinear = rvicvnumULinear, icvnumU; } } } else { if (!defined(rvicvnumULinear)) { intvec rvicvnumULinear = icvnumU; } else { rvicvnumULinear = rvicvnumULinear, icvnumU; } } } //Find first solution matrix matnumUDash = submat(U, 1..nrvnumU, rvicvnumUDash); matrix cvnumX[ncvnumU][1]; cvnumX[rvicvnumUDash, 1] = lusolveInvertible(P, L, matnumUDash, cvnumB); //Derive the remaining affine solution space. int iicvnumULinear; int nicvnumULinear = size(rvicvnumULinear); matrix cvnumBLinear; matrix cvnumXLinear[ncvnumU][1]; matrix matnumLinear[ncvnumU][nicvnumULinear]; for (iicvnumULinear = 1; iicvnumULinear <= nicvnumULinear; iicvnumULinear++) { cvnumBLinear = submat(L * U, 1..nrvnumU, rvicvnumULinear[iicvnumULinear]); cvnumXLinear[rvicvnumUDash, 1] = lusolveInvertible(unitmat(nrvnumP), L, matnumUDash, cvnumBLinear); cvnumXLinear[rvicvnumULinear[iicvnumULinear], 1] = -1; matnumLinear[1..ncvnumU, iicvnumULinear] = cvnumXLinear; } list smatnumResult = cvnumX, matnumLinear; return(smatnumResult); } example { "EXAMPLE:"; echo=2; ring R = 0,x,dp; matrix matnumA[2][3] = 1,1,1,0,-1,-2; matrix cvnumB[2][1] = 1,1; list smatnumPLU = ludecomp(matnumA); list smatnumResult = lusolveUnderdetFullRank(smatnumPLU[1], smatnumPLU[2], smatnumPLU[3], cvnumB); print(smatnumResult); //Should yield a list of two entries: // - A solution for the system of linear equations, e.g. the 3x1-matrix [2,-1,0]; // - Generators of the affine space of solutions, e.g. the 3x1-matrix [-1,2,-1]. kill matnumA; kill cvnumB; kill smatnumPLU; kill smatnumResult; matrix matnumA[2][3] = 0,1,0,1,1,1; matrix cvnumB[2][1] = 0,0; list smatnumPLU = ludecomp(matnumA); list smatnumResult = lusolveUnderdetFullRank(smatnumPLU[1], smatnumPLU[2], smatnumPLU[3], cvnumB); print(smatnumResult); //Should yield a list of two entries: // - A solution for the system of linear equations, e.g. the 3x1-matrix [0,0,0]; // - Generators of the affine space of solutions, e.g. the 3x1-matrix [-1,0,1]. } proc unitMatrix(int n) "USAGE: unitMatrix(int n); with n: int. PURPOSE: Builds an n x n unit matrix with integer entries. RETURN: The n x n unit matrix with integer entries. EXAMPLE: example unitMatrix; shows an example." { int i; intmat result[n][n]; for (i = 1; i <= n; i++) { result[i,i] = 1; } return(result); } example { "EXAMPLE:"; echo=2; intmat result = unitMatrix(3); print(result); //Should return // (1 0 0) // (0 1 0) // (0 0 1) } proc vectorComponentZero(vnum, int nnum) "USAGE: vectorComponentZero(vnum, nnum); with vnum: intvec/vector of numbers, nnum: int. ASSUME: nnum is the number of entries of the passed vector, thus nnum = nrows(vnum) for intvecs. PURPOSE: Marks all those entries of a vector positively which equal zero. RETURN: An intvec of flags whose components are: result[i] == (vnum[i] != 0). EXAMPLE: example vectorComponentZero; shows an example." { //Check whether the passed argument is of appropriate type. if (typeof(vnum) != "intvec" && typeof(vnum) != "vector") { ERROR("The argument passed must be of type \"intvec\" or \"vector\"."); } if (typeof(vnum) == "intvec" && nnum != nrows(vnum)) { ERROR("The passed intvec's size and the passed size differ."); } int inum; //Check which entries of the passed vector are zero. intvec vfResult = 0:nnum; for (inum = 1; inum <= nnum; inum++) { vfResult[inum] = (vnum[inum] == 0); } return(vfResult); } example { "EXAMPLE:"; echo=2; //intvec example intvec vf = 1,0; intvec vfResult = vectorComponentZero(vf, 2); print(vfResult); //Should yield the intvec (0,1). //vector example ring R = (0,a),T,dp; vector vnum = [1/2, 0, a]; intvec vfResult2 = vectorComponentZero(vnum, 4); print(vfResult2); //Should yield the intvec (0,1,0,1). } proc vectorFromIntvec(intvec vz) "USAGE: vectorFromIntvec(vz); with vz: intvec. PURPOSE: Converts an intvec to a vector. RETURN: The vector resulting by converting the passed intvec's components from int to number. EXAMPLE: example vectorFromIntvec; shows an example." { //Make sure a basering is defined if (!defined(basering)) { ERROR("This function requires a basering to be defined before being called."); } int iz; int nz = nrows(vz); vector vnumResult = [0]; for (iz = 1; iz <= nz; iz++) { vnumResult = vnumResult + vz[iz] * gen(iz); } return(vnumResult); } example { "EXAMPLE"; echo=2; ring R = (0,a),T,dp; intvec vz = 1,2,0; vector vnum = vectorFromIntvec(vz); print(vnum); //Should yield the vector [1, 2]. } proc vectorFromMatrix(matrix cvnum) "USAGE: vectorFromMatrix(cvnum); with cvnum: matrix of numbers. ASSUME: cvnum is a matrix consisting of exactly one column. PURPOSE: Converts a matrix consisting of one column to a vector. RETURN: The vector that equalling the passed matrix. EXAMPLE: example vectorFromMatrix; shows an example." { if (ncols(cvnum) != 1) { ERROR("The passed matrix must contain exactly one column."); } int inum; int nnum = nrows(cvnum); vector vResult; for (inum = 1; inum <= nnum; inum++) { vResult = vResult + cvnum[inum, 1] * gen(inum); } return(vResult); } example { "EXAMPLE:"; echo=2; ring R = 0,x,dp; matrix cvnum[3][1] = 1, 2, 3; vector vnum = vectorFromMatrix(cvnum); print(vnum); //Should yield the vector [1, 2, 3]. } proc vectorTake(vector vnum, intvec vfTake) "USAGE: vectorTake(vnum, vfTake); with vnum: vector, vfTake: intvec of flags. PURPOSE: Generates a vector by taking specified entries of an original vector. RETURN: A vector consisting of those entries of the first passed vector that were specified by the second passed vector. EXAMPLE: example vectorTake; shows an example." { //Check whether the sizes of the two passed vectors match. int ifTake; int nfTake = nrows(vfTake); //Make sure at least one index to take is specified. if (vfTake == 0:nfTake) { ERROR("The vector specfifying the indices to take must not be zero."); } //Copy the vector's specified entries to a new intvec. int inumResult = 1; vector vnumResult = [0]; for (ifTake = 1; ifTake <= nfTake; ifTake++) { if (vfTake[ifTake]) { vnumResult = vnumResult + vnum[ifTake] * gen(inumResult); inumResult++; } } return(vnumResult); } example { "EXAMPLE:"; echo=2; ring R = (0,a),T,dp; vector vnum = [1, 2, a, 4]; intvec vf = 1,0,1,0; vector vnumResult = vectorTake(vnum, vf); print(vnumResult); //Should yield the vector [1, a]. } /*********************************************** * Polynomial handling * ***********************************************/ proc containedVariables(poly p) "USAGE: containedVariables(p); with p: poly. PURPOSE: Collects the indices of all variables that are contained within the passed polynomial with a positive exponent. RETURN: An intvec vf with entries vf[i] = 1 if the i-th variable is contained in the passed polynomial and vf[i] = 0 else. EXAMPLE: example containedVariables; shows an example." { //Initialize the resulting vector intvec vfResult = 0:nvars(basering); //Determine the variables used term by term while(p != 0) { vfResult = intvecComponentOr(vfResult, leadexp(p)); p = p - lead(p); } return(vfResult); } example { "EXAMPLE:"; echo=2; ring R = 0,T(1..5),dp; poly p1 = T(1)*T(3)^2 + T(4)^3; intvec vfResult1 = containedVariables(p1); print(vfResult1); //Should return (1,0,1,1,0) poly p2 = 0; intvec vfResult2 = containedVariables(p2); print(vfResult2); //Should return (0,0,0,0,0) } proc evaluatePoly(poly p, vnum) "USAGE: evaluatePoly(p, vnum); with p: poly, vnum: intvec/vector of numbers only. ASSUME: The size of vnum matches the number of ring variables of the currently defined basering. PURPOSE: Calculates p(vnum). RETURN: A number p(vnum). EXAMPLE: example evaluatePoly; shows an example." { //Check whether the second parameter is of an appropriate vector type if (typeof(vnum) != "intvec" && typeof(vnum) != "vector") { ERROR("The second argument passed must be of type \"intvec\" or \"vector\"."); } //Check whether the vector size matches the number of ring variables int imonVars; int nmonVars = nvars(basering); if (typeof(vnum) == "intvec" && nrows(vnum) != nmonVars) { ERROR("The number of the passed point's coordinates and the number of the passed ring's variables differ."); } //Check whether all entries of a vector are numbers if (typeof(vnum) == "vector") { for (imonVars = 1; imonVars <= nmonVars; imonVars++) { if (leadexp(vnum[imonVars]) != 0:nmonVars) { ERROR("The vector passed must consist of numbers only."); } } } //Evaluate for each variable sequentially for (imonVars = 1; imonVars <= nmonVars; imonVars++) { p = subst(p, var(imonVars), vnum[imonVars]); } return(number(p)); } example { "EXAMPLE:"; echo=2; //intvec example ring R = (0,a),T(1..6),dp; poly p = 3/2*T(1)*T(2)^2 + T(3)^2*T(4); intvec vz = 1,2,3,0,4,0; number numResult = evaluatePoly(p, vz); print(numResult); //Should return 3/2*1*2^2+3^2*0 = 6 //vector example vector vnum = [2/3, a, 1/2, 0, 4]; numResult = evaluatePoly(p, vnum); print(numResult); //Should return 3/2*2/3*a^2+(1/2)^2*0 = a^2 } proc exponentMatrixFromPoly(poly p) "USAGE: exponentMatrixFromPoly(p); with p: poly. PURPOSE: Lists the exponents of all polynomial's terms as rows of a matrix. RETURN: An intmat containing the polynomial's terms' exponents as follows: Each row represents one of the polynomial's terms and each coulumn represents one of the ring's variables. The (i,j)-th matrix entry therefore is the exponent of the j-th variable in the i-th term. EXAMPLE: example exponentMatrixFromPoly; shows an example." { //For each term of the passed polynomial, extract its leading exponent as a vector and add it to the resulting matrix as a further row. intmat cvrvzResult = intmat(leadexp(p), 1, nvars(basering)); p = p - lead(p); while (p != 0) { cvrvzResult = intmatAppendRow(cvrvzResult, leadexp(p)); p = p - lead(p); } return(cvrvzResult); } example { "EXAMPLE:"; echo=2; ring R = 0,T(1..4),lp; poly p1 = 2*T(1)*T(2)^3 - 1/2*T(3) + T(1)^2*T(3)*T(4)^3; intmat cvrvzResult1 = exponentMatrixFromPoly(p1); print(cvrvzResult1); //Should yield an intmat with the rows // (2 0 1 3) // (1 3 0 0) // (0 0 1 0) poly p2 = T(1) + 1; intmat cvrvzResult2 = exponentMatrixFromPoly(p2); print(cvrvzResult2); //Should yield the intmat with the rows // (1 0 0 0) // (0 0 0 0) poly p3 = 0; intmat cvrvzResult3 = exponentMatrixFromPoly(p3); print(cvrvzResult3); //Should yield the intmat with the rows // (0 0 0 0) } proc getDegree(intmat rvcvzQ, poly hpQ) "USAGE: getDegree(rvcvzQ, hpQ); with rvcvzQ: intmat, hpQ: poly. ASSUME: hpQ is homogeneous w.r.t. the basering's variables' degrees implied by rvcvzQ. PURPOSE: Calculates the degree of the passed polynomial which is homogeneous w.r.t. to the grading given by the passed matrix. RETURN: An intvec that is the degree of the passed polynomial w.r.t. the grading given by the passed matrix. EXAMPLE: example getDegree; shows an example." { //The matrix's column count must match the number of the current basering's variables (each variable must be assigned a degree). if (nvars(basering) != ncols(rvcvzQ)) { ERROR("The passed matrix is not a valid grading matrix for the current basering."); } //Simply calculate the degree of the passed polynomial's leading term. return(rvcvzQ * leadexp(hpQ)); } example { "EXAMPLE:"; echo=2; ring R = 0,T(1..3),dp; poly hpQ = T(1)^10*T(2) - T(2)^3*T(3); intmat rvcvzQ[2][3] = 1,3,4, 0,-3,6; intvec cvzResult = getDegree(rvcvzQ, hpQ); print(cvzResult); //Should yield the degree vector (13,-3); } proc intersectContainedVariables(list #) "USAGE: intersectContainedVariables(sp); with sp: list of polys. PURPOSE: Collects the indices of all variables that are contained within all passed polynomials with a positive exponent. RETURN: An intvec vf with entries vf[i] = 1 if the i-th variable is contained in all passed polynomials and vf[i] = 0 else. EXAMPLE: example intersectContainedVariables; shows an example." { //Make sure a basering is defined if (!defined(basering)) { ERROR("This function requires a basering to be defined before being called."); } int nmonVars = nvars(basering); list sp = #; int ip; int np = size(sp); //Special case: An empty list was passed, therefore no variable was used. if (np == 0) { return(0:nmonVars); } //Initialize the resulting vector intvec vfResult = 1:nmonVars; //Collect the exponents used by each polynomial and //calculate the intersection with the exponents already used for (ip = 1; ip <= np; ip++) { vfResult = intvecComponentAnd(vfResult, containedVariables(sp[ip])); } return(vfResult); } example { "EXAMPLE:"; echo=2; ring R = 0,T(1..5),dp; poly p1 = T(1)*T(3)^2 + T(4)^3; poly p2 = T(1)^2*T(2)^2 + T(4)^3 + T(5); list sp1 = p1, p2; intvec vzResult1 = intersectContainedVariables(sp1); print(vzResult1); //Should return the vector (1,0,0,1,0). poly p3 = T(5); poly p4 = 1; list sp2 = p3, p4; intvec vzResult2 = intersectContainedVariables(sp2); print(vzResult2); //Should return the vector (0,0,0,0,0). intvec vzResult3 = intersectContainedVariables(p3); print(vzResult3); //Should return the vector (0,0,0,0,1). intvec vzResult4 = intersectContainedVariables(list()); print(vzResult4); //Should return the vector (0,0,0,0,0). } proc isOfTotalDegree(int nDeg, poly p) "USAGE: isOfTotalDegree(nDeg, p); with nDeg: int, p: poly. PURPOSE: Checks whether a polynomial is homogeneous and of a specific degree w.r.t. standard grading. RETURN: 1 if the passed polynomial is homogeneous and of the passed degree w.r.t. standard grading, 0 else. EXAMPLE: example isOfTotalDegree; shows an example." { //First, check whether the passed polynomial is homogeneous at all. if (!isTotalDegreeHomogeneous(p)) { return(0); } int nmonVars = nvars(basering); //The total grading is a special case of a grading given by a matrix's columns (all columns equal 1). intmat rvcvzQ = intmat(1:nmonVars, 1, nmonVars); //An intvec with one component is returned. int nLeadDeg = getDegree(rvcvzQ, p)[1]; //Degree check. return(nLeadDeg == nDeg); } example { "EXAMPLE:"; echo=2; ring R = 0,T(1..3),dp; poly p1 = T(1)*T(2)^3; int fResult1a = isOfTotalDegree(3, p1); print(fResult1a); //Should be false. int fResult1b = isOfTotalDegree(4, p1); print(fResult1b); //Should be true. poly p2 = T(1) + T(3); int fResult2a = isOfTotalDegree(2, p2); print(fResult2a); //Should be false. int fResult2b = isOfTotalDegree(1, p2); print(fResult2b); //Should be true. poly p3 = 1; int fResult3a = isOfTotalDegree(1, p3); print(fResult3a); //Should be false. int fResult3b = isOfTotalDegree(0, p3); print(fResult3b); //Should be true. poly p4 = T(1)*T(2) + T(3); int fResult4 = isOfTotalDegree(1, p4); print(fResult4); //Should be false, as this is not a homogeneous polynomial. } proc isTotalDegreeHomogeneous(poly p) "USAGE: isTotalDegreeHomogeneous(p); with p: poly. PURPOSE: Determines whether a polynomial is homogeneous w.r.t. standard grading. RETURN: 1 if the passed polynomial is homogeneous w.r.t. standard grading, 0 else. EXAMPLE: example isTotalDegreeHomogeneous; shows an example." { //Special case: p = 0 is homogeneous if (p == 0) { return(1); } int nmonVars = nvars(basering); int imonP; int nmonP = size(p); //The total grading is a special case of a grading given by a matrix's columns (all columns equal 1). intmat rvcvzQ = intmat(1:nmonVars, 1, nmonVars); //Calculate the degree of the lead term of the polynomial int nLeadDeg = getDegree(rvcvzQ, lead(p))[1]; //Match the degrees of the further polynomial's terms with the first term's degree. for (imonP = 2; imonP <= nmonP; imonP++) { if (getDegree(rvcvzQ, p[imonP])[1] != nLeadDeg) { //Degree mismatch found -> passed polynomial not homogeneous return(0); } } //No degree mismatch found return(1); } example { "EXAMPLE:"; echo=2; ring R = 0,T(1..3),dp; poly p1 = T(1)*T(2)^3; isTotalDegreeHomogeneous(p1); //Should be true. poly p2 = T(1)*T(2) + T(3)^2; isTotalDegreeHomogeneous(p2); //Should be true. poly p3 = 0; isTotalDegreeHomogeneous(p3); //Should be true. poly p4 = T(1)*T(2) + T(3); isTotalDegreeHomogeneous(p4); //Should be false. } proc veronese(int nDeg, vector vnum) "USAGE: veronese(nDeg, vnum); with nDeg: int, vnum: vector of numbers. PURPOSE: Applies the veronese embedding of a certain degree on a vector. RETURN: A vector of numbers that is the result of the veronese embedding of the passed degree on the passed vector. EXAMPLE: example veronese; shows an example." { int imonVars; int nmonVars = nvars(basering); //Check for positive degree. if (nDeg < 1) { ERROR("The degree passed must be positive."); } //Check whether all entries of the passed vector are numbers. for (imonVars = 1; imonVars <= nmonVars; imonVars++) { if (leadexp(vnum[imonVars]) != 0:nmonVars) { ERROR("The vector passed must consist of numbers only."); } } //Get all monomials of the passed degree ideal smon = maxideal(nDeg); int imon; int nmon = ncols(smon); //Evaluate all monomials in the passed vector, put the results into a vector. vector vnumResult = 0; for (imon = 1; imon <= nmon; imon++) { vnumResult = vnumResult + evaluatePoly(smon[imon], vnum) * gen(nmon - imon + 1); } return(vnumResult); } example { "EXAMPLE:"; echo=2; ring R = 0,T(1..3),dp; vector vnum = [1, 2, 3]; vector vnumResult = veronese(3, vnum); print(vnumResult); //Should yield the vector // [1*1*1, 1*1*2, 1*1*3, 1*2*2, 1*2*3, 1*3*3, 2*2*2, 2*2*3, 2*3*3, 3*3*3] // = [1, 2, 3, 4, 6, 9, 8, 1 2, 18, 27] } /*********************************************** * Geometric functions * ***********************************************/ proc isPointInVariety(ideal spI, vnum) "USAGE: isPointInVariety(spI, vnum); with spI: ideal, vnum: intvec/vector numbers. PURPOSE: Checks whether a point is contained in the variety V(K^n; spI) of the passed ideal spI. RETURN: The int 1 if the passed point is contained in the variety of the passed ideal, 0 else. EXAMPLE: example isPointInVariety; shows an example." { //Check whether the second parameter is of an appropriate vector type if (typeof(vnum) != "intvec" && typeof(vnum) != "vector") { ERROR("The second parameter passed must be of type \"intvec\" or \"vector\""); } //The passed point's component count must match the number of variables in the current basering. int imonVars; int nmonVars = nvars(basering); if (typeof(vnum) == "intvec" && nmonVars != size(vnum)) { ERROR("The passed point's number of entries and the current basering's count of variables differ."); } //Check whether all entries of a vector are numbers if (typeof(vnum) == "vector") { for (imonVars = 1; imonVars <= nmonVars; imonVars++) { if (leadexp(vnum[imonVars]) != 0:nmonVars) { ERROR("The vector passed must consist of numbers only."); } } } int ipI; int npI = ncols(spI); //Check whether the passed point vanishes under each polynomial of the passed ideal for (ipI = 1; ipI <= npI; ipI++) { if (evaluatePoly(spI[ipI], vnum) != 0) { return(0); } } return(1); } example { "EXAMPLE:"; echo=2; //intvec example ring R = 0,T(1..3),dp; poly p1 = T(1) + T(2); poly p2 = T(3) - T(2)^2; ideal sp = p1, p2; intvec vz1 = 1,-1,1; int f1 = isPointInVariety(sp, vz1); print(f1); //Should be true. intvec vz2 = 1,1,1; int f2 = isPointInVariety(sp, vz2); print(f2); //Should be false. //vector example vector vnum1 = [1/2, -1/2, 1/4]; int f3 = isPointInVariety(sp, vnum1); print(f3); //Should be true. vector vnum2 = [1, 1/2]; int f4 = isPointInVariety(sp, vnum2); print(f4); //Should be false. } proc hyperplanes(list svnumPts) "USAGE: hyperplanes(svnumPts); with svnumPts: list of vectors of numbers. ASSUME: The current basering is the polynomial ring in n >= 2 variables. PURPOSE: Computes all hyperplanes through the passed points. RETURN: A list of vectors containing the normal vectors on the resulting hyperplanes. EXAMPLE: example hyperplanes; shows an example." { //Make sure a basering is defined if (!defined(basering)) { ERROR("This function needs a basering to be defined before being called."); } int nmonVars = nvars(basering); if (nmonVars < 2) { ERROR("This function requires a basering with at least two variables."); } int iDim; int nDim = nmonVars - 1; int ivnumPts; int nvnumPts = size(svnumPts); //Check whether enough points were passed. if (nvnumPts < nDim) { ERROR("There were not enough points passed."); } //Compute the set of all linear forms whose zero sets contain at least nDim of the passed points. //Generate the list of all the combinations of n points to iterate through: // - Initialize first combination list ssivnumPts = list(); list sivnumPts = list(); for (ivnumPts = nDim; ivnumPts >= 1; ivnumPts--) { sivnumPts = insert(sivnumPts, ivnumPts); } while (sivnumPts[1] <= nvnumPts - nDim + 1) { // - Add current combination ssivnumPts = insert(ssivnumPts, sivnumPts); // - Modify current combination, obtaining the next one sivnumPts[nDim] = sivnumPts[nDim] + 1; if (sivnumPts[nDim] > nvnumPts) { iDim = nDim; while (iDim > 1 && sivnumPts[iDim] > nvnumPts - nDim + iDim) { iDim--; sivnumPts[iDim] = sivnumPts[iDim] + 1; } for (iDim = iDim + 1; iDim <= nDim; iDim++) { sivnumPts[iDim] = sivnumPts[iDim - 1] + 1; } } } //Compute the normal vectors on the hyperplanes defined by n points. int isivnumPts; int nsivnumPts = size(ssivnumPts); int iivnumPts; matrix rvcvnumCurrentHyperplane; matrix cvrvnumCurrentHyperplane; matrix cvnumZero = matrix(0:nDim, nDim, 1); list smatnumPLU = list(); list smatnumLUResult = list(); vector vnumNormal = 0; list svnumNormal = list(); for (isivnumPts = nsivnumPts; isivnumPts >= 1; isivnumPts--) { // - Given n points, put their coordinates in a row matrix A = p_1, ..., p_n. rvcvnumCurrentHyperplane = matrix(veronese(1, svnumPts[ssivnumPts[isivnumPts][1]]), nmonVars, 1); for (iivnumPts = 2; iivnumPts <= nDim; iivnumPts++) { rvcvnumCurrentHyperplane = concat(rvcvnumCurrentHyperplane, matrix(veronese(1, svnumPts[ssivnumPts[isivnumPts][iivnumPts]]), nmonVars, 1)); } cvrvnumCurrentHyperplane = transpose(rvcvnumCurrentHyperplane); // - Does A have full rank <=> Do p_1, ..., p_n define a hyperplane? smatnumPLU = ludecomp(cvrvnumCurrentHyperplane); if (submat(smatnumPLU[3], nDim, 1..nmonVars) == matrix(0:nmonVars, 1, nmonVars)) { isivnumPts--; continue; } // - The normal vector on the hyperplane defined by these n points is a vector x != 0 with Ax = 0. smatnumLUResult = lusolveUnderdetFullRank(smatnumPLU[1], smatnumPLU[2], smatnumPLU[3], cvnumZero); // - Retrieve normal vector vnumNormal = vectorFromMatrix(smatnumLUResult[2]); if (vnumNormal != 0) { svnumNormal = insert(svnumNormal, vnumNormal); } } return(svnumNormal); } example { "EXAMPLE:"; echo=2; ring R = 0,T(1..3),dp; vector vnumPt1 = [1,0,0]; vector vnumPt2 = [0,1,0]; vector vnumPt3 = [0,0,1]; vector vnumPt4 = [1,1,1]; list svnumPts = vnumPt1, vnumPt2, vnumPt3, vnumPt4; list svnumNormal = hyperplanes(svnumPts); print(svnumNormal); //Should yield a list of normal vectors like //[1,0,0], [0,1,0], [0,0,1], [0,1,-1], [1,0,-1] and [1,-1,0]. } proc quadrics(list svnumPts) "USAGE: quadrics(svnumPts); with svnumPts: list of vectors of numbers. ASSUME: The current basering is the polynomial ring in n >= 2 variables. PURPOSE: Computes all quadrics through the passed points. RETURN: A list of vectors containing the coefficient vectors for the resulting quadrics (for processing with veronesePoly). EXAMPLE: example quadrics; shows an example." { //Make sure a basering is defined if (!defined(basering)) { ERROR("This function needs a basering to be defined before being called."); } int nmonVars = nvars(basering); if (nmonVars < 2) { ERROR("This function requires a basering with at least two variables."); } int nDim = nmonVars - 1; int ivnumPts; int nvnumPts = size(svnumPts); int iDefPts; int nDefPts = ((nDim + 1) * (nDim + 2)) div 2 - 1; //Check whether enough points were passed. if (nvnumPts < nDefPts) { ERROR("There were not enough points passed."); } //Compute the set of all quadrics whose zero sets contain at least nDefPts of the passed points. //Generate the list of all the combinations of nDefPts points to iterate through: // - Initialize first combination list ssivnumPts = list(); list sivnumPts = list(); for (ivnumPts = nDefPts; ivnumPts >= 1; ivnumPts--) { sivnumPts = insert(sivnumPts, ivnumPts); } while (sivnumPts[1] <= nvnumPts - nDefPts + 1) { // - Add current combination ssivnumPts = insert(ssivnumPts, sivnumPts); // - Modify current combination, obtaining the next one sivnumPts[nDefPts] = sivnumPts[nDefPts] + 1; if (sivnumPts[nDefPts] > nvnumPts) { iDefPts = nDefPts; while (iDefPts > 1 && sivnumPts[iDefPts] > nvnumPts - nDefPts + iDefPts) { iDefPts--; sivnumPts[iDefPts] = sivnumPts[iDefPts] + 1; } for (iDefPts = iDefPts + 1; iDefPts <= nDefPts; iDefPts++) { sivnumPts[iDefPts] = sivnumPts[iDefPts - 1] + 1; } } } //Compute the coefficient vectors for the quadrics defined by nDefPts points. int isivnumPts; int nsivnumPts = size(ssivnumPts); int iivnumPts; matrix rvcvnumCurrentQuadric; matrix cvrvnumCurrentQuadric; matrix cvnumZero = matrix(0:nDefPts, nDefPts, 1); list smatnumPLU = list(); list smatnumLUResult = list(); vector vnumCoeffs = 0; list svnumCoeffs = list(); for (isivnumPts = nsivnumPts; isivnumPts >= 1; isivnumPts--) { // - Given nDefPts points, put their coordinates in a row matrix A = p_1, ..., p_nDefPts. rvcvnumCurrentQuadric = matrix(veronese(2, svnumPts[ssivnumPts[isivnumPts][1]]), nDefPts + 1, 1); for (iivnumPts = 2; iivnumPts <= nDefPts; iivnumPts++) { rvcvnumCurrentQuadric = concat(rvcvnumCurrentQuadric, matrix(veronese(2, svnumPts[ssivnumPts[isivnumPts][iivnumPts]]), nDefPts + 1, 1)); } cvrvnumCurrentQuadric = transpose(rvcvnumCurrentQuadric); // - Does A have full rank <=> Do p_1, ..., p_nDefPts define a quadric? smatnumPLU = ludecomp(cvrvnumCurrentQuadric); if (submat(smatnumPLU[3], nDefPts, 1..(nDefPts+1)) == matrix(0:(nDefPts + 1), 1, nDefPts + 1)) { isivnumPts--; continue; } // - The normal vector on the hyperplane defined by these n points is a vector x != 0 with Ax = 0. smatnumLUResult = lusolveUnderdetFullRank(smatnumPLU[1], smatnumPLU[2], smatnumPLU[3], cvnumZero); // - Retrieve coefficient vector vnumCoeffs = vectorFromMatrix(smatnumLUResult[2]); if (vnumCoeffs != 0) { svnumCoeffs = insert(svnumCoeffs, vnumCoeffs); } } return(svnumCoeffs); } example { "EXAMPLE:"; echo=2; ring R = (0,a,b),T(1..3),dp; vector vnumPt1 = [1,0,0]; vector vnumPt2 = [0,1,0]; vector vnumPt3 = [0,0,1]; vector vnumPt4 = [1,1,1]; vector vnumPt5 = [1,a,b]; list svnumPts = vnumPt1, vnumPt2, vnumPt3, vnumPt4, vnumPt5; list svnumNormal = quadrics(svnumPts); print(svnumNormal); //Should yield a list of cofficient vectors with one entry: //[0, (a*b-a-b)/(a-b), (-a*b+2*a)/(a-b), 0, -1, 0]. } /*********************************************** * Cone and fan operations * ***********************************************/ proc containsRay(fan scn, intvec vzRay) "USAGE: containsRay(scn, vzRay); with scn: fan, vzRay: intvec. ASSUME: The passed ray's component count matches the fan's ambient dimension. PURPOSE: Checks whether a ray is contained in a fan. RETURN: The int 1 if the passed fan contains the passed ray, 0 else. EXAMPLE: example containsRay; shows an example." { //Make sure the passed fan is defined if (size(fVector(scn)) == 0) { ERROR("The passed fan is not defined."); } //The passed vector's component count must match the fan's ambient dimension. int nzRay = size(vzRay); if (nzRay != ambientDimension(scn)) { ERROR("The passed vector's component count does not match the passed fan's ambient dimension."); } //Transform the passed vector to an actual ray (as a cone) intmat cvrvzRay = intmat(vzRay, 1, nzRay); cone cnRay = coneViaPoints(cvrvzRay); //For each ray of the passed fan, check whether it contains (and therefore equals) the passed ray int iRays; int nRays = numberOfConesOfDimension(scn, 1, 0, 0); for (iRays = 1; iRays <= nRays; iRays++) { if (containsAsFace(getCone(scn, 1, iRays), cnRay)) { return(1); } } return(0); } example { "EXAMPLE:"; echo=2; intmat cvrvz1[2][2] = 1,0, 0,1; intmat cvrvz2[2][2] = 0,1, -1,0; cone cn1 = coneViaPoints(cvrvz1); cone cn2 = coneViaPoints(cvrvz2); fan scn = fanViaCones(cn1, cn2); intvec vzRay1 = 0,1; int fResult1 = containsRay(scn, vzRay1); print(fResult1); //Should be true. intvec vzRay2 = 0,-1; int fResult2 = containsRay(scn, vzRay2); print(fResult2); //Should be false. intvec vzRay3 = 1,1; int fResult3 = containsRay(scn, vzRay3); print(fResult3); //Should be false. } proc contractRay(fan scn, intvec vzRay) "USAGE: contractRay(scn, vzRay); with scn: fan, vzRay: intvec. ASSUME: The passed ray is indeed contained in the passed fan and contractible. PURPOSE: Computes the fan obtained by contracting one of the rays of the original fan. RETURN: The original fan with the passed ray contracted. EXAMPLE: example contractRay; shows an example." { //Make sure the passed fan is defined if (size(fVector(scn)) == 0) { ERROR("The passed fan is not defined."); } //Make sure the passed ray is contained in the passed fan if (!(containsRay(scn, vzRay))) { ERROR("The passed vector is no ray of the passed fan."); } //Helper variables int iDim; int nDim = dimension(scn); int icnMax; int ncnMax; cone cnMax; list cnMaxList = list(); //Initialize the contraction cone and the list of cones holding all maximal cones to go into the resulting fan cone cnContractionSum = coneViaPoints(intmat(0:nDim, 1, nDim)); list scnResultList = list(); //Iterate over the cone dimension for (iDim = 1; iDim <= nDim; iDim++) { //Iterate over the maximal cones of a dimension ncnMax = numberOfConesOfDimension(scn, iDim, 0, 1); for (icnMax = 1; icnMax <= ncnMax; icnMax++) { cnMax = getCone(scn, iDim, icnMax, 0, 1); //The maximal cones that do not contain the ray to contract simply go into the resulting fan. //All other cones are added to the large contraction cone. if (intmatContainsRow(rays(cnMax), vzRay)) { cnContractionSum = convexHull(cnContractionSum, cnMax); } else { cnMaxList = cnMax; scnResultList = scnResultList + cnMaxList; } } } //Insert the contraction cone to the list of cones that go into the resulting fan. list cnContractionSumList = cnContractionSum; scnResultList = scnResultList + cnContractionSumList; //Finally, build the fan. fan scnResult = fanViaCones(scnResultList); return(scnResult); } example { "EXAMPLE:"; echo=2; //The fan of the blowup of P2 in [1, 0, 0]. We contract the ray (-1, 0) added in the process of blowing up P2, obtaining P2 itself. intmat cvrvz1[2][2] = 1,0, 0,1; intmat cvrvz2[2][2] = 0,1, -1,0; intmat cvrvz3[2][2] = -1,0, -1,-1; intmat cvrvz4[2][2] = -1,-1, 1,0; cone cn1 = coneViaPoints(cvrvz1); cone cn2 = coneViaPoints(cvrvz2); cone cn3 = coneViaPoints(cvrvz3); cone cn4 = coneViaPoints(cvrvz4); fan scn = fanViaCones(cn1, cn2, cn3, cn4); intvec vzRay = -1,0; fan scnContracted = contractRay(scn, vzRay); scnContracted; //Should yield the fan with the maximal cones cn1, cn2 + cn3, cn4. } proc generateOrthantFace(intvec vfIndices) "USAGE: generateOrthantFace(vfIndices); with vfIndices: intvec. PURPOSE: Calculates the orthant face generated by the canonical basis vectors flagged true in the passed vector. RETURN: A cone that is the orthant face generated by the canonical basis vectors flagged true in the passed vector. EXAMPLE: example generateOrthantFace; shows an example." { //Create the rays' matrix for defining the resulting cone. intmat cvrvzRays = unitMatrix(size(vfIndices)); cvrvzRays = intmatTakeRows(cvrvzRays, vfIndices); //Build that cone. cone cnResult = coneViaPoints(cvrvzRays); return(cnResult); } example { "EXAMPLE:"; echo=2; intvec vfIndices = 1,0,1,1,0; cone cnResult = generateOrthantFace(vfIndices); rays(cnResult); //Should yield the cone defined by the rays through (1,0,0,0,0), (0,0,1,0,0) and (0,0,0,1,0). } proc irrelevantIdealFromFan(fan sof) "USAGE: irrelevantIdealFromFan(sof); with sof: fan consisting solely of faces of the positive orthant. PURPOSE: Calculates the irrelevant ideal that is associated with a fan consisting of faces of the positive orthant. RETURN: Returns the ring in which the monomials generating the irrelevant ideal were built in. SIDE EFFECTS: Exports an ideal gsmonResult containing the irrelevant ideal that is associated with the passed fan consisting of faces of the positive orthant. EXAMPLE: example irrelevantIdealFromFan; shows an example." { //Make sure the passed fan is defined if (size(fVector(sof)) == 0) { ERROR("The passed fan is not defined."); } //Set the ring in which the irrelevant ideal will be built. int iAmbientDimension; int nAmbientDimension = ambientDimension(sof); ring R = 0,T(1..nAmbientDimension),dp; int iDimension; int nDimension = dimension(sof); int iofOfDim; int nofOfDim; cone ofOfDim; intvec vfRayInComplement; intvec vnRay; list smonResultList; //Iterate through the fan's cones' dimensions for (iDimension = 1; iDimension <= nDimension; iDimension++) { //Iterate through the maximal orthant faces of the current dimension contained in the fan nofOfDim = numberOfConesOfDimension(sof, iDimension, 0, 1); for (iofOfDim = nofOfDim; iofOfDim >= 1; iofOfDim--) { ofOfDim = getCone(sof, iDimension, iofOfDim, 0, 1); //Determine which of the orthant's rays are not contained in the current face vfRayInComplement = 1:nAmbientDimension; for (iAmbientDimension = 1; iAmbientDimension <= nAmbientDimension; iAmbientDimension++) { vnRay = 0:nAmbientDimension; vnRay[iAmbientDimension] = 1; if (containsInSupport(ofOfDim, vnRay)) { vfRayInComplement[iAmbientDimension] = 0; } } //Create the monomial related to all rays not contained in the current face. smonResultList = insert(smonResultList, monomial(vfRayInComplement)); } } //Compose the resulting ideal ideal gsmonResult = smonResultList[1..size(smonResultList)]; export(gsmonResult); return(basering); } example { "EXAMPLE:"; echo=2; intmat cvrvz1[2][4] = 1,0,0,0, 0,1,0,0; intmat cvrvz2[2][4] = 0,1,0,0, 0,0,0,1; intmat cvrvz3[2][4] = 0,0,0,1, 0,0,1,0; intmat cvrvz4[2][4] = 0,0,1,0, 1,0,0,0; cone of1 = coneViaPoints(cvrvz1); cone of2 = coneViaPoints(cvrvz2); cone of3 = coneViaPoints(cvrvz3); cone of4 = coneViaPoints(cvrvz4); fan sof = fanViaCones(of1, of2, of3, of4); def R = irrelevantIdealFromFan(sof); setring R; ideal smonResult = gsmonResult; print(smonResult); //Should yield the ideal // . } proc isRayContractible(rvcvzP, int icvzRay) "USAGE: isRayContractible(rvcvzP, icvzRay); with rvcvzP: bigintmat/intmat, icvzRay: int. PURPOSE: Checks whether a ray from a list of rays is contractible, i.e. whether its corresponding ray in the gale dual version is extremal. RETURN: The int 1 if the passed ray is contractible w.r.t. the passed other rays, 0 else. EXAMPLE: example isRayContractible; shows an example." { //Check whether a bigintmat or intmat was passed for rvcvzP if (typeof(rvcvzP) != "intmat" && typeof(rvcvzP) != "bigintmat") { ERROR("The first passed argument has to be a matrix of integers."); } //Check whether the ray marked is contractible via the equivalent description as extremal ray of the dual cone intmat cvrvzQ = kernelLattice(rvcvzP); int irvzRay = icvzRay; //Transposition: former column indices are now row indices intvec rvzOmitted = intmatTakeRow(cvrvzQ, irvzRay); intmat cvrvzQContracted = intmatDeleteRow(cvrvzQ, irvzRay); return(!(containsInSupport(coneViaPoints(cvrvzQContracted), rvzOmitted))); } example { "EXAMPLE:"; echo=2; //bigintmat example bigintmat rvcvzP1[2][4] = 1,0,-1,-1, 0,1,-1,0; int fResult1 = isRayContractible(rvcvzP1, 3); print(fResult1); //Should be false. //intmat example intmat rvcvzP2[2][4] = 1,0,-1,-1, 0,1,-1,0; int fResult2 = isRayContractible(rvcvzP2, 4); print(fResult2); //Should be true. } proc mapFan(rvcvzP, fan scn) "USAGE: mapFan(rvcvzP, scn); with rvcvzP: bigintmat/intmat, scn: fan. PURPOSE: Maps a fan under a matrix by mapping each cone under that matrix. RETURN: A fan that is the original fan mapped by rvcvzP. EXAMPLE: example mapFan; shows an example." { //Make sure the passed fan is defined if (size(fVector(scn)) == 0) { ERROR("The passed fan is not defined."); } //Check whether a bigintmat or intmat was passed for rvcvzP if (typeof(rvcvzP) != "intmat" && typeof(rvcvzP) != "bigintmat") { ERROR("The first passed argument has to be a matrix of integers."); } int iDim; int nDim = dimension(scn); //Check whether the passed matrix's dimensions matches the fan's dimensions. if (ncols(rvcvzP) != ambientDimension(scn) || nrows(rvcvzP) != dimension(scn)) { ERROR("The passed matrix is not in relation with the passed fan's dimensions."); } int icnMax; int ncnMax; cone cnMax; cone cnMaxMapped; list cnMaxMappedList; bigintmat cvrvzMappedRays; //Initialize the list of cones to go into the resulting fan list scnList = list(); //Iterate through the cones' dimensions for (iDim = 1; iDim <= nDim; iDim++) { //Iterate through the maximal cones of the current dimension ncnMax = numberOfConesOfDimension(scn, iDim, 0, 1); for (icnMax = 1; icnMax <= ncnMax; icnMax++) { cnMax = getCone(scn, iDim, icnMax, 0, 1); //Extract the cone's rays and map them by rvcvzP cvrvzMappedRays = rays(cnMax) * transpose(rvcvzP); cnMaxMapped = coneViaPoints(cvrvzMappedRays); //Add the cone arising from the mapped rays to the list of cones to go into the resulting fan cnMaxMappedList = cnMaxMapped; scnList = scnList + cnMaxMappedList; } } //Finally, build the fan fan scnResult = fanViaCones(scnList); return(scnResult); } example { "EXAMPLE:"; echo=2; intmat cvrvz1[2][3] = 1,0,0, 0,1,0; intmat cvrvz2[2][3] = 0,1,0, 0,0,1; intmat cvrvz3[2][3] = 0,0,1, 1,0,0; cone cn1 = coneViaPoints(cvrvz1); cone cn2 = coneViaPoints(cvrvz2); cone cn3 = coneViaPoints(cvrvz3); fan scn = fanViaCones(cn1, cn2, cn3); //bigintmat example bigintmat rvcvzP1[2][3] = 1,0,-1, 0,1,-1; fan scnResult1 = mapFan(rvcvzP1, scn); scnResult1; //Should be the fan with the rays (1,0), (0,1), (-1,-1) and full support. //intmat example intmat rvcvzP2[2][3] = 1,0,-1, 0,1,-1; fan scnResult2 = mapFan(rvcvzP2, scn); scnResult2; //Should be the fan with the rays (1,0), (0,1), (-1,-1) and full support. } proc movingConeFromWeights(rvcvzQ) "USAGE: movingConeFromWeights(rvcvzQ); with rvcvzWeights: bigintmat/intmat. PURPOSE: Calculates the moving cone of the cone defined by the column vectors of the passed intmat. RETURN: A cone that is the moving cone of the cone defined by the column vectors of the passed intmat. EXAMPLE: example movingConeFromWeights; shows an example." { //Check whether a bigintmat or intmat was passed for rvcvzQ if (typeof(rvcvzQ) != "intmat" && typeof(rvcvzQ) != "bigintmat") { ERROR("The first passed argument has to be a matrix of integers."); } //Check the column count of the passed matrix int ncvzQ = ncols(rvcvzQ); //No column -> error if (ncvzQ == 0) { ERROR("The passed matrix must have at least one column."); } //Special case: One column -> cone consisting only of the origin if (ncvzQ == 1) { int nzQ = nrows(rvcvzQ); intmat cvrvzZero = intmat(0:nzQ, 1, nzQ); cone cnResult = coneViaPoints(cvrvzZero); return(cnResult); } //Transpose the input data intmat cvrvzQ = transpose(rvcvzQ); int irvzQ; int nrvzQ = ncvzQ; //Initialize the resulting cone cone cnResult = coneViaPoints(cvrvzQ); //Calculate the cones generated by all of the matrix's rows except one and take their intersection. for (irvzQ = 1; irvzQ <= ncvzQ; irvzQ++) { cnResult = convexIntersection(cnResult, coneViaPoints(intmatDeleteRow(cvrvzQ, irvzQ))); } return(cnResult); } example { "EXAMPLE:"; echo=2; //bigintmat example bigintmat rvcvzQ1[2][4] = 1,1,0,1, 0,1,1,0; cone cnMov1 = movingConeFromWeights(rvcvzQ1); rays(cnMov1); //intmat examples //one column intmat rvcvzQ2[2][1] = 1, 0; cone cnMov2 = movingConeFromWeights(rvcvzQ2); cnMov2; //Should yield the cone consisting only of the origin in ambient dimension 2. //more columns intmat rvcvzQ3[2][4] = 1,1,0,1, 0,1,1,0; cone cnMov3 = movingConeFromWeights(rvcvzQ3); rays(cnMov3); //Should yield the cone generated by the rays passing through (1,0) and (1,1). } proc orthantFanFromWeight(rvcvzQ, intvec cvzWeight) "USAGE: orthantFanFromWeight(rvcvzQ, cvzWeight); with rvcvzQ: bigintmat/intmat, cvzWeight: intvec. ASSUME: The weight vector is an element of the relative interior of the moving cone defined by the grading matrix. PURPOSE: Calculates the fan corresponding to a grading matrix and a weight vector. RETURN: A fan that is the fan corresponding to the passed grading matrix and the passed weight vector. EXAMPLE: example orthantFanFromWeight; shows an example." { //Check whether a bigintmat or intmat was passed for rvcvzQ if (typeof(rvcvzQ) != "intmat" && typeof(rvcvzQ) != "bigintmat") { ERROR("The first passed argument has to be a matrix of integers."); } //Initialize the vector of flags describing the cone generated by which rays is considered intvec rvfCurrentRays = 0:ncols(rvcvzQ); intvec rvfNewRays; dbprint(printlevel, ""); dbprint(printlevel, "----- Starting fan computation from a weight, please be patient... -----"); if (printlevel >= 1) { int nStartTime = timer; } int icvzQ = 0; int ncvzQ = ncols(rvcvzQ); list sof = list(); list sofCurrent = list(); list stsRays = list(); list sRay = rvfCurrentRays, icvzQ; int iMaxIndex; stsRays = insert(stsRays, sRay); int nsRays = 1; cone cnComplementary; int fHadResult = 0; while (nsRays != 0) { //Stack handling sRay = stsRays[nsRays]; rvfCurrentRays = sRay[1]; iMaxIndex = sRay[2]; stsRays = delete(stsRays, nsRays); nsRays--; fHadResult = 0; //Iterate through the remaining rays that can additionally be considered. for (icvzQ = iMaxIndex + 1; icvzQ <= ncvzQ; icvzQ++) { //Add a new ray to the already considered ones. rvfNewRays = rvfCurrentRays; rvfNewRays[icvzQ] = 1; //If all rays are considered at once, the passed weight is definitely not contained in the new cone's complementary cone. //Thus, the current cone is maximal. if (rvfNewRays == 1:ncvzQ) { icvzQ++; continue; } //Build the cone complementary to the cone defined by the current rays and //compute that new cone's image under the grading matrix rvcvzQ cnComplementary = coneViaPoints(transpose(intmatTakeCols(rvcvzQ, intvecComponentNot(rvfNewRays)))); //If the complementary cone of degrees does not contain the weight, the new cone does not belong to the resulting fan. //Thus, the current cone is maximal. if (!containsRelatively(cnComplementary, cvzWeight)) { icvzQ++; continue; } //If it does, the originally defined cone belongs to the resulting fan. //But before inserting it into the list of cones for the fan, search for bigger cones containing the current one. sRay = rvfNewRays, icvzQ; stsRays = insert(stsRays, sRay); nsRays++; fHadResult = 1; //Recursive call: Find all cones containing the current one and still belonging to the resulting fan. } if (!fHadResult) { sofCurrent = generateOrthantFace(rvfCurrentRays); sof = sof + sofCurrent; } } //Compose the wanted fan from the returned cones. fan sofResult = fanViaCones(sof); if (printlevel >= 1) { int nTimeElapsed = timer - nStartTime; "Time elapsed while building the fan (in ms):"; print(nTimeElapsed); } dbprint(printlevel, "----- Fan computation done! -----"); return(sofResult); } example { "EXAMPLE:"; echo=2; //bigintmat example bigintmat rvcvzQ1[2][4] = 1,1,1,0, 1,0,0,1; cone cn1 = movingConeFromWeights(rvcvzQ1); intvec cvzWeight = intvec(relativeInteriorPoint(cn1)); fan sof1 = orthantFanFromWeight(rvcvzQ1, cvzWeight); sof1; //Should return four maximal cones in a fan: // - cone generated by rays (1,0,0,0) and (0,1,0,0) // - cone generated by rays (1,0,0,0) and (0,0,1,0) // - cone generated by rays (0,1,0,0) and (0,0,0,1) // - cone generated by rays (0,0,1,0) and (0,0,0,1). //intmat example intmat rvcvzQ2[2][4] = 1,1,1,0, 1,0,0,1; fan sof2 = orthantFanFromWeight(rvcvzQ2, cvzWeight); sof2; //Should return four maximal cones in a fan: // - cone generated by rays (1,0,0,0) and (0,1,0,0) // - cone generated by rays (1,0,0,0) and (0,0,1,0) // - cone generated by rays (0,1,0,0) and (0,0,0,1) // - cone generated by rays (0,0,1,0) and (0,0,0,1). } proc stellarSubdivision(fan scn, intvec vzNewRay) "USAGE: stellarSubdivision(scn, vzNewRay); with scn: fan, vzNewRay: intvec. ASSUME: The passed new ray is contained in the support of the passed fan. PURPOSE: Calculates the stellar subdivision of a fan w.r.t. a new ray to be inserted. RETURN: A fan that is the original fan except that stellar subdivisions were performed on cones containing the new ray to add. EXAMPLE: example stellarSubdivision; shows an example." { //Make sure the passed fan is defined if (size(fVector(scn)) == 0) { ERROR("The passed fan is not defined."); } //Check whether there is something to do. if (containsRay(scn, vzNewRay)) { return(scn); } int iDim; int nDim = dimension(scn); int icnMax; int ncnMax; int icnMaxFacets; int ncnMaxFacets; cone cnMax; fan cnMaxFan; cone cnMaxFacet; list cnHelperList = list(); //Initialize the list of cones to go into the resulting fan list scnResultList = list(); //Iterate over the the cone dimensions. for (iDim = 2; iDim <= nDim; iDim++) { //Iterate over the maximal cones of the current dimension ncnMax = numberOfConesOfDimension(scn, iDim, 0, 1); for (icnMax = 1; icnMax <= ncnMax; icnMax++) { cnMax = getCone(scn, iDim, icnMax, 0, 1); //If the current maximal cone contains the ray to add, perform subdivision. if (containsInSupport(cnMax, vzNewRay)) { //Find all facets of the current maximal cone that do not contain the ray to add cnMaxFan = fanViaCones(cnMax); ncnMaxFacets = numberOfConesOfDimension(cnMaxFan, iDim - 1, 0, 0); for (icnMaxFacets = 1; icnMaxFacets <= ncnMaxFacets; icnMaxFacets++) { cnMaxFacet = getCone(cnMaxFan, iDim - 1, icnMaxFacets); if (!(containsInSupport(cnMaxFacet, vzNewRay))) { //Generate the cone arising from the current facet not containing the ray to add //and the ray to add and insert it into the list of cones to go into the resulting fan cnHelperList = coneViaPoints(intmatAppendRow(rays(cnMaxFacet), vzNewRay)); scnResultList = scnResultList + cnHelperList; } } } //If the current maximal cone does not contain the ray to add, simply add it to the list of cones to go into the resulting fan. else { cnHelperList = cnMax; scnResultList = scnResultList + cnHelperList; } } } //Build the resulting fan fan scnResult = fanViaCones(scnResultList); return(scnResult); } example { "EXAMPLE:"; echo=2; //First example intmat cvrvz1[2][2] = 1,0, 0,1; intmat cvrvz2[2][2] = 0,1, -1,-1; intmat cvrvz3[2][2] = -1,-1, 1,0; cone cn1 = coneViaPoints(cvrvz1); cone cn2 = coneViaPoints(cvrvz2); cone cn3 = coneViaPoints(cvrvz3); fan scn = fanViaCones(cn1, cn2, cn3); intvec vzRay = -1,0; fan scnSubdivision = stellarSubdivision(scn, vzRay); scnSubdivision; //Should yield the fan with the maximal cones cn1, cn3 the two cones resulting of the division of cn2. //Second example intmat cvrvz4[4][3] = 1,0,0, 1,2,0, 1,2,2, 1,0,2; intmat cvrvz5[3][3] = 1,0,0, 1,2,0, 1,1,-1; cone cn4 = coneViaPoints(cvrvz4); cone cn5 = coneViaPoints(cvrvz5); fan scn2 = fanViaCones(cn4, cn5); intvec vzRay2 = 1,1,0; fan scnSubdivision2 = stellarSubdivision(scn2, vzRay2); scnSubdivision2; //Should yield the fan where both cones are divided, the simplicial one in two, the other one in three parts. } /*********************************************** * Ideal and poly generation * ***********************************************/ proc binomialIdeal(intmat rvcvzP) "USAGE: binomialIdeal(rvcvzP); with rvcvzP: intmat. ASSUME: The current basering has as many variables as the passed matrix has columns. PURPOSE: Calculates the lattice ideal generated by the rows of a matrix. RETURN: The ideal that is the lattice ideal associated with the passed matrix. EXAMPLE: example binomialIdeal; shows an example." { //Make sure a basering is defined if (!defined(basering)) { ERROR("This function needs a basering to be defined before being called."); } //Appropriate basering check. int icvzP; int ncvzP = ncols(rvcvzP); if (nvars(basering) != ncvzP) { ERROR("The current basering has not exactly as many variables as the passed matrix has columns."); } int izP; int nzP = nrows(rvcvzP); //Define the ring the binomial ideal will be built in poly monPositive; poly monNegative; list sbinResultList = list(); //Iterate through the rows of the passed matrix for (izP = nzP; izP >= 1; izP--) { monPositive = 1; monNegative = 1; //Iterate through a row's entries for (icvzP = 1; icvzP <= ncvzP; icvzP++) { //Multiply up all variables with the exponents given by thw passed matrix, but separated by the exponents' signs. if (rvcvzP[izP, icvzP] >= 0) { monPositive = monPositive * var(icvzP)^rvcvzP[izP, icvzP]; } else { monNegative = monNegative * var(icvzP)^(-rvcvzP[izP, icvzP]); } } //Compute the difference of the monomials associated with the positive and negative exponents sbinResultList = insert(sbinResultList, monPositive - monNegative); } //Build the resulting ideal ideal sbinResult = sbinResultList[1..size(sbinResultList)]; return(sbinResult); } example { "EXAMPLE:"; echo=2; ring R = 0,X(1..3),dp; intmat rvcvzP[2][3] = 1,0,-1, 0,2,-1; ideal sbinResult = binomialIdeal(rvcvzP); print(sbinResult); //Should yield the ideal . } static proc eliminateVariable(ideal sp, int imonVar) "USAGE: eliminateVariable(sp, imonVar); with sp: ideal, imonVar: int. ASSUME: The ideal sp is already passed in Groebner basis representation and the basering's order is an elimination order w.r.t. the variable to eliminate. PURPOSE: Computes the elimination ideal of an ideal by removal of a variable under suitable circumstances. RETURN: The elimination ideal obtained by removing the variable of passed index from the passed ideal. Under the assumptions made, this is equivalent to simply removing all those polynomials from the passed ideal that contain the variable of the passed index. EXAMPLE: example eliminateVariable; shows an example." { //Check whether the passed variable's index refers to a ring variable if (imonVar < 1 || imonVar > nvars(basering)) { ERROR("The passed index for taking a variable is out of bounds."); } ideal spResult = 0; //Select those polynomials from the ideal's generators that do not contain the variable to eliminate list spResultList = list(); int ip; int np = ncols(sp); poly p; for (ip = np; ip >= 1; ip--) { p = sp[ip]; if (leadexp(p)[imonVar] == 0) { spResultList = insert(spResultList, p); } } //Compose the result int npResultList = size(spResultList); if (npResultList > 0) { spResult = spResultList[1..npResultList]; } return(spResult); } example { "EXAMPLE:"; echo=2; ring R = 0,T(1..4),lp; poly p1 = T(2)*T(3); poly p2 = T(1)*T(3) - T(4); ideal sp = p1, p2; ideal spResult = eliminateVariable(sp, 1); print(spResult); //Should yield the ideal . ideal spResult2 = eliminateVariable(spResult, 2); print(spResult2); //Should yield the zero ideal. } proc minimalizeIdeal(ideal sp) "USAGE: minimalizeIdeal(sp); with sp: ideal. PURPOSE: Computes a minimal representation of an ideal. RETURN: The passed ideal minimally represented. NOTE: This is only a wrapper for SINGULAR's minbase function. EXAMPLE: example minimalizeIdeal; shows an example." { def R = basering; int nmonVars = nvars(basering); //Switch to a local ordering and minimalize using the minbase function execute("ring RMin = (" + charstr(R) + "),(" + varstr(R) + "),ds;"); ideal spMin = minbase(fetch(R, sp)); setring R; sp = fetch(RMin, spMin); return(sp); } example { "EXAMPLE"; echo=2; //Parameterless example ring R = 0,X(1..3),dp; poly p1 = X(1)*X(2) + X(3); poly p2 = 2*X(1)*X(2)^2 + X(3); poly p3 = 5*X(2); ideal sp = p1, p2, p3; ideal spMin = minimalizeIdeal(sp); print(spMin); //Should yield the input ideal, but with a minimal representation, such as //Example with parameter ring S = (0,a),X(1..3),dp; poly p1 = X(1)*X(2) + a*X(3); poly p2 = 2*X(1)*X(2)^2 + a*X(3); poly p3 = 5*X(2); ideal sp = p1, p2, p3; ideal spMin = minimalizeIdeal(sp); print(spMin); //Should yield the input ideal, but with a minimal representation, such as } static proc removeMonomials(list sp) "USAGE: removeMonomials(sp); with sp: list of polynomials. PURPOSE: Filters a list of polynomials: only polynomials with at least two terms are kept. RETURN: A list of polynomials that contains all passed polynomials except those with less than two terms. EXAMPLE: example removeMonomials; shows an example." { //Make sure a basering is defined if (!defined(basering)) { ERROR("This function needs a basering to be defined before being called."); } int ip; int np = size(sp); //Remove all polynomials from the passed list that have less than two terms. list spResult = sp; for (ip = np; ip >= 1; ip--) { if (size(sp[ip]) < 2) { spResult = delete(spResult, ip); } } return(spResult); } example { "EXAMPLE:"; echo=2; ring R = 0,T(1..3),dp; poly p1 = T(1); poly p2 = T(1) + T(2); poly p3 = 2; poly p4 = T(1)*T(2)^2*T(3)^3; list sp = p1,p2,p3,p4; list spFiltered = removeMonomials(sp); print(spFiltered); //Should be the list containing only T(1) + T(2). } proc toricMorphismPullback(intmat matz, ideal sp) "USAGE: toricMorphismPullback(matz, sp); with matz: intmat, sp: ideal. PURPOSE: Pulls back an ideal under a morphism of standard tori determined by a martix. RETURN: The ring (standard polynomial ring) in which the pulled back ideal is defined. SIDE EFFECTS: Exports the ideal gspResult which is the pullback of the passed ideal unter the morphism of standard tori determined by matz. As the returned ring is not a Laurent polynomial ring, care is taken that the generators of gspResult will have positive exponents only. EXAMPLE: example toricMorphismPullback; shows an example." { def R = basering; if (nvars(R) != nrows(matz)) { ERROR("The passed mapping matrix's row count and the basering's variable count differ."); } int icvz; int ncvz = ncols(matz); int ip; int np = ncols(sp); intmat cvrvzExponents; intmat cvrvzPullbackExponents; number numLeadCoef; execute("ring RResult = (" + charstr(R) + "),T(1..ncvz),dp;"); int imonTerms; int nmonTerms; poly pTermAcc; list spResultList = list(); for (ip = np; ip >= 1; ip--) { setring R; cvrvzExponents = exponentMatrixFromPoly(sp[ip]); cvrvzPullbackExponents = cvrvzExponents * matz; cvrvzPullbackExponents = intmatNonNegativeCols(cvrvzPullbackExponents); nmonTerms = size(sp[ip]); setring RResult; pTermAcc = 0; for (imonTerms = 1; imonTerms <= nmonTerms; imonTerms++) { setring R; numLeadCoef = leadcoef(sp[ip][imonTerms]); setring RResult; pTermAcc = pTermAcc + fetch(R, numLeadCoef) * monomial(intmatTakeRow(cvrvzPullbackExponents, imonTerms)); } spResultList = insert(spResultList, pTermAcc); } ideal gspResult = spResultList[1..size(spResultList)]; export(gspResult); return(RResult); } example { "EXAMPLE:"; echo=2; //Parameterless example ring R = 0,X(1..2),dp; intmat matz[2][4] = 0,-1,1,0, -1,-1,3,-3; poly p1 = X(1)*X(2) + X(1) + 3; poly p2 = X(1) + X(2); ideal sp = p1, p2; def RPullback = toricMorphismPullback(matz, sp); setring RPullback; ideal spResult = gspResult; print(spResult); //Should yield the ideal //<3*X(1)*X(2)^2*X(4)^3 + X(1)*X(2)*X(3)*X(4)^3 + X(3)^4, X(1)*X(3)*X(4)^3 + X(3)^3> //in Q[X(1)+-, ..., X(4)+-]. //Example with parameter ring S = (0,a),X(1..2),dp; poly p1 = X(1)*X(2) + X(1) + a; poly p2 = X(1) + X(2); ideal sp = p1, p2; def SPullback = toricMorphismPullback(matz, sp); setring SPullback; ideal spResult = gspResult; print(spResult); //Should yield the ideal // //in Q[X(1)+-, ..., X(4)+-]. } proc varproductOrbitClosureIdeal(intmat rvcvzP, rvnumPointInOrbit) "USAGE: varproductOrbitClosureIdeal(rvcvzP, rvnumPointInOrbit); with rvcvzP: intmat, rvnumPointInOrbit: intvec/vector of numbers. ASSUME: The current basering has as many variables as the passed matrix has columns. PURPOSE: Calculates the ideal of a point's orbit closure under a torus action implied by the rays of the considered ambient toric variety's fan. RETURN: The ideal of the passed point's orbit closure under the torus action defined by the passed matrix. EXAMPLE: example varproductOrbitClosureIdeal; shows an example." { //Make sure a basering is defined if (!defined(basering)) { ERROR("This function needs a basering to be defined before being called."); } int icvzP; int ncvzP = ncols(rvcvzP); //Appropriate basering check. int imonVars; int nmonVars = nvars(basering); if (nmonVars != ncvzP) { ERROR("The current basering has not exactly as many variables as the passed matrix has columns."); } //The number of rays must match the dimension of the vector space the point is defined in. if (typeof(rvnumPointInOrbit) == "intvec" && nrows(rvnumPointInOrbit) != ncvzP) { ERROR("The size of the passed vector differs from the number of columns of the passed matrix."); } //The orbit of the origin always is the origin itself again. Return the appropriate ideal . if (vectorComponentZero(rvnumPointInOrbit, ncvzP) == 1:ncvzP) { ideal spResult = var(1); for (icvzP = 2; icvzP <= ncvzP; icvzP++) { spResult = spResult, var(icvzP); } return(spResult); } //Find the zero components of the passed vector. intvec vfPointCoordinateZero = vectorComponentZero(rvnumPointInOrbit, ncvzP); if (vfPointCoordinateZero != 0:ncvzP) { //If some of the vector's components are zero, take the rays assigned to these components, build a matrix from them, //transpose it and then calculate that matrix's kernel in order to eliminate the zero indexes' rays' span of im(rvcvzP*). //This yields im(rvcvzP*) \cap , where e_i_j denominate those indices of the passed vector that are not zero. //The columns of im(rvcvzP*) \cap describe the lattice basis needed for building the wanted binomial ideal. intmat cvrvzPZero = transpose(intmatTakeCols(rvcvzP, vfPointCoordinateZero)); intmat cvrvzQZero = gale(cvrvzPZero); intmat rvcvzPNonzero = cvrvzQZero * rvcvzP; } else { //Otherwise simply take the passed matrix as is for building the wanted binomial ideal. intmat rvcvzPNonzero = rvcvzP; } //Now calculate the binomial ideal that is the ideal of the toric orbit of the passed point. //First consider every non-zero entry of the passed point to be 1. ideal sbinBinomialIdeal = binomialIdeal(rvcvzPNonzero); //Now consider the actual values of the non-zero entries: map T_i |--> T_i * 1/rvnumPointInOrbit[i] for nonzero entries, T_i |--> 0 else. //Build the map. list spMapList = list(); for (icvzP = ncvzP; icvzP >= 1; icvzP--) { if (!vfPointCoordinateZero[icvzP]) { spMapList = insert(spMapList, var(icvzP) / rvnumPointInOrbit[icvzP]); } else { spMapList = insert(spMapList, 0); } } ideal spMap = spMapList[1..size(spMapList)]; //Apply the map on the binomial ideal obtained above. def R = basering; map mapCoef = R, spMap; sbinBinomialIdeal = mapCoef(sbinBinomialIdeal); //Calculate the binomial ideal's saturation in order to get the toric orbit closure's ideal from the toric orbit's ideal. ideal spResult = minimalizeIdeal(varproductIdealSaturation(sbinBinomialIdeal)); //Finally, add the variables associated with the point's zero entries to the saturated binomial ideal. for (icvzP = 1; icvzP <= ncvzP; icvzP++) { if (vfPointCoordinateZero[icvzP]) { spResult = spResult, var(icvzP); } } spResult = simplify(spResult, 2); kill spMapList, spMap; return(spResult); } example { "EXAMPLE:"; echo=2; //First example: ring R1 = 0,X(1..3),dp; intmat rvcvzP1[2][3] = 1,0,2, 0,3,1; intvec vz1 = 1,1,1; ideal spResult1 = varproductOrbitClosureIdeal(rvcvzP1, vz1); print(spResult1); //Should yield the ideal . //Second example: ring R2 = 0,X(1..6),dp; intmat rvcvzP2[2][6] = 1,0,-1,1,0,-1, 0,1,-1,1,-1,0; intvec vz2 = 2,1,1,1,2,1; ideal spResult2 = varproductOrbitClosureIdeal(rvcvzP2, vz2); print(spResult2); //Should yield the ideal <1/2*X(1)*X(4) - X(3)*X(6), X(2)*X(4) - 1/2*X(3)*X(5), 1/4*X(1)*X(5) - X(2)*X(6)>. //Third example: ring R3 = (0,a),X(1..6),dp; intmat rvcvzP3[2][6] = 1,0,-1,1,0,-1, 0,1,-1,1,-1,0; vector vnum3 = [1, 1, a, 1, 1]; ideal spResult3 = varproductOrbitClosureIdeal(rvcvzP3, vnum3); print(spResult3); //Should yield the ideal . //Fourth example: ring R4 = 0,X(1..3),dp; intmat rvcvzP4[2][3] = 1,0,-1, 0,1,-1; intvec vz4 = 0,0,0; ideal spResult4 = varproductOrbitClosureIdeal(rvcvzP4, vz4); print(spResult4); //Should yield the ideal . } proc varproductIdealSaturation(ideal spI) "USAGE: varproductIdealSaturation(spI); with spI: ideal. PURPOSE: Calculates the saturation of an ideal w.r.t. the monomial var(1)*...*var(nvars). RETURN: The ideal that is the saturation of the passed ideal w.r.t. the monomial var(1)*...*var(nvars). EXAMPLE: example varproductIdealSaturation; shows an example." { //Saturation computation via Rabinovic. def R = basering; int imonVars; int nmonVars = nvars(R); //Switch to a ring with an additional variable and elimination order. execute("ring RElim = (" + charstr(R) + "),(S,T(1..nmonVars)),(lp(1),dp);"); map mapElimFromOrig = R, T(1..nmonVars); ideal spJ = mapElimFromOrig(spI); //Add the additional polynomial for Rabinovic and eliminate the new variable. spJ = spJ, monomial(1:(nmonVars + 1)) - 1; dbprint(printlevel, ""); dbprint(printlevel, "----- Starting elimination, please be patient... -----"); if (printlevel >= 1) { int nStartTime = timer; } spJ = slimgb(spJ); ideal spISat = eliminateVariable(spJ, rvar(S)); if (printlevel >= 1) { int nTimeElapsed = timer - nStartTime; print("Time elapsed while eliminating (in ms):"); print(nTimeElapsed); } dbprint(printlevel, "----- Elimination done! -----"); //Map the saturated ideal back to the original ring. setring R; list smonOrigFromElimMap = list(); for (imonVars = nmonVars; imonVars >= 1; imonVars--) { smonOrigFromElimMap = insert(smonOrigFromElimMap, var(imonVars)); } smonOrigFromElimMap = insert(smonOrigFromElimMap, poly(0)); map mapOrigFromElim = RElim, smonOrigFromElimMap[1..size(smonOrigFromElimMap)]; ideal spResult = mapOrigFromElim(spISat); kill smonOrigFromElimMap; return(spResult); } example { "EXAMPLE:"; echo=2; //Parameterless example ring R = 0,X(1..10),dp; ideal spI = X(7)*X(10) - X(2)*X(4) + X(3)*X(5), X(8)*X(10) - X(1)*X(4) + X(3)*X(6), X(9)*X(10) - X(1)*X(5) + X(2)*X(6); ideal spISat = varproductIdealSaturation(spI); print(spISat); //Should yield the ideal generated by: // X(7)*X(10) - X(2)*X(4) + X(3)*X(5), // X(8)*X(10) - X(1)*X(4) + X(3)*X(6), // X(9)*X(10) - X(1)*X(5) + X(2)*X(6), // X(6)*X(7) - X(5)*X(8) + X(4)*X(9), // X(1)*X(7) - X(2)*X(8) + X(3)*X(9). //Example with parameter ring S = (0,a),X(1..10),dp; ideal spI = a*X(7)*X(10) - X(2)*X(4) + X(3)*X(5), X(8)*X(10) - X(1)*X(4) + X(3)*X(6), X(9)*X(10) - X(1)*X(5) + X(2)*X(6); ideal spISat = varproductIdealSaturation(spI); print(spISat); //Should yield the ideal generated by: // a*X(7)*X(10) - X(2)*X(4) + X(3)*X(5), // X(8)*X(10) - X(1)*X(4) + X(3)*X(6), // X(9)*X(10) - X(1)*X(5) + X(2)*X(6), // a*X(6)*X(7) - X(5)*X(8) + X(4)*X(9), // a*X(1)*X(7) - X(2)*X(8) + X(3)*X(9). } proc xIdealSaturation(ideal spI, list simonLexVars, intvec vfSatVars) "USAGE: xIdealSaturation(spI, simonLexVars, vfSatVars); with spI: ideal, simonLexVars: list of integers (the indices of those variables that shall be ordered lexicographically), vfSatVars: intvec (indices of the variables to consider in saturation). PURPOSE: Calculates the saturation of an ideal w.r.t. the monomial defined by vfSatVars. In contrast to working with the standard elimination order, some variables can be ordered lexicographically. RETURN: The ideal that is the saturation of the passed ideal w.r.t. the monomial defined by vfSatVars. EXAMPLE: example xIdealSaturation; shows an example." { //Saturation computation via Rabinovic. def R = basering; int imonVars; int nmonVars = nvars(R); //Check whether the passed list of variable indices contains valid variable indices solely and each index at most once. Along with this, //extract the used indices from the passed list of variable indices together with their position therein. int imonLexVars; intvec vnLexVars = 0:nmonVars; int iimonLexVars; int nimonLexVars = size(simonLexVars); for (iimonLexVars = nimonLexVars; iimonLexVars >= 1; iimonLexVars--) { imonLexVars = simonLexVars[iimonLexVars]; if (imonLexVars <= 0 || imonLexVars > nmonVars) { ERROR("The passed list of indices contains indices out of bounds."); } if (vnLexVars[imonLexVars] > 0) { ERROR("The passed list of indices contains some indices multiple times.") } vnLexVars[imonLexVars] = iimonLexVars; } poly pSat = monomial(vfSatVars); //Switch to a ring with an additional variable and elimination order extended to a lex order for the variables stated. if (nimonLexVars == nmonVars) { execute("ring RElim = (" + charstr(R) + "),(S,T(1..nmonVars)),(lp(1),lp(nimonLexVars));"); } else { if (nimonLexVars == 0) { execute("ring RElim = (" + charstr(R) + "),(S,T(1..nmonVars)),(lp(1),dp);"); } else { execute("ring RElim = (" + charstr(R) + "),(S,T(1..nmonVars)),(lp(1),lp(nimonLexVars),dp);"); } } int imonNonlex = nmonVars; list smonElimFromOrigMap = list(); for (imonVars = nmonVars; imonVars >= 1; imonVars--) { if (vnLexVars[imonVars] > 0) { smonElimFromOrigMap = insert(smonElimFromOrigMap, T(vnLexVars[imonVars])); } else { smonElimFromOrigMap = insert(smonElimFromOrigMap, T(imonNonlex)); imonNonlex--; } } map mapElimFromOrig = R, smonElimFromOrigMap[1..size(smonElimFromOrigMap)]; ideal spJ = mapElimFromOrig(spI), mapElimFromOrig(pSat) * S - 1; dbprint(printlevel, ""); dbprint(printlevel, "----- Starting elimination, please be patient... -----"); if (printlevel >= 1) { int nStartTime = timer; } spJ = slimgb(spJ); ideal spISat = eliminateVariable(spJ, rvar(S)); if (printlevel >= 1) { int nTimeElapsed = timer - nStartTime; print("Time elapsed while eliminating (in ms):"); print(nTimeElapsed); } dbprint(printlevel, "----- Elimination done! -----"); kill smonElimFromOrigMap; //Map the saturated ideal to the original ring reordering the variables setring R; list smonOrigFromElimMap = list(); imonNonlex = nmonVars; for (imonVars = nmonVars; imonVars > nimonLexVars; imonVars--) { if (imonNonlex > 0) { while (vnLexVars[imonNonlex] > 0) { imonNonlex--; if (imonNonlex == 0) { break; } } if (imonNonlex > 0) { smonOrigFromElimMap = insert(smonOrigFromElimMap, var(imonNonlex)); imonNonlex--; } } } for (iimonLexVars = nimonLexVars; iimonLexVars >= 1; iimonLexVars--) { smonOrigFromElimMap = insert(smonOrigFromElimMap, var(simonLexVars[iimonLexVars])); } smonOrigFromElimMap = insert(smonOrigFromElimMap, poly(0)); map mapOrigFromElim = RElim, smonOrigFromElimMap[1..size(smonOrigFromElimMap)]; ideal spResult = mapOrigFromElim(spISat); kill smonOrigFromElimMap; return(spResult); } example { "EXAMPLE:"; echo=2; //Parameterless example ring R = 0,X(1..10),dp; ideal spI = X(7)*X(10) - X(2)*X(4) + X(3)*X(5), X(8)*X(10) - X(1)*X(4) + X(3)*X(6), X(9)*X(10) - X(1)*X(5) + X(2)*X(6); list simonLexVars = 9,8,7; intvec vfSatVars = 0,0,0,0,0,0,0,0,0,1; ideal spISat = xIdealSaturation(spI, simonLexVars, vfSatVars); print(spISat); //Should yield the ideal generated by: // X(7)*X(10) - X(2)*X(4) + X(3)*X(5), // X(8)*X(10) - X(1)*X(4) + X(3)*X(6), // X(9)*X(10) - X(1)*X(5) + X(2)*X(6), // X(6)*X(7) - X(5)*X(8) + X(4)*X(9), // X(1)*X(7) - X(2)*X(8) + X(3)*X(9). //Example with parameter ring S = (0,a),X(1..10),dp; ideal spI = a*X(7)*X(10) - X(2)*X(4) + X(3)*X(5), X(8)*X(10) - X(1)*X(4) + X(3)*X(6), X(9)*X(10) - X(1)*X(5) + X(2)*X(6); ideal spISat = xIdealSaturation(spI, simonLexVars, vfSatVars); print(spISat); //Should yield the ideal generated by: // a*X(7)*X(10) - X(2)*X(4) + X(3)*X(5), // X(8)*X(10) - X(1)*X(4) + X(3)*X(6), // X(9)*X(10) - X(1)*X(5) + X(2)*X(6), // a*X(6)*X(7) - X(5)*X(8) + X(4)*X(9), // a*X(1)*X(7) - X(2)*X(8) + X(3)*X(9). } proc veronesePoly(int nDeg, vector vnum) "USAGE: veronesePoly(nDeg, vnum); with nDeg: int, vnum: vector of numbers. PURPOSE: Expands a vector of coefficients for the monomials of a certain degree to the sum of those monomials multiplied by their corresponding coefficient. RETURN: A poly that is the sum of all monomials of the passed degree, each of them multiplied by its corresponding passed coefficient. EXAMPLE: example veronesePoly; shows an example." { int imonVars; int nmonVars = nvars(basering); ideal smon = maxideal(nDeg); int imon; int nmon = ncols(smon); poly pResult = 0; for (imon = 1; imon <= nmon; imon++) { pResult = pResult + vnum[imon] * smon[nmon - imon + 1]; } return(pResult); } example { "EXAMPLE"; echo=2; ring R = 0,T(1..3),dp; vector vnum = veronese(3, [1, 2, 3]); print(vnum); //Should yield the vector //[1, 2, 3, 4, 6, 9, 8, 12, 18, 27] poly pResult = veronesePoly(3, vnum); print(pResult); //Should yield the polynomial //T(1)^3 + 2*T(1)^2*T(2) + 3*T(1)^2*T(3) + 4*T(1)*T(2)^2 + 6*T(1)*T(2)*T(3) + 9*T(1)*T(3)^2 + 8*T(2)^3 + 12*T(2)^2*T(3) + 18*T(2)*T(3)^2 + 27*T(3)^3 } /*********************************************** * Algorithms on CEMDSs * ***********************************************/ proc blowupCEMDS(CEMDS cemds, list spC, intvec vnD, list #) "USAGE: blowupCEMDS(cemds, spC, vnD); with cemds: CEMDS, spC: list of polynomials, vnD: intvec consisting of positive integers; blowupCEMDS(cemds, spC, vnD, fComputeFan); with cemds: CEMDS, spC: list of polynomials, vnD: intvec consisting of positive integers; fComputeFan: int. blowupCEMDS(cemds, spC, vnD, fComputeFan, vzWeight); with cemds: CEMDS, C: ideal, fVerify: int, fComputeFan: int, vzWeight: intvec. ASSUME: spC is an ideal in the Cox ring of the given CEMDS and the generators of spC are all K-prime; we assume that the subvariety defined by spC is contained in the smooth locus. PURPOSE: Blows up the given CEMDS in the subvariety defined by the ideal spC. Multiplicities of the generators of spC are specified by vnD. RETURN: The blown up CEMDS. It is automatically checked whether the result is indeed a CEMDS. The CEMDS's fan is computed iff fComputeFan != 0. NOTE: By default: fComputeFan = 0; vzWeight is an element of the relative interior of the moving cone generated by the columns of the grading matrix rvcvzQ = gale(cemds.rvcvzP) extended by the passed polynomials' degrees. SEE ALSO: blowupCEMDSpoints treats the special case of blowing up points on a CEMDS. EXAMPLE: example blowupCEMDS; shows an example." { //Make sure the CEMDS is associated with a ring def R1 = cemds.R; if (!defined(R1)) { ERROR("The passed CEMDS is not associated with any ring."); } setring R1; if (!defined(spC)) { ERROR("The passed polynomials must be elements of the given CEMDS's ring."); } if (size(vnD) != size(spC)) { ERROR("The number of polynomials must equal the number of entries of the given multiplicity vector."); } int nmonVars = nvars(R1); //Search for optional parameters list sParams = #; int nParams = size(sParams); //Default parameters int fComputeFan = 0; //intvec vzWeight: intentionally not defined here //First parameter: compute a fan or not if (nParams > 0) { //Type check for fan computation flag if (typeof(sParams[1]) != "int") { ERROR("The parameter passed for \"fComputeFan\" must be of type \"int\" and equal \"0\" or \"1\"."); } fComputeFan = sParams[1]; //Second parameter: An optional weight for initializing the fan computation. if (nParams > 1) { //Type check for weight if (typeof(sParams[2]) != "intvec") { ERROR("The parameter passed for \"vzWeight\" must be of type \"intvec\"."); } intvec vzWeight = sParams[2]; //There are at most 2 optional parameters if (nParams > 2) { ERROR("Too many optional parameters."); } } } //If there is nothing to blow up, return original CEMDS. if (size(spC) == 0) { return(cemds); } // initialization CEMDS cemdsStretched = stretchCEMDS(cemds, spC, list(), 0); // FAN: To check?! def RStretch = cemdsStretched.R; setring RStretch; dbprint(printlevel - 1, ""); dbprint(printlevel - 1, "The stretched CEMDS:"); dbprint(printlevel - 1, ""); dbprint(printlevel - 1, cemdsStretched); // prepare for blow up intvec rvnVanishingDegrees = 0:nmonVars, vnD; intvec vfFakeVars = 0:nmonVars, 1:size(vnD); CEMDS cemdsModified = modifyCEMDS2(cemdsStretched, rvnVanishingDegrees, 0, vfFakeVars, list(), 0); def RModified = cemdsModified.R; setring RModified; dbprint(printlevel - 1, ""); dbprint(printlevel - 1, "The resulting CEMDS after blowing up the ideal:"); dbprint(printlevel - 1, ""); dbprint(printlevel - 1, cemdsResult); // always verify the result: list spCFetched = fetch(R1, spC); int fVerified = checkDim(cemdsModified, nmonVars, spCFetched); //Compress the resulting CEMDS. //Compress the coordinates of the points to be blown up later along with the CEMDS. CEMDS cemdsCompressed; if (defined(vzWeight)) { cemdsCompressed = compressCEMDS(cemdsModified, 0, list(), fComputeFan, vzWeight); } else { cemdsCompressed = compressCEMDS(cemdsModified, 0, list(), fComputeFan); } def R2 = cemdsCompressed.R; setring R2; if (!fVerified) { dbprint(printlevel + 1, ""); dbprint(printlevel + 1, "CEMDS verification failed on blowing up."); dbprint(printlevel + 1, ""); } else { dbprint(printlevel + 1, ""); dbprint(printlevel + 1, "CEMDS verification successful."); dbprint(printlevel + 1, ""); } dbprint(printlevel - 1, ""); dbprint(printlevel - 1, "The resulting CEMDS after compressing:"); dbprint(printlevel - 1, ""); dbprint(printlevel - 1, cemdsCompressed); return(cemdsCompressed); } example { "EXAMPLE:"; echo=2; //First example: intmat rvcvzQ1[5][8] = 1,0,0,0,2,0,3,-1, 0,1,0,0,1,0,2,-1, 0,0,1,0,-3,0,-2,2, 0,0,0,1,-2,0,-1,1, 0,0,0,0,0,1,1,-1; intmat rvcvzP1[3][8] = gale(rvcvzQ1); fan scnSigma1 = emptyFan(2); ring R1 = 0,T(1..8),dp; ideal spG1 = T(3)^3*T(4)^2*T(5) - T(1)^2*T(2) - T(7)*T(8); CEMDS cemdsX1 = createCEMDS(rvcvzP1, scnSigma1, spG1); list spC = T(1),T(3),T(7); intvec vnD = 1,1,1; CEMDS cemdsX2 = blowupCEMDS(cemdsX1, spC, vnD); print(cemdsX2); //Verification should be unsuccessful kill rvcvzQ1, spG1, cemdsX1, spC, vnD, cemdsX2, scnSigma1, rvcvzP1, R1; // Second example: intmat rvcvzQ1[5][8] = 1,0,0,0,2,0,3,-1, 0,1,0,0,1,0,2,-1, 0,0,1,0,-3,0,-2,2, 0,0,0,1,-2,0,-1,1, 0,0,0,0,0,1,1,-1; intmat rvcvzP1[3][8] = gale(rvcvzQ1); fan scnSigma1 = emptyFan(2); ring R1 = 0,T(1..8),dp; ideal spG1 = T(3)^3*T(4)^2*T(5) - T(1)^2*T(2) - T(7)*T(8); CEMDS cemdsX1 = createCEMDS(rvcvzP1, scnSigma1, spG1); list spC = T(1),T(3),T(7); intvec vnD = 1,1,2; CEMDS cemdsX2 = blowupCEMDS(cemdsX1, spC, vnD); print(cemdsX2); //Verification should be successful } // checks whether dim(I+) > dim(I+ // where T^nu is the product over all T_i with T_i not in C // and the first r1 variables in cemds belong to the ring // before the blow up. static proc checkDim(CEMDS cemds, int r1, list spC) { int b = 0; def RCurrent = cemds.R; ideal C = spC[1..size(spC)]; //setring RCurrent; if (nvars(RCurrent) <= r1) { ERROR("r1 counts the number of variables before the blow up, i.e. must be smaller than the number of variables of the current ring."); } // create T^nu i.e. the product over all T_i with T_i not in C: poly Tnu = 1; ideal I = cemds.spG; ideal G = std(C); int i; for (i = 1; i <= r1; i++) { if (reduce(var(i), G) != 0) { Tnu = Tnu * var(i); } } // Create the ideal J := I+ // where T_r is the exceptional divisor ideal J = I, var(nvars(RCurrent)); ideal Jnu = J, Tnu; if (dim(std(J)) > dim(std(Jnu))) { b = 1; } return(b); } example { "EXAMPLE"; echo=2; // 1st example: // should be successful intmat Q1[5][8] = 1,0,0,0,2,0,3,-1, 0,1,0,0,1,0,2,-1, 0,0,1,0,-3,0,-2,2, 0,0,0,1,-2,0,-1,1, 0,0,0,0,0,1,1,-1; intmat P1[3][8] = gale(Q1); fan Sigma1 = emptyFan(2); ring R1 = 0,T(1..8),dp; ideal G1 = T(3)^3*T(4)^2*T(5) - T(1)^2*T(2) - T(7)*T(8); CEMDS X1 = createCEMDS(P1, Sigma1, G1); ideal C = T(1),T(3),T(7); int b = checkDim(X1,7,C); print(b); // should be successful // 2nd example: // should be unsuccessful G1 = T(3)^3*T(4)^2*T(5) - T(8)^2*T(1)^2*T(2) - T(7)*T(8); X1 = createCEMDS(P1, Sigma1, G1); C = T(1),T(3),T(7); b = checkDim(X1,7,C); print(b); // should be unsuccessful } proc blowupCEMDSpoints(CEMDS cemds, list #) "USAGE: blowupCEMDSpoints(cemds, svnumPts); with cemds: CEMDS, svnumPts: list of vectors of numbers/a vector of numbers; blowupCEMDSpoints(cemds, svnumPts, fVerify); with cemds: CEMDS, svnumPts: list of vectors of numbers, fVerify: int; blowupCEMDSpoints(cemds, svnumPts, fVerify, fComputeFan); with cemds: CEMDS, svnumPts: list of vectors of numbers, fVerify: int, fComputeFan: int. blowupCEMDSpoints(cemds, svnumPts, fVerify, fComputeFan, vzWeight); with cemds: CEMDS, svnumPts: list of vectors of numbers, fVerify: int, fComputeFan: int, vzWeight: intvec. ASSUME: All of the passed points to blow up are contained in the passed CEMDS. PURPOSE: Blows up the ambient toric variety of a CEMDS in some points contained therein. RETURN: The modified CEMDS in the blown up toric ambient variety. Whether the result is indeed a CEMDS is checked iff fVerify != 0. The CEMDS's fan is computed iff fComputeFan != 0. NOTE: By default: svnumPts = list(); fVerify = 0; fComputeFan = 0; vzWeight is an element of the relative interior of the moving cone generated by the columns of the grading matrix rvcvzQ = gale(cemds.rvcvzP) extended by the passed polynomials' degrees. EXAMPLE: example blowupCEMDSpoints; shows an example." { //Make sure the CEMDS is associated with a ring def RInitial = basering; def RCurrent = cemds.R; if (!defined(RCurrent)) { ERROR("The passed CEMDS is not associated with any ring."); } int imonVars; int nmonVars = nvars(RCurrent); //Search for optional parameters list sParams = #; int nParams = size(sParams); //Default parameters //list svnumPts: intentionally not defined here; int fVerify = 0; int fComputeFan = 0; //intvec vzWeight: intentionally not defined here //First parameter: A list of points where the passed CEMDS shall be blown up. if (nParams > 0) { //Type and size check for the list of points to blow up along with the CEMDS. //Either a list of vectors was passed or a list containing a list of vectors as its first entry. if (typeof(sParams[1]) == "vector") { list svnumPts = sParams; } else { if (typeof(sParams[1]) != "list") { ERROR("The parameter passed for \"svnumPts\" must be of type \"list\" and contain vectors only."); } list svnumPts = sParams[1]; } int ivnumPts; int nvnumPts = size(svnumPts); for (ivnumPts = 1; ivnumPts <= nvnumPts; ivnumPts++) { if (typeof(svnumPts[ivnumPts]) != "vector") { ERROR("The passed points to blow up along with the CEMDS must be of type \"vector\" and consist of numbers only."); } for (imonVars = 1; imonVars <= nmonVars; imonVars++) { if (leadexp(svnumPts[ivnumPts][imonVars]) != 0:nmonVars) { ERROR("The vectors passed must consist of numbers only."); } } } //Second parameter: Verify the result or not. if (typeof(sParams[1]) == "list" && nParams > 1) { //Type check for verification flag if (typeof(sParams[2]) != "int") { ERROR("The parameter passed for \"fVerify\" must be of type \"int\" and equal \"0\" or \"1\"."); } fVerify = sParams[2]; //Third parameter: Compute a fan or not. if (nParams > 2) { //Type check for fan computation flag if (typeof(sParams[3]) != "int") { ERROR("The parameter passed for \"fComputeFan\" must be of type \"int\" and equal \"0\" or \"1\"."); } fComputeFan = sParams[3]; //Fourth parameter: An optional weight for initializing the fan computation. if (nParams > 3) { //Type check for weight if (typeof(sParams[4]) != "intvec") { ERROR("The parameter passed for \"vzWeight\" must be of type \"intvec\"."); } intvec vzWeight = sParams[4]; //There are at most 4 optional parameters if (nParams > 4) { ERROR("Too many optional parameters."); } } } } } //If there is no point to blow up, return original CEMDS. if (nvnumPts == 0) { return(cemds); } //Iteration initialization intvec vfFakeVars; CEMDS cemdsResult = cemds; setring RCurrent; int nnumPts; if (!defined(svnumPts)) { list svnumPts = fetch(RInitial, svnumPts); } intvec rvnVanishingDegrees; int fVerified = 1; //Blow up the passed points one by one. for (ivnumPts = 1; ivnumPts <= nvnumPts; ivnumPts++) { //Select the next point to blow up vector vnumPtBlowup = svnumPts[1]; svnumPts = delete(svnumPts, 1); nnumPts = size(svnumPts); dbprint(printlevel - 1, ""); dbprint(printlevel - 1, "The next point to blow up:"); dbprint(printlevel - 1, vnumPtBlowup); //First calculate the point's lattice ideal. ideal spI = varproductOrbitClosureIdeal(cemdsResult.rvcvzP, vnumPtBlowup); list spIList = spI[1..ncols(spI)]; spIList = removeMonomials(spIList); dbprint(printlevel - 1, ""); dbprint(printlevel - 1, "Ideal of the orbit closure of the point to blow up:"); dbprint(printlevel - 1, spI); //Compute the multiplicity vector needed later for blowing up the CEMDS via stellar subdivision in its fan. rvnVanishingDegrees = vanishingDegrees(spIList, vnumPtBlowup); //Then stretch the embedding using this ideal's generators so the orbit of the point considered can be obtained by requiring some variables to be equal to zero. //Stretch the coordinates of the points to be blown up later along with the CEMDS. if (nnumPts > 0) { RCurrent = stretchCEMDS(cemdsResult, spIList, svnumPts, 0); kill vnumPtBlowup; kill svnumPts; kill spI; kill spIList; setring RCurrent; list sStretchResult = gsResult; cemdsResult = sStretchResult[1]; list svnumPts = sStretchResult[2]; } else { cemdsResult = stretchCEMDS(cemdsResult, spIList, list(), 0); kill vnumPtBlowup; kill spI; kill spIList; RCurrent = cemdsResult.R; setring RCurrent; } dbprint(printlevel - 1, ""); dbprint(printlevel - 1, "The stretched CEMDS:"); dbprint(printlevel - 1, ""); dbprint(printlevel - 1, cemdsResult); //Mark the new variables for fake relations vfFakeVars = 0:nmonVars; for (imonVars = nmonVars + 1; imonVars <= nvars(RCurrent); imonVars++) { vfFakeVars = vfFakeVars, 1; } //Blow up the previously selected point, thereby modifying the CEMDS. //Modify the coordinates of the points to be blown up later along with the CEMDS. if (nnumPts > 0) { RCurrent = modifyCEMDS2(cemdsResult, rvnVanishingDegrees, 0, vfFakeVars, svnumPts, 0); kill gsResult; kill sStretchResult; kill svnumPts; setring RCurrent; list sModifyResult = gsResult; cemdsResult = sModifyResult[1]; list svnumPts = sModifyResult[2]; } else { cemdsResult = modifyCEMDS2(cemdsResult, rvnVanishingDegrees, 0, vfFakeVars, list(), 0); RCurrent = cemdsResult.R; setring RCurrent; } dbprint(printlevel - 1, ""); dbprint(printlevel - 1, "The resulting CEMDS after blowing up in the current point:"); dbprint(printlevel - 1, ""); dbprint(printlevel - 1, cemdsResult); //Compress the resulting CEMDS. //Compress the coordinates of the points to be blown up later along with the CEMDS. if (nnumPts > 0) { RCurrent = compressCEMDS(cemdsResult, fVerify, svnumPts, 0); kill gsResult; kill sModifyResult; kill svnumPts; setring RCurrent; list sCompressResult = gsResult; cemdsResult = sCompressResult[1]; if (fVerify) { fVerified = sCompressResult[2]; } kill gsResult; } else { if (defined(vzWeight)) { if (fVerify) { list sCompressResult = compressCEMDS(cemdsResult, fVerify, list(), fComputeFan, vzWeight); cemdsResult = sCompressResult[1]; fVerified = sCompressResult[2]; } else { cemdsResult = compressCEMDS(cemdsResult, fVerify, list(), fComputeFan, vzWeight); } } else { if (fVerify) { list sCompressResult = compressCEMDS(cemdsResult, fVerify, list(), fComputeFan); cemdsResult = sCompressResult[1]; fVerified = sCompressResult[2]; } else { cemdsResult = compressCEMDS(cemdsResult, fVerify, list(), fComputeFan); } } RCurrent = cemdsResult.R; setring RCurrent; } if (fVerify && !fVerified) { dbprint(printlevel, ""); dbprint(printlevel, "CEMDS verification failed on blowing up the " + ordStrFromInt(ivnumPts) + " point."); } dbprint(printlevel - 1, ""); dbprint(printlevel - 1, "The resulting CEMDS after compressing:"); dbprint(printlevel - 1, ""); dbprint(printlevel - 1, cemdsResult); //Extract the points to blow up next or print the verification result in the last iteration. if (nnumPts > 0) { if (fVerify) { list svnumPts = sCompressResult[3]; } else { list svnumPts = sCompressResult[2]; } nmonVars = nvars(RCurrent); kill sCompressResult; } else { dbprint(printlevel + 1, ""); if (!fVerify) { dbprint(printlevel + 1, "NOTE: You have chosen to not perform a verification step; the result may not be a CEMDS."); } else { if (fVerified) { dbprint(printlevel + 1, "The resulting embedded space was successfully verified to be a CEMDS."); } else { dbprint(printlevel + 1, "The resulting embedded space could not be verified to be a CEMDS."); } } } } return(cemdsResult); } example { "EXAMPLE:"; echo=2; //First example: No fan computation. ring R = (0,a),T(1..3),dp; //Define the input data intmat rvcvzP[2][3] = 1,0,-1, 0,1,-1; ideal spG = 0; vector vnum1 = [1, 0, 0]; vector vnum2 = [0, 1, 0]; list svnum = vnum1, vnum2; intmat cvrvz1[2][2] = 1,0, 0,1; intmat cvrvz2[2][2] = 0,1, -1,-1; intmat cvrvz3[2][2] = -1,-1, 1,0; cone cn1 = coneViaPoints(cvrvz1); cone cn2 = coneViaPoints(cvrvz2); cone cn3 = coneViaPoints(cvrvz3); fan scnSigma = fanViaCones(cn1, cn2, cn3); //Compose the CEMDS from the data provided until now. CEMDS cemds = createCEMDS(rvcvzP, scnSigma, spG); //Perform the blowup CEMDS cemdsBlowup = blowupCEMDSpoints(cemds, svnum); //Obtains the CEMDS's ring def RBlowup = cemdsBlowup.R; setring RBlowup; print(cemdsBlowup); //Should yield // - A matrix with the following rays as columns // (1 0 -1 -1 0) // (0 1 -1 0 -1) // or there must exist a linear isomorphism mapping the output rays to those above. // - An empty fan. // - The zero ideal. kill R, rvcvzP, cvrvz1, cvrvz2, cvrvz3, cn1, cn2, cn3, scnSigma, cemds; //Second example: With fan computation and verification ring R = 0,T(1..4),dp; //Define the input data intmat rvcvzP[2][4] = 1,-1,-1,1, 1,1,-1,-1; ideal spG = 0; vector vnum1 = [1, 1, 1, 1]; list svnum = vnum1; intmat cvrvz1[2][2] = 1,1, -1,1; intmat cvrvz2[2][2] = -1,-1, -1,1; intmat cvrvz3[2][2] = -1,-1, 1,-1; intmat cvrvz4[2][2] = 1,1, 1,-1; cone cn1 = coneViaPoints(cvrvz1); cone cn2 = coneViaPoints(cvrvz2); cone cn3 = coneViaPoints(cvrvz3); cone cn4 = coneViaPoints(cvrvz4); fan scnSigma = fanViaCones(cn1, cn2, cn3, cn4); //Compose the CEMDS from the data provided until now. CEMDS cemds = createCEMDS(rvcvzP, scnSigma, spG); //Perform the blowup CEMDS cemdsBlowupWithFan = blowupCEMDSpoints(cemds, svnum, 1, 1); //Obtains the CEMDS's ring def RBlowupWithFan = cemdsBlowupWithFan.R; setring RBlowupWithFan; printCEMDS(cemdsBlowupWithFan); //Should yield // - A matrix with the following rays as columns // ( 0 1 0 -1 0 0 0 0 0) // (-1 0 1 0 0 0 0 0 0) // (-1 1 0 0 1 -1 0 0 0) // ( 0 0 0 0 -1 0 1 0 0) // (-1 1 0 0 -1 0 0 1 0) // (-1 -1 0 0 1 0 0 0 1) // or there must exist a linear isomorphism mapping the output rays to those above. // - Some fan. // - The ideal spanned by // T(5)^2-T(7)^2+T(6)*T(8), // T(4)*T(5)-T(1)*T(6)+T(2)*T(7), // T(3)*T(5)-T(1)*T(7)+T(2)*T(8), // T(2)*T(5)-T(3)*T(6)+T(4)*T(7), // T(1)*T(5)-T(3)*T(7)+T(4)*T(8), // T(2)*T(3)-T(1)*T(4)-T(5)*T(9), // T(2)^2-T(4)^2-T(6)*T(9), // T(1)*T(2)-T(3)*T(4)-T(7)*T(9), // T(1)^2-T(3)^2-T(8)*T(9). } proc compressCEMDS(CEMDS cemds, list #) "USAGE: compressCEMDS(cemds); with cemds: CEMDS; compressCEMDS(cemds, fVerify); with cemds: CEMDS, fVerify: int; compressCEMDS(cemds, fVerify, svnumPts); with cemds: CEMDS, fVerify: int, svnumPts: list of vectors of numbers; compressCEMDS(cemds, fVerify, svnumPts, fComputeFan); with cemds: CEMDS, fVerify: int, svnumPts: list of vectors of numbers, fComputeFan: int; compressCEMDS(cemds, fVerify, svnumPts, fComputeFan, vzWeight); with cemds: CEMDS, fVerify: int, svnumPts: list of vectors of numbers, fComputeFan: int, vzWeight: intvec. ASSUME: (i) All points to compress along with the CEMDS, if there were any passed, are contained in the CEMDS to compress; (ii) The intvec \"vzWeight\", if passed, is element of the ample cone defined by the grading matrix rvcvzQ extended by the passed polynomials' degrees and some true saturated connected collection of facets of the positive orthant. PURPOSE: Embeds a CEMDS in a new ambient toric variety such that abundant relations in the ideal defining the CEMDS (i.e. of the type var(i) - p with p not containing var(i)) are erased. Optionally checks whether the result really is a CEMDS and transfers points from the uncompressed CEMDS to the compressed one. RETURN: If no optional output was requested, the compressed CEMDS only. If only verification was requested, but no points or an empty list of points to compress along with the CEMDS were passed, a list containing the following (in order of appearance): (i) the compressed CEMDS with or without a fan (depending on fComputeFan and vzWeight) (always present), (ii) the verification result (present as fVerify was passed and fVerify != 0). If a non-empty list of points to compress along with the CEMDS was passed, the compressed CEMDS's basering. SIDE EFFECTS: Iff a non-empty list of points to compress along with the CEMDS was passed, exports a list gsResult containing the following (if present and in order of appearance): (i) the compressed CEMDS with or without a fan (depending on fComputeFan and vzWeight) (always present), (ii) the verification result (present iff fVerify != 0), (iii) the list of compressed points (present as there was a list passed). NOTE: By default: fVerify = 0; svnumPts = list(); fComputeFan = 0; vzWeight is an element of the relative interior of the moving cone generated by the columns of the grading matrix rvcvzQ = gale(cemds.rvcvzP) extended by the passed polynomials' degrees. EXAMPLE: example compressCEMDS; shows an example." { //The CEMDS to compress must be associated with a ring def R = cemds.R; if (!defined(R)) { ERROR("The passed CEMDS is not associated with any ring."); } //Switch to the CEMDS's ring setring R; int imonVars; int nmonVars = nvars(R); //Search for optional parameters list sParams = #; int nParams = size(sParams); //Default parameters int fVerify = 0; //list svnumPts: intentionally not defined here; int fComputeFan = 0; //intvec vzWeight: intentionally not defined here //First parameter: Verify whether the resulting ES is a CEMDS or not. if (nParams > 0) { //Type check for verification flag if (typeof(sParams[1]) != "int") { ERROR("The parameter passed for \"fVerify\" must be of type \"int\" and equal \"0\" or \"1\"."); } fVerify = sParams[1]; //Second parameter: A list of points whose coordinates shall be transferred to the compressed CEMDS. if (nParams > 1) { //Type and size check for the list of points to compress along with the CEMDS if (typeof(sParams[2]) != "list") { ERROR("The parameter passed for \"svnumPts\" must be of type \"list\" and contain vectors only."); } list svnumPts = sParams[2]; int ivnumPts; int nvnumPts = size(svnumPts); for (ivnumPts = 1; ivnumPts <= nvnumPts; ivnumPts++) { if (typeof(svnumPts[ivnumPts]) != "vector") { ERROR("The passed points to compress along with the CEMDS must be of type \"vector\" and consist of numbers only."); } for (imonVars = 1; imonVars <= nmonVars; imonVars++) { if (leadexp(svnumPts[ivnumPts][imonVars]) != 0:nmonVars) { ERROR("The vectors passed must consist of numbers only."); } } } //Third parameter: Compute a fan or not. if (nParams > 2) { //Type check for fan computation flag if (typeof(sParams[3]) != "int") { ERROR("The parameter passed for \"fComputeFan\" must be of type \"int\" and equal \"0\" or \"1\"."); } fComputeFan = sParams[3]; //Fourth parameter: An optional weight for initializing the fan computation. if (nParams > 3) { //Type check for weight if (typeof(sParams[4]) != "intvec") { ERROR("The parameter passed for \"vzWeight\" must be of type \"intvec\"."); } intvec vzWeight = sParams[4]; //There are at most 4 optional parameters if (nParams > 4) { ERROR("Too many optional parameters."); } } } } } //Initialize the loop compressing the CEMDS's embedding ideal int ipG; int npG; ideal spG = cemds.spG; poly pGi; int imonGi; int nmonGi; poly monGi; intvec rvfRemainingVars = 1:nmonVars; int fFoundFake = 1; //As long as fake polynomials are found in the embedding ideal, the single variables therein are eliminated from the other generators of the ideal. while (fFoundFake) { fFoundFake = 0; //Erase zero entries (if spG != 0) and redundant generators. spG = simplify(spG, 10); //Search the ideal's polynomials for fake polynomials, i.e., //polynomials containing a term being a ring's variable not contained in the same polynomial's other terms. ipG = 1; npG = ncols(spG); //Iterate through the ideal's polynomials while (ipG <= npG && !fFoundFake) { pGi = spG[ipG]; imonGi = 1; nmonGi = size(pGi); //Search the current polynomial for a term being a ring's variable solely. while (imonGi <= nmonGi && !fFoundFake) { if (isOfTotalDegree(1, pGi[imonGi])) { monGi = pGi[imonGi]; //Check whether the found ring's varible appears in other terms of the polynomial, too. If it does not, the current polynomial is fake. if (intersectContainedVariables(monGi, (monGi - pGi)) == 0:nmonVars) { fFoundFake = 1; } } imonGi++; } ipG++; } //If a fake polynomial was found, compress the ideal. if (fFoundFake) { //Replace the fake variable in all other polynomials of the passed ideal for (ipG = 1; ipG <= npG; ipG++) { spG[ipG] = subst(spG[ipG], leadmonom(monGi), (monGi - pGi) / leadcoef(monGi) ); } rvfRemainingVars[rvar(leadmonom(monGi))] = 0; } } //Transfer the compressed result to a new ring with as few variables as possible. //Generate ring int nmonNewVars = intvecSumFlags(rvfRemainingVars); execute("ring RCompressed = (" + charstr(R) + "),T(1..nmonNewVars),dp;"); //Build map for transferring list spMapList = list(); int imonNewVars = nmonNewVars; for (imonVars = nmonVars; imonVars >= 1; imonVars--) { if (rvfRemainingVars[imonVars]) { spMapList = insert(spMapList, rvfRemainingVars[imonVars] * var(imonNewVars)); imonNewVars--; } else { spMapList = insert(spMapList, 0); } } ideal spMap = spMapList[1..nmonVars]; map mapCompress = R, spMap; //Transfer the resulting ideal and compute a minimal representation for it ideal spGCompressed = mapCompress(spG); //Extract the matrix of remaining rays intmat rvcvzQ = gale(cemds.rvcvzP); intmat rvcvzQDash = intmatTakeCols(rvcvzQ, rvfRemainingVars); intmat rvcvzPDash = gale(rvcvzQDash); //Calculate the new CEMDS's fan, if wanted fan scnSigmaDash = emptyFan(nrows(rvcvzPDash)); if (fComputeFan) { //If no weight was passed, find one in the new grading matrix's moving cone's relative interior. if (!defined(vzWeight)) { cone cnMov = movingConeFromWeights(rvcvzQDash); intvec vzWeight = intvec(relativeInteriorPoint(cnMov)); } fan sof = orthantFanFromWeight(rvcvzQDash, vzWeight); scnSigmaDash = mapFan(rvcvzPDash, sof); } kill spMapList, spMap; //Compose and return result list gsResult = list(); //The compressed points setring R; if (defined(svnumPts)) { if (nvnumPts > 0) { setring RCompressed; list svnumFetchedPts = fetch(R, svnumPts); list svnumPtsCompressed = list(); for (ivnumPts = nvnumPts; ivnumPts >= 1; ivnumPts--) { svnumPtsCompressed = insert(svnumPtsCompressed, vectorTake(svnumFetchedPts[ivnumPts], rvfRemainingVars)); } gsResult = insert(gsResult, svnumPtsCompressed); dbprint(printlevel - 1, ""); dbprint(printlevel - 1, "Compressed coordinates of the points passed:"); dbprint(printlevel - 1, svnumPtsCompressed); } } setring RCompressed; //Whether the verification was positive if (fVerify) { dbprint(printlevel, ""); dbprint(printlevel, "----- Starting verification, please be patient... -----"); if (printlevel >= 1) { int nStartTime = timer; } int fVerified = 1; //Test for being almost free fVerified = verifyAlmostFree(rvcvzPDash); if (fVerified) { ideal spGCompressedStd = std(spGCompressed); //Dimension test fVerified = verifyDimension(spGCompressedStd); if (fVerified) { //Primality test for the variables' classes in the CEMDS's ring fVerified = verifyVarPrimality(spGCompressedStd); if (!fVerified) { dbprint(printlevel - 1, "Primality test on the ring's variables failed."); } } else { dbprint(printlevel - 1, "Dimension test failed."); } } else { dbprint(printlevel - 1, "Test for being almost free failed."); } if (printlevel >= 1) { int nTimeElapsed = timer - nStartTime; print("Time elapsed while verifying (in ms):"); print(nTimeElapsed); } dbprint(printlevel, "----- Verification done! -----"); gsResult = insert(gsResult, fVerified); } //The compressed CEMDS CEMDS cemdsResult = createCEMDS(rvcvzPDash, scnSigmaDash, minimalizeIdeal(spGCompressed)); if (size(gsResult) > 0) { list scemdsResult = cemdsResult; gsResult = scemdsResult + gsResult; if (defined(svnumFetchedPts)) { export(gsResult); return(RCompressed); } else { return(gsResult); } } else { return(cemdsResult); } } example { "EXAMPLE:"; echo=2; //First example: No optional arguments ring R = (0,a),T(1..5),dp; //Blowup of P2 (stretched in order to blow up [1,1,1]) intmat rvcvzP[4][5] = 1,0,-1,0,0, 0,1,-1,0,0, 0,0,-1,1,0, 0,0,-1,0,1; intmat cvrvz1[4][4] = -1,-1,-1,-1, 0,0,0,1, 0,0,1,0, 0,1,0,0; intmat cvrvz2[4][4] = -1,-1,-1,-1, 0,0,0,1, 0,0,1,0, 1,0,0,0; intmat cvrvz3[4][4] = -1,-1,-1,-1, 0,0,0,1, 0,1,0,0, 1,0,0,0; intmat cvrvz4[4][4] = -1,-1,-1,-1, 0,0,1,0, 0,1,0,0, 1,0,0,0; intmat cvrvz5[4][4] = 0,0,0,1, 0,0,1,0, 0,1,0,0, 1,0,0,0; cone cn1 = coneViaPoints(cvrvz1); cone cn2 = coneViaPoints(cvrvz2); cone cn3 = coneViaPoints(cvrvz3); cone cn4 = coneViaPoints(cvrvz4); cone cn5 = coneViaPoints(cvrvz5); fan scnSigma = fanViaCones(cn1, cn2, cn3, cn4, cn5); poly p1 = T(1)-T(2)+T(4); poly p2 = T(2)-T(3)+T(5); ideal spG = p1, p2; CEMDS cemds = createCEMDS(rvcvzP, scnSigma, spG); CEMDS cemdsResult = compressCEMDS(cemds); //Obtains the CEMDS's ring def RCompressed = cemdsResult.R; setring RCompressed; print(cemdsResult); //Should yield the CEMDS consisting of // - A rays' matrix with the columns // (1 0 -1) // (0 1 -1) // or a matrix with its rows linearly dependent on that // - The empty fan of ambient dimension 2 // - The zero ideal //Second example: Verify that the result indeed is a CEMDS. setring R; list sCompressResult2 = compressCEMDS(cemds, 1); print(sCompressResult2); //Should yield a list with two entries: // - The same CEMDS as in the first example // - The verification result (1) //Third example: Points are transferred to the compressed CEMDS. setring R; vector vnumPt1 = [1, 2, 3, 1, 1]; vector vnumPt2 = [1, 2, a, 1, 1]; list svnumPts = vnumPt1, vnumPt2; def RCompressed3 = compressCEMDS(cemds, 0, svnumPts); setring RCompressed3; list sCompressResult3 = gsResult; print(sCompressResult3); //Should yield a list with two entries: // - The same CEMDS as in the first example // - A list containing the transferred points ([3, 1, 1] and [a, 1, 1], if the first variables are compressed first). //Fourth example: The CEMDS's fan shall also be computed. setring R; CEMDS cemdsResult4 = compressCEMDS(cemds, 0, list(), 1); print(cemdsResult4); //Should yield the CEMDS consisting of // - The rays' matrix // (1 0 -1) // (0 1 -1) // or a matrix with its rows linearly dependent on that // - The fan consisting of the maximal cones generated by the following rays: // * (1,0), (0,1); // * (0,1), (-1,-1); // * (-1,-1), (1,0); // - The zero ideal //Fifth example: Verify that the result indeed is a CEMDS, transfer a point to the compressed CEMDS and compute the CEMDS's fan with a user defined start vector. setring R; intvec vzWeight = 1; def RCompressed5 = compressCEMDS(cemds, 1, svnumPts, 1, vzWeight); setring RCompressed5; list sCompressResult5 = gsResult; print(sCompressResult5); //Should yield a list with three entries: // - The same CEMDS as in the fourth example // - The verification result (1) // - A list containing the transferred points ([3, 1, 1] and [a, 1, 1], if the first variables are compressed first). } proc contractCEMDS(CEMDS cemds, intvec vfKeepRay, list #) "USAGE: contractCEMDS(cemds, vfKeepRay); with cemds: CEMDS, vfKeepRay: intvec; contractCEMDS(cemds, vfKeepRay, fComputeFan); with cemds: CEMDS, vfKeepRay: intvec, fComputeFan: int; ASSUME: The set of rays to contract is contractible. PURPOSE: Blows down the ambient toric variety of a CEMDS; all rays not to be kept are deleted. RETURN: The original CEMDS with the specified rays contracted including the contracted fan, if not specified otherwise. NOTE: By default: fComputeFan = 1. EXAMPLE: example contractCEMDS; shows an example." { //The CEMDS to contract must be associated with a ring. def R = cemds.R; if (!defined(R)) { ERROR("The passed CEMDS is not associated with any ring."); } //The intvec describing which rays to keep must have as many entries as the CEMDS has rays in its fan. int ifKeepRay; int nfKeepRay = size(vfKeepRay); if (nfKeepRay != ncols(cemds.rvcvzP)) { ERROR("The passed vector and the CEMDS's column count differ."); } //Search for optional arguments list sParams = #; int nParams = size(sParams); //Default arguments int fComputeFan = 1; //First argument: Whether to compute a fan or not. if (nParams > 0) { //Type check for fan computation flag if (typeof(sParams[1]) != "int") { ERROR("The parameter passed for \"fComputeFan\" must be of type \"int\" and equal \"0\" or \"1\"."); } fComputeFan = sParams[1]; //There is at most one optional argument. if (nParams > 1) { ERROR("Too many optional arguments."); } } //Switch to the CEMDS's ring setring R; //Contraction of the rays' matrix intmat rvcvzPResult = intmatTakeCols(cemds.rvcvzP, vfKeepRay); int ncvzPResult = ncols(rvcvzPResult); //Contraction of the CEMDS's fan and the embedding ideal fan scnResult = cemds.scnSigma; ideal spG = cemds.spG; //Keep only the ring variables needed int nmonNewVars = intvecSumFlags(vfKeepRay); execute("ring RContracted = (" + charstr(R) + "),T(1..nmonNewVars),dp;"); int imonNewVars = nmonNewVars; list smonMappingList = list(); for (ifKeepRay = nfKeepRay; ifKeepRay >= 1; ifKeepRay--) { if (!vfKeepRay[ifKeepRay]) { if (fComputeFan) { //Contract one ray in the fan scnResult = contractRay(scnResult, intmatTakeCol(cemds.rvcvzP, ifKeepRay)); } //Replace the corresponding variable by 1 smonMappingList = insert(smonMappingList, 1); } else { smonMappingList = insert(smonMappingList, var(imonNewVars)); imonNewVars--; } } ideal smonMapping = smonMappingList[1..nfKeepRay]; map mapIdeal = R, smonMapping; ideal spGResult = mapIdeal(spG); //Compose new CEMDS, compress it and return it CEMDS cemdsContracted = createCEMDS(rvcvzPResult, scnResult, spGResult); CEMDS cemdsResult = compressCEMDS(cemdsContracted, 0, list(), fComputeFan); kill smonMappingList, smonMapping; return(cemdsResult); } example { "EXAMPLE:"; echo=2; //First example: Standard mode ring R = (0,a),T(1..6),dp; //Blowup of P2 in [1,1,1] (already stretched) intmat rvcvzP[4][6] = 1,0,-1,0,0,0, 0,1,-1,0,0,0, 0,0,-1,1,0,1, 0,0,-1,0,1,1; intmat cvrvz1[4][4] = -1,-1,-1,-1, 0,0,0,1, 0,0,1,1, 0,1,0,0; intmat cvrvz2[4][4] = -1,-1,-1,-1, 0,0,1,0, 0,0,1,1, 0,1,0,0; intmat cvrvz3[4][4] = -1,-1,-1,-1, 0,0,0,1, 0,0,1,1, 1,0,0,0; intmat cvrvz4[4][4] = -1,-1,-1,-1, 0,0,1,0, 0,0,1,1, 1,0,0,0; intmat cvrvz5[4][4] = -1,-1,-1,-1, 0,0,0,1, 0,1,0,0, 1,0,0,0; intmat cvrvz6[4][4] = -1,-1,-1,-1, 0,0,1,0, 0,1,0,0, 1,0,0,0; intmat cvrvz7[4][4] = 0,0,0,1, 0,0,1,1, 0,1,0,0, 1,0,0,0; intmat cvrvz8[4][4] = 0,0,1,0, 0,0,1,1, 0,1,0,0, 1,0,0,0; cone cn1 = coneViaPoints(cvrvz1); cone cn2 = coneViaPoints(cvrvz2); cone cn3 = coneViaPoints(cvrvz3); cone cn4 = coneViaPoints(cvrvz4); cone cn5 = coneViaPoints(cvrvz5); cone cn6 = coneViaPoints(cvrvz6); cone cn7 = coneViaPoints(cvrvz7); cone cn8 = coneViaPoints(cvrvz8); fan scnSigma = fanViaCones(cn1, cn2, cn3, cn4, cn5, cn6, cn7, cn8); poly p1 = a*T(1)-T(2)+T(4)*T(6); poly p2 = T(2)-T(3)+T(5)*T(6); ideal spG = p1, p2; intvec vfKeepRay = 1,1,1,1,1,0; CEMDS cemds = createCEMDS(rvcvzP, scnSigma, spG); CEMDS cemdsContracted = contractCEMDS(cemds, vfKeepRay); //Obtains the CEMDS's ring def RContracted = cemdsContracted.R; setring RContracted; print(cemdsContracted); //Should yield // - A matrix with the following rays as columns // (1 0 -1) // (0 1 -1) //or there must exist a linear isomorphism mapping the output rays to those above. // - The fan having the three maximal cones generated by the rays // [(1,0), (0,1)], // [(0,1), (-1,-1)] and // [(-1,-1), (1,0)] // - The zero ideal //Second example: Do not compute a fan. setring R; CEMDS cemdsContracted2 = contractCEMDS(cemds, vfKeepRay, 0); //Obtains the CEMDS's ring def RContracted2 = cemdsContracted2.R; setring RContracted2; print(cemdsContracted2); //Should yield // - A matrix with the following rays as columns // (1 0 -1) // (0 1 -1) //or there must exist a linear isomorphism mapping the output rays to those above. // - The empty fan of dimension two. // - The zero ideal } proc linearBlowup(list svnumPts, list #) "USAGE: linearBlowup(svnumPts); with svnumPts: list of vectors of numbers; linearBlowup(svnumPts, fVerify); with svnumPts: list of vectors of numbers, fVerify: int; linearBlowup(svnumPts, fVerify, fComputeFan); with svnumPts: list of vectors of numbers, fVerify: int, fComputeFan: int; linearBlowup(svnumPts, fVerify, fComputeFan, vzWeight); with svnumPts: list of vectors of numbers, fVerify: int, fComputeFan: int, vzWeight: intvec. ASSUME: The current basering contains at least n >= 2 variables. There are at least n points to blow up containing a subset of n points in general position. PURPOSE: Blows up the projective space of dimension (n-1) in some of its points using linear relations for the stretching part only. RETURN: If starting with a basering with n variables, the projective space of dimension (n-1) blown up in the points specified. Whether the result is indeed a CEMDS is checked iff fVerify != 0. The CEMDS's fan is computed iff fComputeFan != 0. NOTE: By default: fVerify = 0; fComputeFan = 0; vzWeight is an element of the relative interior of the moving cone generated by the columns of the grading matrix rvcvzQ = gale(cemds.rvcvzP) extended by the passed polynomials' degrees. EXAMPLE: example linearBlowup; shows an example." { if (!defined(basering)) { ERROR("This function requires a basering to be defined before being called."); } int imonVars; int nmonVars = nvars(basering); if (nmonVars < 2) { ERROR("This function requires a basering with at least two variables."); } int iDim; int nDim = nmonVars - 1; //Search for optional parameters list sParams = #; int nParams = size(sParams); //Default parameters //list svnumPts: intentionally not defined here; int fVerify = 0; int fComputeFan = 0; //intvec vzWeight: intentionally not defined here //Size check for the list of points to blow up int ivnumPts; int nvnumPts = size(svnumPts); //Check whether enough points were passed. if (nvnumPts < nDim) { ERROR("There were not enough points passed."); } for (ivnumPts = 1; ivnumPts <= nvnumPts; ivnumPts++) { if (typeof(svnumPts[ivnumPts]) != "vector") { ERROR("The passed points to modify along with the CEMDS must be of type \"vector\" and consist of numbers only."); } for (imonVars = 1; imonVars <= nmonVars; imonVars++) { if (leadexp(svnumPts[ivnumPts][imonVars]) != 0:nmonVars) { ERROR("The vectors passed must consist of numbers only."); } } } if (nParams > 0) { //First parameter: Verify the result or not. //Type check for verification flag if (typeof(sParams[1]) != "int") { ERROR("The parameter passed for \"fVerify\" must be of type \"int\" and equal \"0\" or \"1\"."); } fVerify = sParams[1]; //Second parameter: Compute a fan or not. if (nParams > 1) { //Type check for fan computation flag if (typeof(sParams[2]) != "int") { ERROR("The parameter passed for \"fComputeFan\" must be of type \"int\" and equal \"0\" or \"1\"."); } fComputeFan = sParams[2]; //Third parameter: An optional weight for initializing the fan computation. if (nParams > 2) { //Type check for weight if (typeof(sParams[3]) != "intvec") { ERROR("The parameter passed for \"vzWeight\" must be of type \"intvec\"."); } intvec vzWeight = sParams[3]; //There are at most 3 optional parameters if (nParams > 3) { ERROR("Too many optional parameters."); } } } } //Create the initial CEMDS for P_n intmat rvcvzP = intmatAppendCols(intmat(-1:nDim, nDim, 1), unitMatrix(nDim)); fan scnSigma = emptyFan(nDim); ideal spG = 0; CEMDS cemdsPnDim = createCEMDS(rvcvzP, scnSigma, spG); //Stretch the CEMDS in all hyperplanes containing at least n of the points to blow up. poly pNormal; ideal shpQPrimesIdeal; //Find the normal vectors to these hyperplanes. list svnumNormal = hyperplanes(svnumPts); int ivnumNormal; int nvnumNormal = size(svnumNormal); for (ivnumNormal = 1; ivnumNormal <= nvnumNormal; ivnumNormal++) { //From the normal vectors, generate the linear forms under which the associated hyperplanes vanish. pNormal = veronesePoly(1, svnumNormal[ivnumNormal]); shpQPrimesIdeal = shpQPrimesIdeal, pNormal; } //Remove redundant linear forms from the list of linear forms. shpQPrimesIdeal = simplify(shpQPrimesIdeal, 10); //Remove monomials from the list of linear forms. int nhpQPrimesIdeal = ncols(shpQPrimesIdeal); list shpQPrimes = shpQPrimesIdeal[1..nhpQPrimesIdeal]; shpQPrimes = removeMonomials(shpQPrimes); dbprint(printlevel - 1, ""); dbprint(printlevel - 1, "linearBlowup: The following linear relations were found:"); dbprint(printlevel - 1, shpQPrimes); //Stretch the initial CEMDS by the linear forms obtained above. def RStretched = stretchCEMDS(cemdsPnDim, shpQPrimes, svnumPts); setring RStretched; CEMDS cemdsStretched = gsResult[1]; list svnumPtsStretched = gsResult[2]; int nmonNewVars = nvars(RStretched); ideal spGDash = cemdsStretched.spG; //Compute the vanishing degrees of the stretched points to blow up. Add them to the fan of the stretched CEMDS as new rays. intmat rvcvzPStretched = cemdsStretched.rvcvzP; intmat rvcvzPDash = cemdsStretched.rvcvzP; int ivnumPtsStretched; int nvnumPtsStretched = size(svnumPtsStretched); //Compute the new ray associated with the first point to blow up and add it to the matrix of rays. intvec cvnVanishingDegrees = vanishingDegrees(list(), svnumPtsStretched[1]); intmat rvcvnVanishingDegrees = intmat(cvnVanishingDegrees, nmonNewVars, 1); rvcvzPDash = intmatAppendCol(rvcvzPDash, rvcvzPStretched * cvnVanishingDegrees); for (ivnumPtsStretched = 2; ivnumPtsStretched <= nvnumPtsStretched; ivnumPtsStretched++) { //Compute the new ray associated with the point to blow up in the modified CEMDS's fan and add it to the matrix of rays. cvnVanishingDegrees = vanishingDegrees(list(), svnumPtsStretched[ivnumPtsStretched]); rvcvnVanishingDegrees = intmatAppendCol(rvcvnVanishingDegrees, cvnVanishingDegrees); rvcvzPDash = intmatAppendCol(rvcvzPDash, rvcvzPStretched * cvnVanishingDegrees); } dbprint(printlevel - 1, ""); dbprint(printlevel - 1, "linearBlowup: The new matrix of rays:"); dbprint(printlevel - 1, rvcvzPDash); intvec rvfFakeVars = 0:nmonVars, 1:(nmonNewVars-nmonVars); //Perform the CEMDS modification defined by the rays added. CEMDS cemdsModified = modifyCEMDSCore(rvcvzPDash, scnSigma, spGDash, rvcvnVanishingDegrees, 0, rvfFakeVars); CEMDS cemdsResult; int fVerified = 0; if (defined(vzWeight)) { if (fVerify) { list sCompressResult = compressCEMDS(cemdsModified, fVerify, list(), fComputeFan, vzWeight); cemdsResult = sCompressResult[1]; fVerified = sCompressResult[2]; } else { cemdsResult = compressCEMDS(cemdsModified, fVerify, list(), fComputeFan, vzWeight); } } else { if (fVerify) { list sCompressResult = compressCEMDS(cemdsModified, fVerify, list(), fComputeFan); cemdsResult = sCompressResult[1]; fVerified = sCompressResult[2]; } else { cemdsResult = compressCEMDS(cemdsModified, fVerify, list(), fComputeFan); } } dbprint(printlevel + 1, ""); if (fVerified) { dbprint(printlevel + 1, "The resulting embedded space was successfully verified to be a CEMDS."); } else { dbprint(printlevel + 1, "The resulting embedded space could not be verified to be a CEMDS."); } dbprint(printlevel - 1, ""); dbprint(printlevel - 1, "The resulting CEMDS:"); dbprint(printlevel - 1, ""); dbprint(printlevel - 1, cemdsResult); return(cemdsResult); } example { "EXAMPLE:"; echo=2; //First example: No verification or fan computation ring R = 0,T(1..3),dp; vector vnumPt1 = [1,0,0]; vector vnumPt2 = [0,1,0]; vector vnumPt3 = [0,0,1]; vector vnumPt4 = [1,1,1]; list svnumPts = vnumPt1, vnumPt2, vnumPt3, vnumPt4; CEMDS cemdsBlowup = linearBlowup(svnumPts); def RBlowup = cemdsBlowup.R; setring RBlowup; print(cemdsBlowup); //Should yield // - A matrix with the following rays as columns // (1 0 -1 -1 0 1 0 0 0 0) // (0 -1 1 1 -1 0 0 0 0 0) // (1 -1 0 0 0 0 -1 1 0 0) // (1 0 -1 0 0 0 -1 0 1 0) // (0 -1 0 1 0 0 -1 0 0 1) // or there must exist a linear isomorphism mapping the output rays to those above. // - The empty fan of dimension 5; // - The ideal generated by // T(3)*T(8)-T(2)*T(9)+T(6)*T(10), // T(6)*T(7)-T(5)*T(8)+T(4)*T(9), // T(3)*T(7)-T(1)*T(9)+T(5)*T(10), // T(2)*T(7)-T(1)*T(8)+T(4)*T(10), // T(3)*T(4)-T(2)*T(5)+T(1)*T(6). setring R; //Second example: Verification and fan computation CEMDS cemdsBlowup2 = linearBlowup(svnumPts, 1, 1); def RBlowup2 = cemdsBlowup2.R; setring RBlowup2; print(cemdsBlowup2); //Should yield the same as above, but with an appropriate fan. } proc linearQuadricBlowupP2(list svnumPts, list #) "USAGE: linearQuadricBlowupP2(svnumPts); with svnumPts: list of vectors of numbers; linearQuadricBlowupP2(svnumPts, fVerify); with svnumPts: list of vectors of numbers, fVerify: int; linearQuadricBlowupP2(svnumPts, fVerify, fComputeFan); with svnumPts: list of vectors of numbers, fVerify: int, fComputeFan: int; linearQuadricBlowupP2(svnumPts, fVerify, fComputeFan, vzWeight); with svnumPts: list of vectors of numbers, fVerify: int, fComputeFan: int, vzWeight: intvec. ASSUME: The current basering contains at exactly three variables (thus representing the Cox ring of P2). There are at least 3 points to blow up containing a subset of 3 points in general position. PURPOSE: Blows up the projective space of dimension 2 in some of its points using relations of degree 1 and 2 for the stretching part only. RETURN: The projective space of dimension 2 blown up in the points specified. Whether the result is indeed a CEMDS is checked iff fVerify != 0. The CEMDS's fan is computed iff fComputeFan != 0. NOTE: By default: fVerify = 0; fComputeFan = 0; vzWeight is an element of the relative interior of the moving cone generated by the columns of the grading matrix rvcvzQ = gale(cemds.rvcvzP) extended by the passed polynomials' degrees. EXAMPLE: example linearQuadricBlowupP2; shows an example." { if (!defined(basering)) { ERROR("This function requires a basering to be defined before being called."); } int imonVars; int nmonVars = nvars(basering); if (nmonVars != 3) { ERROR("This function requires a basering with exactly three variables."); } int iDim; int nDim = nmonVars - 1; //Search for optional parameters list sParams = #; int nParams = size(sParams); //Default parameters //list svnumPts: intentionally not defined here; int fVerify = 0; int fComputeFan = 0; //intvec vzWeight: intentionally not defined here //Size check for the list of points to blow up int ivnumPts; int nvnumPts = size(svnumPts); //Check whether enough points were passed. if (nvnumPts < nDim) { ERROR("There were not enough points passed."); } for (ivnumPts = 1; ivnumPts <= nvnumPts; ivnumPts++) { if (typeof(svnumPts[ivnumPts]) != "vector") { ERROR("The passed points to modify along with the CEMDS must be of type \"vector\" and consist of numbers only."); } for (imonVars = 1; imonVars <= nmonVars; imonVars++) { if (leadexp(svnumPts[ivnumPts][imonVars]) != 0:nmonVars) { ERROR("The vectors passed must consist of numbers only."); } } } if (nParams > 0) { //First parameter: Verify the result or not. //Type check for verification flag if (typeof(sParams[1]) != "int") { ERROR("The parameter passed for \"fVerify\" must be of type \"int\" and equal \"0\" or \"1\"."); } fVerify = sParams[1]; //Second parameter: Compute a fan or not. if (nParams > 1) { //Type check for fan computation flag if (typeof(sParams[2]) != "int") { ERROR("The parameter passed for \"fComputeFan\" must be of type \"int\" and equal \"0\" or \"1\"."); } fComputeFan = sParams[2]; //Third parameter: An optional weight for initializing the fan computation. if (nParams > 2) { //Type check for weight if (typeof(sParams[3]) != "intvec") { ERROR("The parameter passed for \"vzWeight\" must be of type \"intvec\"."); } intvec vzWeight = sParams[3]; //There are at most 3 optional parameters if (nParams > 3) { ERROR("Too many optional parameters."); } } } } //Create the initial CEMDS for P_2 intmat rvcvzP = intmatAppendCols(intmat(-1:nDim, nDim, 1), unitMatrix(nDim)); fan scnSigma = emptyFan(nDim); ideal spG = 0; CEMDS cemdsPnDim = createCEMDS(rvcvzP, scnSigma, spG); //Stretch the CEMDS in all hyperplanes containing at least 2 of the points to blow up //and in all quadrics containing at least 5 of the points to blow up. poly pNormal; ideal shpQPrimesIdeal; //Find the normal vectors to the hyperplanes. list svnumNormal = hyperplanes(svnumPts); int ivnumNormal; int nvnumNormal = size(svnumNormal); for (ivnumNormal = 1; ivnumNormal <= nvnumNormal; ivnumNormal++) { //From the normal vectors, generate the linear forms under which the associated hyperplanes vanish. pNormal = veronesePoly(1, svnumNormal[ivnumNormal]); shpQPrimesIdeal = shpQPrimesIdeal, pNormal; } //Find the coefficient vectors to the quadrics. list svnumCoeffs = quadrics(svnumPts); int ivnumCoeffs; int nvnumCoeffs = size(svnumCoeffs); for (ivnumCoeffs = 1; ivnumCoeffs <= nvnumCoeffs; ivnumCoeffs++) { //From the coefficient vectors, generate the quadric forms under which the associated quadrics vanish. pNormal = veronesePoly(2, svnumCoeffs[ivnumCoeffs]); shpQPrimesIdeal = shpQPrimesIdeal, pNormal; } //Remove redundant linear forms from the list of linear forms. shpQPrimesIdeal = simplify(shpQPrimesIdeal, 10); //Remove monomials from the list of linear forms. int nhpQPrimesIdeal = ncols(shpQPrimesIdeal); list shpQPrimes = shpQPrimesIdeal[1..nhpQPrimesIdeal]; shpQPrimes = removeMonomials(shpQPrimes); dbprint(printlevel - 1, ""); dbprint(printlevel - 1, "linearQuadricBlowupP2: The following linear relations were found:"); dbprint(printlevel - 1, shpQPrimes); //Stretch the initial CEMDS by the linear forms obtained above. def RStretched = stretchCEMDS(cemdsPnDim, shpQPrimes, svnumPts); setring RStretched; CEMDS cemdsStretched = gsResult[1]; list svnumPtsStretched = gsResult[2]; dbprint(printlevel - 1, ""); dbprint(printlevel - 1, "linearQuadricBlowupP2: The stretched CEMDS:"); dbprint(printlevel - 1, cemdsStretched); int nmonNewVars = nvars(RStretched); ideal spGDash = cemdsStretched.spG; //Compute the vanishing degrees of the stretched points to blow up. Add them to the fan of the stretched CEMDS as new rays. intmat rvcvzPStretched = cemdsStretched.rvcvzP; intmat rvcvzPDash = cemdsStretched.rvcvzP; int ivnumPtsStretched; int nvnumPtsStretched = size(svnumPtsStretched); //Compute the new ray associated with the first point to blow up and add it to the matrix of rays. intvec cvnVanishingDegrees = vanishingDegrees(list(), svnumPtsStretched[1]); intmat rvcvnVanishingDegrees = intmat(cvnVanishingDegrees, nmonNewVars, 1); rvcvzPDash = intmatAppendCol(rvcvzPDash, rvcvzPStretched * cvnVanishingDegrees); for (ivnumPtsStretched = 2; ivnumPtsStretched <= nvnumPtsStretched; ivnumPtsStretched++) { //Compute the new ray associated with the point to blow up in the modified CEMDS's fan and add it to the matrix of rays. cvnVanishingDegrees = vanishingDegrees(list(), svnumPtsStretched[ivnumPtsStretched]); rvcvnVanishingDegrees = intmatAppendCol(rvcvnVanishingDegrees, cvnVanishingDegrees); rvcvzPDash = intmatAppendCol(rvcvzPDash, rvcvzPStretched * cvnVanishingDegrees); } dbprint(printlevel - 1, ""); dbprint(printlevel - 1, "linearQuadricBlowupP2: The new matrix of rays:"); dbprint(printlevel - 1, rvcvzPDash); intvec rvfFakeVars = 1:(nmonNewVars); //Perform the CEMDS modification defined by the rays added. CEMDS cemdsModified = modifyCEMDSCore(rvcvzPDash, scnSigma, spGDash, rvcvnVanishingDegrees, 0, rvfFakeVars); dbprint(printlevel - 1, ""); dbprint(printlevel - 1, "linearQuadricBlowupP2: The modified CEMDS:"); dbprint(printlevel - 1, cemdsModified); CEMDS cemdsResult; int fVerified = 0; if (defined(vzWeight)) { if (fVerify) { list sCompressResult = compressCEMDS(cemdsModified, fVerify, list(), fComputeFan, vzWeight); cemdsResult = sCompressResult[1]; fVerified = sCompressResult[2]; } else { cemdsResult = compressCEMDS(cemdsModified, fVerify, list(), fComputeFan, vzWeight); } } else { if (fVerify) { list sCompressResult = compressCEMDS(cemdsModified, fVerify, list(), fComputeFan); cemdsResult = sCompressResult[1]; fVerified = sCompressResult[2]; } else { cemdsResult = compressCEMDS(cemdsModified, fVerify, list(), fComputeFan); } } dbprint(printlevel + 1, ""); if (fVerified) { dbprint(printlevel + 1, "The resulting embedded space was successfully verified to be a CEMDS."); } else { dbprint(printlevel + 1, "The resulting embedded space could not be verified to be a CEMDS."); } return(cemdsResult); } example { "EXAMPLE:"; echo=2; ring R = (0,a,b),T(1..3),dp; vector vnumPt1 = [1,0,0]; vector vnumPt2 = [0,1,0]; vector vnumPt3 = [0,0,1]; vector vnumPt4 = [1,1,1]; vector vnumPt5 = [1,1,0]; list svnumPts = vnumPt1, vnumPt2, vnumPt3, vnumPt4, vnumPt5; CEMDS cemdsBlowup = linearQuadricBlowupP2(svnumPts); def RBlowup = cemdsBlowup.R; setring RBlowup; print(cemdsBlowup); //Should yield // - A matrix with the following rays as columns // (-1 1 0 0 1 -1 0 0 0 0 0) // ( 1 -1 0 0 0 0 -1 1 0 0 0) // ( 1 0 -1 1 0 -1 -2 0 2 0 0) // ( 1 1 -1 1 0 0 -2 0 1 1 0) // (-1 0 1 0 0 0 1 0 -1 0 1) // or there must exist a linear isomorphism mapping the output rays to those above. // - The empty fan of dimension 5; // - The ideal generated by // -T(3)*T(8)*T(11)+T(2)*T(9)-T(6)*T(10), // -T(3)*T(7)*T(11)+T(1)*T(9)-T(5)*T(10), // T(4)*T(9)*T(11)+T(6)*T(7)-T(5)*T(8), // T(4)*T(10)*T(11)+T(2)*T(7)-T(1)*T(8), // -T(3)*T(4)*T(11)^2+T(2)*T(5)-T(1)*T(6). setring R; //Second example: Verification and fan computation CEMDS cemdsBlowup2 = linearQuadricBlowupP2(svnumPts, 1, 1); def RBlowup2 = cemdsBlowup2.R; setring RBlowup2; print(cemdsBlowup2); //Should yield the same as above, but with an appropriate fan. } proc modifyCEMDS(CEMDS cemds, list sRaysFan, list #) "USAGE: modifyCEMDS(cemds, sRaysFan); with cemds: CEMDS, sRaysFan: list of an intmat and a fan; modifyCEMDS(cemds, sRaysFan, fVerify); with cemds: CEMDS, sRaysFan: list of an intmat and a fan, fVerify: int; modifyCEMDS(cemds, sRaysFan, fVerify, rvfFakeVars); with cemds: CEMDS, sRaysFan: list of an intmat and a fan, fVerify: int, rvfFakeVars: intvec; modifyCEMDS(cemds, sRaysFan, fVerify, rvfFakeVars, svnumPts); with cemds: CEMDS, sRaysFan: list of an intmat and a fan, fVerify: int, rvfFakeVars: intvec, svnumPts: list of vectors of numbers. ASSUME: (i) The matrix passed in \"sRaysFan\" is the matrix \"cemds.rvcvzP\" extended by further columns that are integer sums of columns of this matrix; (ii) The fan passed in \"sRaysFan\" has the columns of the matrix passed in \"sRaysFan\" as its rays; (iii) rvfFakeVars marks only indices of fake variables as true (i.e. variables for which there exists a generator of the type var(i) - p with p not containing var(i) in the CEMDS's ideal); (iv) All points to modify along with the CEMDS, if there were any passed, are contained in the CEMDS to modify. PURPOSE: Blows up the ambient toric variety of a CEMDS by a refining of its fan. Optionally checks whether the result really is a CEMDS and transfers points from the unmodified CEMDS to the modified one. RETURN: If no optional output was requested, the modified CEMDS only. If only verification was requested, but no points or an empty list of points to modify along with the CEMDS were passed, a list containing the following (in order of appearance): (i) the modified CEMDS that is obtained by blowing up the passed CEMDS's ambient toric variety by adding some rays to its fan and performing stellar subdivision (with or without a fan, depending on fComputeFan) (always present); (ii) the verification result (present as fVerify was passed and fVerify != 0). If a non-empty list of points to modify along with the CEMDS was passed, the modified CEMDS's basering. SIDE EFFECTS: Iff a non-empty list of points to modify along with the CEMDS was passed, exports a list gsResult containing the following (if present and in order of appearance): (i) the modified CEMDS that is obtained by blowing up the passed CEMDS's ambient toric variety by adding some rays to its fan and performing stellar subdivision (with or without a fan, depending on fComputeFan) (always present); (ii) the verification result (present iff fVerify was passed and fVerify != 0), (iii) the list of modified points (present as there was a list passed). NOTE: By default: fVerify = 0; rvfFakeVars = 0:nvars; svnumPts = list(). EXAMPLE: example modifyCEMDS; shows an example." { //The CEMDS to modify must be assigned to a ring. def R = cemds.R; if (!defined(R)) { ERROR("The passed CEMDS is not assigned to any ring."); } //Switch to the CEMDS's ring setring R; int imonVars; int nmonVars = nvars(R); //The passed list must contain the matrix and the fan to use for the new CEMDS. if (typeof(sRaysFan[1]) != "intmat") { ERROR("The first parameter in \"sRaysFan\" must be of type \"intmat\"."); } if (ncols(sRaysFan[1]) <= ncols(cemds.rvcvzP)) { ERROR("The matrix passed in \"sRaysFan\" must contain the matrix of the original CEMDS."); } if (typeof(sRaysFan[2]) != "fan") { ERROR("The second parameter in \"sRaysFan\" must be of type \"fan\"."); } intmat rvcvzPDash = sRaysFan[1]; int izPDash; int nzPDash = nrows(rvcvzPDash); fan scnSigmaDash = sRaysFan[2]; int icvzNewRays; int ncvzNewRays = ncols(rvcvzPDash) - ncols(cemds.rvcvzP); //Search for optional arguments list sParams = #; int nParams = size(sParams); //Default arguments int fVerify = 0; intvec rvfFakeVars = 0:nmonVars; list svnumPts = list(); int ivzPts; int nvzPts = 0; //First argument: Verify whether the resulting ES is a CEMDS or not. if (nParams > 0) { //Type check for verification flag if (typeof(sParams[1]) != "int") { ERROR("The argument passed for \"fVerify\" must be of type \"int\" and equal \"0\" or \"1\"."); } fVerify = sParams[1]; //Second argument: An intvec specifying which of the variables are fake. if (nParams > 1) { //Type check for the fake variable flags if (typeof(sParams[2]) != "intvec") { ERROR("The argument passed for \"rvfFakeVars\" must be of type \"intvec\"."); } if (nmonVars != size(rvfFakeVars)) { ERROR("The CEMDS's ring's variable count and the number of entries of \"rvfFakeVars\" differ."); } rvfFakeVars = sParams[2]; //Third argument: A list of points whose coordinates shall be transferred to the modified CEMDS. if (nParams > 2) { //Type and size check for the list of points to modify along with the CEMDS if (typeof(sParams[3]) != "list") { ERROR("The argument passed for \"svnumPts\" must be of type \"list\" and contain and vectors of numbers only."); } svnumPts = sParams[3]; int ivnumPts; int nvnumPts = size(svnumPts); for (ivnumPts = 1; ivnumPts <= nvnumPts; ivnumPts++) { if (typeof(svnumPts[ivnumPts]) != "vector") { ERROR("The passed points to modify along with the CEMDS must be of type \"vector\" and consist of numbers only."); } for (imonVars = 1; imonVars <= nmonVars; imonVars++) { if (leadexp(svnumPts[ivnumPts][imonVars]) != 0:nmonVars) { ERROR("The vectors passed must consist of numbers only."); } } } //There are at most 3 optional arguments. if (nParams > 3) { ERROR("Too many optional parameters."); } } } } ideal spG = cemds.spG; //Idea for getting the vanishing degrees: Solve Px = r (with P: the old rays' matrix cemds.rvcvzP, r: the new ray) and x >= 0 (component-by-component) //Approach by finding an arbitrary edge via the simplex algorithm, each column having cost 0. //NOTE: Implementation in Q or Q(a,b,...) may be more accurate as coefficients should be integers only. ring RReal = (real,10),T(1..nmonVars),dp; intvec cvzNewRay; intmat rvcvzPExtended = intmatAppendRows(intmat(0:nmonVars, 1, nmonVars), cemds.rvcvzP); matrix matnumSimplex; list sSimplexResult; intvec cvnVanishingDegree = 0:ncols(cemds.rvcvzP); intmat rvcvnVanishingDegrees; // for (icvzNewRays = 1; icvzNewRays <= ncvzNewRays; icvzNewRays++) { //Initialize the matrix for the simplex algorithm. cvzNewRay = 0,intmatTakeCol(rvcvzPDash, nmonVars + icvzNewRays); matnumSimplex = matrix(cvzNewRay, nzPDash + 1, 1); matnumSimplex = concat(matnumSimplex, matrix(-rvcvzPExtended)); for (izPDash = 2; izPDash <= nzPDash + 1; izPDash++) { if (matnumSimplex[izPDash, 1] < 0) { matnumSimplex = multrow(matnumSimplex, izPDash, -1); } } //Find coefficients via simplex algorithm. These should be integers if the extended CEMDS's matrix allows it. sSimplexResult = simplex(matnumSimplex, nzPDash, nmonVars, 0, 0, nzPDash); cvnVanishingDegree = 0:ncols(cemds.rvcvzP); for (izPDash = 1; izPDash <= nrows(sSimplexResult[1]); izPDash++) { if (izPDash == 1) { } else { if ( sSimplexResult[3][izPDash-1] <= sSimplexResult[6] ) { cvnVanishingDegree[sSimplexResult[3][izPDash-1]] = int(sSimplexResult[1][izPDash, 1]); if (cvnVanishingDegree[sSimplexResult[3][izPDash-1]] - sSimplexResult[1][izPDash, 1] != 0.0) { ERROR("Rounding error: The simplex algorithm did not return integers only."); } } } } //Compose vanishing degree matrix (of integer coefficients) if (icvzNewRays == 1) { rvcvnVanishingDegrees = intmat(cvnVanishingDegree, ncols(cemds.rvcvzP), 1); } else { rvcvnVanishingDegrees = intmatAppendCol(rvcvnVanishingDegrees, cvnVanishingDegree); } } dbprint(printlevel - 1, ""); dbprint(printlevel - 1, "modifyCEMDS: Found the following weights for the added rays:"); dbprint(printlevel - 1, rvcvnVanishingDegrees); //Having obtained the vanishing degrees, we can transfer the CEMDS's ideal in the core function. setring R; return(modifyCEMDSCore(rvcvzPDash, scnSigmaDash, spG, rvcvnVanishingDegrees, fVerify, rvfFakeVars, svnumPts)); } example { "EXAMPLE:"; echo=2; // //First example: No optional arguments. ring R = (0,a,b),(T(1),T(2),T(3),T(4),T(5)),(dp(5),C); intmat rvcvzP = intmat(intvec( 1,0,-1,0,0, 0,1,-1,0,0, -1,0,0,1,0, -1,0,0,0,1 ), 4, 5); intmat cvrvz; list scn = list(); list scnList = list(); cvrvz = intmat(intvec( -1, -1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1 ), 4, 4); scn = coneViaPoints(cvrvz); scnList = scnList + scn; cvrvz = intmat(intvec( 1, 0, -1, -1, -1, -1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1 ), 4, 4); scn = coneViaPoints(cvrvz); scnList = scnList + scn; cvrvz = intmat(intvec( -1, -1, 0, 0, 0, 1, 0, 0, 1, 0, -1, -1, 0, 0, 0, 1 ), 4, 4); scn = coneViaPoints(cvrvz); scnList = scnList + scn; cvrvz = intmat(intvec( -1, -1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, -1, -1 ), 4, 4); scn = coneViaPoints(cvrvz); scnList = scnList + scn; cvrvz = intmat(intvec( 1, 0, -1, -1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1 ), 4, 4); scn = coneViaPoints(cvrvz); scnList = scnList + scn; fan scnSigma = fanViaCones(scnList); ideal spG = (-a)*T(1)+T(2)+T(4), (-b)*T(1)+T(2)+T(5); CEMDS cemds = createCEMDS(rvcvzP, scnSigma, spG); intmat P2 = cemds.rvcvzP; fan Sigma2 = cemds.scnSigma; intvec newRay1 = intmatTakeCol(P2, 2) + intmatTakeCol(P2, 3); P2 = intmatAppendCol(P2, newRay1); Sigma2 = stellarSubdivision(Sigma2, newRay1); intvec newRay2 = intmatTakeCol(P2, 1) + intmatTakeCol(P2, 3); P2 = intmatAppendCol(P2, newRay2); Sigma2 = stellarSubdivision(Sigma2, newRay2); intvec newRay3 = intmatTakeCol(P2, 3) + intmatTakeCol(P2, 4); P2 = intmatAppendCol(P2, newRay3); Sigma2 = stellarSubdivision(Sigma2, newRay3); intvec newRay4 = intmatTakeCol(P2, 3) + intmatTakeCol(P2, 5); P2 = intmatAppendCol(P2, newRay4); Sigma2 = stellarSubdivision(Sigma2, newRay4); list sRaysFan = P2, Sigma2; intvec rvfFakeVars = 0,0,0,1,1; CEMDS cemdsResult = modifyCEMDS(cemds, sRaysFan); print(cemdsResult); //Should yield the CEMDS consisting of // - The matrix rvcvzPDash; // - The fan scnSigmaDash; // - The ideal generated by // # T(2)*T(6)+(-a)*T(1)*T(7)+T(4)*T(8), // # (a-b)*T(1)*T(7)-T(4)*T(8)+T(5)*T(9). //Second example: With verification list smodifyResult = modifyCEMDS(cemds, sRaysFan, 1); print(smodifyResult); //Should yield a list containing // - The CEMDS consisting of // * The matrix rvcvzPDash; // * The fan scnSigmaDash; // * The ideal generated by // # T(2)*T(6)+(-a)*T(1)*T(7)+T(4)*T(8), // # (a-b)*T(1)*T(7)-T(4)*T(8)+T(5)*T(9). // - The verification result (1). //Third example: Without verification, with fake vars. // intvec rvfFakeVars = 0,0,0,1,1,1; CEMDS cemdsResult2 = modifyCEMDS(cemds, sRaysFan, 0, rvfFakeVars); print(cemdsResult2); //Should yield the CEMDS consisting of // - The matrix rvcvzPDash; // - The fan scnSigmaDash; // - The ideal generated by // # T(2)*T(6)+(-a)*T(1)*T(7)+T(4)*T(8), // # (a-b)*T(1)*T(7)-T(4)*T(8)+T(5)*T(9). //Fourth example: With point to modify along the CEMDS. setring R; vector vnumPt1 = [1,1,0,0,0]; list svnumPts = vnumPt1; def RModified = modifyCEMDS(cemds, sRaysFan, 1, rvfFakeVars, svnumPts); setring RModified; print(gsResult); //Should yield a list containing // - The CEMDS consisting of // * The matrix rvcvzPDash; // * The fan scnSigmaDash; // * The ideal generated by // # T(2)*T(6)+(-a)*T(1)*T(7)+T(4)*T(8), // # (a-b)*T(1)*T(7)-T(4)*T(8)+T(5)*T(9). // - The verification result (1). // - A list containing the modified point [1,1,0,0,0,1,1,1,1]. } proc modifyCEMDS2(CEMDS cemds, intvec rvnVanishingDegrees, list #) "USAGE: modifyCEMDS2(cemds, rvnVanishingDegrees); with cemds: CEMDS, rvnVanishingDegrees: intvec; modifyCEMDS2(cemds, rvnVanishingDegrees, fVerify); with cemds: CEMDS, rvnVanishingDegrees: intvec, fVerify: int; modifyCEMDS2(cemds, rvnVanishingDegrees, fVerify, rvfFakeVars); with cemds: CEMDS, rvnVanishingDegrees: intvec, fVerify: int, rvfFakeVars: intvec; modifyCEMDS2(cemds, rvnVanishingDegrees, fVerify, rvfFakeVars, svnumPts); with cemds: CEMDS, rvnVanishingDegrees: intvec, fVerify: int, rvfFakeVars: intvec, svnumPts: list of vectors of numbers; modifyCEMDS2(cemds, rvnVanishingDegrees, fVerify, rvfFakeVars, svnumPts, fComputeFan); with cemds: CEMDS, rvnVanishingDegrees: intvec, fVerify: int, rvfFakeVars: intvec, svnumPts: list of vectors of numbers, fComputeFan: int. ASSUME: (i) rvfFakeVars marks only indices of fake variables as true (i.e. variables for which there exists a generator of the type var(i) - p with p not containing var(i) in the CEMDS's ideal) (ii) All points to modify along with the CEMDS, if there were any passed, are contained in the CEMDS to modify. PURPOSE: Blows up the ambient toric variety of a CEMDS by adding some rays to its fan and performing stellar subdivision. Optionally checks whether the result really is a CEMDS and transfers points from the unmodified CEMDS to the modified one. RETURN: If no optional output was requested, the modified CEMDS only. If only verification was requested, but no points or an empty list of points to modify along with the CEMDS were passed, a list containing the following (in order of appearance): (i) the modified CEMDS that is obtained by blowing up the passed CEMDS's ambient toric variety by adding some rays to its fan and performing stellar subdivision (with or without a fan, depending on fComputeFan) (always present); (ii) the verification result (present as fVerify was passed and fVerify != 0). If a non-empty list of points to modify along with the CEMDS was passed, the modified CEMDS's basering. SIDE EFFECTS: Iff a non-empty list of points to modify along with the CEMDS was passed, exports a list gsResult containing the following (if present and in order of appearance): (i) the modified CEMDS that is obtained by blowing up the passed CEMDS's ambient toric variety by adding some rays to its fan and performing stellar subdivision (with or without a fan, depending on fComputeFan) (always present); (ii) the verification result (present iff fVerify was passed and fVerify != 0), (iii) the list of modified points (present as there was a list passed). NOTE: By default: fVerify = 0; rvfFakeVars = 0:nvars; svnumPts = list(); fComputeFan = 0. EXAMPLE: example modifyCEMDS2; shows an example." { //The CEMDS to modify must be assigned to a ring. def R = cemds.R; if (!defined(R)) { ERROR("The passed CEMDS is not assigned to any ring."); } //Switch to the CEMDS's ring setring R; int imonVars; int nmonVars = nvars(R); //Search for parameters list sParams = #; int nParams = size(sParams); //Default parameters int fVerify = 0; intvec rvfFakeVars = 0:nmonVars; list svnumPts = list(); int fComputeFan = 0; //First parameter: Verify whether the resulting ES is a CEMDS or not. if (nParams > 0) { //Type check for verification flag if (typeof(sParams[1]) != "int") { ERROR("The parameter passed for \"fVerify\" must be of type \"int\" and equal \"0\" or \"1\"."); } fVerify = sParams[1]; //Second parameter: An intvec specifying which of the variables are fake. if (nParams > 1) { //Type check for the fake variable flags if (typeof(sParams[2]) != "intvec") { ERROR("The parameter passed for \"rvfFakeVars\" must be of type \"intvec\"."); } rvfFakeVars = sParams[2]; if (nmonVars != size(rvfFakeVars)) { ERROR("The CEMDS's ring's variable count and the number of entries of \"rvfFakeVars\" differ."); } //Third parameter: A list of points whose coordinates shall be transferred to the modified CEMDS. if (nParams > 2) { //Type and size check for the list of points to modify along with the CEMDS if (typeof(sParams[3]) != "list") { ERROR("The parameter passed for \"svnumPts\" must be of type \"list\" and contain vectors of numbers only."); } svnumPts = sParams[3]; int ivnumPts; int nvnumPts = size(svnumPts); for (ivnumPts = 1; ivnumPts <= nvnumPts; ivnumPts++) { if (typeof(svnumPts[ivnumPts]) != "vector") { ERROR("The passed points to modify along with the CEMDS must be of type \"vector\" and consist of numbers only."); } for (imonVars = 1; imonVars <= nmonVars; imonVars++) { if (leadexp(svnumPts[ivnumPts][imonVars]) != 0:nmonVars) { ERROR("The vectors passed must consist of numbers only."); } } } //Fourth parameter: Whether to compute a fan or not. if (nParams > 3) { //Type check for fan computation flag if (typeof(sParams[4]) != "int") { ERROR("The parameter passed for \"fComputeFan\" must be of type \"int\" and equal \"0\" or \"1\"."); } fComputeFan = sParams[4]; //There are at most 4 optional parameters. if (nParams > 4) { ERROR("Too many optional parameters."); } } } } } int ncvzP = ncols(cemds.rvcvzP); if (ncvzP != size(rvnVanishingDegrees)) { ERROR("\"rvnVanishingDegrees\" must have as many entries as the CEMDS's matrix has columns."); } //Compute the new ray in the modified CEMDS's fan intvec cvzNewRay = cemds.rvcvzP * rvnVanishingDegrees; //Add the ray to the matrix of rays intmat rvcvzPDash = intmatAppendCol(cemds.rvcvzP, cvzNewRay); //Add the ray to the fan itself by performing a stellar subdivision fan scnSigmaDash = emptyFan(nrows(rvcvzPDash)); if (fComputeFan) { scnSigmaDash = cemds.scnSigma; scnSigmaDash = stellarSubdivision(scnSigmaDash, cvzNewRay); } ideal spG = cemds.spG; return(modifyCEMDSCore(rvcvzPDash, scnSigmaDash, spG, intmat(rvnVanishingDegrees, ncvzP, 1), fVerify, rvfFakeVars, svnumPts)); } example { "EXAMPLE:"; echo=2; //First example: No optional arguments. ring R = 0,T(1..3),dp; //The CEMDS's rays' matrix intmat rvcvzP[2][3] = 1,0,-1, 0,1,-1; //The CEMDS's fan intmat cvrvz1[2][2] = 1,0, 0,1; intmat cvrvz2[2][2] = 0,1, -1,-1; intmat cvrvz3[2][2] = -1,-1, 1,0; cone cn1 = coneViaPoints(cvrvz1); cone cn2 = coneViaPoints(cvrvz2); cone cn3 = coneViaPoints(cvrvz3); fan scnSigma = fanViaCones(cn1, cn2, cn3); //The CEMDS's embedding ideal ideal spG = 0; CEMDS cemds = createCEMDS(rvcvzP, scnSigma, spG); intvec rvnVanishingDegrees = 0,1,1; CEMDS cemdsResult = modifyCEMDS2(cemds, rvnVanishingDegrees); //Obtains the CEMDS's ring def RModified = cemdsResult.R; setring RModified; print(cemdsResult); //Should yield the CEMDS consisting of // - A matrix with the entries // (1 0 -1 -1) // (0 1 -1 0) // - The empty fan of ambient dimension 2. // - The zero ideal. //Second example: Verify that the result indeed is a CEMDS setring R; list sModifyResult = modifyCEMDS2(cemds, rvnVanishingDegrees, 1); print(sModifyResult); //Should yield a list with two entries: // - the same CEMDS as in the first example // - the verification result (1). kill cemds; kill rvcvzP; kill cvrvz1, cvrvz2, cvrvz3; kill cn1, cn2, cn3; kill scnSigma; kill spG; kill rvnVanishingDegrees; kill cemdsResult; kill sModifyResult; kill RModified; kill R; //Third example: Fake variables are marked. ring R = (0,a),T(1..9),dp; //The CEMDS's rays' matrix intmat rvcvzP[5][9] = 1, 0, -1, 1, 0, -1, 0, 0, 0, 0, 1, -1, 1, -1, 0, 0, 0, 0, 0, -1, 0, -1, 0, 0, 1, 0, 0, -1, 0, 0, -1, 0, 0, 0, 1, 0, -1, 0, 0, 0, -1, 0, 0, 0, 1; //The CEMDS's fan intmat cvrvz; list scn = list(); list scnList = list(); cvrvz = intmat(intvec( 0, -1, 0, 0, 0, -1, -1, 0, 0, 1, 0, 1, -1, 0, -1, 1, 1, -1, -1, -1 ), 4, 5); scn = coneViaPoints(cvrvz); scnList = scnList + scn; cvrvz = intmat(intvec( -1, 0, 0, 0, 0, 1, 1, -1, -1, -1, -1, -1, 0, 0, 1, 1, 0, 0, -1, -1 ), 4, 5); scn = coneViaPoints(cvrvz); scnList = scnList + scn; cvrvz = intmat(intvec( -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, -1, -1, 0, 1, -1, 0, -1 ), 4, 5); scn = coneViaPoints(cvrvz); scnList = scnList + scn; cvrvz = intmat(intvec( -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, -1, -1, 0, 0, 1 ), 5, 5); scn = coneViaPoints(cvrvz); scnList = scnList + scn; cvrvz = intmat(intvec( -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1, -1, 0, -1, 0, 0, 0, 1, 0, -1, -1, 0, 0, 1 ), 5, 5); scn = coneViaPoints(cvrvz); scnList = scnList + scn; cvrvz = intmat(intvec( -1, 0, 0, 0, 0, -1, -1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1 ), 5, 5); scn = coneViaPoints(cvrvz); scnList = scnList + scn; cvrvz = intmat(intvec( -1, 0, 0, 0, 0, -1, -1, 0, 0, 1, 0, 1, -1, 0, -1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1 ), 5, 5); scn = coneViaPoints(cvrvz); scnList = scnList + scn; cvrvz = intmat(intvec( -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, -1, -1, -1, -1, 0, 0, 1 ), 5, 5); scn = coneViaPoints(cvrvz); scnList = scnList + scn; cvrvz = intmat(intvec( -1, -1, 0, 0, 1, 0, -1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1 ), 5, 5); scn = coneViaPoints(cvrvz); scnList = scnList + scn; cvrvz = intmat(intvec( -1, 0, 0, 0, 0, -1, -1, 0, 0, 1, 0, 1, -1, 0, -1, 1, 1, -1, -1, -1, 0, 0, 0, 0, 1 ), 5, 5); scn = coneViaPoints(cvrvz); scnList = scnList + scn; cvrvz = intmat(intvec( -1, 0, 0, 0, 0, 0, 1, -1, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1 ), 5, 5); scn = coneViaPoints(cvrvz); scnList = scnList + scn; cvrvz = intmat(intvec( -1, -1, 0, 0, 1, 0, -1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, -1, -1, 0, 0, 0, 0, 1 ), 5, 5); scn = coneViaPoints(cvrvz); scnList = scnList + scn; cvrvz = intmat(intvec( -1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, -1, -1, -1, 0, 1, -1, 0, -1 ), 5, 5); scn = coneViaPoints(cvrvz); scnList = scnList + scn; cvrvz = intmat(intvec( -1, -1, 0, 0, 1, 0, -1, 0, 0, 0, 1, 1, -1, -1, -1, 1, 0, 0, -1, -1, 0, 0, 0, 0, 1 ), 5, 5); scn = coneViaPoints(cvrvz); scnList = scnList + scn; cvrvz = intmat(intvec( 1, 0, 0, -1, -1, 0, -1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1 ), 5, 5); scn = coneViaPoints(cvrvz); scnList = scnList + scn; cvrvz = intmat(intvec( -1, 0, 0, 0, 0, 1, 0, 0, -1, -1, 0, 0, 1, 0, 0, 0, 1, -1, 0, -1, 1, 1, -1, -1, -1 ), 5, 5); scn = coneViaPoints(cvrvz); scnList = scnList + scn; cvrvz = intmat(intvec( 0, 1, -1, 0, -1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, -1, -1, -1 ), 5, 5); scn = coneViaPoints(cvrvz); scnList = scnList + scn; cvrvz = intmat(intvec( 0, 0, 0, 0, 1, 0, -1, 0, 0, 0, 1, 1, -1, -1, -1, 0, 0, 0, 1, 0, 1, 0, 0, -1, -1 ), 5, 5); scn = coneViaPoints(cvrvz); scnList = scnList + scn; cvrvz = intmat(intvec( 0, 1, -1, 0, -1, 0, -1, 0, 0, 0, 1, 0, 0, -1, -1, 0, 0, 0, 1, 0, 1, 1, -1, -1, -1 ), 5, 5); scn = coneViaPoints(cvrvz); scnList = scnList + scn; cvrvz = intmat(intvec( 1, 0, 0, -1, -1, 1, 1, -1, -1, -1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1 ), 5, 5); scn = coneViaPoints(cvrvz); scnList = scnList + scn; cvrvz = intmat(intvec( 0, 1, -1, 0, -1, 1, 0, 0, -1, -1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, -1, -1, -1 ), 5, 5); scn = coneViaPoints(cvrvz); scnList = scnList + scn; fan scnSigma = fanViaCones(scnList); //The CEMDS's embedding ideal poly p1 = T(7) - T(2)*T(4) + a*T(3)*T(5); poly p2 = T(8) - T(1)*T(4) + T(3)*T(6); poly p3 = T(9) - a*T(1)*T(5) + T(2)*T(6); ideal spG = p1, p2, p3; CEMDS cemds = createCEMDS(rvcvzP, scnSigma, spG); intvec rvnVanishingDegrees = 0,0,0,0,0,0,1,1,1; intvec rvfFakeVars = 0,0,0,0,0,0,1,1,1; CEMDS cemdsResult = modifyCEMDS2(cemds, rvnVanishingDegrees, 0, rvfFakeVars); print(cemdsResult); //Should yield the CEMDS consisting of: // - A matrix with the entries // ( 1 0 -1 1 0 -1 0 0 0 0) // ( 0 1 -1 1 -1 0 0 0 0 0) // ( 0 -1 0 -1 0 0 1 0 0 1) // (-1 0 0 -1 0 0 0 1 0 1) // (-1 0 0 0 -1 0 0 0 1 1) // - The empty fan of ambient dimension 5 // - The ideal generated by // T(7)*T(10) - T(2)*T(4) + a*T(3)*T(5), // T(8)*T(10) - T(1)*T(4) + T(3)*T(6), // T(9)*T(10) - a*T(1)*T(5) + T(2)*T(6), // T(6)*T(7) - a*T(5)*T(8) + T(4)*T(9), // T(1)*T(7) - T(2)*T(8) + T(3)*T(9). //Fourth example: Fake variables are marked and points are transferred to the modified CEMDS. vector vnumPt = [1, 2, 1, 1, 1, 2, 1, -1, -3]; list svnumPts = vnumPt; def RModified = modifyCEMDS2(cemds, rvnVanishingDegrees, 0, rvfFakeVars, svnumPts); setring RModified; list sModifyResult = gsResult; print(sModifyResult); //Should yield a list with two entries: // - the same CEMDS as in the third example // - a list containing the transferred point [1, 2, 1, 1, 1, 2, 1, -1, -3, 1] kill cemds; kill rvcvzP; kill scnSigma; kill rvnVanishingDegrees; kill cemdsResult; kill sModifyResult; kill RModified; kill R; //Fifth example: Verification of the resulting CEMDS, fake variables are marked, points are transferred to the modified CEMDS and its fan is computed. ring R = 0,T(1..3),dp; //The CEMDS's rays' matrix intmat rvcvzP[2][3] = 1,0,-1, 0,1,-1; //The CEMDS's fan intmat cvrvz1[2][2] = 1,0, 0,1; intmat cvrvz2[2][2] = 0,1, -1,-1; intmat cvrvz3[2][2] = -1,-1, 1,0; cone cn1 = coneViaPoints(cvrvz1); cone cn2 = coneViaPoints(cvrvz2); cone cn3 = coneViaPoints(cvrvz3); fan scnSigma = fanViaCones(cn1, cn2, cn3); //The CEMDS's embedding ideal ideal spG = 0; CEMDS cemds = createCEMDS(rvcvzP, scnSigma, spG); intvec rvnVanishingDegrees = 0,1,1; list sModifyResult = modifyCEMDS2(cemds, rvnVanishingDegrees, 1, 0:3, list(), 1); print(sModifyResult); //Should yield a list with two entries: // - The CEMDS containing: // * The matrix with the entries // (1 0 -1 -1) // (0 1 -1 0) // * The fan consisting of the maximal cones generated by the following rays: // . (1,0), (0,1); // . (0,1), (-1,0); // . (-1,0, (-1,-1); // . (-1,-1), (1,0); // * The zero ideal // - the verification result (1) } static proc modifyCEMDSCore(intmat rvcvzPDash, fan scnSigmaDash, ideal spG, intmat rvcvnVanishingDegrees, int fVerify, intvec rvfFakeVars, list #) "USAGE: modifyCEMDSCore(rvcvzPDash, scnSigmaDash, spG, rvcvnVanishingDegrees, fVerify, rvfFakeVars); modifyCEMDSCore(rvcvzPDash, scnSigmaDash, spG, rvcvnVanishingDegrees, fVerify, rvfFakeVars, svnumPts); ASSUME: All input data was already checked for compliance with the requirements of this function. PURPOSE: Main functionality of modifyCEMDS and modifyCEMDS2. RETURN: Same as modifyCEMDS and modifyCEMDS2. EXAMPLE: none." { def R = basering; list sParams = #; int nParams = size(#); //First parameter: A list of points whose coordinates shall be transferred to the modified CEMDS. if (nParams > 0) { //Type and size check for the list of points to modify along with the CEMDS if (typeof(sParams[1]) == "vector") { list svnumPts = sParams; nParams = 1; } else { if (typeof(sParams[1]) == "list") { list svnumPts = sParams[1]; } else { ERROR("\"modifyCEMDSCore\": Something went terribly wrong in passing the points to modify along with the CEMDS."); } } int ivnumPts; int nvnumPts = size(svnumPts); } int imonVars; int nmonVars = nvars(R); //The modification of the CEMDS by a refinement of its fan is //represented by the morphism of fans of ambient toric varieties given by this matrix intmat matzRefinement = unitMatrix(nmonVars); matzRefinement = intmatAppendCols(matzRefinement, rvcvnVanishingDegrees); dbprint(printlevel - 1, ""); dbprint(printlevel - 1, "Matrix refining the CEMDS's fan:"); dbprint(printlevel - 1, matzRefinement); //The blown up CEMDS's ideal is given by the saturation w.r.t var(1)*...*var(nvars) //of the inverse image of the original CEMDS's ideal under the morphism just defined. def RModified = toricMorphismPullback(matzRefinement, spG); setring RModified; ideal spGDash = gspResult; //Order the variables to enhance the saturation computation int imonNewVars; int nmonNewVars = nvars(RModified) - nmonVars; list simonLexVars = list(); intvec vfSatVars = 0:nvars(RModified); for (imonNewVars = 1; imonNewVars <= nmonNewVars; imonNewVars++) { //simonLexVars = insert(simonLexVars, nmonVars + imonNewVars); vfSatVars[nmonVars + imonNewVars] = 1; } for (imonVars = 1; imonVars <= nmonVars; imonVars++) { if (rvfFakeVars[imonVars]) { simonLexVars = insert(simonLexVars, imonVars); } } spGDash = xIdealSaturation(spGDash, simonLexVars, vfSatVars); spGDash = simplify(spGDash, 2); //Compose and return the result list gsResult = list(); //The modified points setring R; if (defined(svnumPts)) { if (nvnumPts > 0) { setring RModified; list svnumFetchedPts = fetch(R, svnumPts); list svnumPtsModified = list(); vector vnumFetchedPts; for (ivnumPts = nvnumPts; ivnumPts >= 1; ivnumPts--) { vnumFetchedPts = svnumFetchedPts[ivnumPts]; for (imonNewVars = 1; imonNewVars <= nmonNewVars; imonNewVars++) { vnumFetchedPts = vnumFetchedPts + gen(nmonVars + imonNewVars); } svnumPtsModified = insert(svnumPtsModified, vnumFetchedPts); } gsResult = insert(gsResult, svnumPtsModified); dbprint(printlevel - 1, ""); dbprint(printlevel - 1, "Modified coordinates of the points passed:"); dbprint(printlevel - 1, svnumPtsModified); } } setring RModified; //Whether the verification was positive if (fVerify) { dbprint(printlevel, ""); dbprint(printlevel, "----- Starting verification, please be patient... -----"); if (printlevel >= 1) { int nStartTime = timer; } ideal spGDashStd = std(spGDash); int fVerified = 1; //Dimension test fVerified = verifyDimension(spGDashStd); if (fVerified) { //Primality test for the variables' classes in the CEMDS's ring fVerified = verifyVarPrimality(spGDashStd); if (!fVerified) { dbprint(printlevel - 1, "Primality test on the ring's variables failed."); } } else { dbprint(printlevel - 1, "Dimension test failed."); } dbprint(printlevel + 1, "WARNING: The verification tests do not include the test whether the CEMDS's ring is normal. Please check that manually."); if (printlevel >= 1) { int nTimeElapsed = timer - nStartTime; print("Time elapsed while verifying (in ms):"); print(nTimeElapsed); } dbprint(printlevel, "----- Verification done! -----"); gsResult = insert(gsResult, fVerified); } //The modified CEMDS CEMDS cemdsResult = createCEMDS(rvcvzPDash, scnSigmaDash, spGDash); if (size(gsResult) > 0) { list scemdsResult = cemdsResult; gsResult = scemdsResult + gsResult; if (defined(svnumFetchedPts)) { export(gsResult); return(RModified); } else { return(gsResult); } } else { return(cemdsResult); } } proc stretchCEMDS(CEMDS cemds, list shpQPrimes, list #) "USAGE: stretchCEMDS(cemds, shpQPrimes); with cemds: CEMDS, shpQPrimes: list of poly; stretchCEMDS(cemds, shpQPrimes, svnumPts); with cemds: CEMDS, shpQPrimes: list of poly, svnumPts: list of vectors of numbers; stretchCEMDS(cemds, shpQPrimes, svnumPts, fComputeFan); with cemds: CEMDS, shpQPrimes: list of poly, svnumPts: list of vectors of numbers, fComputeFan: int; stretchCEMDS(cemds, shpQPrimes, svnumPts, fComputeFan, vzWeight); with cemds: CEMDS, shpQPrimes: list of poly, svnumPts: list of vectors of numbers, fComputeFan: int, vzWeight: intvec. ASSUME: (i) The passed list of polynomials shpQPrimes is homogeneous w.r.t. the grading given by rvcvzQ = gale(cemds.rvcvzP); (ii) all points to stretch along with the CEMDS, if there were any passed, are contained in the CEMDS to stretch; (iii) the intvec \"vzWeight\", if passed, is element of the ample cone defined by the grading matrix rvcvzQ extended by the passed polynomials' degrees and some true saturated connected collection of facets of the positive orthant. PURPOSE: Embeds a CEMDS in a new ambient toric variety such that each passed polynomial is set in relation to a new variable. Optionally transfers points from the unstretched CEMDS to the stretched one. RETURN: The stretched CEMDS only or, iff a non-empty list of points to stretch along with the CEMDS was passed, the basering the stretched CEMDS is defined in. SIDE EFFECTS: Iff a non-empty list of points to stretch along with the CEMDS were passed, a list gsResult containing the following (in order of appearance): (i) the stretched CEMDS with or without a fan (depending on fComputeFan and vzWeight) (always present), (ii) the list of stretched points (present as there was a non-empty list list passed). NOTE: By default: fComputeFan = 0; vzWeight is an element of the relative interior of the moving cone generatwd by the columns of the grading matrix rvcvzQ = gale(cemds.rvcvzP) extended by the passed polynomials' degrees. EXAMPLE: example stretchCEMDS; shows an example." { //The CEMDS to stretch must be associated with a ring def R = cemds.R; if (!defined(R)) { ERROR("The passed CEMDS is not associated with any ring."); } //Switch to the CEMDS's ring setring R; int imonVars; int nmonVars = nvars(R); //Search for optional parameters list sParams = #; int nParams = size(sParams); //Default parameters //list svnumPts: intentionally not defined here. int fComputeFan = 0; //intvec vzWeight: intentionally not defined here. //First parameter: A list of points whose coordinates shall be transferred to the stretched CEMDS. if (nParams > 0) { //Type and size check for the list of points to stretch along with the CEMDS if (typeof(sParams[1]) == "vector") { list svnumPts = sParams; nParams = 1; } else { if (typeof(sParams[1]) == "list") { list svnumPts = sParams[1]; } else { ERROR("The parameter passed for \"svnumPts\" must be of type \"list\" and contain vectors of numbers only."); } } int ivnumPts; int nvnumPts = size(svnumPts); for (ivnumPts = 1; ivnumPts <= nvnumPts; ivnumPts++) { if (typeof(svnumPts[ivnumPts]) != "vector") { ERROR("The passed points to modify along with the CEMDS must be of type \"vector\" and consist of numbers only."); } for (imonVars = 1; imonVars <= nmonVars; imonVars++) { if (leadexp(svnumPts[ivnumPts][imonVars]) != 0:nmonVars) { ERROR("The vectors passed must consist of numbers only."); } } } //Second parameter: Whether to compute a fan or not. if (nParams > 1) { //Type check for fan computation flag if (typeof(sParams[2]) != "int") { ERROR("The parameter passed for \"fComputeFan\" must be of type \"int\" and equal \"0\" or \"1\"."); } fComputeFan = sParams[2]; //Third parameter: An optional weight for initializing the fan computation. if (nParams > 2) { //Type check for weight if (typeof(sParams[3]) != "intvec") { ERROR("The parameter passed for \"vzWeight\" must be of type \"intvec\"."); } intvec vzWeight = sParams[3]; //There are at most 3 optional parameters if (nParams > 3) { ERROR("Too many optional parameters."); } } } } //Check whether there is something to do int ihpQPrimes; int nhpQPrimes = size(shpQPrimes); if (nhpQPrimes == 0) { if (defined(svnumPts)) { if (nvnumPts > 0) { list gsResult = list(); gsResult = insert(gsResult, svnumPts); list scemdsResult = cemds; gsResult = scemdsResult + gsResult; export(gsResult); return(R); } else { return(cemds); } } else { return(cemds); } } //Clear the the passed polynomials' denominators list spStretchedPolys = list(); for (ihpQPrimes = nhpQPrimes; ihpQPrimes >= 1; ihpQPrimes--) { spStretchedPolys = insert(spStretchedPolys, cleardenom(shpQPrimes[ihpQPrimes])); } //Generate the grading matrix intmat rvcvzQ = gale(cemds.rvcvzP); //Calculate the homogeneous polynomials' degrees and add them to the grading matrix intmat rvcvzQDash = rvcvzQ; for (ihpQPrimes = 1; ihpQPrimes <= nhpQPrimes; ihpQPrimes++) { rvcvzQDash = intmatAppendCol(rvcvzQDash, getDegree(rvcvzQ, spStretchedPolys[ihpQPrimes])); } //Extend the original rays' matrix so it is gale dual to the new grading matrix intmat rvcvzPDash = galeExtension(rvcvzQDash, cemds.rvcvzP); //Calculate the new CEMDS's fan, if wanted fan scnSigmaDash = emptyFan(nrows(rvcvzPDash)); if (fComputeFan) { //If no weight was passed, find one in the new grading matrix's moving cone's relative interior. if (!defined(vzWeight)) { cone cnMov = movingConeFromWeights(rvcvzQDash); intvec vzWeight = intvec(relativeInteriorPoint(cnMov)); } fan sof = orthantFanFromWeight(rvcvzQDash, vzWeight); scnSigmaDash = mapFan(rvcvzPDash, sof); } //Generate new relations for the ideal describing how the CEMDS is embedded in its ambient toric variety. //First, switch to an appropriate ring. ideal spG = cemds.spG; execute("ring RExtended = (" + charstr(R) + "),T(1..(nmonVars + nhpQPrimes)),dp;"); ideal spResult = fetch(R, spG); list spStretchedFetchedPolys = fetch(R, spStretchedPolys); //Add the new relations to the ideal for (ihpQPrimes = 1; ihpQPrimes <= nhpQPrimes; ihpQPrimes++) { spResult = spResult, var(nmonVars + ihpQPrimes) - spStretchedFetchedPolys[ihpQPrimes]; } spResult = simplify(spResult, 2); CEMDS cemdsResult = createCEMDS(rvcvzPDash, scnSigmaDash, spResult); //Compose and return the result list gsResult = list(); setring R; //The points to be stretched along with the CEMDS if (defined(svnumPts)) { if (nvnumPts > 0) { number numEval; setring RExtended; number numEvalFetched; list svnumFetchedPts = fetch(R, svnumPts); list svnumPtsStretched = list(); vector vnumFetchedPts; for (ivnumPts = nvnumPts; ivnumPts >= 1; ivnumPts--) { vnumFetchedPts = svnumFetchedPts[ivnumPts]; for (ihpQPrimes = 1; ihpQPrimes <= nhpQPrimes; ihpQPrimes++) { setring R; numEval = evaluatePoly(spStretchedPolys[ihpQPrimes], svnumPts[ivnumPts]); setring RExtended; numEvalFetched = fetch(R, numEval); vnumFetchedPts = vnumFetchedPts + numEvalFetched * gen(nmonVars + ihpQPrimes); } svnumPtsStretched = insert(svnumPtsStretched, vnumFetchedPts); } gsResult = insert(gsResult, svnumPtsStretched); dbprint(printlevel - 1, ""); dbprint(printlevel - 1, "Stretched coordinates of the points passed:"); dbprint(printlevel - 1, svnumPtsStretched); } } setring RExtended; if (size(gsResult) > 0) { list scemdsResult = cemdsResult; gsResult = scemdsResult + gsResult; export(gsResult); return(RExtended); } else { return(cemdsResult); } } example { "EXAMPLE"; echo=2; //First example: Parameterized basering with no optional arguments ring R = (0,a),T(1..6),dp; //The CEMDS's rays' matrix intmat rvcvzP[2][6] = ( 1,0,-1,1,0,-1, 0,1,-1,1,-1,0); //The CEMDS's fan intmat cvrvz1[2][2] = 1,0, 1,1; intmat cvrvz2[2][2] = 1,1, 0,1; intmat cvrvz3[2][2] = 0,1, -1,0; intmat cvrvz4[2][2] = -1,0, -1,-1; intmat cvrvz5[2][2] = -1,-1, 0,-1; intmat cvrvz6[2][2] = 0,-1, 1,0; cone cn1 = coneViaPoints(cvrvz1); cone cn2 = coneViaPoints(cvrvz2); cone cn3 = coneViaPoints(cvrvz3); cone cn4 = coneViaPoints(cvrvz4); cone cn5 = coneViaPoints(cvrvz5); cone cn6 = coneViaPoints(cvrvz6); fan scnSigma = fanViaCones(cn1, cn2, cn3, cn4, cn5, cn6); //The CEMDS's embedding ideal ideal spG = 0; CEMDS cemds = createCEMDS(rvcvzP, scnSigma, spG); //The orbit ideal generators of (1,a,1,1,1,1); poly p1 = 1/a*T(2)*T(4) - T(3)*T(5); poly p2 = T(1)*T(4) - T(3)*T(6); poly p3 = T(1)*T(5) - 1/a*T(2)*T(6); list shpQPrimes = p1,p2,p3; CEMDS cemdsResult = stretchCEMDS(cemds, shpQPrimes); //Obtains the CEMDS's ring def RResult = cemdsResult.R; setring RResult; print(cemdsResult); //Should yield // - The rays' matrix // ( 1 0 -1 1 0 -1 0 0 0) // ( 0 1 -1 1 -1 0 0 0 0) // ( 0 -1 0 -1 0 0 1 0 0) // (-1 0 0 -1 0 0 0 1 0) // (-1 0 0 0 -1 0 0 0 1) // or some other matrix whose first two rows do not differ and whose further rows each linearly depend on the rows of this matrix // - An empty fan of ambient dimension 5. // - The ideal generated by // T(7) - T(2)*T(4) + a*T(3)*T(5), // T(8) - T(1)*T(4) + T(3)*T(6), // T(9) - a*T(1)*T(5) + T(2)*T(6). //Second example with one optional argument: Points whose coordinates shall be stretched along with the CEMDS. setring R; vector vnumPt1 = [1, 1, 1, 1, 1, 1]; vector vnumPt2 = [0, 1, 2, 3, 4, 5]; list svnumPts = vnumPt1, vnumPt2; def RResult2 = stretchCEMDS(cemds, shpQPrimes, svnumPts); setring RResult2; list sStretchResult2 = gsResult; print(sStretchResult2); //Should yield a list with two entries: // - The same CEMDS as in the first example. // - A list of two vectors: [1,1,1,1,1,1,-a+1,0,a-1], [0,1,2,3,4,5,-8a+3,-10,-5]. //Third example with two optional arguments: Fan computation enabled. setring R; CEMDS cemdsResult3 = stretchCEMDS(cemds, list(), list(), 1); print(cemdsResult3); //Should yield the original CEMDS. //Fourth example with three optional arguments: Points whose coordinates shall be stretched along with the CEMDS, fan computation enabled and starting with a user defined weight. setring R; intvec vzWeight = 3,-1,2,2; def RResult4 = stretchCEMDS(cemds, shpQPrimes, svnumPts, 1, vzWeight); setring RResult4; list sStretchResult4 = gsResult; print(sStretchResult4); //Should yield a list with two entries: // - The same CEMDS as in the first example, but with a fan computed // - A list of two vectors: [1,1,1,1,1,1,-a+1,0,a-1], [0,1,2,3,4,5,-8a+3,-10,-5]. } proc verifyAlmostFree(intmat rvcvzP) "USAGE: verifyAlmostFree(rvcvzP); with rvcvzP: intmat. PURPOSE: Checks whether the grading defined by a gale dual matrix to the passed matrix is an almost free grading. RETURN: The int 1 if the grading defined by a gale dual matrix to the passed matrix is an almost free grading, the int 0 else. EXAMPLE: example verifyAlmostFree; shows an example." { //Test for being almost free int izP; int icvzP2; int nzP = nrows(rvcvzP); int ncvzP = ncols(rvcvzP); int nGcd; intvec cvzP; intvec cvzP2; int fVerified = 1; int icvzP = 1; while (fVerified && icvzP <= ncvzP) { cvzP = intmatTakeCol(rvcvzP, icvzP); icvzP2 = icvzP + 1; //Each pair of columns of P must differ while (fVerified && icvzP2 <= ncvzP) { cvzP2 = intmatTakeCol(rvcvzP, icvzP2); fVerified = (cvzP != cvzP2); icvzP2++; } //Each column has to be primitive nGcd = rvcvzP[1, icvzP]; for (izP = 2; izP <= nzP; izP++) { nGcd = gcd(nGcd, rvcvzP[izP, icvzP]); } fVerified = (nGcd == 1); icvzP++; } return(fVerified); } example { "EXAMPLE:"; echo=2; intmat rvcvzP = intmat(intvec( -1,0,1,1,0,-1,0,0,0,0, 1,-1,0,-1,1,0,0,0,0,0, 1,-1,0,0,0,0,1,-1,0,0, -1,0,1,0,0,0,-1,0,1,0, 1,-1,-1,-1,0,0,1,0,0,1 ), 5, 10); int fVerified = verifyAlmostFree(rvcvzP); print(fVerified); //Should be true. } proc verifyDimension(ideal spG) "USAGE: verifyDimension(spG); with spG: ideal. PURPOSE: Checks whether dim(spG) - dim(spG + + ) >= 2 for all i != j. RETURN: The int 1 if dim(spG) - dim(spG + + ) >= 2 for all i != j, the int 0 else. EXAMPLE: example verifyDimension; shows an example." { int imonVars; int nmonVars = nvars(basering); ideal spGStd = std(spG); int nDimG = dim(spGStd); ideal spGStdVars; int imonVars2; int fVerified = 1; imonVars = 1; while (fVerified && imonVars <= nmonVars) { imonVars2 = imonVars + 1; while (fVerified && imonVars2 <= nmonVars) { spGStdVars = spGStd, var(imonVars), var(imonVars2); spGStdVars = std(spGStdVars); fVerified = fVerified && (nDimG - dim(spGStdVars) >= 2); imonVars2++; } imonVars++; } return(fVerified); } example { "EXAMPLE:"; echo=2; ring R = 0,T(1..10),dp; ideal spG = T(4)*T(7)-T(5)*T(8)+T(6)*T(9), T(1)*T(7)-T(2)*T(8)+T(3)*T(9), T(3)*T(5)-T(2)*T(6)-T(7)*T(10), T(3)*T(4)-T(1)*T(6)-T(8)*T(10), T(2)*T(4)-T(1)*T(5)-T(9)*T(10); int fVerified = verifyDimension(spG); print(fVerified); //Should be true. } proc verifyVarPrimality(ideal spG) "USAGE: verifyVarPrimality(spG); with spG: ideal. PURPOSE: Checks whether each ring variable's class in the basering modulo spG defines a prime ideal. RETURN: The int 1 if each ring variable's class in the basering modulo spG defines a prime ideal, the int 0 else. EXAMPLE: example verifyVarPrimality; shows an example." { int imonVars; int nmonVars = nvars(basering); //Primality test for the variables' classes in the CEMDS's ring int fVerified = 1; def RBase = basering; ideal spVar; int iPrimeCompGen; int nPrimeCompGen; imonVars = 1; while (fVerified && imonVars <= nmonVars) { spVar = spG, var(imonVars); def RPrimes = absPrimdecGTZ(spVar); setring RPrimes; list sAbsPrimes = absolute_primes; list sPrimDecomp = primary_decomp; //More than one irreducible component --> ideal non-prime. if (size(sAbsPrimes) != 1) { fVerified = 0; break; } //More than one irreducible component (by conjugates) --> ideal non-prime. if (sAbsPrimes[1][2] != 1) { fVerified = 0; break; } //Is the primary component prime, i.e., is its prime radical in [1][2] subset of [1][1]? ideal spPrimary = std(primary_decomp[1][1]); nPrimeCompGen = size(primary_decomp[1][2]); for (iPrimeCompGen = 1; iPrimeCompGen <= nPrimeCompGen; iPrimeCompGen++) { if (reduce(primary_decomp[1][2][iPrimeCompGen], spPrimary) != 0) { fVerified = 0; break; } } setring RBase; kill RPrimes; imonVars++; } return(fVerified); } example { "EXAMPLE:"; echo=2; ring R = 0,T(1..10),dp; ideal spG = T(4)*T(7)-T(5)*T(8)+T(6)*T(9), T(1)*T(7)-T(2)*T(8)+T(3)*T(9), T(3)*T(5)-T(2)*T(6)-T(7)*T(10), T(3)*T(4)-T(1)*T(6)-T(8)*T(10), T(2)*T(4)-T(1)*T(5)-T(9)*T(10); int fVerified = verifyVarPrimality(spG); print(fVerified); //Should be true. } static proc vanishingDegrees(list spI, vector vnum) "USAGE: vanishingDegrees(spI, vnum); with spI: list of polynomials, vnum: vector. PURPOSE: Computes a vector pointing out whether the point to blow up vanishes under existing variables and variables to be defined later. RETURN: an intvec pointing out whether the point to blow up vanishes under existing variables and variables to be defined later." { int imonVars; int nmonVars = nvars(basering); int ipI; int npI = size(spI); int nvnResult = npI + nmonVars; intvec rvnResult = 0:nvnResult; //For the first nmonVars variables, simply determine whether the passed point is vanishing or not. for (imonVars = 1; imonVars <= nmonVars; imonVars++) { if (vnum[imonVars] == 0) { rvnResult[imonVars] = 1; } } //For the new variables, simply determine whether the passed point is vanishing or not. for (ipI = 1; ipI <= npI; ipI++) { rvnResult[nmonVars + ipI] = 1; } return(rvnResult); }