////////////////////////////////////////////////////////////////////////////// version="version arr.lib 4.1.1.0 Dec_2017 "; // $Id$ category="Miscellaneous"; /* ** TOPICS ** (Ctrl+f to search) #1 NEWSTRUCTS & OVERLOADS #2 CONSTRUCTORS #3 ACCESS & ASSIGNEMENT #4 DELETION #5 COMPERATORS #6 TYPE CASTING #7 INHERITED FUNCTIONS #8 PRINTING #9 MANIPULATING VARIABLES #10 CENTER COMPUTATIONS #11 GEOMETRIC CONSTRUCTIONS #11 EXAMPLES OF ARRANGEMENTS #13 ORLIK SOLOMON AND POINCARE POLYNOMIAL #14 FREENESS #15 MULTI-ARRANGEMENTS #16 COMBINATORICS */ info=" LIBRARY: arr.lib a library of algorithms for arrangements of hyperplanes AUTHORS: Randolf Scholz (rscholz@rhrk.uni-kl.de), Patrick Serwene (serwene@mathematik.uni-kl.de), Lukas Kuehne (lf.kuehne@gmail.com) OVERLOADS: // OPERATORS = arrAdd assignment + arrAdd union of two arrs [ arrGet access to a single/multiple hyperplane(s) - arrMinus deletes given hyperplanes from the arr <= arrLEQ comparison >= arrGEQ comparison == arrEQ comparison != arrNEQ comparison < arrLNEQ comparison > arrGNEQ comparison // TYPECASTING matrix arr2mat coeff matrix poly arr2poly defining polynomial // OTHER variables arrVariables ideal generated by the variables the arr depends on nvars arrNvars number of variables the arr depends on delete arrDelete deletes hyperplanes by indices print arrPrint prints the arr on the screen // IDEAL INHERITED FUNCTIONS homog arrHomog checks if arrangement is homogeneous simplify arrSimplify simplifies arrangement size arrSize number of planes subst arrSubst substitute variables // MULTI-ARRANGEMENTS = multarrAdd assignement of multarr + multarrAdd union of multarr poly multarr2poly defining polynomial size multarrSize number of hyperplanes with mult. print multarrPrint displays multiarr delete multarrDelete deletes hyperplane PROCEDURES: arrSet(arr A, int k, poly p) replaces the k-th Hyperplane with poly p type2arr(#) converts general input to 'arr' using arrAdd. mat2arr( matrix M) affine arrangement from coeff matrix mat2carr(matrix M) central arrangement from coeff matrix arrPrintMatrix(arr A) readable output as a coeff matrix varMat(intvec v) matrix of the corresponding ring_variables varNum(def u) number of given variable (enh. version of varNum in dmod.lib) arrSwapVar(arr A, i, j) swaps two variables in the arrangement arrLastVar(arr A) ring_variable of largest index used in arrangement arrCenter(arr A) computes center of an arrangement arrCentral(arr A) checks if arrangement is central arrCentered(arr A) checks if arrangement is centered arrCentralize(arr A) makes centered arrangement central arrCoordChange(A, T, #) performs coordinate change arrCoordNormalize(A, v) performs projection onto coordinate hyperplane arrCone(arr A, var) coned arrangement arrDecone(arr A, int k) deconed arrangement arrLocalize(arr A, intvec v) localization of an arrangement onto a flat arrRestrict(arr A, intvec v) restricted arrangement onto a flat arrIsEssential(arr A) checks if arrangement is essential arrEssentialize(arr A) essentialized arragnement arrBoolean(int v) boolean arrangement arrBraid(int v) braid arrangement arrTypeB(int v) type B arrangement arrTypeD(int v) type D arrangement arrRandom(d,m,n) random (affine) arrangement arrRandomCentral(d,m,n) random central arrangement arrEdelmanReiner() Edelman-Reiner arrangement arrOrlikSolomon(arr A) Orlik-Solomon algebra of the arrangement arrDer(A) module of derivation arrIsFree(A) checks if arrangement is free arrExponents(A) exponents of a (free) arrangement arr2multarr(arr A, intvec v) converts normal arrangement to multiarrangement multarr2arr(multarr A) converts multiarrangement to normal arrangement multarrRestrict(arr A, v) restriction of A (as arr) to a flat with multiplicities multarrMultRestrict(A, int k) restriction of A (as multarr) to a hyperplane with multiplicities arrFlats(arr A) intersection lattice arrLattice(arr A) computes the intersection lattice / poset moebius(arrposet P) computes moebius values arrCharPoly(arr A) characteristic polynomial arrPoincare(arr A) poincare polynomial of the arrangement arrChambers(arr A) number of chambers of the arrangement arrBoundedChambers(arr A) number of bounded chambers of the arrangement printMoebius(arr A) displays the moebius values of all the flats in the poset "; //============================================================================// //-------------------------- #1 NEWSTRUCTS & OVERLOADS -----------------------// //============================================================================// // initialization of the library static proc mod_init() { // NEWSTRUCTS newstruct("arr","ideal l"); newstruct("flat", "intvec REL, int moebius, intvec parents, int flag"); newstruct("arrposet","arr A, list r"); newstruct("arrflats","arr A, list r"); newstruct("multarr","ideal l, intvec m"); // intvec: multiplicities of hyperplanes // OPERATORS system("install","arr","=" ,arrAdd ,1); // assignment system("install","arr","+" ,arrAdd ,2); // union of arrs system("install","arr","[" ,arrGet ,2); // access system("install","arr","-" ,arrMinus ,2); // delete plane system("install","arr","<=" ,arrLEQ ,2); // comparison system("install","arr",">=" ,arrGEQ ,2); // comparison system("install","arr","==" ,arrEQ ,2); // comparison system("install","arr","!=" ,arrNEQ ,2); // comparison system("install","arr","<" ,arrLNEQ ,2); // comparison system("install","arr",">" ,arrGNEQ ,2); // comparison // TYPECASTING system("install","arr","matrix" ,arr2mat ,1); // coeff matrix system("install","arr","poly" ,arr2poly ,1); // defining polynomial system("install","arr","list" ,arr2list ,1); // list of defining polynomials system("install","arr","ideal" ,arr2ideal ,1); // list of defining polynomials // OTHER system("install","arr","variables" ,arrVariables ,1); system("install","arr","nvars" ,arrNvars ,1); system("install","arr","delete" ,arrDelete ,2); system("install","arr","print" ,arrPrint ,1); // IDEAL INHERITED FUNCTIONS system("install","arr","homog" ,arrHomog ,1); // checks if homogeneous system("install","arr","homog" ,arrHomog ,2); // checks if homogeneous system("install","arr","simplify" ,arrSimplify ,2); // simplifies arrangement system("install","arr","size" ,arrSize ,1); // number of planes system("install","arr","subst" ,arrSubst ,4); // substitute variables // MULTI-ARRANGEMENTS system("install","multarr","=" ,multarrAdd ,1); // assignement of multarr system("install","multarr","+" ,multarrAdd ,2); // union of multarr system("install","multarr","poly" ,multarr2poly ,1); // defining polynomial system("install","multarr","size" ,multarrSize ,1); // number of hyperplanes with mult. system("install","multarr","print" ,multarrPrint ,1); // displays multiarr system("install","multarr","delete" ,multarrDelete ,2); // deletes hyperplane // COMBINATORICS system("install","arr","rank" ,arrRank ,1); system("install","arrflats","print" ,arrPrintFlats ,1); system("install","arrposet","print" ,arrPrintPoset ,1); // NEEDED LIBRARIES LIB "general.lib"; LIB "monomialideal.lib"; } //============================================================================// //--------------------------- #2 CONSTRUCTORS --------------------------------// //============================================================================// // general method for creating arrangements static proc arrAdd "USAGE: A = #; A +#; # list containing arr/ideal/list/matrix/poly RETURN: [arr] Arrangement constructed by input parameters. REMARKS: The algorithm splits up the list # and uses the appropiate procedure to handle the input. NOTE: arrAdd simplifies certain inputs, for example A = (x,y,2x); gives the same arrangement as A = (x,y); SEE ALSO: arrAdd, type2arr KEYWORDS: arrangement; equal; constructor; operator EXAMPLE: example arrAdd; shows an example" { arr A; for(int k=1; k<=size(#); k++){ while(1){ //simulates switch, which singular doesn't offer if(typeof(#[k]) == "arr" ){A = arrAddArr (A, #[k]);break;} if(typeof(#[k]) == "poly" ){A = arrAddPoly (A, #[k]);break;} if(typeof(#[k]) == "ideal" ){A = arrAddIdeal (A, #[k]);break;} if(typeof(#[k]) == "matrix"){A = arrAddMatrix(A, #[k]);break;} if(typeof(#[k]) == "intmat"){A = arrAddMatrix(A, #[k]);break;} if(typeof(#[k]) == "list" ){A = arrAdd ( #[k]);break;} ERROR("bad input type"); } } return (A); } example { "EXAMPLE: Creating a few arrangements"; echo = 2; ring R = 0,(x,y,z),dp; arr A= ideal(x,y,z); A; arr B = A + ideal(x+1, x-1); B; arr C = list(A, x+1, x-1); C; arr D = x2 - y2; D; } // union of two arrangements static proc arrAddArr(arr A, arr B){ return (arrAddIdeal(A, B.l)); } // adds a single poly to the arrangement // if the poly is linear, it is just added, otherwise Singular factorizes static proc arrAddPoly(arr A, poly p){ if(deg(p) == 0){ ERROR("Given poly is not linear or Singluar is not able to factorize it");} else{ if(deg(p) == 1){ A.l = A.l + p; return (A); } else{ ideal I = factorize(p,1); if(size(I) == 1){ERROR("Given poly is not a hyperplane");} else{ return (arrAdd(A,I)); } } } return(A); } // adds defining polys to the arrangement static proc arrAddIdeal(arr A, ideal I){ for(int k=1; k<=size(I); k++){ A = arrAddPoly(A,I[k]); } return (A); } // adds defining polys to the arrangement static proc arrAddMatrix(arr A, matrix M){ return ( arrAddArr(A,mat2arr(M)) ); } //============================================================================// //--------------------------- #3 ACCESS & ASSIGNEMENT ------------------------// //============================================================================// // access to hyperplanes, overloads [] operator static proc arrGet(arr A, intvec v) "USAGE: A[v]; v int/intvec RETURN: [poly] if v is [int] The defining poly of the the v-th hyperplane+ [arr] if v is [intvec] The corresponing subarrangement SEE ALSO: arrGet, arrSet KEYWORDS: arrangement; get; operator EXAMPLE: example arrGet; shows an example" { if(size(v) == 1){ return (A.l[v[1]]); } //returns poly if v is integer ideal I = A.l; A = ideal(I[v]); return (A); } example { "EXAMPLE: "; echo = 2; ring R = 0,(x,y,z),dp; arr A = ideal(x,y,z); A[2]; intvec v = 1,3; A[v]; } // replaces the k-th plane with poly p proc arrSet(arr A, int k, poly p) "USAGE: arrSet(A, k, p); arr A, int k, poly p; RETURN: [arr] Arrangement where the k-th hyperplane is replaced by p. NOTE: p must be linear KEYWORDS: arrangement; hyperplane; assign; set EXAMPLE: example arrSet; shows an example" { if(deg(p) != 1){ERROR("Given poly is not a hyperplane");} else{ // looks akward but needs to be done this way ideal I = A.l; I[k] = p; A = I; } return (A); } example { "EXAMPLE: "; echo = 2; ring R = 0,(x,y,z),dp; arr A = ideal(x,y,z); arrSet(A,1,x+1); } //============================================================================// //------------------------------- #4 DELETION --------------------------------// //============================================================================// // deletes all hyperplanes of the given indices static proc arrDelete(arr A, intvec v) "USAGE: delete(A, v); v integer/intvec RETURN: [arr] Arrangement A without the hyperplanes given by v; NOTE: for deleting hyperplanes via polynomials, use arrMinus instead SEE ALSO: arrDelete, arrMinus KEYWORDS: arrangement; delete EXAMPLE: example arrDelete; shows an example" { int i = 0; int k; int n = size(A); intvec u = 1..n; // puts 0 in u for every element that needs to be deleted // is done this way to deal with the case that the user gives the same index multiple times. for(k=1; k<=size(v); k++){ if(u[v[k]] != 0){ u[v[k]] = 0; i++; } // i = #elts that need to be deleted } if( i == n){return (emptyArr); } // create intvec of the remaining indices v = 1..(n-i); i=1; for(k=1; k<=n; k++){ if(u[k] != 0) { v[i] = u[k]; i++; } } arr A' = arrGet(A,v); return (A'); } example { "EXAMPLE: "; echo = 2; ring R = 0,(x,y,z),dp; arr A = ideal(x,y,z); delete(A,2); intvec v = (1,3); delete(A,v); } // deletes hyperplanes, overloads - operator static proc arrMinus "USAGE: A - #; # list containing arr/ideal/list/matrix/poly RETURN: [arr] arrangement A without the hyperplanes of the arrangement defined by #. REMARKS: algorithm creates an arrangement B from # using arrAdd and then deletes hyperplanes which occur in both A and B. NOTE: The alorithm does not simplify by scalars, i.e. some hyperplanes might not be deleted. See example. SEE ALSO: arrDelete, arrMinus KEYWORDS: arrangement; delete; operator EXAMPLE: example arrMinus; shows an example" { if(typeof(#[1]) != "arr"){ERROR("First input must be arr!");} arr A = #[1]; arr B = #[2..size(#)]; // collects hyperplanes to be deleted list L; int k, l; for(k=1; k<=size(A); k++){ // create list of hyperplanes to be deleted for(l=1; l<=size(B); l++){ if(A[k] == B[l]){L = insert(L,k);} } } l = size(L); // transforms list to intvec if(l != 0){ intvec v = 1..l; for(k=1; k<=l; k++){v[k] = L[k];} A = delete(A,v); } return (A); } example { "EXAMPLE: "; echo = 2; ring R = 0,(x,y,z),dp; arr A = ideal(x,y,z); A - ideal(x,y); A - poly(2x); } //============================================================================// //------------------------------- #5 COMPERATORS -----------------------------// //============================================================================// // returns logical value of 'A<=B' static proc arrLEQ(arr A, arr B) "USAGE: A<=B; A,B arr RETURN: [0,1] true if A is a subarrangement of B NOTE: algorithm is based on arrMinus and does not simplify by scalars, hence some technically equal hyperplanes might not be detected. See example. SEE ALSO: arrLEQ, arrLNEQ, arrGEQ, arrGNEQ, arrEQ, arrNEQ KEYWORDS: comparison; subarrangement; operator EXAMPLE: example arrLEQ; shows an example" { arr C = A - B; if(C[1] == 0){ return (1); } return (0); } example{example arrEQ;} // returns logical value of 'A>=B' static proc arrGEQ(arr A, arr B) "USAGE: A>=B; A,B arr RETURN: [0,1] true if B is a subarrangement of A NOTE: algorithm is based on arrMinus and does not simplify by scalars, hence some technically equal hyperplanes might not be detected. See example. SEE ALSO: arrLEQ, arrLNEQ, arrGEQ, arrGNEQ, arrEQ, arrNEQ KEYWORDS: comparison; subarrangement; operator EXAMPLE: example arrLEQ; shows an example" { arr C = B - A; if(C[1] == 0){ return (1); } return (0); } example{example arrEQ;} // returns logical value of 'A==B' static proc arrEQ(arr A, arr B) "USAGE: A==B; A,B arr RETURN: [0,1] true if A equals B NOTE: algorithm is based on arrMinus and does not simplify by scalars, hence some technically equal hyperplanes might not be detected. See example. SEE ALSO: arrLEQ, arrLNEQ, arrGEQ, arrGNEQ, arrEQ, arrNEQ KEYWORDS: comparison; equal; operator EXAMPLE: example arrLEQ; shows an example" { return ((A<=B) & (A>=B)); } example { "EXAMPLE: Relationships between a few arrangements."; echo = 2; ring R = 0,(x,y,z),dp; arr A = ideal(x,y,z); arr B = ideal(x,y); A<=B; A>=B; A==B; A!=B; AB; } // returns logical value of 'A!=B' static proc arrNEQ(arr A, arr B) "USAGE: A!=B; A,B arr RETURN: [0,1] true if A is not equal to B NOTE: algorithm is based on arrMinus and does not simplify by scalars, hence some technically equal hyperplanes might not be detected. See example. SEE ALSO: arrLEQ, arrLNEQ, arrGEQ, arrGNEQ, arrEQ, arrNEQ KEYWORDS: comparison; equal; operator EXAMPLE: example arrLEQ; shows an example" { return (!(A==B)); } example{example arrEQ;} // returns logical value of 'AB' static proc arrGNEQ(arr A, arr B) "USAGE: A>B; A,B arr RETURN: [0,1] true if B is a proper subarrangement of A NOTE: algorithm is based on arrMinus and does not simplify by scalars, hence some technically equal hyperplanes might not be detected. See example. SEE ALSO: arrLEQ, arrLNEQ, arrGEQ, arrGNEQ, arrEQ, arrNEQ KEYWORDS: comparison; subarrangement; operator EXAMPLE: example arrLEQ; shows an example" { return ((A>=B) & (A!=B)); } example{example arrEQ;} //============================================================================// //------------------------------ #6 TYPE CASTING -----------------------------// //============================================================================// // TYPE => ARRANGEMENT proc type2arr(list #) "USAGE: type2arr(#); # def RETURN: [arr] Arrangement defined by the input NOTE: The procedure tries to cast the input to [arr] using arrAdd KEYWORDS: typecasting; arrangement EXAMPLE: example tye2arr; shows an example" { return (arrAdd(#)); } example { "EXAMPLE:"; echo = 2; ring R = 0,(x,y,z),dp; ideal I = x,y,z; typeof(type2arr(I)); } // ARRANGEMENT => POLY static proc arr2poly(arr A) "USAGE: poly(A); arr A RETURN: [poly] The defining polynomial which is the product of polynomials occuring in arr NOTE: The procedure will automatically simplify the polynomial by scalar multiplication. SEE ALSO: arrAdd, arr2poly, arr2mat, arr2list, arr2ideal, type2arr KEYWORDS: typecasting; defining polynomial EXAMPLE: example arr2poly; shows an example" { poly f = simplify(product(A.l),1); if(f == 0){ return (1); } return (f); } example { "EXAMPLE:"; echo = 2; ring R = 0,(x,y,z),dp; arr A = ideal(x+1, 2x-2, y); arr2poly(A); } // ARRANGEMENT => LIST static proc arr2list(arr A) "USAGE: arr2list(A); A arr RETURN: [list] containing all generators of A SEE ALSO: arrAdd, arr2poly, arr2mat, arr2list, arr2ideal, type2arr KEYWORDS: typecasting; list EXAMPLE: example arr2list; shows an example" { int n = size(A); list L; for(int k=1; k<=n; k++){L[k] = A[k];} return (L); } example { "EXAMPLE:"; echo = 2; ring R = 0,(x,y,z),dp; arr A= ideal(x,y,z); arr2list(A); } // ARRANGEMENT => IDEAL static proc arr2ideal(arr A) "USAGE: arr2ideal(A); A arr RETURN: [ideal], the internal ideal of A NOTE: arr2ideal(A); is the same as A.l - which may become private SEE ALSO: arrAdd, arr2poly, arr2mat, arr2list, arr2ideal, type2arr KEYWORDS: typecasting; ideal EXAMPLE: example arr2ideal; shows an example" { return (A.l); } example { "EXAMPLE:"; echo = 2; ring R = 0,(x,y,z),dp; arr A = ideal(x,y,z); arr2ideal(A); } // ARRANGEMENT => MATRIX static proc arr2mat(arr A, list #) "USAGE: matrix(A); A arr matrix(A, 'c') for central arrangements (shorter matrix) RETURN: [matrix] M = [T|b] representing the arrangement NOTE: If the arrangement is central or one is not interested in the const terms one can use "matrix(A, 'c')" instead to get the same matrix without the last column. SEE ALSO: arr2mat, mat2arr, mat2carr KEYWORDS: typecasting; matrix; coefficient EXAMPLE: example arr2mat;" { matrix M = jacob(A.l); if(size(#) == 0){return ( concat(M, transpose(jet(A.l,0))) );} if(#[1] == 'c'){return (M);} ERROR("Bad optional input parameter!") } example { "EXAMPLE:"; echo = 2; ring R = 0,(x,y,z),dp; arr A = ideal(x,y,z); arr D = ideal(x+1, y-2, z, x+y+4); print(arr2mat(A)); print(arr2mat(D)); } // MATRIX => ARRANGEMENT proc mat2arr(matrix M) "USAGE: mat2arr(M); matrix (M|b) RETURN: [arr] interprets the rows of the matrix as the defining polynomial equations of the arrangement where the last column will be considered as the constant terms, i.e. if M is an m*(n+1) matrix we have H_i = Ker( M_i1*x_1 +...+ M_in*x_n + M_i(n+1) ) for i=1...m and A = {H_1,...,H_m} the resulting arrangement. SEE ALSO: mat2carr KEYWORDS: typecasting; matrix; coefficient EXAMPLE: example mat2arr;" { if(ncols(M) > nvars(basering)+1) { ERROR("Matrix too big! Please add variables to basering."); } int n = ncols(M)-1; matrix X[n+1][1]; X[1..n,1] = varMat(1..n); X[n+1 ,1] = 1; arr A = ideal(M*X); return(A); } example { "EXAMPLE:"; echo = 2; ring R = 0,(x,y,z),dp; matrix M[4][4] = 1,0,1,1,1,1,0,2,0,1,1,3,2,1,1,4; print(M); mat2arr(M); } // MATRIX => ARRANGEMENT (central) proc mat2carr(matrix M) "USAGE: mat2carr(M); matrix M RETURN: [arr] interprets the rows of the matrix as the defining polynomial equations of the arrangement. I.e. if M is an m*n matrix we have H_i = Ker( M_i1*x_1 +...+ M_in*x_n) for i=1...m and A = {H_1,...,H_m} the resulting arrangement. SEE ALSO: mat2arr KEYWORDS: typecasting; matrix; coefficient; central EXAMPLE: example mat2carr;" { if( ncols(M) > nvars(basering) ){ ERROR("Error! not enough variables in the basering."); } arr A = ideal(M*varMat(1..ncols(M))); return(A); } example { "EXAMPLE:"; echo = 2; ring R = 0,(x,y,z),dp; matrix M[4][3] = 1,0,1,1,1,0,0,1,1,2,1,1; print(M); mat2carr(M); } //============================================================================// //---------------------------- #7 IDEAL FUNCTIONS ----------------------------// //============================================================================// // checks if A is central, homogenizes static proc arrHomog(arr A, list #) "USAGE: homog(A); arr A homog(A, p); arr A, ring_variable p RETURN: [0,1] homog(A) is the same as arrCentral(A) [arr] homog(A,p) homogenizes A with respect to p NOTE: homog(A,p) is not the same as arrCone(A,p) as it does not add the additional hyperplane SEE ALSO: arrHomog, arrCentral, arrCone KEYWORDS: central; homogenize EXAMPLE: example arrHomog; shows an example" { if(size(#) == 0){ return (homog(A.l)); } if(size(#) == 1){ A = homog(A.l, #[1]); return (A); } ERROR("Too many innput arguments!"); } example { "EXAMPLE: "; echo = 2; ring R = 0,(x,y,z),dp; arr A = ideal(x,y); arr B = ideal(x,y,x+y+1); homog(A); homog(B); homog(B,z); homog(_); } // simplified arrangement static proc arrSimplify(arr A, list #) "USAGE: arrSimplify(A); simplify(A, n); arr A, int n RETURN: [arr] simplified arrangement. NOTE: arrSimplify(A) is the same as simplify(A, 1), simplify with higher ints is not needed SEE ALSO: arrSimplify KEYWORDS: simplify EXAMPLE: example arrSimplify; shows an example" { if(size(#) == 0){ A = simplify(A.l, 1); return (A); } if(size(#) == 1){ A = simplify(A.l, #[1]); return (A); } ERROR("Too many input arguments!"); } example { "EXAMPLE: "; echo = 2; ring R = 0,(x,y,z),dp; arr A = ideal(2x,2y-1,2z-2); arrSimplify(A); } // number of hyperplanes static proc arrSize(arr A) "USAGE: size(A); arr A; RETURN: [int] number of hyperplanes in the arangement NOTE: size(A) also useable for multi-arrangements SEE ALSO: arrSize, multarrSize KEYWORDS: hyperplanes; size; number EXAMPLE: example arrSize; shows an example" { return (size(A.l)); } example { "EXAMPLE: "; echo = 2; ring R = 0,(x,y,z),dp; arr A = ideal(x,y,z); size(A); } // substitute variables static proc arrSubst(arr A, list #) "USAGE: arrSubst(A, #); arr A, ring_variables/integers i,j; RETURN: [arr] with the corresponding substitutions made NOTE: applies 'subst' on the arrangement SEE ALSO: arrSubst KEYWORDS: variables; ring_variable; substitute EXAMPLE: example arrSubst; shows an example" { if(size(#)%2 != 0){ ERROR("Odd number of parameter inputs!"); } for(int i=1; ia+b)=> (a+b|b) =(b->b-a)=> (b|b-a) =(a->b-a)=> (b|a) A = subst(A.l, u, u+v); A = subst(A.l, v, v-u); A = subst(A.l, u, v-u); return (A); } example { "EXAMPLE: "; echo = 2; ring R = 0,(x,y,z),lp; arr A = ideal(x+1,x+y,z); arrSwapVar(A,x,z); } //ring_variable of largest index used in arrangement proc arrLastVar(arr A) "USAGE: arrLastVar(A); arr A RETURN: [int] number of the last variable A uses NOTE: useful if you want a list containing all variables x_1 ... x_k used in A, but you do not want to skip any like variables(A) does. SEE ALSO: varMat, varNum, arrSwapVar, arrLastVar KEYWORDS: variables; ring_variable EXAMPLE: example arrLastVar; shows an example" { return ( rvar(variables(A)[arrNvars(A)]) ); } example { "EXAMPLE: "; echo = 2; ring R = 0,x(1..10),dp; arr A = ideal(x(1), x(2), x(3), x(6)); int n = arrLastVar(A); varMat(1..n); variables(A); } //============================================================================// //--------------------------- #10 CENTER COMPUTATIONS ------------------------// //============================================================================// // checks if arr is central proc arrCentral(arr A) "USAGE: arrCentral(A); arr A RETURN: [0,1] true if arr is central(i.e. all planes intersect in 0) NOTE: This is the same as homog(A) SEE ALSO: arrCentered, arrCentral, arrCenter, arrCentralize KEYWORDS: center; central EXAMPLE: example arrCentral;" { return (homog(A)); } example { "EXAMPLE:"; echo = 2; ring R = 0,(x,y,z),dp; // centered and central arr A = ideal(x,y,z); arrCentered(A); arrCentral(A); // centered but not central (center: (-1,-1/2, 1)) arr B = ideal(x+1,2y+1,-z+1); arrCentered(B); arrCentral(B); } // checks wether arrangement has a center proc arrCentered(arr A) "USAGE: arrCentered(A); arr A RETURN: [0,1] true if A is centered(i.e. intersection of all planes not empty) NOTE: The algorithm uses the rank of matrix: Ax=b has a solution iff rank(A) = rank(A|b) SEE ALSO: arrCentered, arrCentral, arrCenter, arrCentralize KEYWORDS: center EXAMPLE: example arrCentered;" { matrix M = matrix(A); // classic test: Ax=b has a solution if and only if rank(A|b) == rank(A) if( rank(M) == rank(submat(M,1..nrows(M), 1..ncols(M)-1))){ return (1); } return (0); } example { "EXAMPLE:"; echo = 2; ring R = 0,(x,y,z),dp; arr A= ideal(x,y,x-y+1); // centerless arrCentral(A); arr B= ideal(x,y,z); // central with center being the origin arrCentral(B); arr C= ideal(x+1,2y+1,-z+1); // central with center (-1,-1/2, 1) arrCentral(C); } // computes center of an arrangement proc arrCenter(arr A) "USAGE: arrCenter(A); arr A RETURN: [list] L entry 0 if A not centered or entries 1, x, H, where x is any particular point of the center and H is a matrix consisting of vectors which spanning linear intersection space. If there is exactly one solution, then H = 0. NOTE: The intersection of all hyperplanes can be expressed in forms of a linear system Ax=b, where (A|b) is the coeff. matrix of the arrange- ment, which is then solved using L-U decomposition SEE ALSO: arrCentered, arrCentral, arrCenter, arrCentralize KEYWORDS: center EXAMPLE: example arrCenter;" { matrix M = matrix(A); //return matrix (T|b) matrix T = submat(M, 1..nrows(M), 1..ncols(M)-1); matrix b = submat(M, 1..nrows(M), ncols(M))*(-1); list L = ludecomp(T); list Q = lusolve(L[1], L[2], L[3], b); return (Q); } example { "EXAMPLE:"; echo = 2; ring R = 0,(x,y,z),dp; arr A= ideal(x,y,x-y+1); // centerless arrCenter(A); arr B= ideal(x,y,z); // center is a single point arrCenter(B); arr C= ideal(x,z,x+z); // center is a line // here we get a wrong result because the matrix is simplified since A doesn't // contain any "y" the matrix (A|b) will be 3x3 only. arrCenter(C); } // makes centered arrangement central proc arrCentralize(arr A) "USAGE: arrCenteralize(A); arr A RETURN: [arr] A after centralization via coordinate change NOTE: The coordinate change only does translation, vector of translation is the second output of arrCenter SEE ALSO: arrCentered, arrCentral, arrCenter, arrCentralize KEYWORDS: central; center; coordinate change EXAMPLE: example arrCenter;" { if(arrCentral(A)){ print("The arrangement is already central!"); return (); } list L = arrCenter(A); if(L[1] == 0){ print("The arrangement has no center and therefor cannot be centralized!"); return (); } int n = nvars(basering); matrix T = diag(1,n); matrix c = L[2]; A = arrCoordChange(A, T, c); return (A); } example { "EXAMPLE:"; echo = 2; ring R = 0,(x,y,z),dp; arr A = ideal(x-1,y,x-z-1,x-z-1); arrCentralize(A); } //============================================================================// //------------------------ #11 GEOMETRIC CONSTRUCTIONS -----------------------// //============================================================================// // performs coordinate change proc arrCoordChange(arr A, matrix T, list #) "USAGE: arrCoordChange(A, T); arr A, (m*n mat) n*n or n*n+1 matrix T arrCoordChange(A, T , c); arr A, n*n matrix T, n*1 matrix/vector RETURN: [arr]: Arrangement A [A|b] after a coordinate change f: x -> Tx + c with T invertible i.e. [A|b] => [AT^-1|b+AT^-1c] since we have f(H) = f(ker(a1*x1 + ... + an*xn - b)) = {f(x) : a'x -b = 0} = {y : a'f^-1(y) -b = 0)} = {y : a'(T^-1(y-c)) - b = 0} = {y : a'T^-1y -(b + a'T^-1c) = 0} NOTE: There are 3 options how you can give the input (in each case n <= nvars(basering)) 1. Just a nxn matrix with => Will automatically complete T by a unit matrix and perform x -> Tx 2. A nxn matrix T and a nx1 vector/matrix c with => Will automatically complete T and c and perform x -> Tx +c 3. A nxn+1 matrix T with => will use last column as translation vector c SEE ALSO: arrCoordChange, arrCoordNormalize KEYWORDS: coordinate change EXAMPLE: example arrCoordChange; shows an example" { int n = nvars(basering); int k = nrows(T); int l = ncols(T); matrix c; if( k>n || l-1>n ){ ERROR("Matrix too big! (It has more cols/rows than there are variables.)"); } // const vector integrated in matrix => split [T|c] into T and c if( l == k+1 ){ if(size(#) > 0){ ERROR("Bad input. If given a constant vector the matrix must be square!") } c = submat(T, 1..k, l); T = submat(T, 1..k, 1..k); l=k; } if( l != k || rank(T) != k){ ERROR("Given matrix is not a base change matrix.") } if(size(#) > 0){ if(size(#) >= 2){ERROR("Too many input arguments!");} c = matrix(#[1]); if(nrows(c) > n || ncols(c) > 1){ ERROR("Constant vector maldefined. Dimension too big.") } } matrix M[n][n] = diag(1,n); M[1..k,1..k] = T; T = M; T = luinverse(T)[2]; M = jacob(A.l); matrix b[n][1]; b[1..nrows(c),1] = c; c = b; // gives c the right length. b = transpose(jet(A.l, 0)); // constants b = b - M*T*c; M = M*T; A = mat2arr(concat(M,b)); return(A); } example { "EXAMPLE: "; echo =2; ring r = 0,(x,y,z),lp; arr A = x,y,z; arrCoordChange(A,1,[0,0,1]); //lifts z-hyperplane by 1 unit matrix T[2][2] = [0,1,1,0]; // swaps x and y arrCoordChange(A,T); matrix c[2][1] = [1,0]; T = concat(T,c); // now swap x and y and add 1 to x afterwards arrCoordChange(A,T); // Note how T doesn't even need to be a full 3x3 base change matrix. } // performs projection onto coordinate hyperplanes proc arrCoordNormalize(arr A, intvec v) "USAGE: arrCoordChange(A, v); RETURN: [arr]: Arrangement after a coordinate change that transforms the arrangement such that after a tranformation x -> Tx + c we have the arrangement has the matrix representation [AT^-1|b+AT^-1c] such that [AT^-1]_v = I and [b+AT^-1c]_v = 0; NOTE: algorithm performs a base change if H_k is homogenous (i.e. has no) constant term and an affine transformation otherwise Ax+b = 0, Transformation x = Ty+c: AT^-1y + AT^-1c + b = 0 Now we want to have (AT^-1)_v = I and (AT^-1c +b)_v = AT^-1_v*c + b_v = 0 SEE ALSO: arrCoordChange, arrCoordNormalize KEYWORDS: coordinate change EXAMPLE: example arrCoordChange; shows an example" { int n = nvars(basering); if(size(v) > n){ ERROR("Too many rows chosen, at max you can choose "+string(n)); } matrix M = matrix(A); matrix Av[n][n]; matrix bv[n][1]; Av[1..size(v), 1..n] = submat(M,v,1..n); bv[1..size(v),1] = submat(M,v, n+1); if(rank(Av) != n){ if(rank(Av) != size(v)){ ERROR("Normal vectors of the given hyperplanes are not linearily " + "independent. Cannot perform coordinate change!"); } // Adds linearly independent lines to Av in order to make it invertible module F = freemodule(n); module Add = reduce(F,std(transpose(Av))); Add = transpose(simplify(transpose(Add),2)); Av[(size(v)+1)..n, 1..n] = matrix(Add); } if(rank(Av) != n){ERROR("Av not invertible!");} A = arrCoordChange(A,Av,bv); return (A); } example { "EXAMPLE: "; echo=2; ring r = 0,(x,y,z),lp; arr A = ideal(x,y,z,x+z+4); intvec v = 1,2,4; arrCoordNormalize(A,v); } // coned arrangement proc arrCone(arr A, list #) "USAGE: arrCone(A); arrCone(A, ring_variable); arr A arrangement in variables x_1...x_n; RETURN: arr, the coned hyperplane Arrangement cA with respect to the given ring_variable, or the last ring_variable if none was given. NOTE: The hyperplanes are homogenized w.r.t. v and a new hyperplane H = ker(x_n+1) is added. SEE ALSO: arrCone, arrDecone, arrRestrict, arrIsEssential, arrEssentialize KEYWORDS: cone; decone EXAMPLE: example arrCone; shows an example" { poly p; int i; // case 1: no ring_variable given if(size(#) == 0){ p = var(nvars(basering)); } // case 2: a ring_variable is given else{ p = #[1]; } // computation A = homog(A, p); A = A + p; return(A); } example { "EXAMPLE: Coning the arrangement A = (x+1, y) alongside x, y and z:"; echo=2; ring R = 0,(x,y,z),dp; arr A = ideal(x+1, x,x-2,x-1); arrCone(A, y); arr B= ideal(x,y,x+y-1); arrCone(B); } // deconed arrangement proc arrDecone(arr A, int k) "USAGE: arrDecone(A, k); arrangement A, integer k; RETURN: arr: the deconed hyperplane Arrangement dA NOTE: A has to be non-empty and central. arrDecone is an inverse operation to arrCone since A == arrDecone(arrCone(A),size(A)+1) for any A. One can also decone a central arrangement with respect to any hyper- plane k, but than a coordinate change is necessary to make H_k = ker(x_k). Since such a coordinate change is not unique, use arrCoordchange to do so. SEE ALSO: arrCone, arrDecone, arrRestrict, arrIsEssential, arrEssentialize KEYWORDS: cone; decone EXAMPLE: example arrDecone; shows an example" { if( !arrCentral(A) ){ ERROR("Non-central arrangement can not be deconed!");} if( size(A) == 0){ ERROR("Empty arrangement can not be deconed!");} if( size(A) < k){ ERROR("There is no k-th hyperplane");} poly p = A[k]; int n = rvar(p); if( n == 0 ){ ERROR("H_" + string(k) + " = " + string(p) + " is not of the form ker(x_i). Please do a coordinate change first." + "You can use arrCoordinateChange to transform the arrangement accordingly."); } A = A - p; A = arrSubst(A, p, 1); return(A); } example { "EXAMPLE: We decone arr consisting of (x,y,x+z) with respect to y"; echo = 2; ring R = 0,(x,y,z),dp; arr A= ideal(x,y,z,x+y-z); arrDecone(A,3); } proc arrLocalize(arr A, intvec v) "USAGE: arrLocalize(A, v); arrangement A, intvec v; RETURN: arr: the localized arrangement A_X, i.e. A_X only contains the hyperplanes which contain the flat X, which is defined by the equations A[v] SEE ALSO: arrCone, arrDecone, arrRestrict, arrIsEssential, arrEssentialize KEYWORDS: localization EXAMPLE: example arrLocalize; shows an example" { ideal I=std(arr2ideal(A[v])); arr L; int k; for(k=1; k nvars(basering)){ return("Error, too few variables or too high dimension"); } intmat M = random(d,m,n+1); arr A = mat2arr(M); return(A); } example { "EXAMPLE:"; echo = 2; ring R = 0,x(1..20),dp; arrRandom(7,3,15); } // random central arrangement proc arrRandomCentral(int d, int m, int n) "USAGE: arrRandomCentral(d,m,n); int d,m,n RETURN: Random central arrangement, where m is the number of hyperplanes, n the dimension, d the upper bound for absolute value of coefficients. SEE ALSO: arrBoolean, arrBraid, arrTypeB, arrTypeD, arrRandom, arrEdelmanReiner KEYWORDS: example; random; central EXAMPLE: example arrRandomCentral;" { if(n > nvars(basering)){ return("Error, too few variables or too high dimension"); } intmat M = random(d,m,n); arr A = mat2carr(M); return(A); } example { "EXAMPLE:"; echo = 2; ring R = 0,x(1..20),dp; arrRandomCentral(7,3,15); } // Edelman-Reiner arrangement proc arrEdelmanReiner() "USAGE: arrEdelmanReiner(); RETURN: the Edelman-Reiner arrangement, which is a free arrangement but the restriction to the 6-th hyperplane is nonfree. (i.e. counterexample for Orlik-Conjecture) NOTE: the active ring must have at least five variables SEE ALSO: arrBoolean, arrBraid, arrTypeB, arrTypeD, arrRandom, arrEdelmanReiner KEYWORDS: example; Edelman-Reiner EXAMPLE: example arrEdelmanReiner;" { if(nvars(basering) < 5){ ERROR("not enough variables"); } arr arrER = arrBoolean(5); int a,b,c,d; for(a=1; a<=2; a++){ for(b=1; b<=2; b++){ for(c=1; c<=2; c++){ for(d=1; d<=2; d++){ arrER = arrER + (var(1)+(-1)^a*var(2)+(-1)^b*var(3)+(-1)^c*var(4)+(-1)^d*var(5)); }}}} return(arrER); } example { "EXAMPLE:"; echo = 2; ring r=0,x(1..5),dp; arrEdelmanReiner(); } //============================================================================// //--------------------- #13 Orlik Solomon and Poincare Poly ------------------// //============================================================================// // Orlik-Solomon algebra of the arrangement proc arrOrlikSolomon(arr A) "USAGE: arrOrlikSolomon(A); arr A RETURN: [ring] exterior Algebra E as ring with Orlik-Solomon ideal as attribute I. The Orlik-Solomon ideal is generated by the differentials of dependent tuples of hyperplanes. For a complex arrangement the quotient E/I is isomorphic to the cohomology ring of the complement of the arrangement. NOTE: In order to access this ideal I activate this exterior algebra with setring. SEE ALSO: arrOrlikSolomon KEYWORDS: Orlik-Solomon EXAMPLE: example arrOrlikSolomon;" { int central = arrCentral(A); // 0: non-central, 1: central if(central == 0){ A = arrCone(A); } module M = syz(A.l); M = jet(M,0); //Only use linear syzygies M = simplify(M,2); // drop zeros int n = ncols(A.l); def startRing = basering; ring R = 0,e(1..n),dp; def ER = Exterior(); setring ER; // defines the Exterioralgebra ideal I; //final orlik solomon ideal matrix X = transpose(varMat(1..n)); module M = fetch(startRing,M); //brings the module M to the Exterior Algebra ideal OSI = ideal(X*M); if(size(OSI) == 0){ // no relations among hyperplanes, I==0 export(I); return(basering); } // monomial subideal procedure due to Sorin Popescu ideal K = OSI; ideal J; while(isMonomial(K) == 0){ J = lead(std(K)); K = intersect(J,OSI); K = interred(K); } OSI = K; poly diffOp; //now applying the differential operator ideal vars; int j; for(int i=1; i<=ncols(OSI); i++) { vars = variables(OSI[i]); diffOp = 0; for(j=1; j<=ncols(vars); j++) { diffOp = diffOp+(-1)^(j+1)*vars[j]; } I = I + ideal( diff(diffOp,ideal(OSI[i])) ); } if(central ==0 ){//project back in the non-central case n = nvars(basering)-1; ring R = 0,e(1..n),dp; def ERDC = Exterior(); setring ERDC; ideal I = fetch(ER,I); } export(I); return (basering); } example { "EXAMPLE: Computing the Orlik-Solomon-Ideal for the D3-Arrangement"; echo = 2; ring R = 0,(x,y,z),dp; arr A = arrTypeB(3); def E = arrOrlikSolomon(A); setring E; //The generators of the Orlik-Solomon-Ideal are: I; } // The following 3 procedures were replaced by their faster combinatorical counterparts /* // poincare polynomial of the arrangement proc arrPoincare(arr A) "USAGE: arrPoincare(A); arr A RETURN: [intvec] The Poincare polynomial as integer vector of the arrangement, which is equal to the second kind Poincare-Series of the Orlik-Solomon Algebra. SEE ALSO: arrPoincare, arrChambers, arrBoundedChambers KEYWORDS: Poincare polynomial EXAMPLE: example arrPoincare;" { def startRing = basering; def Ext = arrOrlikSolomon(A); setring Ext; ideal OSI = lead(I); OSI = OSI + ideal(Ext); int n = nvars(Ext); ring suppRing = 0,x(1..n),dp; ideal I = fetch(Ext,OSI); intvec HP = hilb(std(I),2); return(HP); } example { "EXAMPLE: Computing the Poincare polynomial as intvec for the D3-Arrangement"; echo = 2; ring R = 0,(x,y,z),dp; arr A = arrTypeB(3); //The coefficients of the Poincare polynomial are: arrPoincare(A); } // number of chambers of the arrangement proc arrChambers(arr A) "USAGE: arrChambers(A); arr A RETURN: [int] The number of chambers of an arrangement, which is equal to the evaluation of the Poincare polynomial at 1. SEE ALSO: arrPoincare, arrChambers, arrBoundedChambers KEYWORDS: chambers EXAMPLE: example arrChambers;" { intvec HP = arrPoincare(A); intvec ones = 1:(size(HP)); return ( transpose(ones)*HP ); } example { "EXAMPLE: Computing the number of chambers for the D3-Arrangement"; echo = 2; ring R = 0,(x,y,z),dp; arr A = arrTypeD(3); //The number of chambers of the D3-Arrangement is: arrChambers(A); } // number of bounded chambers of the arrangement proc arrBoundedChambers(arr A) "USAGE: arrBoundedChambers(A); arr A RETURN: [int] The number of bounded chambers of an arrangement, which is equal to the evaluation of the Poincare polynomial at -1. SEE ALSO: arrPoincare, arrChambers, arrBoundedChambers KEYWORDS: chambers EXAMPLE: example arrBoundedChambers;" { intvec HP = arrPoincare(A); intvec altOnes; for(int i=1; i<=size(HP); i++){ altOnes[i] = (-1)^(i-1); } return ( transpose(altOnes)*HP ); } example { "EXAMPLE: Computing the number of chambers for the D3-Arrangement"; echo = 2; ring R = 0,(x,y,z),dp; arr A = arrTypeD(3); // The number of bounded chambers of the D3-Arrangement is: arrBoundedChambers(A); } */ //============================================================================// //------------------------------- #14 Freeness -------------------------------// //============================================================================// // module of derivation proc arrDer "USAGE: arrDer(A); arr A , multarr A RETURN: [module] The module Der(A) of derivations of the (multi-)arrangement A, i.e. the derivations tangent to each hyperplane of A (resp. with multiplicities) NOTE: This is only defined for central (multi-)arrangements SEE ALSO: arrDer, arrIsFree, arrExponents KEYWORDS: derivation; multiarrangement EXAMPLE: example arrDer;" { def A = #[1]; if( (typeof(A) != "arr") && ( typeof(A) != "multarr") ){ ERROR("bad input type!"); } if(homog(A.l) == 0){ ERROR("Arrangement not central!"); } module fA = jacob(A.l); ideal J; if( typeof(A) == "arr" ){ J = arr2ideal(A); } if(typeof(A) == "multarr"){ J = multIdeal(A); } module tA = diag(matrix(J)); module K = modulo(fA,tA); module derivations = lessGenerators(K); return (derivations); } example { "EXAMPLE: Computing the derivation module of the boolean and braid arrangement"; echo = 2; ring R = 0,(x,y,z),dp; arr A3 = arrBoolean(3); arr B3 = arrTypeB(3); arr G = ideal(x,y,z,x+y+z); //The derivation module of the Boolean 3-arrangement: arrDer(A3); //The derivation module of the Braid 3-arrangement: arrDer(B3); //The derivation module of the generic arrangement: arrDer(G); } // checks if arrangement is free proc arrIsFree(list #) "USAGE: arrIsFree(A); arr A, multarr A RETURN: [0,1] 1 if the (multi-)arrangement is free, i.e. Der(A) is a free module NOTE: only defined for central arrangements SEE ALSO: arrDer, arrIsFree, arrExponents KEYWORDS: free; multiarrangement EXAMPLE: example arrIsFree;" { def A = #[1]; module derivations = arrDer(A); return ( nvars(basering) == ncols(derivations) ); } example { "EXAMPLE: checking freeness of the Edelman-Reiner arrangement and its restriction: "; echo = 2; ring R = 0,(x,y,z),dp; arr A3 = arrBoolean(3); arr B3 = arrTypeB(3); arr G = ideal(x,y,z,x+y+z); arrIsFree(A3); arrIsFree(B3); arrIsFree(G); } // exponents of a (free) arrangement proc arrExponents "USAGE: arrExponents(A); arr A, multarr A RETURN: [intvec] The exponents of a free (multi-) arrangement, i.e. the degrees of a basis of D(A) the derivation module. NOTE: only defined for central arrangements SEE ALSO: arrDer, arrIsFree, arrExponents KEYWORDS: free; exponents; multiarrangement EXAMPLE: example arrExponents;" { def A = #[1]; if(arrIsFree(A) != 1){ ERROR("Arrangement is not free!"); } module der = arrDer(A); intvec exp; for(int i =1; i<=size(der); i++){ exp[i] = deg(der[i]); } return(exp); } example { "EXAMPLE: computing the exponents of the Edelman-Reiner arrangement and its restriction: "; echo = 2; ring R = 0,(x,y,z),dp; arr A3 = arrBoolean(3); arr B3 = arrTypeB(3); arr G = ideal(x,y,z,x+y+z); arrExponents(A3); arrExponents(B3); } // Tries to reduce the number of generators of a generating set for a module static proc lessGenerators(module X){ module Z = syz(X); matrix K = getColumnIndependentUnitPositions(Z); if( K == unitmat(nrows(Z)) ){ return(X); } module Xnew = X*K; Xnew = simplify(Xnew,2); return(lessGenerators(Xnew)); } // Looks for a generator which is redundant static proc getColumnIndependentUnitPositions(module Z){ int n = nrows(Z); // number of generators of D matrix K = unitmat(n); int i; for(int j=1; j<=ncols(Z); j++){ for(i=1; i<=nrows(Z); i++){ if(deg(Z[i,j]) == 0){ K[i,i] = 0; return(K); } } } return(K); } // Outputs the ideal of powers of hyperplanes of a multiarrangement. // Needed for the computation of arrDer. static proc multIdeal(multarr A){ ideal I; for(int i=1; i<=size(A.l); i++){ I[i] = A.l[i]^A.m[i]; } return(I); } //============================================================================// //-------------------------- #15 MULTI-ARRANGEMENTS --------------------------// //============================================================================// //============================================================================// //------------------------------- CONSTRUCTORS -------------------------------// //============================================================================// // general method for creating multiarrangements static proc multarrAdd "USAGE: A = #; A +#; # list containing arr/ideal/list/matrix/poly RETURN: [multarr] multiarrangement constructed by the input parameters. NOTE: algorithm splits up the list # and uses appropiate procedure to handle the input KEYWORDS: multiarrangement; equal; constructor; operator EXAMPLE: example multarrAdd; shows an example" { multarr A; for(int k=1; k<=size(#); k++){ while(1){ //simulates switch, which singular doesn't offer if(typeof(#[k]) == "poly" ) {A = multarrAddPoly (A, #[k]);break;} if(typeof(#[k]) == "ideal") {A = multarrAddIdeal(A, #[k]);break;} if(typeof(#[k]) == "multarr"){A = multarrAddArr (A, #[k]);break;} if(typeof(#[k]) == "list" ) {A = multarrAdd ( #[k]);break;} ERROR("bad input type"); } } return (A); } example { "EXAMPLE: Creating a few multiarrangements"; echo = 2; ring R = 0,(x,y,z),dp; multarr A = ideal(x,y,z,x); A; multarr B = A + ideal(x+1, x-1,y); B; multarr C = list(B, x+1, x-1); C; multarr D = (x2 - y2)^2; D; } // adds a single poly to the arrangement // if the poly is linear, it is just added. If not Singular tries factorization static proc multarrAddPoly(multarr A, poly p) { p=simplify(p,1); if(deg(p) == 0){ ERROR("Given poly is not linear or Singluar is not able to factorize it");} else{ if(deg(p) == 1){ int k; int b = 0; for(k=1; k<=size(A.l); k++){ if(A.l[k] == p){ A.m[k] = A.m[k]+1; b = 1; } } if(b == 0){ A.l[size(A.l)+1] = p; A.m[size(A.l)] = 1; } return(A); } else{ list I = factorize(p,2); if((size(I[1]) == 1) && (deg(I[1][1] > 1))){ERROR("Given poly is not a hyperplane");} else{ int j,i; ideal J; for(i=1; i<=size(I[1]); i++) { for(j=1; j<=I[2][i]; j++){ J[size(J)+1] = I[1][i]; } } return (multarrAdd(A,J)); } } } return(A); } // adds defining polys to the arrangement static proc multarrAddIdeal(multarr A, ideal I){ for(int k=1; k<=size(I); k++){ A = multarrAddPoly(A,I[k]); } return (A); } // union of two multarrangements static proc multarrAddArr(multarr A, multarr B){ return (multarrAddIdeal(A, multIdeal(B))); } //============================================================================// //------------------------------- TYPE CASTING -------------------------------// //============================================================================// // computes the defining polynomial static proc multarr2poly(multarr A) "USAGE: multArrQPoly(multarr A); RETURN: [poly] q: defining polynomial of a multiarrangement with multiplicities SEE ALSO: arr2poly, multarr2poly KEYWORDS: multiarrangement; defining polynomial EXAMPLE: example multArrQPolys;" { poly q=1; for(int i=1;i<=size(A.l);i=i+1) { q=q*A.l[i]^A.m[i]; } return(q); } example { "EXAMPLE: Computing the Q-Poly for a multiarrangement"; echo = 2; ring R = 0,(x,y,z),dp; multarr A = ideal(x,y,z,x,y,x-y,x-z,x-z,y-z); A; poly q=multarr2poly(A); q; } // converts simple arrangement to multiarrangement proc arr2multarr(arr A, intvec v) "USAGE: multArrFromIntvec(arr A, intvec v); RETURN: [multarr] multiarrangement MA, which is the arrangement A with multiplicities v NOTE: the size of v must match the number of hyperplanes of the arrangement A SEE ALSO: arr2multarr, multarr2arr KEYWORDS: multiarrangement; arrangement; constructor EXAMPLE: example arr2multarr" { if(ncols(A.l)!=size(v)){ ERROR("Vector's size does not match the hyperplanes."); } multarr B; B.l=A.l; B.m=v; return(B); } example { "EXAMPLE:"; echo = 2; ring R = 0,(x,y,z),dp; arr A = arrTypeB(3); A; intvec v=2:9; v; multarr MA=arr2multarr(A,v); MA; } // converts multiarrangement to simple arrangement proc multarr2arr(multarr A) "USAGE: multarr2arr(multarr A, intvec v); RETURN: [arr] arrangement A, with all multiplicities removed SEE ALSO: arr2multarr, multarr2arr KEYWORDS: multiarrangement; arrangement; constructor EXAMPLE: example multarr2arr" { arr B; B.l=A.l; return(B); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),dp; multarr A=x2y3z5; A; arr AS = multarr2arr(A); AS; } //============================================================================// //------------------------------- PRINTING -----------------------------------// //============================================================================// // prints arrangement in the console static proc multarrPrint(multarr A) "USAGE: A; A arr RETURN: [] better readable output in the console as newstruct print. SEE ALSO: arrPrint, multarrPrint KEYWORDS: print EXAMPLE: example multarrPrint;" { for(int j=1;j<=ncols(A.l);j++){ print("_["+string(j)+"]=("+string(A.l[j])+")^"+ string(A.m[j])); } } example { "EXAMPLE:"; echo = 2; ring R = 0,(x,y,z),dp; multarr A = ideal(x2,y3,z); A; } // number of hyperplanes with multiplicities static proc multarrSize(multarr A) "USAGE: size(A); A multarr RETURN: [int] Number of hyperplanes with multiplicities SEE ALSO: arrSize KEYWORDS: multiarrangement; size; number; hyperplanes EXAMPLE: example multarrSize;" { return (sum(A.m)); } example { "EXAMPLE:"; echo = 2; ring R = 0,(x,y,z),dp; multarr A = ideal(x2,y3,z); A; } //============================================================================// //------------------------ RESTRICTION & DELETION ----------------------------// //============================================================================// // decrements the multiplicity of a hyperplane by one static proc multarrDelete(multarr A, int k) "USAGE: multarrDelete(A, k); arrangement A, integer k; RETURN: [multarr] the hyperplane Multiarrangement A', i.e. the multiarrangement with multiplicity of H_k decremented by one. If m(H_k)=1, then the hyperplane H_k is deleted SEE ALSO: arrDelete, multarrDelete, KEYWORDS: multiarrangement; delete EXAMPLE: example multarrDelete;" { if(k>ncols(A.l)){ ERROR("There is no k-th hyperplane"); } poly q = multarr2poly(A); q = q/(A.l[k]); multarr MA = q; return (MA); } example { "EXAMPLE:"; echo = 2; ring R = 0,(x,y,z),dp; multarr A =ideal(x2,y3,z); multarr AD = multarrDelete(A,2); AD; } // restriction of A (as arr) to a hyperplane with multiplicities proc multarrRestrict(arr A,intvec v, list #) "USAGE: multarrRestrict(A, v); arrangement A, int/intvec v, optional argument "CC"; RETURN: [multarr] the restricted hyperplane Multi-Arrangement (A^X) with multiplicities i.e. counting how often one element of the restricted arrangement occurs as intersetion of hyperplane of the first arrangement. This definition is due to Guenter M. Ziegler. NOTE: A has to be non-empty. REMARKS: We restrict A to the flat X, defined by the equations in A[v]. The restriction will only be performed, if the ideal defining the flat X is monomial (i.e. X is an intersection of coordinate planes). If the optional argument CC is given, the arrangement is transformed in such a way that X has the above form. SEE ALSO: multarrRestrict, multarrMultRestrict, arrRestrict KEYWORDS: multiarrangement; restriction EXAMPLE: example multarrRestrict; " { ideal I=A.l; I=I[v]; option(redSB); I=std(I); // defining equations for flat X //option(none); ideal Al=arr2ideal(A); int i; if(isMonomial(I)==1){ for(i=1;i<= size(I);i++){ Al=subst(Al,I[i],0); } multarr AR=simplify(Al,3); return(AR); } if(size(#) == 0){ ERROR("The flat X is not defined by a monomial ideal. " + "Please do a coordinate change first. You can use arrCoordNormalize to " + "transform the arrangement accordingly by adding the argument CC."); } if(#[1]!="CC"){ ERROR("The flat X is not defined by a monomial ideal. " + "Please do a coordinate change first. You can use arrCoordNormalize to " + "transform the arrangement accordingly by adding the argument CC."); } if(size(v)==0){return(A);} intvec w=v[1]; intvec tmp; for(i=2;i<=size(v);i++){ tmp=w,v[i]; if(rank(matrix(A[w]))==size(tmp)){ w=tmp; } } arr B=arrCoordNormalize(A,w); ideal Bl=arr2ideal(B); for(i=1;i<= size(w);i++){ Bl=subst(Bl,Bl[w[i]],0); } multarr BR=simplify(Bl,3); return(BR); } example { "EXAMPLE:"; echo = 2; ring R = 0,x(1..5),dp; arr A = arrEdelmanReiner(); A; multarr AR = multarrRestrict(A,6,"CC"); AR; } // restriction of A (as multarr) to a hyperplane with multiplicities proc multarrMultRestrict(multarr A,int k) "USAGE: multarrMultRestrict(A, k); multiarrangement A, integer k; RETURN: [multarr] the restricted hyperplane Multi-Arrangement (A^H_k) with multiplicities, i.e. counting with multiplicities how often one element of the restricted arrangement occurs as intersetion of hyperplane of the first multiarrangement. This definition is due to Guenter M. Ziegler. NOTE: A has to be non-empty. REMARKS: The restriction will only be performed, if H_k = ker(x_i) for some i. One can also restrict an arrangement with respect to any hyper- plane k, but than a coordinate change is necessary first to make H_k = ker(x_k). Since such a coordinate change is not unique, please use arrCoordchange to do so. SEE ALSO: multarrRestrict, multarrMultRestrict, arrRestrict KEYWORDS: multiarrangement; restriction EXAMPLE: example multarrMultRestrict; " { ideal I = variables(A.l); ideal J; int j,i; for(i=1;i<=size(A.l);i=i+1) { for(j=1;j<=A.m[i];j++){ J[ncols(J)+1]=A.l[i]; } } poly p = simplify(A.l[k],1); int n = rvar(p); if( n == 0 ){ ERROR("H_" + string(k) + " = " + string(p) + " is not of the form ker(x_i). Please do a coordinate change first." + "You can use arrCoordinateChange to transform the arrangement accordingly."); } J = subst(J, p, 0); multarr MA = simplify(J,3); return(MA); } example { "EXAMPLE:"; echo = 2; ring R = 0,(x,y,z),dp; multarr A =ideal(x2,y2,z2,(x-y)^3,(x-z)^2,(y-z)); A; //The restriction of the multiarrangement is: multarr AR = multarrMultRestrict(A,1); AR; } //============================================================================// //----------------------------- #16 COMBINATORICS ---------------------------// //============================================================================// /* RELATIONS can have 3 values: -1 : flat&hyperplane are parallel 0 : hyperplane intersects the flat, but only further on in the lattice. +1 : hyperplane is part of the flat. */ // intvec: multiplicities of hyperplanes proc arrLattice(arr ARR) "USAGE: arrLattice(arr ARR) RETURN: [arrposet] intersection poset of the arrangement NOTE: The algorithm works by a bottom up approach, i.e. it calculates the SEE ALSO: arrLattice, arrFlats KEYWORDS: intersection lattice; poset; lattice EXAMPLE: example arrFlats;" { //Performance test system("--ticks-per-sec",1); timer = 0; int start = timer; // start dbprint(3-voice,newline); dbprint(3-voice,"=== Computing poset ==="); dbprint(3-voice,newline); // Initialization list L,L2; intvec u,v, NFLAT; int i,j,l,m,n,rk, rk2,rk3, tic, toc, numberOfFlats,numberOfHplanes, counter, ntests; matrix NewRel, M; // Initialization of matrix ntests = 0; M = matrix(ARR); // m x n+1 matrix m = size(ARR); // m = # hplanes n = nvars(basering); // n = # variables numberOfHplanes = m; // Initialization of poset and flat flat F; arrposet P; P.A = ARR; F.REL = 0:m; F.moebius = -1; for(i=1; i<=n; i++){ P.r = insert(P.r,list());} for(i=1; i<=m; i++){ F.REL[i] = 1; L[i] = F; F.REL[i] = 0;} F.moebius = 0; P.r[1] = L; // calculating r[i] step by step for(rk=2; rk<=rank(ARR); rk++){ // maximal flat possible has rank A.v ofc counter = 0; tic = timer; // Initialization of current round L = P.r[rk-1]; // flats from the last iteration give rise to current ones. numberOfFlats = size(L); L2 = list(); // new list for elts of rk= NewRel = getOldRel(L,m); // contains all the information of the old arrangement // find the next leftmost child by looking for intersetion of F_j with hyperplanes for(j=1; j<=numberOfFlats; j++){ // goes over the elts of L.r[i-1] NFLAT = L[j].REL; for(i=1; i<=numberOfHplanes; i++){ if (NewRel[i,j] == 0){ //test this hyperplane NFLAT[i] = 1; v = collectHplanes(NFLAT); rk2 = rank(submat(M,v,1..n)); ntests = ntests+2; // CHILD FOUND BY INTERSECTING WITH H_i if(rk2 == rank(submat(M,v,1..n+1))){ NewRel[i,j] = 1; counter++; // potentially there exist more hyperplanes which intersect in the same space for(l=i+1; l<=m; l++){ if(NewRel[l,j] == 0){ NFLAT[l] = 1; u = collectHplanes(NFLAT); rk3 = rank(submat(M, u, 1..n)); ntests = ntests +2; // FLAT EXISTS if(rk3 == rank(submat(M,u,1..n+1)) ){ if(rk2 == rk3){NewRel[l,j] = 1;} else{NFLAT[l] = 0;} //itersecting with both yields higher rank flat } // FLAT DOES NOT EXIST else{NFLAT[l] = -1;} } } // Add new flat to list, reset NFLAT. F.REL = NFLAT; NFLAT = L[j].REL; // refresh. // Find Parents: (most expensive part it seems) F.parents = j:1; //parent in any case for(l=j+1; l<=numberOfFlats; l++){ if(isParent(L[l].REL, F.REL)){ F.parents = F.parents,l; NewRel[1..m,l] = updateRel(F.REL, submat(NewRel,1..m,l)); } } L2 = insert(L2,F,size(L2)); // means that the flat is parallel } else{ L[j].REL[i] = -1; NewRel[i,j] = -1; NFLAT[i] = -1; } } } } toc = timer; dbprint(3-voice,"rank "+string(rk)+": found "+ string(counter) + " flats in "+string(toc - tic)+"s"); counter = 0; P.r[rk] = L2; } dbprint(3-voice,newline); //dbprint(3-voice,"Computation time: "+string(timer - start)+"s"); dbprint(3-voice,"Matrix tests: "+string(ntests)); return (P); } example { "EXAMPLE: Intersection lattice of the braid arrangement in 3 dimensions "; echo = 2; ring r; arrLattice(arrTypeB(3)); } static proc arrRank(arr A){ int n = rank(matrix(A)); int m = nvars(A); if(n < m) {return (n);} return (m); } static proc getOldRel(list L, int m){ int n = size(L); matrix NewRel[m][n]; for(int i=1; i<=n; i++){ NewRel[1..m,i] = L[i].REL; } return (NewRel); } static proc updateRel(intvec REL, matrix NewRel){ for(int i=1; i<=size(REL); i++){ if(NewRel[i,1] == 0 && REL[i] == 1){NewRel[i,1] = 1;} } return (NewRel); } // Flat F is a parent of Flat G if F[i] == 1 => G[i] == 1 for all i static proc isParent(intvec parent, intvec child){ int counter = 0; for(int i=1; i<= size(parent); i++){ if(parent[i]==1){ if(child[i] == 1){counter++;} else{ return (0); } } } if(counter == 0){ return (0); } return (1); } // transforms the "rel" intvec into an intvec which contains the indices of the hyperplanes. static proc collectHplanes(intvec RELATIONS){ intvec result; for(int i=1; i<=size(RELATIONS); i++){ if(RELATIONS[i] > 0){result = result,i;} } result = result[2..size(result)]; return (result); } static proc getFlat(arrposet P, int i, int j){ list L = P.r[i]; return (L[j]); } static proc setFlat(arrposet P, int i, int j, flat F){ list L = P.r[i]; L[j] = F; P.r[i] = L; return (P); } static proc getFlag(arrposet P, int i, int j){ flat F = getFlat(P,i,j); return (F.flag); } static proc setFlag(arrposet P, int i, int j, int flag){ flat F = getFlat(P,i,j); F.flag = flag; return (setFlat(P,i,j,F)); } //============================================================================// //------------------------- #16 INTERSECTION-LATTICE -------------------------// //============================================================================// //============================================================================// //-------------------------- MOEBIUS -------------------------------------// //============================================================================// // calculates the moebius values of the poset proc moebius(arrposet P) "USAGE: moebius(arrposet P) RETURN: [arrposet] fills in the moebius values of the flats in the poset SEE ALSO: moebius KEYWORDS: moebius function EXAMPLE: example arrCharPoly; shows an example" { list L; int i,j; for(i=2; i<=rank(P.A); i++){ L = P.r[i]; for(j=1; j<=size(L); j++){ arrposet Q = P; export Q; L[j].moebius = -(1 + moebiusRecursion(i, j)); kill Q; } P.r[i] = L; } return (P); } example { "EXAMPLE: We look at the moebius values of the braid arrangement in 4 dimensions: "; echo = 2; ring R = 0,(x,y,z,t),dp; arr A = arrBraid(4); arrposet P = arrLattice(A); P; //As you can see the values are not calculated yet: printMoebius(P); P = moebius(P); //Now all entries are initialized: printMoebius(P); } // moebiusrecursion(i,j) is the sum moebius(X) over {X | V>X>Flat_ij} // Vss: all moebius values of rank < i are known and correctly saved in P static proc moebiusRecursion(int i, int j){ flat F = getFlat(Q, i, j); int m = F.moebius; if(getFlag(Q,i,j) == 1){return (0);} else { Q = setFlag(Q,i,j,1); if(i > 1){ for(int k=1; k<=size(F.parents); k++){ m = m + moebiusRecursion(i-1, F.parents[k]); } } } return (m); } // pi(A,t) = sum[X in L(A); moebius(X)*(-t)^r(X)] proc arrPoincare(def input) "USAGE: arrPoincare(A); arr A RETURN: [intvec] The Poincare polynomial as integer vector of the arrangement, which is equal to the second kind Poincare-Series of the Orlik-Solomon Algebra. SEE ALSO: arrPoincare, arrChambers, arrBoundedChambers KEYWORDS: Poincare polynomial EXAMPLE: example arrPoincare;" { while(1){ if(typeof(input) == "arr"){ arr A = input; arrposet P = arrLattice(A); break; } if(typeof(input) == "arrposet"){ arr A = input.A; arrposet P = input; break; } ERROR("Bad input type"); } P = moebius(P); int i, j, coeff, sign; list L; intvec v = 1; sign = 1; for(i=1; i<= rank(A); i++){ L = P.r[i]; sign = (-1)*sign; coeff = 0; for(j=1;j<=size(L);j++){ coeff = coeff + L[j].moebius; } coeff = sign*coeff; v = v,coeff; } return (v); } example { "EXAMPLE: The poincare polynomial of the braid arrangement in k dimensions is given as: " + "pi(A,t) = (1 + t)*...*(1 + (k-1)*t)"; echo = 2; ring R = 0,(x,y,z,u,v),dp; arr A = arrBraid(5); intvec v = arrPoincare(A); (1+x)*(1+2x)*(1+3x)*(1+4x); v; } // X(A,t) = t^l * pi(A,-t^-1) proc arrCharPoly(def input) "USAGE: arrCharPoly(arr A) RETURN: [intvec] coefficients of the characteristic polynomial of A in incresing order REMARKS: The algorithm only returns the coefficients of the characteristic polynomial since they are whole numbers but the basering could be something different. SEE ALSO: arrCharPoly, arrPoincare KEYWORDS: characteristic polynomial EXAMPLE: example arrCharPoly; shows an example" { intvec v = arrPoincare(input); int l = size(v); v = v[l..1]; //reverses order for(int i=2; i<=l; i=i+2){ v[i] = -v[i]; } return (v); } example { "EXAMPLE: The characteristic polynomial of the Braid arrangement in k dimensions is given as: " + "X(A,t) = t*(t-1)*...*(t-(k-1)))"; echo = 2; ring R = 0,(x,y,z,u,v),dp; arr A = arrBraid(5); intvec v = arrCharPoly(A); x*(x-1)*(x-2)*(x-3)*(x-4); v; } // number of chambers of the arrangement proc arrChambers(arr A) "USAGE: arrChambers(A); arr A RETURN: [int] The number of chambers of an arrangement, which is equal to the evaluation of the Poincare polynomial at 1. SEE ALSO: arrPoincare, arrChambers, arrBoundedChambers KEYWORDS: chambers EXAMPLE: example arrChambers;" { intvec v = arrPoincare(A); return(sum(v)); } example { "EXAMPLE: Computing the number of number of chambers in a simple arrangement"; echo = 2; ring R = 0,(x,y),dp; arr A = ideal(x,y,x+y-1); arrChambers(A); } // number of bounded chambers of the arrangement proc arrBoundedChambers(arr A) "USAGE: arrBoundedChambers(A); arr A RETURN: [int] The number of bounded chambers of an arrangement, which is equal to the evaluation of the Poincare polynomial at -1. SEE ALSO: arrPoincare, arrChambers, arrBoundedChambers KEYWORDS: chambers EXAMPLE: example arrBoundedChambers;" { intvec v = arrPoincare(A); for(int i=2; i<=size(v); i=i+2){ v[i] = -v[i]; } return (sum(v)); } example { "EXAMPLE: Computing the number of bounded chambers in a simple arrangement"; echo = 2; ring R = 0,(x,y),dp; arr A = ideal(x,y,x+y-1); arrBoundedChambers(A); } //============================================================================// //------------------------------- OUTPUT ------------------------------------// //============================================================================// static proc arrPrintPoset(arrposet P){ print("Given Arrangement:"); P.A; print("Corresponding poset:"); string s; int i,j; list L; for(i=1; i<= size(P.r); i++){ print("====== rank "+string(i)+": "+string(size(P.r[i]))+" flats ======"); s = ""; L = P.r[i]; for(j=1; j<=size(L); j++){ s = s + " (" + string(collectHplanes(L[j].REL)) + "), "; } print(s); if(size(P.r[i]) == 0){break;} } } proc printMoebius(arrposet P) "USAGE: printMoebius(A); arr A RETURN: [] displays the moebius values of all the flats in the poset REMARKS: Mainly used for debugging. EXAMPLE: example printMoebius;" { print("Moebius values: "); string s; int i,j; list L; for(i=1; i<= size(P.r); i++){ print("====== rank "+string(i)+": "+string(size(P.r[i]))+" flats ======"); s = ""; L = P.r[i]; for(j=1; j<=size(L); j++){ s = s + " (" + string(L[j].moebius) + "), "; } print(s); if(size(P.r[i]) == 0){break;} } } example { "EXAMPLE: We look at the moebius values of the braid arrangement in 4 dimensions: "; echo = 2; ring R = 0,(x,y,z,t),dp; arr A = arrBraid(4); arrposet P = arrLattice(A); P; //As you can see the values are not calculated yet: printMoebius(P); P = moebius(P); //Now all entries are initialized: printMoebius(P); } //============================================================================// //-------------------------------- arrFlats Stuff --------------------------------// //============================================================================// // test via hyperplanes. proc arrFlats(arr ARR) "USAGE: size(A); A arr RETURN: [arrposet] Intersection lattice SEE ALSO: arrFlats KEYWORDS: intersection lattice EXAMPLE: example arrFlats;" { print(newline); print("=== Computing poset ==="); print(newline); //Performance test system("--ticks-per-sec",1000); timer = 0; int time = timer; // Initialization list L,L2; intvec u,v,w,src; int i,j,k,l,m,n,d, rk,rk2; int counter, ntests; ntests = 0; matrix S,T; intmat REL; // Initialization of matrix matrix M = matrix(ARR); // m x n+1 matrix m = size(ARR); n = nvars(basering); // Initialization of P arrflats P; P.A = ARR; for(i=1; i<=n; i++){ P.r = insert(P.r,list());} for(i=1; i<=m; i++){ L[i] = i;} P.r[1] = L; // calculating r[i] step by step for(i=2; i<=n; i++){ // maximal flat possible has rank A.v ofc counter = 0; L = P.r[i-1]; // flats from the last iteration give rise to current ones. if(size(L) <= 1){break;} // finish early if no more intersections possible. L2 = list(); // new list for elts of rk=i // REL is an intmat that contains information about the relations between // the i-1 flats and the Hyperplanes. // REL[i,j] = 1 means that H_i intersects F_j // REL[i,j] = 0 means that it has not been tested yet // REL[i,j] = -1 means that H_i does NOT intersect F_j, i.e. they are parallel // resetting the REL-matrix REL = intmat(intvec(0),m,size(L)); // filles in the trivial entries, i.e. each flat intersects all hplanes it is composed of. for(j=1; j<=size(L); j++){ u = L[j]; for(k=1; k<= size(u); k++){ REL[u[k],j] = 1; } } //find new flat of rank i for(j=1; j<=size(L); j++){ // goes over the elts of L.r[i-1] // Looking for intersetion of F_j with hyperplanes u = L[j]; for(k=1; k<=m; k++){ if (REL[k,j] == 0 ){ v = insertVal(u,k); counter++; ntests = ntests +2; rk = rank(submat(M,v,1..n)); if(rk == rank(submat(M,v,1..n+1))){ //INTERSECTION FOUND // potentially there exist more hyperplanes which intersect in the same // space. For example (x,y,x+y) all intersect in (0,0) REL[k,j] = 1; for(l=k+1; l<=m; l++){ if(REL[l,j] == 0){ u = insertVal(v,l); ntests = ntests +2; rk2 = rank(submat(M, u, 1..n)); if( rk == rk2 && rk == rank(submat(M,u,1..n+1)) ){ v = u; REL[l,j] = 1; } }} // Add new flat to list, reset u. L2 = insert(L2,v,size(L2)); u = L[j]; //if k = m this is not necessary.... } }} } print("rank: "+string(i)+", expected OPS: "+string(binomial(size(L),2)) +", excecuted OPS: " + string(counter)); counter = 0; //cleanup => doing this during computation increases speed!! for(j=1; ju[size(u)]){ u = u,k; return (u); } int i = 2; while(u[i] v[n]){res[k] = v[n]; n++; break;} if(u[m] == v[n]){res[k] = u[m]; m++; n++; break;} ERROR("Something went wrong in proc mergeIV"); } } while(m <= size(u)){res[k] = u[m]; m++; k++;} while(n <= size(v)){res[k] = v[n]; n++; k++;} return (res); } // checks if v[i] < u[i] for all i static proc isSmaller(intvec u, intvec v){ for(int i= 1; i<=size(u); i++){ if(u[i] > v[i]){ return (0); } } return (1); } static proc new2old(arrposet P){ arrflats Q; Q.A = P.A; int i,j; list L = P.r; for(i=1; i<=size(L); i++){ list L2 = L[i]; for(j=1; j<=size(L2); j++){ L2[j] = collectHplanes(L2[j].REL); } L[i] = L2; } Q.r = L; return (Q); } /********************************************************************* newstruct("arr","ideal l"); Defines a hyperplane arrangement by a list of linear polynomials such that the hyperplanes are the varieties of those polynomials. supported operators: A = input defines an arrangement by the input which may consist of the types arr/ideal/poly/matrix/list A + input adds arrangement defined by the right hand side to the left hand side A - input deletes hyperplanes from the arrangement <, >, <=, >=, ==, != set theoretical comparison of two arrangements A[int] access to a single hyperplane A[intvec] subarrangement defined by the indexed hyperplanes A; prints the arrangement on the screen supported type conversions: matrix(A) returns coefficient matrix of the defining polynomials poly(A) returns the defining polynomial which is the product of all the single polynomials list(A) returns the defining polynomials as a list ideal(A) returns the defining polynomials as an ideal type2arr(input) returns arrangement defined by the input which may consist of the types arr/ideal/poly/matrix/list supported inherited functions: delete(A, intvec) deletes indexed hyperplanes from the arrangement homog(A) checks wether the defining polys are all homogenous <=> arr is central homog(A, rvar) homogenizes the def. polys with respect to the given ring variable size(A) number of hyperplanes in the arrangement subst(A,rvar,poly,...) substitutes ringvariables with given polynomials, though they need to remain linear variables(A) ideal generated by the ring variables that A depends on nvars(A) number of ring variables that A depends on newstruct("multarr","ideal l, intvec m"); Defines a hyperplane arrangement with multiplicities by a list of linear polynomials such that the hyperplanes are the varieties of those polynomials and an intvec in which the multiplicities are saved. supported operators: M = input defines an multarr by the input which may consist of the types multarr/ideal/poly/list M + input adds arrangement defined by the right hand side to the left hand side M; prints the multarr on the screen supported type conversions: poly(M) returns defining polynomial with multiplicities arr2multarr(A,intvec) returns multarr with multiplicities defined by the intvec. multarr2arr(M) returns the internal arrangement (all multiplicities set to 1) delete(M,int) decrements the multiplicity of the hyperplane defined by the index by 1 size(M) returns number of hyperplanes counting multiplicities. *-----------------------------------------------------------------------*/