1 | ////////////////////////////////////////////////////////////////////////////// |
---|
2 | version="version arr.lib 4.1.2.0 Feb_2019 "; // $Id$ |
---|
3 | category="Miscellaneous"; |
---|
4 | |
---|
5 | /* |
---|
6 | |
---|
7 | ** TOPICS ** (Ctrl+f to search) |
---|
8 | #1 NEWSTRUCTS & OVERLOADS |
---|
9 | #2 CONSTRUCTORS |
---|
10 | #3 ACCESS & ASSIGNMENT |
---|
11 | #4 DELETION |
---|
12 | #5 COMPERATORS |
---|
13 | #6 TYPE CASTING |
---|
14 | #7 INHERITED FUNCTIONS |
---|
15 | #8 PRINTING |
---|
16 | #9 MANIPULATING VARIABLES |
---|
17 | #10 CENTER COMPUTATIONS |
---|
18 | #11 GEOMETRIC CONSTRUCTIONS |
---|
19 | #11 EXAMPLES OF ARRANGEMENTS |
---|
20 | #13 ORLIK SOLOMON AND POINCARE POLYNOMIAL |
---|
21 | #14 FREENESS |
---|
22 | #15 MULTI-ARRANGEMENTS |
---|
23 | #16 COMBINATORICS |
---|
24 | |
---|
25 | */ |
---|
26 | |
---|
27 | info=" |
---|
28 | LIBRARY: arr.lib a library of algorithms for arrangements of hyperplanes |
---|
29 | |
---|
30 | AUTHORS: Randolf Scholz (rscholz@rhrk.uni-kl.de), |
---|
31 | Patrick Serwene (serwene@mathematik.uni-kl.de), |
---|
32 | Lukas Kuehne (lf.kuehne@gmail.com) |
---|
33 | |
---|
34 | OVERLOADS: |
---|
35 | |
---|
36 | // OPERATORS |
---|
37 | = arrAdd assignment |
---|
38 | + arrAdd union of two arrs |
---|
39 | [ arrGet access to a single/multiple hyperplane(s) |
---|
40 | - arrMinus deletes given hyperplanes from the arr |
---|
41 | <= arrLEQ comparison |
---|
42 | >= arrGEQ comparison |
---|
43 | == arrEQ comparison |
---|
44 | != arrNEQ comparison |
---|
45 | < arrLNEQ comparison |
---|
46 | > arrGNEQ comparison |
---|
47 | |
---|
48 | // TYPECASTING |
---|
49 | matrix arr2mat coeff matrix |
---|
50 | poly arr2poly defining polynomial |
---|
51 | |
---|
52 | // OTHER |
---|
53 | variables arrVariables ideal generated by the variables the arr depends on |
---|
54 | nvars arrNvars number of variables the arr depends on |
---|
55 | delete arrDelete deletes hyperplanes by indices |
---|
56 | print arrPrint prints the arr on the screen |
---|
57 | |
---|
58 | // IDEAL INHERITED FUNCTIONS |
---|
59 | homog arrHomog checks if arrangement is homogeneous |
---|
60 | simplify arrSimplify simplifies arrangement |
---|
61 | size arrSize number of planes |
---|
62 | subst arrSubst substitute variables |
---|
63 | |
---|
64 | // MULTI-ARRANGEMENTS |
---|
65 | = multarrAdd assignment of multarr |
---|
66 | + multarrAdd union of multarr |
---|
67 | poly multarr2poly defining polynomial |
---|
68 | size multarrSize number of hyperplanes with mult. |
---|
69 | print multarrPrint displays multiarr |
---|
70 | delete multarrDelete deletes hyperplane |
---|
71 | |
---|
72 | PROCEDURES: |
---|
73 | arrSet(arr A, int k, poly p) replaces the k-th Hyperplane with poly p |
---|
74 | |
---|
75 | type2arr(#) converts general input to 'arr' using arrAdd. |
---|
76 | mat2arr( matrix M) affine arrangement from coeff matrix |
---|
77 | mat2carr(matrix M) central arrangement from coeff matrix |
---|
78 | arrPrintMatrix(arr A) readable output as a coeff matrix |
---|
79 | varMat(intvec v) matrix of the corresponding ring_variables |
---|
80 | varNum(def u) number of given variable (enh. version of varNum in dmod.lib) |
---|
81 | arrSwapVar(arr A, i, j) swaps two variables in the arrangement |
---|
82 | arrLastVar(arr A) ring_variable of largest index used in arrangement |
---|
83 | arrCenter(arr A) computes center of an arrangement |
---|
84 | arrCentral(arr A) checks if arrangement is central |
---|
85 | arrCentered(arr A) checks if arrangement is centered |
---|
86 | arrCentralize(arr A) makes centered arrangement central |
---|
87 | arrCoordChange(A, T, #) performs coordinate change |
---|
88 | arrCoordNormalize(A, v) performs projection onto coordinate hyperplane |
---|
89 | arrCone(arr A, var) coned arrangement |
---|
90 | arrDecone(arr A, int k) deconed arrangement |
---|
91 | arrLocalize(arr A, intvec v) localization of an arrangement onto a flat |
---|
92 | arrRestrict(arr A, intvec v) restricted arrangement onto a flat |
---|
93 | arrIsEssential(arr A) checks if arrangement is essential |
---|
94 | arrEssentialize(arr A) essentialized arragnement |
---|
95 | arrBoolean(int v) boolean arrangement |
---|
96 | arrBraid(int v) braid arrangement |
---|
97 | arrTypeB(int v) type B arrangement |
---|
98 | arrTypeD(int v) type D arrangement |
---|
99 | arrRandom(d,m,n) random (affine) arrangement |
---|
100 | arrRandomCentral(d,m,n) random central arrangement |
---|
101 | arrEdelmanReiner() Edelman-Reiner arrangement |
---|
102 | arrOrlikSolomon(arr A) Orlik-Solomon algebra of the arrangement |
---|
103 | arrDer(A) module of derivation |
---|
104 | arrIsFree(A) checks if arrangement is free |
---|
105 | arrExponents(A) exponents of a (free) arrangement |
---|
106 | arr2multarr(arr A, intvec v) converts normal arrangement to multiarrangement |
---|
107 | multarr2arr(multarr A) converts multiarrangement to normal arrangement |
---|
108 | multarrRestrict(arr A, v) restriction of A (as arr) to a flat with multiplicities |
---|
109 | multarrMultRestrict(A, int k) restriction of A (as multarr) to a hyperplane with multiplicities |
---|
110 | arrFlats(arr A) intersection lattice |
---|
111 | arrLattice(arr A) computes the intersection lattice / poset |
---|
112 | moebius(arrposet P) computes moebius values |
---|
113 | arrCharPoly(arr A) characteristic polynomial |
---|
114 | arrPoincare(arr A) poincare polynomial of the arrangement |
---|
115 | arrChambers(arr A) number of chambers of the arrangement |
---|
116 | arrBoundedChambers(arr A) number of bounded chambers of the arrangement |
---|
117 | printMoebius(arr A) displays the moebius values of all the flats in the poset |
---|
118 | "; |
---|
119 | |
---|
120 | |
---|
121 | //============================================================================// |
---|
122 | //-------------------------- #1 NEWSTRUCTS & OVERLOADS -----------------------// |
---|
123 | //============================================================================// |
---|
124 | |
---|
125 | // initialization of the library |
---|
126 | static proc mod_init() |
---|
127 | { |
---|
128 | // NEWSTRUCTS |
---|
129 | newstruct("arr","ideal l"); |
---|
130 | newstruct("flat", "intvec REL, int moebius, intvec parents, int flag"); |
---|
131 | newstruct("arrposet","arr A, list r"); |
---|
132 | newstruct("arrflats","arr A, list r"); |
---|
133 | newstruct("multarr","ideal l, intvec m"); // intvec: multiplicities of hyperplanes |
---|
134 | |
---|
135 | // OPERATORS |
---|
136 | system("install","arr","=" ,arrAdd ,1); // assignment |
---|
137 | system("install","arr","+" ,arrAdd ,2); // union of arrs |
---|
138 | system("install","arr","[" ,arrGet ,2); // access |
---|
139 | system("install","arr","-" ,arrMinus ,2); // delete plane |
---|
140 | system("install","arr","<=" ,arrLEQ ,2); // comparison |
---|
141 | system("install","arr",">=" ,arrGEQ ,2); // comparison |
---|
142 | system("install","arr","==" ,arrEQ ,2); // comparison |
---|
143 | system("install","arr","!=" ,arrNEQ ,2); // comparison |
---|
144 | system("install","arr","<" ,arrLNEQ ,2); // comparison |
---|
145 | system("install","arr",">" ,arrGNEQ ,2); // comparison |
---|
146 | |
---|
147 | // TYPECASTING |
---|
148 | system("install","arr","matrix" ,arr2mat ,1); // coeff matrix |
---|
149 | system("install","arr","poly" ,arr2poly ,1); // defining polynomial |
---|
150 | system("install","arr","list" ,arr2list ,4); // list of defining polynomials |
---|
151 | system("install","arr","ideal" ,arr2ideal ,4); // list of defining polynomials |
---|
152 | |
---|
153 | // OTHER |
---|
154 | system("install","arr","variables" ,arrVariables ,1); |
---|
155 | system("install","arr","nvars" ,arrNvars ,1); |
---|
156 | system("install","arr","delete" ,arrDelete ,2); |
---|
157 | system("install","arr","print" ,arrPrint ,1); |
---|
158 | |
---|
159 | // IDEAL INHERITED FUNCTIONS |
---|
160 | system("install","arr","homog" ,arrHomog ,1); // checks if homogeneous |
---|
161 | system("install","arr","homog" ,arrHomog ,2); // checks if homogeneous |
---|
162 | system("install","arr","simplify" ,arrSimplify ,2); // simplifies arrangement |
---|
163 | system("install","arr","size" ,arrSize ,1); // number of planes |
---|
164 | system("install","arr","subst" ,arrSubst ,4); // substitute variables |
---|
165 | |
---|
166 | // MULTI-ARRANGEMENTS |
---|
167 | system("install","multarr","=" ,multarrAdd ,1); // assignment of multarr |
---|
168 | system("install","multarr","+" ,multarrAdd ,2); // union of multarr |
---|
169 | system("install","multarr","poly" ,multarr2poly ,1); // defining polynomial |
---|
170 | system("install","multarr","size" ,multarrSize ,1); // number of hyperplanes with mult. |
---|
171 | system("install","multarr","print" ,multarrPrint ,1); // displays multiarr |
---|
172 | system("install","multarr","delete" ,multarrDelete ,2); // deletes hyperplane |
---|
173 | |
---|
174 | // COMBINATORICS |
---|
175 | system("install","arr","rank" ,arrRank ,1); |
---|
176 | system("install","arrflats","print" ,arrPrintFlats ,1); |
---|
177 | system("install","arrposet","print" ,arrPrintPoset ,1); |
---|
178 | |
---|
179 | // NEEDED LIBRARIES |
---|
180 | LIB "general.lib"; |
---|
181 | LIB "monomialideal.lib"; |
---|
182 | |
---|
183 | } |
---|
184 | |
---|
185 | |
---|
186 | //============================================================================// |
---|
187 | //--------------------------- #2 CONSTRUCTORS --------------------------------// |
---|
188 | //============================================================================// |
---|
189 | |
---|
190 | |
---|
191 | // general method for creating arrangements |
---|
192 | static proc arrAdd |
---|
193 | "USAGE: A = #; A +#; # list containing arr/ideal/list/matrix/poly |
---|
194 | RETURN: [arr] Arrangement constructed by input parameters. |
---|
195 | REMARKS: The algorithm splits up the list # and uses the appropriate procedure |
---|
196 | to handle the input. |
---|
197 | NOTE: arrAdd simplifies certain inputs, for example A = (x,y,2x); gives the same arrangement |
---|
198 | as A = (x,y); |
---|
199 | SEE ALSO: arrAdd, type2arr |
---|
200 | KEYWORDS: arrangement; equal; constructor; operator |
---|
201 | EXAMPLE: example arrAdd; shows an example" |
---|
202 | |
---|
203 | { |
---|
204 | arr A; |
---|
205 | for(int k=1; k<=size(#); k++){ |
---|
206 | while(1){ //simulates switch, which singular doesn't offer |
---|
207 | if(typeof(#[k]) == "arr" ){A = arrAddArr (A, #[k]);break;} |
---|
208 | if(typeof(#[k]) == "poly" ){A = arrAddPoly (A, #[k]);break;} |
---|
209 | if(typeof(#[k]) == "ideal" ){A = arrAddIdeal (A, #[k]);break;} |
---|
210 | if(typeof(#[k]) == "matrix"){A = arrAddMatrix(A, #[k]);break;} |
---|
211 | if(typeof(#[k]) == "intmat"){A = arrAddMatrix(A, #[k]);break;} |
---|
212 | if(typeof(#[k]) == "list" ){A = arrAdd ( #[k]);break;} |
---|
213 | |
---|
214 | ERROR("bad input type"); |
---|
215 | } |
---|
216 | } |
---|
217 | |
---|
218 | return (A); |
---|
219 | } |
---|
220 | |
---|
221 | example |
---|
222 | { |
---|
223 | "EXAMPLE: Creating a few arrangements"; echo = 2; |
---|
224 | ring R = 0,(x,y,z),dp; |
---|
225 | arr A= ideal(x,y,z); A; |
---|
226 | arr B = A + ideal(x+1, x-1); B; |
---|
227 | arr C = list(A, x+1, x-1); C; |
---|
228 | arr D = x2 - y2; D; |
---|
229 | } |
---|
230 | |
---|
231 | // union of two arrangements |
---|
232 | static proc arrAddArr(arr A, arr B){ |
---|
233 | return (arrAddIdeal(A, B.l)); |
---|
234 | } |
---|
235 | |
---|
236 | // adds a single poly to the arrangement |
---|
237 | // if the poly is linear, it is just added, otherwise Singular factorizes |
---|
238 | static proc arrAddPoly(arr A, poly p){ |
---|
239 | if(deg(p) == 0){ |
---|
240 | ERROR("Given poly is not linear or Singluar is not able to factorize it");} |
---|
241 | else{ |
---|
242 | if(deg(p) == 1){ |
---|
243 | A.l = A.l + p; |
---|
244 | return (A); |
---|
245 | } |
---|
246 | else{ |
---|
247 | ideal I = factorize(p,1); |
---|
248 | if(size(I) == 1){ERROR("Given poly is not a hyperplane");} |
---|
249 | else{ return (arrAdd(A,I)); } |
---|
250 | } |
---|
251 | } |
---|
252 | return(A); |
---|
253 | } |
---|
254 | |
---|
255 | // adds defining polys to the arrangement |
---|
256 | static proc arrAddIdeal(arr A, ideal I){ |
---|
257 | for(int k=1; k<=size(I); k++){ |
---|
258 | A = arrAddPoly(A,I[k]); |
---|
259 | } |
---|
260 | |
---|
261 | return (A); |
---|
262 | } |
---|
263 | |
---|
264 | // adds defining polys to the arrangement |
---|
265 | static proc arrAddMatrix(arr A, matrix M){ |
---|
266 | return ( arrAddArr(A,mat2arr(M)) ); |
---|
267 | } |
---|
268 | |
---|
269 | |
---|
270 | //============================================================================// |
---|
271 | //--------------------------- #3 ACCESS & ASSIGNMENT -------------------------// |
---|
272 | //============================================================================// |
---|
273 | |
---|
274 | |
---|
275 | // access to hyperplanes, overloads [] operator |
---|
276 | static proc arrGet(arr A, intvec v) |
---|
277 | "USAGE: A[v]; v int/intvec |
---|
278 | RETURN: [poly] if v is [int] The defining poly of the the v-th hyperplane+ |
---|
279 | [arr] if v is [intvec] The corresponding subarrangement |
---|
280 | SEE ALSO: arrGet, arrSet |
---|
281 | KEYWORDS: arrangement; get; operator |
---|
282 | EXAMPLE: example arrGet; shows an example" |
---|
283 | |
---|
284 | { |
---|
285 | if(size(v) == 1){ return (A.l[v[1]]); } //returns poly if v is integer |
---|
286 | ideal I = A.l; |
---|
287 | A = ideal(I[v]); |
---|
288 | return (A); |
---|
289 | } |
---|
290 | |
---|
291 | example |
---|
292 | { |
---|
293 | "EXAMPLE: "; echo = 2; |
---|
294 | ring R = 0,(x,y,z),dp; |
---|
295 | arr A = ideal(x,y,z); |
---|
296 | A[2]; |
---|
297 | intvec v = 1,3; |
---|
298 | A[v]; |
---|
299 | } |
---|
300 | |
---|
301 | // replaces the k-th plane with poly p |
---|
302 | proc arrSet(arr A, int k, poly p) |
---|
303 | "USAGE: arrSet(A, k, p); arr A, int k, poly p; |
---|
304 | RETURN: [arr] Arrangement where the k-th hyperplane is replaced by p. |
---|
305 | NOTE: p must be linear |
---|
306 | KEYWORDS: arrangement; hyperplane; assign; set |
---|
307 | EXAMPLE: example arrSet; shows an example" |
---|
308 | |
---|
309 | { |
---|
310 | if(deg(p) != 1){ERROR("Given poly is not a hyperplane");} |
---|
311 | else{ // looks akward but needs to be done this way |
---|
312 | ideal I = A.l; |
---|
313 | I[k] = p; |
---|
314 | A = I; |
---|
315 | } |
---|
316 | return (A); |
---|
317 | } |
---|
318 | |
---|
319 | example |
---|
320 | { |
---|
321 | "EXAMPLE: "; echo = 2; |
---|
322 | ring R = 0,(x,y,z),dp; |
---|
323 | arr A = ideal(x,y,z); |
---|
324 | arrSet(A,1,x+1); |
---|
325 | } |
---|
326 | |
---|
327 | |
---|
328 | //============================================================================// |
---|
329 | //------------------------------- #4 DELETION --------------------------------// |
---|
330 | //============================================================================// |
---|
331 | |
---|
332 | |
---|
333 | // deletes all hyperplanes of the given indices |
---|
334 | static proc arrDelete(arr A, intvec v) |
---|
335 | "USAGE: delete(A, v); v integer/intvec |
---|
336 | RETURN: [arr] Arrangement A without the hyperplanes given by v; |
---|
337 | NOTE: for deleting hyperplanes via polynomials, use arrMinus instead |
---|
338 | SEE ALSO: arrDelete, arrMinus |
---|
339 | KEYWORDS: arrangement; delete |
---|
340 | EXAMPLE: example arrDelete; shows an example" |
---|
341 | |
---|
342 | { |
---|
343 | int i = 0; int k; |
---|
344 | int n = size(A); |
---|
345 | intvec u = 1..n; |
---|
346 | |
---|
347 | // puts 0 in u for every element that needs to be deleted |
---|
348 | // is done this way to deal with the case that the user gives the same index multiple times. |
---|
349 | for(k=1; k<=size(v); k++){ |
---|
350 | if(u[v[k]] != 0){ u[v[k]] = 0; i++; } // i = #elts that need to be deleted |
---|
351 | } |
---|
352 | |
---|
353 | if( i == n){arr emptyArr; return (emptyArr); } |
---|
354 | |
---|
355 | // create intvec of the remaining indices |
---|
356 | v = 1..(n-i); i=1; |
---|
357 | for(k=1; k<=n; k++){ |
---|
358 | if(u[k] != 0) { v[i] = u[k]; i++; } |
---|
359 | } |
---|
360 | |
---|
361 | arr A' = arrGet(A,v); |
---|
362 | return (A'); |
---|
363 | } |
---|
364 | |
---|
365 | example |
---|
366 | { |
---|
367 | "EXAMPLE: "; echo = 2; |
---|
368 | ring R = 0,(x,y,z),dp; |
---|
369 | arr A = ideal(x,y,z); |
---|
370 | delete(A,2); |
---|
371 | intvec v = (1,3); |
---|
372 | delete(A,v); |
---|
373 | } |
---|
374 | |
---|
375 | |
---|
376 | // deletes hyperplanes, overloads - operator |
---|
377 | static proc arrMinus |
---|
378 | "USAGE: A - #; # list containing arr/ideal/list/matrix/poly |
---|
379 | RETURN: [arr] arrangement A without the hyperplanes of the arrangement defined by #. |
---|
380 | REMARKS: algorithm creates an arrangement B from # using arrAdd and |
---|
381 | then deletes hyperplanes which occur in both A and B. |
---|
382 | NOTE: The algorithm does not simplify by scalars, i.e. some hyperplanes might |
---|
383 | not be deleted. See example. |
---|
384 | SEE ALSO: arrDelete, arrMinus |
---|
385 | KEYWORDS: arrangement; delete; operator |
---|
386 | EXAMPLE: example arrMinus; shows an example" |
---|
387 | |
---|
388 | { |
---|
389 | if(typeof(#[1]) != "arr"){ERROR("First input must be arr!");} |
---|
390 | arr A = #[1]; |
---|
391 | arr B = #[2..size(#)]; // collects hyperplanes to be deleted |
---|
392 | list L; |
---|
393 | int k, l; |
---|
394 | for(k=1; k<=size(A); k++){ // create list of hyperplanes to be deleted |
---|
395 | for(l=1; l<=size(B); l++){ |
---|
396 | if(A[k] == B[l]){L = insert(L,k);} |
---|
397 | } |
---|
398 | } |
---|
399 | |
---|
400 | l = size(L); // transforms list to intvec |
---|
401 | if(l != 0){ |
---|
402 | intvec v = 1..l; |
---|
403 | for(k=1; k<=l; k++){v[k] = L[k];} |
---|
404 | A = delete(A,v); |
---|
405 | } |
---|
406 | |
---|
407 | return (A); |
---|
408 | } |
---|
409 | |
---|
410 | example |
---|
411 | { |
---|
412 | "EXAMPLE: "; echo = 2; |
---|
413 | ring R = 0,(x,y,z),dp; |
---|
414 | arr A = ideal(x,y,z); |
---|
415 | A - ideal(x,y); |
---|
416 | A - poly(2x); |
---|
417 | } |
---|
418 | |
---|
419 | |
---|
420 | //============================================================================// |
---|
421 | //------------------------------- #5 COMPERATORS -----------------------------// |
---|
422 | //============================================================================// |
---|
423 | |
---|
424 | |
---|
425 | // returns logical value of 'A<=B' |
---|
426 | static proc arrLEQ(arr A, arr B) |
---|
427 | "USAGE: A<=B; A,B arr |
---|
428 | RETURN: [0,1] true if A is a subarrangement of B |
---|
429 | NOTE: algorithm is based on arrMinus and does not simplify by scalars, hence some |
---|
430 | technically equal hyperplanes might not be detected. See example. |
---|
431 | SEE ALSO: arrLEQ, arrLNEQ, arrGEQ, arrGNEQ, arrEQ, arrNEQ |
---|
432 | KEYWORDS: comparison; subarrangement; operator |
---|
433 | EXAMPLE: example arrLEQ; shows an example" |
---|
434 | |
---|
435 | { |
---|
436 | arr C = A - B; |
---|
437 | if(C[1] == 0){ return (1); } |
---|
438 | return (0); |
---|
439 | } |
---|
440 | |
---|
441 | example{example arrEQ;} |
---|
442 | |
---|
443 | |
---|
444 | // returns logical value of 'A>=B' |
---|
445 | static proc arrGEQ(arr A, arr B) |
---|
446 | "USAGE: A>=B; A,B arr |
---|
447 | RETURN: [0,1] true if B is a subarrangement of A |
---|
448 | NOTE: algorithm is based on arrMinus and does not simplify by scalars, hence some |
---|
449 | technically equal hyperplanes might not be detected. See example. |
---|
450 | SEE ALSO: arrLEQ, arrLNEQ, arrGEQ, arrGNEQ, arrEQ, arrNEQ |
---|
451 | KEYWORDS: comparison; subarrangement; operator |
---|
452 | EXAMPLE: example arrLEQ; shows an example" |
---|
453 | |
---|
454 | { |
---|
455 | arr C = B - A; |
---|
456 | if(C[1] == 0){ return (1); } |
---|
457 | return (0); |
---|
458 | } |
---|
459 | |
---|
460 | example{example arrEQ;} |
---|
461 | |
---|
462 | |
---|
463 | // returns logical value of 'A==B' |
---|
464 | static proc arrEQ(arr A, arr B) |
---|
465 | "USAGE: A==B; A,B arr |
---|
466 | RETURN: [0,1] true if A equals B |
---|
467 | NOTE: algorithm is based on arrMinus and does not simplify by scalars, hence some |
---|
468 | technically equal hyperplanes might not be detected. See example. |
---|
469 | SEE ALSO: arrLEQ, arrLNEQ, arrGEQ, arrGNEQ, arrEQ, arrNEQ |
---|
470 | KEYWORDS: comparison; equal; operator |
---|
471 | EXAMPLE: example arrLEQ; shows an example" |
---|
472 | |
---|
473 | { |
---|
474 | return ((A<=B) & (A>=B)); |
---|
475 | } |
---|
476 | |
---|
477 | example |
---|
478 | { |
---|
479 | "EXAMPLE: Relationships between a few arrangements."; echo = 2; |
---|
480 | ring R = 0,(x,y,z),dp; |
---|
481 | arr A = ideal(x,y,z); |
---|
482 | arr B = ideal(x,y); |
---|
483 | A<=B; |
---|
484 | A>=B; |
---|
485 | A==B; |
---|
486 | A!=B; |
---|
487 | A<B; |
---|
488 | A>B; |
---|
489 | } |
---|
490 | |
---|
491 | |
---|
492 | // returns logical value of 'A!=B' |
---|
493 | static proc arrNEQ(arr A, arr B) |
---|
494 | "USAGE: A!=B; A,B arr |
---|
495 | RETURN: [0,1] true if A is not equal to B |
---|
496 | NOTE: algorithm is based on arrMinus and does not simplify by scalars, hence some |
---|
497 | technically equal hyperplanes might not be detected. See example. |
---|
498 | SEE ALSO: arrLEQ, arrLNEQ, arrGEQ, arrGNEQ, arrEQ, arrNEQ |
---|
499 | KEYWORDS: comparison; equal; operator |
---|
500 | EXAMPLE: example arrLEQ; shows an example" |
---|
501 | |
---|
502 | { |
---|
503 | return (!(A==B)); |
---|
504 | } |
---|
505 | |
---|
506 | example{example arrEQ;} |
---|
507 | |
---|
508 | |
---|
509 | // returns logical value of 'A<B' |
---|
510 | static proc arrLNEQ(arr A, arr B) |
---|
511 | "USAGE: A<B; A,B arr |
---|
512 | RETURN: [0,1] true if A is a proper subarrangement of B |
---|
513 | NOTE: algorithm is based on arrMinus and does not simplify by scalars, hence some |
---|
514 | technically equal hyperplanes might not be detected. See example. |
---|
515 | SEE ALSO: arrLEQ, arrLNEQ, arrGEQ, arrGNEQ, arrEQ, arrNEQ |
---|
516 | KEYWORDS: comparison; subarrangement; operator |
---|
517 | EXAMPLE: example arrLEQ; shows an example" |
---|
518 | |
---|
519 | { |
---|
520 | return ((A<=B) & (A!=B)); |
---|
521 | } |
---|
522 | |
---|
523 | example{example arrEQ;} |
---|
524 | |
---|
525 | |
---|
526 | // returns logical value of 'A>B' |
---|
527 | static proc arrGNEQ(arr A, arr B) |
---|
528 | "USAGE: A>B; A,B arr |
---|
529 | RETURN: [0,1] true if B is a proper subarrangement of A |
---|
530 | NOTE: algorithm is based on arrMinus and does not simplify by scalars, hence some |
---|
531 | technically equal hyperplanes might not be detected. See example. |
---|
532 | SEE ALSO: arrLEQ, arrLNEQ, arrGEQ, arrGNEQ, arrEQ, arrNEQ |
---|
533 | KEYWORDS: comparison; subarrangement; operator |
---|
534 | EXAMPLE: example arrLEQ; shows an example" |
---|
535 | |
---|
536 | { |
---|
537 | return ((A>=B) & (A!=B)); |
---|
538 | } |
---|
539 | |
---|
540 | example{example arrEQ;} |
---|
541 | |
---|
542 | |
---|
543 | //============================================================================// |
---|
544 | //------------------------------ #6 TYPE CASTING -----------------------------// |
---|
545 | //============================================================================// |
---|
546 | |
---|
547 | |
---|
548 | // TYPE => ARRANGEMENT |
---|
549 | proc type2arr(list #) |
---|
550 | "USAGE: type2arr(#); # def |
---|
551 | RETURN: [arr] Arrangement defined by the input |
---|
552 | NOTE: The procedure tries to cast the input to [arr] using arrAdd |
---|
553 | KEYWORDS: typecasting; arrangement |
---|
554 | EXAMPLE: example tye2arr; shows an example" |
---|
555 | |
---|
556 | { |
---|
557 | return (arrAdd(#)); |
---|
558 | } |
---|
559 | |
---|
560 | example |
---|
561 | { |
---|
562 | "EXAMPLE:"; echo = 2; |
---|
563 | ring R = 0,(x,y,z),dp; |
---|
564 | ideal I = x,y,z; |
---|
565 | typeof(type2arr(I)); |
---|
566 | } |
---|
567 | |
---|
568 | |
---|
569 | // ARRANGEMENT => POLY |
---|
570 | static proc arr2poly(arr A) |
---|
571 | "USAGE: poly(A); arr A |
---|
572 | RETURN: [poly] The defining polynomial which is the product of polynomials occurring in arr |
---|
573 | NOTE: The procedure will automatically simplify the polynomial by scalar multiplication. |
---|
574 | SEE ALSO: arrAdd, arr2poly, arr2mat, arr2list, arr2ideal, type2arr |
---|
575 | KEYWORDS: typecasting; defining polynomial |
---|
576 | EXAMPLE: example arr2poly; shows an example" |
---|
577 | |
---|
578 | { |
---|
579 | poly f = simplify(product(A.l),1); |
---|
580 | if(f == 0){ return (1); } |
---|
581 | return (f); |
---|
582 | } |
---|
583 | |
---|
584 | example |
---|
585 | { |
---|
586 | "EXAMPLE:"; echo = 2; |
---|
587 | ring R = 0,(x,y,z),dp; |
---|
588 | arr A = ideal(x+1, 2x-2, y); |
---|
589 | arr2poly(A); |
---|
590 | } |
---|
591 | |
---|
592 | |
---|
593 | // ARRANGEMENT => LIST |
---|
594 | static proc arr2list(arr A) |
---|
595 | "USAGE: arr2list(A); A arr |
---|
596 | RETURN: [list] containing all generators of A |
---|
597 | SEE ALSO: arrAdd, arr2poly, arr2mat, arr2list, arr2ideal, type2arr |
---|
598 | KEYWORDS: typecasting; list |
---|
599 | EXAMPLE: example arr2list; shows an example" |
---|
600 | |
---|
601 | { |
---|
602 | int n = size(A); |
---|
603 | list L; |
---|
604 | for(int k=1; k<=n; k++){L[k] = A[k];} |
---|
605 | |
---|
606 | return (L); |
---|
607 | } |
---|
608 | |
---|
609 | example |
---|
610 | { |
---|
611 | "EXAMPLE:"; echo = 2; |
---|
612 | ring R = 0,(x,y,z),dp; |
---|
613 | arr A= ideal(x,y,z); |
---|
614 | arr2list(A); |
---|
615 | } |
---|
616 | |
---|
617 | |
---|
618 | // ARRANGEMENT => IDEAL |
---|
619 | static proc arr2ideal(arr A) |
---|
620 | "USAGE: arr2ideal(A); A arr |
---|
621 | RETURN: [ideal], the internal ideal of A |
---|
622 | NOTE: arr2ideal(A); is the same as A.l - which may become private |
---|
623 | SEE ALSO: arrAdd, arr2poly, arr2mat, arr2list, arr2ideal, type2arr |
---|
624 | KEYWORDS: typecasting; ideal |
---|
625 | EXAMPLE: example arr2ideal; shows an example" |
---|
626 | |
---|
627 | { |
---|
628 | return (A.l); |
---|
629 | } |
---|
630 | |
---|
631 | example |
---|
632 | { |
---|
633 | "EXAMPLE:"; echo = 2; |
---|
634 | ring R = 0,(x,y,z),dp; |
---|
635 | arr A = ideal(x,y,z); |
---|
636 | arr2ideal(A); |
---|
637 | } |
---|
638 | |
---|
639 | |
---|
640 | // ARRANGEMENT => MATRIX |
---|
641 | static proc arr2mat(arr A, list #) |
---|
642 | "USAGE: matrix(A); A arr |
---|
643 | matrix(A, 'c') for central arrangements (shorter matrix) |
---|
644 | RETURN: [matrix] M = [T|b] representing the arrangement |
---|
645 | NOTE: If the arrangement is central or one is not interested in the const |
---|
646 | terms one can use "matrix(A, 'c')" instead to get the same matrix without |
---|
647 | the last column. |
---|
648 | SEE ALSO: arr2mat, mat2arr, mat2carr |
---|
649 | KEYWORDS: typecasting; matrix; coefficient |
---|
650 | EXAMPLE: example arr2mat;" |
---|
651 | |
---|
652 | { |
---|
653 | matrix M = jacob(A.l); |
---|
654 | if(size(#) == 0){return ( concat(M, transpose(jet(A.l,0))) );} |
---|
655 | if(#[1] == 'c'){return (M);} |
---|
656 | ERROR("Bad optional input parameter!") |
---|
657 | } |
---|
658 | |
---|
659 | example |
---|
660 | { |
---|
661 | "EXAMPLE:"; echo = 2; |
---|
662 | ring R = 0,(x,y,z),dp; |
---|
663 | arr A = ideal(x,y,z); |
---|
664 | arr D = ideal(x+1, y-2, z, x+y+4); |
---|
665 | print(arr2mat(A)); |
---|
666 | print(arr2mat(D)); |
---|
667 | } |
---|
668 | |
---|
669 | |
---|
670 | // MATRIX => ARRANGEMENT |
---|
671 | proc mat2arr(matrix M) |
---|
672 | "USAGE: mat2arr(M); matrix (M|b) |
---|
673 | RETURN: [arr] interprets the rows of the matrix as the defining polynomial equations |
---|
674 | of the arrangement where the last column will be considered as the constant |
---|
675 | terms, i.e. if M is an m*(n+1) matrix we have |
---|
676 | H_i = Ker( M_i1*x_1 +...+ M_in*x_n + M_i(n+1) ) for i=1...m and |
---|
677 | A = {H_1,...,H_m} the resulting arrangement. |
---|
678 | SEE ALSO: mat2carr |
---|
679 | KEYWORDS: typecasting; matrix; coefficient |
---|
680 | EXAMPLE: example mat2arr;" |
---|
681 | |
---|
682 | { |
---|
683 | if(ncols(M) > nvars(basering)+1) |
---|
684 | { |
---|
685 | ERROR("Matrix too big! Please add variables to basering."); |
---|
686 | } |
---|
687 | |
---|
688 | int n = ncols(M)-1; |
---|
689 | matrix X[n+1][1]; |
---|
690 | X[1..n,1] = varMat(1..n); |
---|
691 | X[n+1 ,1] = 1; |
---|
692 | arr A = ideal(M*X); |
---|
693 | |
---|
694 | return(A); |
---|
695 | } |
---|
696 | |
---|
697 | example |
---|
698 | { |
---|
699 | "EXAMPLE:"; echo = 2; |
---|
700 | ring R = 0,(x,y,z),dp; |
---|
701 | matrix M[4][4] = 1,0,1,1,1,1,0,2,0,1,1,3,2,1,1,4; |
---|
702 | print(M); |
---|
703 | mat2arr(M); |
---|
704 | } |
---|
705 | |
---|
706 | |
---|
707 | // MATRIX => ARRANGEMENT (central) |
---|
708 | proc mat2carr(matrix M) |
---|
709 | "USAGE: mat2carr(M); matrix M |
---|
710 | RETURN: [arr] interprets the rows of the matrix as the defining polynomial equations |
---|
711 | of the arrangement. I.e. if M is an m*n matrix we have |
---|
712 | H_i = Ker( M_i1*x_1 +...+ M_in*x_n) for i=1...m and |
---|
713 | A = {H_1,...,H_m} the resulting arrangement. |
---|
714 | SEE ALSO: mat2arr |
---|
715 | KEYWORDS: typecasting; matrix; coefficient; central |
---|
716 | EXAMPLE: example mat2carr;" |
---|
717 | |
---|
718 | { |
---|
719 | if( ncols(M) > nvars(basering) ){ |
---|
720 | ERROR("Error! not enough variables in the basering."); |
---|
721 | } |
---|
722 | |
---|
723 | arr A = ideal(M*varMat(1..ncols(M))); |
---|
724 | |
---|
725 | return(A); |
---|
726 | } |
---|
727 | |
---|
728 | example |
---|
729 | { |
---|
730 | "EXAMPLE:"; echo = 2; |
---|
731 | ring R = 0,(x,y,z),dp; |
---|
732 | matrix M[4][3] = 1,0,1,1,1,0,0,1,1,2,1,1; |
---|
733 | print(M); |
---|
734 | mat2carr(M); |
---|
735 | } |
---|
736 | |
---|
737 | |
---|
738 | //============================================================================// |
---|
739 | //---------------------------- #7 IDEAL FUNCTIONS ----------------------------// |
---|
740 | //============================================================================// |
---|
741 | |
---|
742 | |
---|
743 | // checks if A is central, homogenizes |
---|
744 | static proc arrHomog(arr A, list #) |
---|
745 | "USAGE: homog(A); arr A |
---|
746 | homog(A, p); arr A, ring_variable p |
---|
747 | RETURN: [0,1] homog(A) is the same as arrCentral(A) |
---|
748 | [arr] homog(A,p) homogenizes A with respect to p |
---|
749 | NOTE: homog(A,p) is not the same as arrCone(A,p) as it does not add the additional hyperplane |
---|
750 | SEE ALSO: arrHomog, arrCentral, arrCone |
---|
751 | KEYWORDS: central; homogenize |
---|
752 | EXAMPLE: example arrHomog; shows an example" |
---|
753 | |
---|
754 | { |
---|
755 | if(size(#) == 0){ |
---|
756 | return (homog(A.l)); |
---|
757 | } |
---|
758 | |
---|
759 | if(size(#) == 1){ |
---|
760 | A = homog(A.l, #[1]); |
---|
761 | return (A); |
---|
762 | } |
---|
763 | |
---|
764 | ERROR("Too many innput arguments!"); |
---|
765 | |
---|
766 | } |
---|
767 | |
---|
768 | example |
---|
769 | { |
---|
770 | "EXAMPLE: "; echo = 2; |
---|
771 | ring R = 0,(x,y,z),dp; |
---|
772 | arr A = ideal(x,y); |
---|
773 | arr B = ideal(x,y,x+y+1); |
---|
774 | homog(A); |
---|
775 | homog(B); |
---|
776 | homog(B,z); |
---|
777 | homog(_); |
---|
778 | } |
---|
779 | |
---|
780 | |
---|
781 | |
---|
782 | // simplified arrangement |
---|
783 | static proc arrSimplify(arr A, list #) |
---|
784 | "USAGE: arrSimplify(A); |
---|
785 | simplify(A, n); arr A, int n |
---|
786 | RETURN: [arr] simplified arrangement. |
---|
787 | NOTE: arrSimplify(A) is the same as simplify(A, 1), simplify with higher ints |
---|
788 | is not needed |
---|
789 | SEE ALSO: arrSimplify |
---|
790 | KEYWORDS: simplify |
---|
791 | EXAMPLE: example arrSimplify; shows an example" |
---|
792 | |
---|
793 | { |
---|
794 | if(size(#) == 0){ |
---|
795 | A = simplify(A.l, 1); |
---|
796 | return (A); |
---|
797 | } |
---|
798 | |
---|
799 | if(size(#) == 1){ |
---|
800 | A = simplify(A.l, #[1]); |
---|
801 | return (A); |
---|
802 | } |
---|
803 | |
---|
804 | ERROR("Too many input arguments!"); |
---|
805 | } |
---|
806 | |
---|
807 | example |
---|
808 | { |
---|
809 | "EXAMPLE: "; echo = 2; |
---|
810 | ring R = 0,(x,y,z),dp; |
---|
811 | arr A = ideal(2x,2y-1,2z-2); |
---|
812 | arrSimplify(A); |
---|
813 | } |
---|
814 | |
---|
815 | |
---|
816 | // number of hyperplanes |
---|
817 | static proc arrSize(arr A) |
---|
818 | "USAGE: size(A); arr A; |
---|
819 | RETURN: [int] number of hyperplanes in the arrangement |
---|
820 | NOTE: size(A) also useable for multi-arrangements |
---|
821 | SEE ALSO: arrSize, multarrSize |
---|
822 | KEYWORDS: hyperplanes; size; number |
---|
823 | EXAMPLE: example arrSize; shows an example" |
---|
824 | |
---|
825 | { |
---|
826 | return (size(A.l)); |
---|
827 | } |
---|
828 | |
---|
829 | example |
---|
830 | { |
---|
831 | "EXAMPLE: "; echo = 2; |
---|
832 | ring R = 0,(x,y,z),dp; |
---|
833 | arr A = ideal(x,y,z); |
---|
834 | size(A); |
---|
835 | } |
---|
836 | |
---|
837 | |
---|
838 | // substitute variables |
---|
839 | static proc arrSubst(arr A, list #) |
---|
840 | "USAGE: arrSubst(A, #); arr A, ring_variables/integers i,j; |
---|
841 | RETURN: [arr] with the corresponding substitutions made |
---|
842 | NOTE: applies 'subst' on the arrangement |
---|
843 | SEE ALSO: arrSubst |
---|
844 | KEYWORDS: variables; ring_variable; substitute |
---|
845 | EXAMPLE: example arrSubst; shows an example" |
---|
846 | |
---|
847 | { |
---|
848 | if(size(#)%2 != 0){ |
---|
849 | ERROR("Odd number of parameter inputs!"); |
---|
850 | } |
---|
851 | |
---|
852 | for(int i=1; i<size(#); i=i+2){ |
---|
853 | A = subst(A.l, #[i], #[i+1]); |
---|
854 | } |
---|
855 | |
---|
856 | return (A); |
---|
857 | } |
---|
858 | |
---|
859 | example |
---|
860 | { |
---|
861 | "EXAMPLE: "; echo = 2; |
---|
862 | ring R = 0,(x,y,z),dp; |
---|
863 | arr A = ideal(x,y,z,x+z); |
---|
864 | arrSubst(A,x,y); |
---|
865 | arrSubst(A,x,y,y,z); |
---|
866 | } |
---|
867 | |
---|
868 | |
---|
869 | //============================================================================// |
---|
870 | //------------------------------- #8 PRINTING --------------------------------// |
---|
871 | //============================================================================// |
---|
872 | |
---|
873 | |
---|
874 | // prints arrangement in the console |
---|
875 | static proc arrPrint(arr A) |
---|
876 | "USAGE: A; A arr |
---|
877 | RETURN: [] better readable output in the console as the newstruct print |
---|
878 | SEE ALSO: arrPrint, arrPrintMatrix |
---|
879 | KEYWORDS: print |
---|
880 | EXAMPLE: example arrPrint;" |
---|
881 | |
---|
882 | { |
---|
883 | A.l; |
---|
884 | } |
---|
885 | |
---|
886 | example |
---|
887 | { |
---|
888 | "EXAMPLE:"; echo = 2; |
---|
889 | ring R = 0,(x,y,z),dp; |
---|
890 | arr A = ideal(x,y,z); |
---|
891 | A; |
---|
892 | } |
---|
893 | |
---|
894 | |
---|
895 | // readable as a coeff matrix |
---|
896 | proc arrPrintMatrix(arr A) |
---|
897 | "USAGE: arPrintMatrix(arr A) |
---|
898 | RETURN: [] prints arr in matrix form |
---|
899 | NOTE: differs print(matrix(arr A)) since variables included |
---|
900 | KEYWORDS: print; matrix; coefficient |
---|
901 | EXAMPLE: example arrPrintMatrix;" |
---|
902 | |
---|
903 | { |
---|
904 | matrix M = matrix(A); |
---|
905 | ideal I = variables(A); |
---|
906 | for(int k=1; k<=size(I); k++){ |
---|
907 | M[1..nrows(M),k] = (M[1..nrows(M),k])*(I[k]); |
---|
908 | } |
---|
909 | |
---|
910 | print(M); |
---|
911 | } |
---|
912 | |
---|
913 | example |
---|
914 | { |
---|
915 | "EXAMPLE:"; echo = 2; |
---|
916 | ring R = 0,(x,y,z),dp; |
---|
917 | arr A = mat2arr(random(20,5,4)); |
---|
918 | A; |
---|
919 | arrPrintMatrix(A); |
---|
920 | } |
---|
921 | |
---|
922 | |
---|
923 | //============================================================================// |
---|
924 | //------------------------- #9 MANIPULATING VARIABLES ------------------------// |
---|
925 | //============================================================================// |
---|
926 | |
---|
927 | |
---|
928 | // matrix of the corresponding ring_variables |
---|
929 | proc varMat(intvec v) |
---|
930 | "USAGE: varMat(v); v intvec |
---|
931 | RETURN: [matrix] M containing the corresponding ring_variables |
---|
932 | SEE ALSO: varMat, varNum, arrSwapVar, arrLastVar |
---|
933 | KEYWORDS: variables; ring_variable |
---|
934 | EXAMPLE: example varMat; shows an example" |
---|
935 | |
---|
936 | { |
---|
937 | matrix M[size(v)][1]; |
---|
938 | for(int k=1; k<=size(v); k++){ |
---|
939 | M[k,1] = var(v[k]); |
---|
940 | } |
---|
941 | return (M); |
---|
942 | } |
---|
943 | |
---|
944 | example |
---|
945 | { |
---|
946 | "EXAMPLE: 'even' ringvariables"; echo = 2; |
---|
947 | ring R = 0,(x(1..6)),dp; |
---|
948 | intvec v = 2,4,6; |
---|
949 | varMat(v); |
---|
950 | } |
---|
951 | |
---|
952 | |
---|
953 | // number of given variable (enh. version of varNum in dmod.lib) |
---|
954 | proc varNum(def u) |
---|
955 | "USAGE: varnum(string s); |
---|
956 | varnum(ring_variable); |
---|
957 | RETURN: [int] number of given ring variable, or 0 if it does not appear |
---|
958 | NOTE: This procedure has the same functionality as varNum from the dmod.lib |
---|
959 | package, but also accepts polys as input. |
---|
960 | SEE ALSO: varMat, varNum, arrSwapVar, arrLastVar |
---|
961 | KEYWORDS: variables; ring_variable; number |
---|
962 | EXAMPLE: example varNum; shows an example" |
---|
963 | |
---|
964 | { |
---|
965 | if(typeof(u) == "string"){ |
---|
966 | for(int i=1; i<=nvars(basering); i++){ |
---|
967 | if(string(var(i)) == u){return (i);} |
---|
968 | } |
---|
969 | return (0); |
---|
970 | } |
---|
971 | |
---|
972 | if(typeof(u) == "poly"){ |
---|
973 | for(int i=1; i<=nvars(basering); i++){ |
---|
974 | if(var(i) == u){return (i);} |
---|
975 | } |
---|
976 | return (0); |
---|
977 | } |
---|
978 | |
---|
979 | ERROR("Wrong input type, expected string or ring_variable (poly)!"); |
---|
980 | } |
---|
981 | |
---|
982 | example |
---|
983 | { |
---|
984 | "EXAMPLE: "; echo = 2; |
---|
985 | ring R = 0,(x,y,z),dp; |
---|
986 | varNum(y); |
---|
987 | ring S = 0,(x(1..5),y(1..5)),dp; |
---|
988 | varNum("y(3)"); |
---|
989 | } |
---|
990 | |
---|
991 | |
---|
992 | // ideal of all variables the arrangement depends on |
---|
993 | static proc arrVariables(arr A) |
---|
994 | "USAGE: variables(A); A arr |
---|
995 | RETURN: [ideal] whereas generators are variables A uses |
---|
996 | NOTE: inherited from the ideal class |
---|
997 | SEE ALSO: varMat, varNum, arrVariables ,arrNvars, arrSwapVar, arrLastVar |
---|
998 | KEYWORDS: variables; ring_variable |
---|
999 | EXAMPLE: example arrVariables; shows an example" |
---|
1000 | |
---|
1001 | { |
---|
1002 | return (variables(A.l)); |
---|
1003 | } |
---|
1004 | |
---|
1005 | example |
---|
1006 | { |
---|
1007 | "EXAMPLE: "; echo = 2; |
---|
1008 | ring R = 0,(x,y,z),lp; |
---|
1009 | arr A = ideal(x,y,z); |
---|
1010 | variables(A); |
---|
1011 | variables(A-y); |
---|
1012 | } |
---|
1013 | |
---|
1014 | |
---|
1015 | // number of variables the arrangement uses |
---|
1016 | static proc arrNvars(arr A) |
---|
1017 | "USAGE: arrNvars(A); A arr |
---|
1018 | RETURN: [int] number of variables A uses |
---|
1019 | NOTE: inherited from the ideal class |
---|
1020 | SEE ALSO: varMat, varNum, arrVariables ,arrNvars, arrSwapVar, arrLastVar |
---|
1021 | KEYWORDS: variables; variables; ring_variable; number |
---|
1022 | EXAMPLE: example arrNvars; shows an example" |
---|
1023 | |
---|
1024 | { |
---|
1025 | return (size(variables(A))); |
---|
1026 | } |
---|
1027 | |
---|
1028 | example |
---|
1029 | { |
---|
1030 | "EXAMPLE: "; echo = 2; |
---|
1031 | ring R = 0,(x,y,z),lp; |
---|
1032 | arr A= ideal(x,y,z); |
---|
1033 | arrNvars(A); |
---|
1034 | arrNvars(A-y); |
---|
1035 | } |
---|
1036 | |
---|
1037 | |
---|
1038 | // swaps two variables in the arrangement |
---|
1039 | proc arrSwapVar(arr A, def i, def j) |
---|
1040 | "USAGE: arrSwapVar(A, i, j); arr A, ring_variables/integers i,j |
---|
1041 | RETURN: [arr] A where variables i and j are swapped. |
---|
1042 | NOTE: if i and/or j are integers the algorithm considers the variables |
---|
1043 | variables(A)[i] and/or variables(A)[j] |
---|
1044 | SEE ALSO: varMat, varNum, arrSwapVar, arrLastVar |
---|
1045 | KEYWORDS: swap; variables; ring_variable |
---|
1046 | EXAMPLE: example arrSwapVar; shows an example" |
---|
1047 | |
---|
1048 | { |
---|
1049 | ideal I = variables(A); |
---|
1050 | poly u,v; |
---|
1051 | if(typeof(i) == "int"){ u = I[i]; } else{ u = i; } |
---|
1052 | if(typeof(j) == "int"){ v = I[j]; } else{ v = j; } |
---|
1053 | |
---|
1054 | if(u == v){ return (A); } // special case which messes up the rest |
---|
1055 | |
---|
1056 | // using the old trick on how to swap 2 cells without needing a third: |
---|
1057 | // (a|b) =(a->a+b)=> (a+b|b) =(b->b-a)=> (b|b-a) =(a->b-a)=> (b|a) |
---|
1058 | A = subst(A.l, u, u+v); |
---|
1059 | A = subst(A.l, v, v-u); |
---|
1060 | A = subst(A.l, u, v-u); |
---|
1061 | |
---|
1062 | return (A); |
---|
1063 | } |
---|
1064 | |
---|
1065 | example |
---|
1066 | { |
---|
1067 | "EXAMPLE: "; echo = 2; |
---|
1068 | ring R = 0,(x,y,z),lp; |
---|
1069 | arr A = ideal(x+1,x+y,z); |
---|
1070 | arrSwapVar(A,x,z); |
---|
1071 | } |
---|
1072 | |
---|
1073 | |
---|
1074 | //ring_variable of largest index used in arrangement |
---|
1075 | proc arrLastVar(arr A) |
---|
1076 | "USAGE: arrLastVar(A); arr A |
---|
1077 | RETURN: [int] number of the last variable A uses |
---|
1078 | NOTE: useful if you want a list containing all variables x_1 ... x_k used in A, |
---|
1079 | but you do not want to skip any like variables(A) does. |
---|
1080 | SEE ALSO: varMat, varNum, arrSwapVar, arrLastVar |
---|
1081 | KEYWORDS: variables; ring_variable |
---|
1082 | EXAMPLE: example arrLastVar; shows an example" |
---|
1083 | |
---|
1084 | { |
---|
1085 | return ( rvar(variables(A)[arrNvars(A)]) ); |
---|
1086 | } |
---|
1087 | |
---|
1088 | example |
---|
1089 | { |
---|
1090 | "EXAMPLE: "; echo = 2; |
---|
1091 | ring R = 0,x(1..10),dp; |
---|
1092 | arr A = ideal(x(1), x(2), x(3), x(6)); |
---|
1093 | int n = arrLastVar(A); |
---|
1094 | varMat(1..n); |
---|
1095 | variables(A); |
---|
1096 | } |
---|
1097 | |
---|
1098 | |
---|
1099 | //============================================================================// |
---|
1100 | //--------------------------- #10 CENTER COMPUTATIONS ------------------------// |
---|
1101 | //============================================================================// |
---|
1102 | |
---|
1103 | |
---|
1104 | // checks if arr is central |
---|
1105 | proc arrCentral(arr A) |
---|
1106 | "USAGE: arrCentral(A); arr A |
---|
1107 | RETURN: [0,1] true if arr is central(i.e. all planes intersect in 0) |
---|
1108 | NOTE: This is the same as homog(A) |
---|
1109 | SEE ALSO: arrCentered, arrCentral, arrCenter, arrCentralize |
---|
1110 | KEYWORDS: center; central |
---|
1111 | EXAMPLE: example arrCentral;" |
---|
1112 | |
---|
1113 | { |
---|
1114 | return (homog(A)); |
---|
1115 | } |
---|
1116 | |
---|
1117 | example |
---|
1118 | { |
---|
1119 | "EXAMPLE:"; echo = 2; |
---|
1120 | ring R = 0,(x,y,z),dp; |
---|
1121 | |
---|
1122 | // centered and central |
---|
1123 | arr A = ideal(x,y,z); |
---|
1124 | arrCentered(A); |
---|
1125 | arrCentral(A); |
---|
1126 | |
---|
1127 | // centered but not central (center: (-1,-1/2, 1)) |
---|
1128 | arr B = ideal(x+1,2y+1,-z+1); |
---|
1129 | arrCentered(B); |
---|
1130 | arrCentral(B); |
---|
1131 | } |
---|
1132 | |
---|
1133 | |
---|
1134 | // checks whether arrangement has a center |
---|
1135 | proc arrCentered(arr A) |
---|
1136 | "USAGE: arrCentered(A); arr A |
---|
1137 | RETURN: [0,1] true if A is centered(i.e. intersection of all planes not empty) |
---|
1138 | NOTE: The algorithm uses the rank of matrix: Ax=b has a solution iff |
---|
1139 | rank(A) = rank(A|b) |
---|
1140 | SEE ALSO: arrCentered, arrCentral, arrCenter, arrCentralize |
---|
1141 | KEYWORDS: center |
---|
1142 | EXAMPLE: example arrCentered;" |
---|
1143 | |
---|
1144 | { |
---|
1145 | matrix M = matrix(A); |
---|
1146 | // classic test: Ax=b has a solution if and only if rank(A|b) == rank(A) |
---|
1147 | if( rank(M) == rank(submat(M,1..nrows(M), 1..ncols(M)-1))){ return (1); } |
---|
1148 | return (0); |
---|
1149 | } |
---|
1150 | |
---|
1151 | example |
---|
1152 | { |
---|
1153 | "EXAMPLE:"; echo = 2; |
---|
1154 | ring R = 0,(x,y,z),dp; |
---|
1155 | arr A= ideal(x,y,x-y+1); // centerless |
---|
1156 | arrCentral(A); |
---|
1157 | arr B= ideal(x,y,z); // central with center being the origin |
---|
1158 | arrCentral(B); |
---|
1159 | arr C= ideal(x+1,2y+1,-z+1); // central with center (-1,-1/2, 1) |
---|
1160 | arrCentral(C); |
---|
1161 | } |
---|
1162 | |
---|
1163 | |
---|
1164 | // computes center of an arrangement |
---|
1165 | proc arrCenter(arr A) |
---|
1166 | "USAGE: arrCenter(A); arr A |
---|
1167 | RETURN: [list] L entry 0 if A not centered or entries 1, x, H, where x is |
---|
1168 | any particular point of the center and H is a matrix consisting of |
---|
1169 | vectors which spanning linear intersection space. |
---|
1170 | If there is exactly one solution, then H = 0. |
---|
1171 | NOTE: The intersection of all hyperplanes can be expressed in forms of a |
---|
1172 | linear system Ax=b, where (A|b) is the coeff. matrix of the arrange- |
---|
1173 | ment, which is then solved using L-U decomposition |
---|
1174 | SEE ALSO: arrCentered, arrCentral, arrCenter, arrCentralize |
---|
1175 | KEYWORDS: center |
---|
1176 | EXAMPLE: example arrCenter;" |
---|
1177 | |
---|
1178 | { |
---|
1179 | matrix M = matrix(A); //return matrix (T|b) |
---|
1180 | matrix T = submat(M, 1..nrows(M), 1..ncols(M)-1); |
---|
1181 | matrix b = submat(M, 1..nrows(M), ncols(M))*(-1); |
---|
1182 | list L = ludecomp(T); |
---|
1183 | list Q = lusolve(L[1], L[2], L[3], b); |
---|
1184 | return (Q); |
---|
1185 | } |
---|
1186 | |
---|
1187 | example |
---|
1188 | { |
---|
1189 | "EXAMPLE:"; echo = 2; |
---|
1190 | ring R = 0,(x,y,z),dp; |
---|
1191 | |
---|
1192 | arr A= ideal(x,y,x-y+1); // centerless |
---|
1193 | arrCenter(A); |
---|
1194 | |
---|
1195 | arr B= ideal(x,y,z); // center is a single point |
---|
1196 | arrCenter(B); |
---|
1197 | |
---|
1198 | arr C= ideal(x,z,x+z); // center is a line |
---|
1199 | // here we get a wrong result because the matrix is simplified since A doesn't |
---|
1200 | // contain any "y" the matrix (A|b) will be 3x3 only. |
---|
1201 | arrCenter(C); |
---|
1202 | } |
---|
1203 | |
---|
1204 | |
---|
1205 | // makes centered arrangement central |
---|
1206 | proc arrCentralize(arr A) |
---|
1207 | "USAGE: arrCenteralize(A); arr A |
---|
1208 | RETURN: [arr] A after centralization via coordinate change |
---|
1209 | NOTE: The coordinate change only does translation, vector of translation is the second |
---|
1210 | output of arrCenter |
---|
1211 | SEE ALSO: arrCentered, arrCentral, arrCenter, arrCentralize |
---|
1212 | KEYWORDS: central; center; coordinate change |
---|
1213 | EXAMPLE: example arrCenter;" |
---|
1214 | |
---|
1215 | { |
---|
1216 | if(arrCentral(A)){ |
---|
1217 | print("The arrangement is already central!"); |
---|
1218 | return (); |
---|
1219 | } |
---|
1220 | |
---|
1221 | list L = arrCenter(A); |
---|
1222 | |
---|
1223 | if(L[1] == 0){ |
---|
1224 | print("The arrangement has no center and therefor cannot be centralized!"); |
---|
1225 | return (); |
---|
1226 | } |
---|
1227 | |
---|
1228 | int n = nvars(basering); |
---|
1229 | |
---|
1230 | matrix T = diag(1,n); |
---|
1231 | matrix c = L[2]; |
---|
1232 | |
---|
1233 | A = arrCoordChange(A, T, c); |
---|
1234 | |
---|
1235 | return (A); |
---|
1236 | } |
---|
1237 | |
---|
1238 | example |
---|
1239 | { |
---|
1240 | "EXAMPLE:"; echo = 2; |
---|
1241 | ring R = 0,(x,y,z),dp; |
---|
1242 | arr A = ideal(x-1,y,x-z-1,x-z-1); |
---|
1243 | arrCentralize(A); |
---|
1244 | } |
---|
1245 | |
---|
1246 | |
---|
1247 | //============================================================================// |
---|
1248 | //------------------------ #11 GEOMETRIC CONSTRUCTIONS -----------------------// |
---|
1249 | //============================================================================// |
---|
1250 | |
---|
1251 | // performs coordinate change |
---|
1252 | proc arrCoordChange(arr A, matrix T, list #) |
---|
1253 | "USAGE: arrCoordChange(A, T); arr A, (m*n mat) n*n or n*n+1 matrix T |
---|
1254 | arrCoordChange(A, T , c); arr A, n*n matrix T, n*1 matrix/vector |
---|
1255 | RETURN: [arr]: Arrangement A [A|b] after a coordinate change f: x -> Tx + c with T invertible |
---|
1256 | i.e. [A|b] => [AT^-1|b+AT^-1c] since we have |
---|
1257 | f(H) = f(ker(a1*x1 + ... + an*xn - b)) = {f(x) : a'x -b = 0} |
---|
1258 | = {y : a'f^-1(y) -b = 0)} |
---|
1259 | = {y : a'(T^-1(y-c)) - b = 0} |
---|
1260 | = {y : a'T^-1y -(b + a'T^-1c) = 0} |
---|
1261 | NOTE: There are 3 options how you can give the input (in each case n <= nvars(basering)) |
---|
1262 | 1. Just a nxn matrix with |
---|
1263 | => Will automatically complete T by a unit matrix and perform x -> Tx |
---|
1264 | 2. A nxn matrix T and a nx1 vector/matrix c with |
---|
1265 | => Will automatically complete T and c and perform x -> Tx +c |
---|
1266 | 3. A nxn+1 matrix T with |
---|
1267 | => will use last column as translation vector c |
---|
1268 | SEE ALSO: arrCoordChange, arrCoordNormalize |
---|
1269 | KEYWORDS: coordinate change |
---|
1270 | EXAMPLE: example arrCoordChange; shows an example" |
---|
1271 | |
---|
1272 | { |
---|
1273 | int n = nvars(basering); |
---|
1274 | int k = nrows(T); |
---|
1275 | int l = ncols(T); |
---|
1276 | matrix c; |
---|
1277 | |
---|
1278 | if( k>n || l-1>n ){ |
---|
1279 | ERROR("Matrix too big! (It has more cols/rows than there are variables.)"); |
---|
1280 | } |
---|
1281 | |
---|
1282 | // const vector integrated in matrix => split [T|c] into T and c |
---|
1283 | if( l == k+1 ){ |
---|
1284 | if(size(#) > 0){ |
---|
1285 | ERROR("Bad input. If given a constant vector the matrix must be square!") |
---|
1286 | } |
---|
1287 | c = submat(T, 1..k, l); |
---|
1288 | T = submat(T, 1..k, 1..k); l=k; |
---|
1289 | } |
---|
1290 | |
---|
1291 | if( l != k || rank(T) != k){ |
---|
1292 | ERROR("Given matrix is not a base change matrix.") |
---|
1293 | } |
---|
1294 | |
---|
1295 | if(size(#) > 0){ |
---|
1296 | if(size(#) >= 2){ERROR("Too many input arguments!");} |
---|
1297 | c = matrix(#[1]); |
---|
1298 | if(nrows(c) > n || ncols(c) > 1){ |
---|
1299 | ERROR("Constant vector maldefined. Dimension too big.") |
---|
1300 | } |
---|
1301 | |
---|
1302 | } |
---|
1303 | |
---|
1304 | matrix M[n][n] = diag(1,n); |
---|
1305 | M[1..k,1..k] = T; |
---|
1306 | T = M; |
---|
1307 | T = luinverse(T)[2]; |
---|
1308 | M = jacob(A.l); |
---|
1309 | |
---|
1310 | matrix b[n][1]; |
---|
1311 | b[1..nrows(c),1] = c; |
---|
1312 | c = b; // gives c the right length. |
---|
1313 | b = transpose(jet(A.l, 0)); // constants |
---|
1314 | |
---|
1315 | b = b - M*T*c; |
---|
1316 | M = M*T; |
---|
1317 | A = mat2arr(concat(M,b)); |
---|
1318 | |
---|
1319 | return(A); |
---|
1320 | } |
---|
1321 | |
---|
1322 | example |
---|
1323 | { |
---|
1324 | "EXAMPLE: "; echo =2; |
---|
1325 | ring r = 0,(x,y,z),lp; |
---|
1326 | arr A = x,y,z; |
---|
1327 | arrCoordChange(A,1,[0,0,1]); //lifts z-hyperplane by 1 unit |
---|
1328 | matrix T[2][2] = [0,1,1,0]; // swaps x and y |
---|
1329 | arrCoordChange(A,T); |
---|
1330 | matrix c[2][1] = [1,0]; |
---|
1331 | T = concat(T,c); // now swap x and y and add 1 to x afterwards |
---|
1332 | arrCoordChange(A,T); |
---|
1333 | // Note how T doesn't even need to be a full 3x3 base change matrix. |
---|
1334 | } |
---|
1335 | |
---|
1336 | // performs projection onto coordinate hyperplanes |
---|
1337 | proc arrCoordNormalize(arr A, intvec v) |
---|
1338 | "USAGE: arrCoordChange(A, v); |
---|
1339 | RETURN: [arr]: Arrangement after a coordinate change that transforms the arrangement such that |
---|
1340 | after a transformation x -> Tx + c we have the arrangement has the matrix representation |
---|
1341 | [AT^-1|b+AT^-1c] such that [AT^-1]_v = I and [b+AT^-1c]_v = 0; |
---|
1342 | NOTE: algorithm performs a base change if H_k is homogeneous (i.e. has no) |
---|
1343 | constant term and an affine transformation otherwise |
---|
1344 | Ax+b = 0, Transformation x = Ty+c: AT^-1y + AT^-1c + b = 0 |
---|
1345 | Now we want to have (AT^-1)_v = I and (AT^-1c +b)_v = AT^-1_v*c + b_v = 0 |
---|
1346 | SEE ALSO: arrCoordChange, arrCoordNormalize |
---|
1347 | KEYWORDS: coordinate change |
---|
1348 | EXAMPLE: example arrCoordChange; shows an example" |
---|
1349 | |
---|
1350 | { |
---|
1351 | int n = nvars(basering); |
---|
1352 | if(size(v) > n){ |
---|
1353 | ERROR("Too many rows chosen, at max you can choose "+string(n)); |
---|
1354 | } |
---|
1355 | |
---|
1356 | matrix M = matrix(A); |
---|
1357 | matrix Av[n][n]; |
---|
1358 | matrix bv[n][1]; |
---|
1359 | Av[1..size(v), 1..n] = submat(M,v,1..n); |
---|
1360 | bv[1..size(v),1] = submat(M,v, n+1); |
---|
1361 | |
---|
1362 | if(rank(Av) != n){ |
---|
1363 | if(rank(Av) != size(v)){ |
---|
1364 | ERROR("Normal vectors of the given hyperplanes are not linearly " + |
---|
1365 | "independent. Cannot perform coordinate change!"); |
---|
1366 | } |
---|
1367 | // Adds linearly independent lines to Av in order to make it invertible |
---|
1368 | module F = freemodule(n); |
---|
1369 | module Add = reduce(F,std(transpose(Av))); |
---|
1370 | Add = transpose(simplify(transpose(Add),2)); |
---|
1371 | Av[(size(v)+1)..n, 1..n] = matrix(Add); |
---|
1372 | } |
---|
1373 | if(rank(Av) != n){ERROR("Av not invertible!");} |
---|
1374 | |
---|
1375 | A = arrCoordChange(A,Av,bv); |
---|
1376 | |
---|
1377 | return (A); |
---|
1378 | } |
---|
1379 | example |
---|
1380 | { |
---|
1381 | "EXAMPLE: "; echo=2; |
---|
1382 | ring r = 0,(x,y,z),lp; |
---|
1383 | arr A = ideal(x,y,z,x+z+4); |
---|
1384 | intvec v = 1,2,4; |
---|
1385 | arrCoordNormalize(A,v); |
---|
1386 | } |
---|
1387 | |
---|
1388 | // coned arrangement |
---|
1389 | proc arrCone(arr A, list #) |
---|
1390 | "USAGE: arrCone(A); |
---|
1391 | arrCone(A, ring_variable); arr A arrangement in variables x_1...x_n; |
---|
1392 | RETURN: arr, the coned hyperplane Arrangement cA with respect to the given |
---|
1393 | ring_variable, or the last ring_variable if none was given. |
---|
1394 | NOTE: The hyperplanes are homogenized w.r.t. v and a new hyperplane |
---|
1395 | H = ker(x_n+1) is added. |
---|
1396 | SEE ALSO: arrCone, arrDecone, arrRestrict, arrIsEssential, arrEssentialize |
---|
1397 | KEYWORDS: cone; decone |
---|
1398 | EXAMPLE: example arrCone; shows an example" |
---|
1399 | |
---|
1400 | { |
---|
1401 | poly p; int i; |
---|
1402 | |
---|
1403 | // case 1: no ring_variable given |
---|
1404 | if(size(#) == 0){ |
---|
1405 | p = var(nvars(basering)); |
---|
1406 | } |
---|
1407 | |
---|
1408 | // case 2: a ring_variable is given |
---|
1409 | else{ p = #[1]; } |
---|
1410 | |
---|
1411 | // computation |
---|
1412 | A = homog(A, p); |
---|
1413 | A = A + p; |
---|
1414 | |
---|
1415 | return(A); |
---|
1416 | } |
---|
1417 | |
---|
1418 | example |
---|
1419 | { |
---|
1420 | "EXAMPLE: Coning the arrangement A = (x+1, y) alongside x, y and z:"; echo=2; |
---|
1421 | ring R = 0,(x,y,z),dp; |
---|
1422 | arr A = ideal(x+1, x,x-2,x-1); |
---|
1423 | arrCone(A, y); |
---|
1424 | arr B= ideal(x,y,x+y-1); |
---|
1425 | arrCone(B); |
---|
1426 | } |
---|
1427 | |
---|
1428 | |
---|
1429 | // deconed arrangement |
---|
1430 | proc arrDecone(arr A, int k) |
---|
1431 | "USAGE: arrDecone(A, k); arrangement A, integer k; |
---|
1432 | RETURN: arr: the deconed hyperplane Arrangement dA |
---|
1433 | NOTE: A has to be non-empty and central. arrDecone is an inverse operation |
---|
1434 | to arrCone since A == arrDecone(arrCone(A),size(A)+1) for any A. |
---|
1435 | One can also decone a central arrangement with respect to any hyper- |
---|
1436 | plane k, but than a coordinate change is necessary to make |
---|
1437 | H_k = ker(x_k). Since such a coordinate change is not unique, |
---|
1438 | use arrCoordchange to do so. |
---|
1439 | SEE ALSO: arrCone, arrDecone, arrRestrict, arrIsEssential, arrEssentialize |
---|
1440 | KEYWORDS: cone; decone |
---|
1441 | EXAMPLE: example arrDecone; shows an example" |
---|
1442 | |
---|
1443 | { |
---|
1444 | if( !arrCentral(A) ){ |
---|
1445 | ERROR("Non-central arrangement can not be deconed!");} |
---|
1446 | if( size(A) == 0){ |
---|
1447 | ERROR("Empty arrangement can not be deconed!");} |
---|
1448 | if( size(A) < k){ |
---|
1449 | ERROR("There is no k-th hyperplane");} |
---|
1450 | |
---|
1451 | poly p = A[k]; |
---|
1452 | int n = rvar(p); |
---|
1453 | |
---|
1454 | if( n == 0 ){ ERROR("H_" + string(k) + " = " + string(p) + |
---|
1455 | " is not of the form ker(x_i). Please do a coordinate change first." + |
---|
1456 | "You can use arrCoordinateChange to transform the arrangement accordingly."); |
---|
1457 | } |
---|
1458 | |
---|
1459 | A = A - p; |
---|
1460 | A = arrSubst(A, p, 1); |
---|
1461 | |
---|
1462 | return(A); |
---|
1463 | } |
---|
1464 | |
---|
1465 | example |
---|
1466 | { |
---|
1467 | "EXAMPLE: We decone arr consisting of (x,y,x+z) with respect to y"; |
---|
1468 | echo = 2; |
---|
1469 | ring R = 0,(x,y,z),dp; |
---|
1470 | arr A= ideal(x,y,z,x+y-z); |
---|
1471 | arrDecone(A,3); |
---|
1472 | } |
---|
1473 | |
---|
1474 | proc arrLocalize(arr A, intvec v) |
---|
1475 | "USAGE: arrLocalize(A, v); arrangement A, intvec v; |
---|
1476 | RETURN: arr: the localized arrangement A_X, i.e. A_X only contains the hyperplanes |
---|
1477 | which contain the flat X, which is defined by the equations A[v] |
---|
1478 | SEE ALSO: arrCone, arrDecone, arrRestrict, arrIsEssential, arrEssentialize |
---|
1479 | KEYWORDS: localization |
---|
1480 | EXAMPLE: example arrLocalize; shows an example" |
---|
1481 | { |
---|
1482 | ideal I=std(arr2ideal(A[v])); |
---|
1483 | arr L; |
---|
1484 | int k; |
---|
1485 | for(k=1; k<size(A); k++){ |
---|
1486 | if(NF(A[k],I)==0){ |
---|
1487 | L=L+A[k]; |
---|
1488 | } |
---|
1489 | } |
---|
1490 | return(L); |
---|
1491 | } |
---|
1492 | |
---|
1493 | example |
---|
1494 | { |
---|
1495 | "EXAMPLE: We localize the Braid 3 arrangement at x and y"; |
---|
1496 | echo = 2; |
---|
1497 | ring R = 0,(x,y,z),dp; |
---|
1498 | arr A=arrTypeB(3); |
---|
1499 | intvec v=5,8; |
---|
1500 | arr B=arrLocalize(A,v); |
---|
1501 | B; |
---|
1502 | } |
---|
1503 | |
---|
1504 | // restricted arrangement onto a hyperplane |
---|
1505 | proc arrRestrict(arr A, intvec v, list #) |
---|
1506 | "USAGE: arrRestrict(A, v); arrangement A, int/intvec v, optional argument "CC"; |
---|
1507 | RETURN: arr: the restricted hyperplane Arrangement (A^X) |
---|
1508 | NOTE: A has to be non-empty. |
---|
1509 | REMARKS: We restrict A to the flat X, defined by the equations in A[v]. |
---|
1510 | The restriction will only be performed, if the ideal defining |
---|
1511 | the flat X is monomial (i.e. X is an intersection of coordinate planes). |
---|
1512 | If the optional argument CC is given, the arrangement is transformed |
---|
1513 | in such a way that X has the above form. |
---|
1514 | SEE ALSO: arrCone, arrDecone, arrRestrict, arrIsEssential, arrEssentialize |
---|
1515 | KEYWORDS: restriction |
---|
1516 | EXAMPLE: example arrRestrict; shows an example" |
---|
1517 | |
---|
1518 | { |
---|
1519 | ideal I=A.l; |
---|
1520 | I=I[v]; |
---|
1521 | option(redSB); |
---|
1522 | I=std(I); // defining equations for flat X |
---|
1523 | //option(none); |
---|
1524 | ideal Al=arr2ideal(A); |
---|
1525 | int i; |
---|
1526 | |
---|
1527 | if(isMonomial(I)==1){ |
---|
1528 | for(i=1;i<= size(I);i++){ |
---|
1529 | Al=subst(Al,I[i],0); |
---|
1530 | } |
---|
1531 | arr AR=Al; |
---|
1532 | return(AR); |
---|
1533 | } |
---|
1534 | |
---|
1535 | if(size(#) == 0){ |
---|
1536 | ERROR("The flat X is not defined by a monomial ideal. " + |
---|
1537 | "Please do a coordinate change first. You can use arrCoordNormalize to " + |
---|
1538 | "transform the arrangement accordingly by adding the argument CC."); |
---|
1539 | } |
---|
1540 | if(#[1]!="CC"){ |
---|
1541 | ERROR("The flat X is not defined by a monomial ideal. " + |
---|
1542 | "Please do a coordinate change first. You can use arrCoordNormalize to " + |
---|
1543 | "transform the arrangement accordingly by adding the argument CC."); |
---|
1544 | } |
---|
1545 | if(size(v)==0){return(A);} |
---|
1546 | intvec w=v[1]; |
---|
1547 | intvec tmp; |
---|
1548 | |
---|
1549 | for(i=2;i<=size(v);i++){ |
---|
1550 | tmp=w,v[i]; |
---|
1551 | if(rank(matrix(A[w]))==size(tmp)){ |
---|
1552 | w=tmp; |
---|
1553 | } |
---|
1554 | } |
---|
1555 | |
---|
1556 | arr B=arrCoordNormalize(A,w); |
---|
1557 | ideal Bl=arr2ideal(B); |
---|
1558 | |
---|
1559 | for(i=1;i<= size(w);i++){ |
---|
1560 | Bl=subst(Bl,Bl[w[i]],0); |
---|
1561 | } |
---|
1562 | arr BR=Bl; |
---|
1563 | return(BR); |
---|
1564 | } |
---|
1565 | |
---|
1566 | example |
---|
1567 | { |
---|
1568 | "EXAMPLE: We consider the possible restrictions of the type B3 arrangement "; |
---|
1569 | echo = 2; |
---|
1570 | ring S = 0,(x,y,z),dp; |
---|
1571 | arr A = arrTypeB(3); |
---|
1572 | A; |
---|
1573 | arrRestrict(A,9); |
---|
1574 | arrRestrict(A,4,"CC"); |
---|
1575 | intvec v=5,8; |
---|
1576 | arrRestrict(A,v); |
---|
1577 | } |
---|
1578 | |
---|
1579 | |
---|
1580 | // checks if arrangement is essential |
---|
1581 | proc arrIsEssential(arr A) |
---|
1582 | "USAGE: arrIsEssential(A); arrangement A; |
---|
1583 | RETURN: boolean: 1 if arr is essential, i.e. rank of maximal element of |
---|
1584 | poset is dimension |
---|
1585 | NOTE: A has to be non-empty. |
---|
1586 | SEE ALSO: arrCone, arrDecone, arrRestrict, arrIsEssential, arrEssentialize |
---|
1587 | KEYWORDS: essential |
---|
1588 | EXAMPLE: example arrIsEssential; shows an example" |
---|
1589 | |
---|
1590 | { |
---|
1591 | ideal I = variables(A); |
---|
1592 | |
---|
1593 | if( var(size(I)) != I[size(I)] ){ |
---|
1594 | return(0); |
---|
1595 | } if( rank(jacob(A.l)) == size(I) ){ |
---|
1596 | return(1); |
---|
1597 | } |
---|
1598 | |
---|
1599 | return(0); |
---|
1600 | } |
---|
1601 | |
---|
1602 | example |
---|
1603 | { |
---|
1604 | "EXAMPLE: We check whether these are essential and a non-essential arrangement "; |
---|
1605 | echo = 2; |
---|
1606 | ring S = 0,(x,y,z),lp; |
---|
1607 | arr A = ideal(x,y,z); |
---|
1608 | arr B = ideal(x+y+z,x,y+z); |
---|
1609 | arrIsEssential(A); |
---|
1610 | arrIsEssential(B); |
---|
1611 | } |
---|
1612 | |
---|
1613 | |
---|
1614 | // essentialized arragnement |
---|
1615 | proc arrEssentialize(arr A) |
---|
1616 | "USAGE: arrEssentialize(A); arrangement A; |
---|
1617 | RETURN: essential arrangement by transformation |
---|
1618 | NOTE: A has to be non-empty. |
---|
1619 | SEE ALSO: arrCone, arrDecone, arrRestrict, arrIsEssential, arrEssentialize |
---|
1620 | KEYWORDS: essential |
---|
1621 | EXAMPLE: example arrEssentialize; shows an example" |
---|
1622 | |
---|
1623 | { |
---|
1624 | ideal I = variables(A); |
---|
1625 | |
---|
1626 | if( var(size(I)) != I[size(I)] ){ |
---|
1627 | for(int k=1; k<=size(I); k++){ |
---|
1628 | if( NF(var(k),std(I)) != 0 ){ |
---|
1629 | A = subst(A.l, I[size(I)], var(k)); |
---|
1630 | return (arrEssentialize(A)); |
---|
1631 | } |
---|
1632 | } |
---|
1633 | } |
---|
1634 | |
---|
1635 | while(arrIsEssential(A) == 0){ |
---|
1636 | A = subst(A.l, I[size(I)], 0); |
---|
1637 | } |
---|
1638 | |
---|
1639 | return (A); |
---|
1640 | } |
---|
1641 | |
---|
1642 | example |
---|
1643 | { |
---|
1644 | "EXAMPLE: We essentialize a non essential arrangement "; echo = 2; |
---|
1645 | ring S = 0,(x,y,z),dp; |
---|
1646 | arr A=arrBraid(3); |
---|
1647 | arrEssentialize(A); |
---|
1648 | arr B = ideal(x+y+z,x,y+z); |
---|
1649 | arrEssentialize(B); |
---|
1650 | } |
---|
1651 | |
---|
1652 | //============================================================================// |
---|
1653 | //------------------------ #12 EXAMPLES OF ARRANGEMENTS ----------------------// |
---|
1654 | //============================================================================// |
---|
1655 | |
---|
1656 | |
---|
1657 | // boolean arrangement |
---|
1658 | proc arrBoolean(int v) |
---|
1659 | "USAGE: arrBoolean(v); int v |
---|
1660 | RETURN: arr, which uses the first v variables of ring for boolean arrangement |
---|
1661 | SEE ALSO: arrBoolean, arrBraid, arrTypeB, arrTypeD, arrRandom, arrEdelmanReiner |
---|
1662 | KEYWORDS: example; boolean |
---|
1663 | EXAMPLE: example arrBoolean;" |
---|
1664 | |
---|
1665 | { |
---|
1666 | if(nvars(basering)<v){ return("not enough variables"); } |
---|
1667 | |
---|
1668 | arr A = mat2carr(unitmat(v)); |
---|
1669 | return(A); |
---|
1670 | } |
---|
1671 | |
---|
1672 | example |
---|
1673 | { |
---|
1674 | "EXAMPLE:Boolean Arrangement, uses the seven first coordinate hyperplanes"; |
---|
1675 | echo = 2; |
---|
1676 | ring R = 0,x(1..10),dp; |
---|
1677 | arrBoolean(7); |
---|
1678 | } |
---|
1679 | |
---|
1680 | |
---|
1681 | // braid arrangement |
---|
1682 | proc arrBraid(int v) |
---|
1683 | "USAGE: arrBraid(v); int v |
---|
1684 | RETURN: Type A (braid) arrangement of dimension v |
---|
1685 | SEE ALSO: arrBoolean, arrBraid, arrTypeB, arrTypeD, arrRandom, arrEdelmanReiner |
---|
1686 | KEYWORDS: example; braid |
---|
1687 | EXAMPLE: example arrBraid;" |
---|
1688 | |
---|
1689 | { |
---|
1690 | if(nvars(basering)<v){ |
---|
1691 | return("not enough variables"); |
---|
1692 | } |
---|
1693 | |
---|
1694 | arr A; |
---|
1695 | int i,j; |
---|
1696 | for(i=1; i<=v; i++){ |
---|
1697 | for(j=i+1; j<=v; j++){ |
---|
1698 | A = A + (var(i) - var(j)); |
---|
1699 | } |
---|
1700 | } |
---|
1701 | return(A); |
---|
1702 | } |
---|
1703 | |
---|
1704 | example |
---|
1705 | { |
---|
1706 | "EXAMPLE: The braid arrangement consists of the hyperplanes x(i)=x(j)"; |
---|
1707 | echo = 2; |
---|
1708 | ring R = 0,x(1..10),dp; |
---|
1709 | arrBraid(7); |
---|
1710 | } |
---|
1711 | |
---|
1712 | |
---|
1713 | // type B arrangement |
---|
1714 | proc arrTypeB(int v) |
---|
1715 | "USAGE: arrTypeB(v); int v |
---|
1716 | RETURN: arrangement, which uses first v variables of ring for reflection |
---|
1717 | arrangement of type B |
---|
1718 | SEE ALSO: arrBoolean, arrBraid, arrTypeB, arrTypeD, arrRandom, arrEdelmanReiner |
---|
1719 | KEYWORDS: example; type B |
---|
1720 | EXAMPLE: example arrTypeB;" |
---|
1721 | |
---|
1722 | { |
---|
1723 | if(nvars(basering)<v){ |
---|
1724 | return("not enough variables"); |
---|
1725 | } |
---|
1726 | |
---|
1727 | arr A; |
---|
1728 | int i,j; |
---|
1729 | for(i=1; i<=v; i++){ |
---|
1730 | for(j=i+1; j<=v; j++){ |
---|
1731 | A = A + (var(i) - var(j)); |
---|
1732 | A = A + (var(i) + var(j)); |
---|
1733 | |
---|
1734 | } |
---|
1735 | A = A + var(i); |
---|
1736 | } |
---|
1737 | return(A); |
---|
1738 | } |
---|
1739 | |
---|
1740 | example |
---|
1741 | { |
---|
1742 | "EXAMPLE: The type B reflection arrangement consists of the hyperplanes " + |
---|
1743 | "x(i)=x(j), x(i)=-x(j) and the coordinate hyperplanes"; |
---|
1744 | echo = 2; |
---|
1745 | ring R = 0,x(1..10),dp; |
---|
1746 | arrTypeB(5); |
---|
1747 | } |
---|
1748 | |
---|
1749 | |
---|
1750 | // type D arrangement |
---|
1751 | proc arrTypeD(int v) |
---|
1752 | "USAGE: arrTypeD(v); int v |
---|
1753 | RETURN: arrangement, which uses first v variables of ring for reflection |
---|
1754 | arrangement of type D |
---|
1755 | SEE ALSO: arrBoolean, arrBraid, arrTypeB, arrTypeD, arrRandom, arrEdelmanReiner |
---|
1756 | KEYWORDS: example; type D |
---|
1757 | EXAMPLE: example arrTypeD;" |
---|
1758 | |
---|
1759 | { |
---|
1760 | if(nvars(basering)<v){ |
---|
1761 | return("not enough variables"); |
---|
1762 | } |
---|
1763 | |
---|
1764 | arr A; |
---|
1765 | int i,j; |
---|
1766 | for(i=1; i<=v; i++){ |
---|
1767 | for(j=i+1; j<=v; j++){ |
---|
1768 | A = A + (var(i) - var(j)); |
---|
1769 | A = A + (var(i) + var(j)); |
---|
1770 | } |
---|
1771 | } |
---|
1772 | return(A); |
---|
1773 | } |
---|
1774 | |
---|
1775 | example |
---|
1776 | { |
---|
1777 | "EXAMPLE: The type D reflection arrangement consists of the hyperplanes " + |
---|
1778 | "x(i)=x(j) and x(i)=-x(j)"; |
---|
1779 | echo = 2; |
---|
1780 | ring R = 0,x(1..10),dp; |
---|
1781 | arrTypeD(5); |
---|
1782 | } |
---|
1783 | |
---|
1784 | |
---|
1785 | // random (affine) arrangement |
---|
1786 | proc arrRandom(int d, int m, int n) |
---|
1787 | "USAGE: arrRandom(n,v,N); int n,v,N |
---|
1788 | RETURN: Random arrangement, where m is the number of hyperplanes, n the |
---|
1789 | dimension, d the upper bound for absolute value of coefficients. |
---|
1790 | NOTE: You can also write arr = random(d,m,n) to create random arrangements |
---|
1791 | SEE ALSO: arrBoolean, arrBraid, arrTypeB, arrTypeD, arrRandom, arrEdelmanReiner |
---|
1792 | KEYWORDS: example; random |
---|
1793 | EXAMPLE: example arrRandom;" |
---|
1794 | |
---|
1795 | { |
---|
1796 | if(n > nvars(basering)){ |
---|
1797 | return("Error, too few variables or too high dimension"); |
---|
1798 | } |
---|
1799 | |
---|
1800 | intmat M = random(d,m,n+1); |
---|
1801 | arr A = mat2arr(M); |
---|
1802 | return(A); |
---|
1803 | } |
---|
1804 | |
---|
1805 | example |
---|
1806 | { |
---|
1807 | "EXAMPLE:"; echo = 2; |
---|
1808 | ring R = 0,x(1..20),dp; |
---|
1809 | arrRandom(7,3,15); |
---|
1810 | } |
---|
1811 | |
---|
1812 | |
---|
1813 | // random central arrangement |
---|
1814 | proc arrRandomCentral(int d, int m, int n) |
---|
1815 | "USAGE: arrRandomCentral(d,m,n); int d,m,n |
---|
1816 | RETURN: Random central arrangement, where m is the number of hyperplanes, n |
---|
1817 | the dimension, d the upper bound for absolute value of coefficients. |
---|
1818 | SEE ALSO: arrBoolean, arrBraid, arrTypeB, arrTypeD, arrRandom, arrEdelmanReiner |
---|
1819 | KEYWORDS: example; random; central |
---|
1820 | EXAMPLE: example arrRandomCentral;" |
---|
1821 | |
---|
1822 | { |
---|
1823 | if(n > nvars(basering)){ |
---|
1824 | return("Error, too few variables or too high dimension"); |
---|
1825 | } |
---|
1826 | |
---|
1827 | intmat M = random(d,m,n); |
---|
1828 | arr A = mat2carr(M); |
---|
1829 | return(A); |
---|
1830 | } |
---|
1831 | |
---|
1832 | example |
---|
1833 | { |
---|
1834 | "EXAMPLE:"; echo = 2; |
---|
1835 | ring R = 0,x(1..20),dp; |
---|
1836 | arrRandomCentral(7,3,15); |
---|
1837 | } |
---|
1838 | |
---|
1839 | |
---|
1840 | // Edelman-Reiner arrangement |
---|
1841 | proc arrEdelmanReiner() |
---|
1842 | "USAGE: arrEdelmanReiner(); |
---|
1843 | RETURN: the Edelman-Reiner arrangement, which is a free arrangement but the |
---|
1844 | restriction to the 6-th hyperplane is nonfree. |
---|
1845 | (i.e. counterexample for Orlik-Conjecture) |
---|
1846 | NOTE: the active ring must have at least five variables |
---|
1847 | SEE ALSO: arrBoolean, arrBraid, arrTypeB, arrTypeD, arrRandom, arrEdelmanReiner |
---|
1848 | KEYWORDS: example; Edelman-Reiner |
---|
1849 | EXAMPLE: example arrEdelmanReiner;" |
---|
1850 | |
---|
1851 | { |
---|
1852 | if(nvars(basering) < 5){ ERROR("not enough variables"); } |
---|
1853 | |
---|
1854 | arr arrER = arrBoolean(5); |
---|
1855 | int a,b,c,d; |
---|
1856 | for(a=1; a<=2; a++){ |
---|
1857 | for(b=1; b<=2; b++){ |
---|
1858 | for(c=1; c<=2; c++){ |
---|
1859 | for(d=1; d<=2; d++){ |
---|
1860 | arrER = arrER + (var(1)+(-1)^a*var(2)+(-1)^b*var(3)+(-1)^c*var(4)+(-1)^d*var(5)); |
---|
1861 | }}}} |
---|
1862 | |
---|
1863 | return(arrER); |
---|
1864 | } |
---|
1865 | |
---|
1866 | example |
---|
1867 | { |
---|
1868 | "EXAMPLE:"; echo = 2; |
---|
1869 | ring r=0,x(1..5),dp; |
---|
1870 | arrEdelmanReiner(); |
---|
1871 | } |
---|
1872 | |
---|
1873 | |
---|
1874 | //============================================================================// |
---|
1875 | //--------------------- #13 Orlik Solomon and Poincare Poly ------------------// |
---|
1876 | //============================================================================// |
---|
1877 | |
---|
1878 | |
---|
1879 | // Orlik-Solomon algebra of the arrangement |
---|
1880 | proc arrOrlikSolomon(arr A) |
---|
1881 | "USAGE: arrOrlikSolomon(A); arr A |
---|
1882 | RETURN: [ring] exterior Algebra E as ring with Orlik-Solomon ideal as attribute I. |
---|
1883 | The Orlik-Solomon ideal is generated by the differentials of dependent |
---|
1884 | tuples of hyperplanes. For a complex arrangement the quotient E/I is |
---|
1885 | isomorphic to the cohomology ring of the complement of the arrangement. |
---|
1886 | NOTE: In order to access this ideal I activate this exterior algebra with setring. |
---|
1887 | SEE ALSO: arrOrlikSolomon |
---|
1888 | KEYWORDS: Orlik-Solomon |
---|
1889 | EXAMPLE: example arrOrlikSolomon;" |
---|
1890 | |
---|
1891 | { |
---|
1892 | int central = arrCentral(A); // 0: non-central, 1: central |
---|
1893 | if(central == 0){ |
---|
1894 | A = arrCone(A); |
---|
1895 | } |
---|
1896 | |
---|
1897 | module M = syz(A.l); |
---|
1898 | M = jet(M,0); //Only use linear syzygies |
---|
1899 | M = simplify(M,2); // drop zeros |
---|
1900 | |
---|
1901 | int n = ncols(A.l); |
---|
1902 | def startRing = basering; |
---|
1903 | ring R = 0,e(1..n),dp; |
---|
1904 | def ER = Exterior(); |
---|
1905 | setring ER; // defines the Exterioralgebra |
---|
1906 | ideal I; //final orlik solomon ideal |
---|
1907 | |
---|
1908 | matrix X = transpose(varMat(1..n)); |
---|
1909 | module M = fetch(startRing,M); //brings the module M to the Exterior Algebra |
---|
1910 | ideal OSI = ideal(X*M); |
---|
1911 | |
---|
1912 | if(size(OSI) == 0){ // no relations among hyperplanes, I==0 |
---|
1913 | export(I); |
---|
1914 | return(basering); |
---|
1915 | } |
---|
1916 | |
---|
1917 | // monomial subideal procedure due to Sorin Popescu |
---|
1918 | ideal K = OSI; |
---|
1919 | ideal J; |
---|
1920 | while(isMonomial(K) == 0){ |
---|
1921 | J = lead(std(K)); |
---|
1922 | K = intersect(J,OSI); |
---|
1923 | K = interred(K); |
---|
1924 | } |
---|
1925 | |
---|
1926 | OSI = K; |
---|
1927 | poly diffOp; //now applying the differential operator |
---|
1928 | ideal vars; |
---|
1929 | int j; |
---|
1930 | |
---|
1931 | for(int i=1; i<=ncols(OSI); i++) |
---|
1932 | { |
---|
1933 | vars = variables(OSI[i]); |
---|
1934 | diffOp = 0; |
---|
1935 | for(j=1; j<=ncols(vars); j++) |
---|
1936 | { |
---|
1937 | diffOp = diffOp+(-1)^(j+1)*vars[j]; |
---|
1938 | } |
---|
1939 | |
---|
1940 | I = I + ideal( diff(diffOp,ideal(OSI[i])) ); |
---|
1941 | } |
---|
1942 | |
---|
1943 | if(central ==0 ){//project back in the non-central case |
---|
1944 | n = nvars(basering)-1; |
---|
1945 | ring R = 0,e(1..n),dp; |
---|
1946 | def ERDC = Exterior(); |
---|
1947 | setring ERDC; |
---|
1948 | ideal I = fetch(ER,I); |
---|
1949 | } |
---|
1950 | |
---|
1951 | export(I); |
---|
1952 | return (basering); |
---|
1953 | } |
---|
1954 | |
---|
1955 | example |
---|
1956 | { |
---|
1957 | "EXAMPLE: Computing the Orlik-Solomon-Ideal for the D3-Arrangement"; echo = 2; |
---|
1958 | ring R = 0,(x,y,z),dp; |
---|
1959 | arr A = arrTypeB(3); |
---|
1960 | def E = arrOrlikSolomon(A); |
---|
1961 | setring E; |
---|
1962 | //The generators of the Orlik-Solomon-Ideal are: |
---|
1963 | I; |
---|
1964 | } |
---|
1965 | |
---|
1966 | |
---|
1967 | // The following 3 procedures were replaced by their faster combinatorial counterparts |
---|
1968 | |
---|
1969 | /* |
---|
1970 | |
---|
1971 | // poincare polynomial of the arrangement |
---|
1972 | proc arrPoincare(arr A) |
---|
1973 | "USAGE: arrPoincare(A); arr A |
---|
1974 | RETURN: [intvec] The Poincare polynomial as integer vector of the arrangement, which |
---|
1975 | is equal to the second kind Poincare-Series of the Orlik-Solomon Algebra. |
---|
1976 | SEE ALSO: arrPoincare, arrChambers, arrBoundedChambers |
---|
1977 | KEYWORDS: Poincare polynomial |
---|
1978 | EXAMPLE: example arrPoincare;" |
---|
1979 | |
---|
1980 | { |
---|
1981 | def startRing = basering; |
---|
1982 | def Ext = arrOrlikSolomon(A); |
---|
1983 | setring Ext; |
---|
1984 | ideal OSI = lead(I); |
---|
1985 | OSI = OSI + ideal(Ext); |
---|
1986 | int n = nvars(Ext); |
---|
1987 | ring suppRing = 0,x(1..n),dp; |
---|
1988 | ideal I = fetch(Ext,OSI); |
---|
1989 | intvec HP = hilb(std(I),2); |
---|
1990 | |
---|
1991 | return(HP); |
---|
1992 | } |
---|
1993 | |
---|
1994 | example |
---|
1995 | { |
---|
1996 | "EXAMPLE: Computing the Poincare polynomial as intvec for the D3-Arrangement"; echo = 2; |
---|
1997 | ring R = 0,(x,y,z),dp; |
---|
1998 | arr A = arrTypeB(3); |
---|
1999 | //The coefficients of the Poincare polynomial are: |
---|
2000 | arrPoincare(A); |
---|
2001 | } |
---|
2002 | |
---|
2003 | |
---|
2004 | // number of chambers of the arrangement |
---|
2005 | proc arrChambers(arr A) |
---|
2006 | "USAGE: arrChambers(A); arr A |
---|
2007 | RETURN: [int] The number of chambers of an arrangement, which is equal to the |
---|
2008 | evaluation of the Poincare polynomial at 1. |
---|
2009 | SEE ALSO: arrPoincare, arrChambers, arrBoundedChambers |
---|
2010 | KEYWORDS: chambers |
---|
2011 | EXAMPLE: example arrChambers;" |
---|
2012 | |
---|
2013 | { |
---|
2014 | intvec HP = arrPoincare(A); |
---|
2015 | intvec ones = 1:(size(HP)); |
---|
2016 | return ( transpose(ones)*HP ); |
---|
2017 | } |
---|
2018 | |
---|
2019 | example |
---|
2020 | { |
---|
2021 | "EXAMPLE: Computing the number of chambers for the D3-Arrangement"; echo = 2; |
---|
2022 | ring R = 0,(x,y,z),dp; |
---|
2023 | arr A = arrTypeD(3); |
---|
2024 | //The number of chambers of the D3-Arrangement is: |
---|
2025 | arrChambers(A); |
---|
2026 | } |
---|
2027 | |
---|
2028 | |
---|
2029 | // number of bounded chambers of the arrangement |
---|
2030 | proc arrBoundedChambers(arr A) |
---|
2031 | "USAGE: arrBoundedChambers(A); arr A |
---|
2032 | RETURN: [int] The number of bounded chambers of an arrangement, which is equal to |
---|
2033 | the evaluation of the Poincare polynomial at -1. |
---|
2034 | SEE ALSO: arrPoincare, arrChambers, arrBoundedChambers |
---|
2035 | KEYWORDS: chambers |
---|
2036 | EXAMPLE: example arrBoundedChambers;" |
---|
2037 | |
---|
2038 | { |
---|
2039 | intvec HP = arrPoincare(A); |
---|
2040 | intvec altOnes; |
---|
2041 | for(int i=1; i<=size(HP); i++){ |
---|
2042 | altOnes[i] = (-1)^(i-1); |
---|
2043 | } |
---|
2044 | |
---|
2045 | return ( transpose(altOnes)*HP ); |
---|
2046 | } |
---|
2047 | |
---|
2048 | example |
---|
2049 | { |
---|
2050 | "EXAMPLE: Computing the number of chambers for the D3-Arrangement"; echo = 2; |
---|
2051 | ring R = 0,(x,y,z),dp; |
---|
2052 | arr A = arrTypeD(3); |
---|
2053 | // The number of bounded chambers of the D3-Arrangement is: |
---|
2054 | arrBoundedChambers(A); |
---|
2055 | } |
---|
2056 | |
---|
2057 | */ |
---|
2058 | |
---|
2059 | //============================================================================// |
---|
2060 | //------------------------------- #14 Freeness -------------------------------// |
---|
2061 | //============================================================================// |
---|
2062 | |
---|
2063 | |
---|
2064 | // module of derivation |
---|
2065 | proc arrDer |
---|
2066 | "USAGE: arrDer(A); arr A , multarr A |
---|
2067 | RETURN: [module] The module Der(A) of derivations of the (multi-)arrangement A, i.e. |
---|
2068 | the derivations tangent to each hyperplane of A (resp. with multiplicities) |
---|
2069 | NOTE: This is only defined for central (multi-)arrangements |
---|
2070 | SEE ALSO: arrDer, arrIsFree, arrExponents |
---|
2071 | KEYWORDS: derivation; multiarrangement |
---|
2072 | EXAMPLE: example arrDer;" |
---|
2073 | |
---|
2074 | { |
---|
2075 | def A = #[1]; |
---|
2076 | if( (typeof(A) != "arr") && ( typeof(A) != "multarr") ){ |
---|
2077 | ERROR("bad input type!"); |
---|
2078 | } if(homog(A.l) == 0){ |
---|
2079 | ERROR("Arrangement not central!"); |
---|
2080 | } |
---|
2081 | |
---|
2082 | module fA = jacob(A.l); |
---|
2083 | ideal J; |
---|
2084 | |
---|
2085 | if( typeof(A) == "arr" ){ |
---|
2086 | J = arr2ideal(A); |
---|
2087 | } if(typeof(A) == "multarr"){ |
---|
2088 | J = multIdeal(A); |
---|
2089 | } |
---|
2090 | |
---|
2091 | module tA = diag(matrix(J)); |
---|
2092 | module K = modulo(fA,tA); |
---|
2093 | module derivations = lessGenerators(K); |
---|
2094 | |
---|
2095 | return (derivations); |
---|
2096 | } |
---|
2097 | |
---|
2098 | example |
---|
2099 | { |
---|
2100 | "EXAMPLE: Computing the derivation module of the boolean and braid arrangement"; |
---|
2101 | echo = 2; |
---|
2102 | ring R = 0,(x,y,z),dp; |
---|
2103 | arr A3 = arrBoolean(3); |
---|
2104 | arr B3 = arrTypeB(3); |
---|
2105 | arr G = ideal(x,y,z,x+y+z); |
---|
2106 | //The derivation module of the Boolean 3-arrangement: |
---|
2107 | arrDer(A3); |
---|
2108 | //The derivation module of the Braid 3-arrangement: |
---|
2109 | arrDer(B3); |
---|
2110 | //The derivation module of the generic arrangement: |
---|
2111 | arrDer(G); |
---|
2112 | } |
---|
2113 | |
---|
2114 | |
---|
2115 | // checks if arrangement is free |
---|
2116 | proc arrIsFree(list #) |
---|
2117 | "USAGE: arrIsFree(A); arr A, multarr A |
---|
2118 | RETURN: [0,1] 1 if the (multi-)arrangement is free, i.e. Der(A) is a free module |
---|
2119 | NOTE: only defined for central arrangements |
---|
2120 | SEE ALSO: arrDer, arrIsFree, arrExponents |
---|
2121 | KEYWORDS: free; multiarrangement |
---|
2122 | EXAMPLE: example arrIsFree;" |
---|
2123 | |
---|
2124 | { |
---|
2125 | def A = #[1]; |
---|
2126 | module derivations = arrDer(A); |
---|
2127 | return ( nvars(basering) == ncols(derivations) ); |
---|
2128 | } |
---|
2129 | |
---|
2130 | example |
---|
2131 | { |
---|
2132 | "EXAMPLE: checking freeness of the Edelman-Reiner arrangement and its restriction: "; |
---|
2133 | echo = 2; |
---|
2134 | ring R = 0,(x,y,z),dp; |
---|
2135 | arr A3 = arrBoolean(3); |
---|
2136 | arr B3 = arrTypeB(3); |
---|
2137 | arr G = ideal(x,y,z,x+y+z); |
---|
2138 | |
---|
2139 | arrIsFree(A3); |
---|
2140 | |
---|
2141 | arrIsFree(B3); |
---|
2142 | |
---|
2143 | arrIsFree(G); |
---|
2144 | } |
---|
2145 | |
---|
2146 | |
---|
2147 | // exponents of a (free) arrangement |
---|
2148 | proc arrExponents |
---|
2149 | "USAGE: arrExponents(A); arr A, multarr A |
---|
2150 | RETURN: [intvec] The exponents of a free (multi-) arrangement, i.e. the degrees of a |
---|
2151 | basis of D(A) the derivation module. |
---|
2152 | NOTE: only defined for central arrangements |
---|
2153 | SEE ALSO: arrDer, arrIsFree, arrExponents |
---|
2154 | KEYWORDS: free; exponents; multiarrangement |
---|
2155 | EXAMPLE: example arrExponents;" |
---|
2156 | |
---|
2157 | { |
---|
2158 | def A = #[1]; |
---|
2159 | |
---|
2160 | if(arrIsFree(A) != 1){ ERROR("Arrangement is not free!"); } |
---|
2161 | |
---|
2162 | module der = arrDer(A); |
---|
2163 | intvec exp; |
---|
2164 | |
---|
2165 | for(int i =1; i<=size(der); i++){ |
---|
2166 | exp[i] = deg(der[i]); |
---|
2167 | } |
---|
2168 | |
---|
2169 | return(exp); |
---|
2170 | } |
---|
2171 | |
---|
2172 | example |
---|
2173 | { |
---|
2174 | "EXAMPLE: computing the exponents of the Edelman-Reiner arrangement and its restriction: "; |
---|
2175 | echo = 2; |
---|
2176 | ring R = 0,(x,y,z),dp; |
---|
2177 | arr A3 = arrBoolean(3); |
---|
2178 | arr B3 = arrTypeB(3); |
---|
2179 | arr G = ideal(x,y,z,x+y+z); |
---|
2180 | |
---|
2181 | arrExponents(A3); |
---|
2182 | |
---|
2183 | arrExponents(B3); |
---|
2184 | |
---|
2185 | } |
---|
2186 | |
---|
2187 | |
---|
2188 | // Tries to reduce the number of generators of a generating set for a module |
---|
2189 | static proc lessGenerators(module X){ |
---|
2190 | module Z = syz(X); |
---|
2191 | matrix K = getColumnIndependentUnitPositions(Z); |
---|
2192 | |
---|
2193 | if( K == unitmat(nrows(Z)) ){ |
---|
2194 | return(X); |
---|
2195 | } |
---|
2196 | |
---|
2197 | module Xnew = X*K; |
---|
2198 | Xnew = simplify(Xnew,2); |
---|
2199 | |
---|
2200 | return(lessGenerators(Xnew)); |
---|
2201 | } |
---|
2202 | |
---|
2203 | |
---|
2204 | // Looks for a generator which is redundant |
---|
2205 | static proc getColumnIndependentUnitPositions(module Z){ |
---|
2206 | int n = nrows(Z); // number of generators of D |
---|
2207 | matrix K = unitmat(n); |
---|
2208 | int i; |
---|
2209 | for(int j=1; j<=ncols(Z); j++){ |
---|
2210 | for(i=1; i<=nrows(Z); i++){ |
---|
2211 | if(deg(Z[i,j]) == 0){ |
---|
2212 | K[i,i] = 0; |
---|
2213 | return(K); |
---|
2214 | } |
---|
2215 | } |
---|
2216 | } |
---|
2217 | |
---|
2218 | return(K); |
---|
2219 | } |
---|
2220 | |
---|
2221 | |
---|
2222 | // Outputs the ideal of powers of hyperplanes of a multiarrangement. |
---|
2223 | // Needed for the computation of arrDer. |
---|
2224 | static proc multIdeal(multarr A){ |
---|
2225 | ideal I; |
---|
2226 | |
---|
2227 | for(int i=1; i<=size(A.l); i++){ |
---|
2228 | I[i] = A.l[i]^A.m[i]; |
---|
2229 | } |
---|
2230 | |
---|
2231 | return(I); |
---|
2232 | } |
---|
2233 | |
---|
2234 | |
---|
2235 | //============================================================================// |
---|
2236 | //-------------------------- #15 MULTI-ARRANGEMENTS --------------------------// |
---|
2237 | //============================================================================// |
---|
2238 | |
---|
2239 | //============================================================================// |
---|
2240 | //------------------------------- CONSTRUCTORS -------------------------------// |
---|
2241 | //============================================================================// |
---|
2242 | |
---|
2243 | |
---|
2244 | // general method for creating multiarrangements |
---|
2245 | static proc multarrAdd |
---|
2246 | "USAGE: A = #; A +#; # list containing arr/ideal/list/matrix/poly |
---|
2247 | RETURN: [multarr] multiarrangement constructed by the input parameters. |
---|
2248 | NOTE: algorithm splits up the list # and uses appropriate procedure |
---|
2249 | to handle the input |
---|
2250 | KEYWORDS: multiarrangement; equal; constructor; operator |
---|
2251 | EXAMPLE: example multarrAdd; shows an example" |
---|
2252 | |
---|
2253 | { |
---|
2254 | multarr A; |
---|
2255 | for(int k=1; k<=size(#); k++){ |
---|
2256 | while(1){ //simulates switch, which singular doesn't offer |
---|
2257 | if(typeof(#[k]) == "poly" ) {A = multarrAddPoly (A, #[k]);break;} |
---|
2258 | if(typeof(#[k]) == "ideal") {A = multarrAddIdeal(A, #[k]);break;} |
---|
2259 | if(typeof(#[k]) == "multarr"){A = multarrAddArr (A, #[k]);break;} |
---|
2260 | if(typeof(#[k]) == "list" ) {A = multarrAdd ( #[k]);break;} |
---|
2261 | ERROR("bad input type"); |
---|
2262 | } |
---|
2263 | } |
---|
2264 | return (A); |
---|
2265 | } |
---|
2266 | |
---|
2267 | example |
---|
2268 | { |
---|
2269 | "EXAMPLE: Creating a few multiarrangements"; echo = 2; |
---|
2270 | ring R = 0,(x,y,z),dp; |
---|
2271 | multarr A = ideal(x,y,z,x); A; |
---|
2272 | multarr B = A + ideal(x+1, x-1,y); B; |
---|
2273 | multarr C = list(B, x+1, x-1); C; |
---|
2274 | multarr D = (x2 - y2)^2; D; |
---|
2275 | } |
---|
2276 | |
---|
2277 | |
---|
2278 | // adds a single poly to the arrangement |
---|
2279 | // if the poly is linear, it is just added. If not Singular tries factorization |
---|
2280 | static proc multarrAddPoly(multarr A, poly p) |
---|
2281 | { |
---|
2282 | p=simplify(p,1); |
---|
2283 | if(deg(p) == 0){ |
---|
2284 | ERROR("Given poly is not linear or Singluar is not able to factorize it");} |
---|
2285 | else{ |
---|
2286 | if(deg(p) == 1){ |
---|
2287 | int k; |
---|
2288 | int b = 0; |
---|
2289 | for(k=1; k<=size(A.l); k++){ |
---|
2290 | if(A.l[k] == p){ |
---|
2291 | A.m[k] = A.m[k]+1; |
---|
2292 | b = 1; |
---|
2293 | } |
---|
2294 | } |
---|
2295 | if(b == 0){ |
---|
2296 | A.l[size(A.l)+1] = p; |
---|
2297 | A.m[size(A.l)] = 1; |
---|
2298 | } |
---|
2299 | return(A); |
---|
2300 | } |
---|
2301 | else{ |
---|
2302 | list I = factorize(p,2); |
---|
2303 | if((size(I[1]) == 1) && (deg(I[1][1] > 1))){ERROR("Given poly is not a hyperplane");} |
---|
2304 | else{ |
---|
2305 | int j,i; |
---|
2306 | ideal J; |
---|
2307 | for(i=1; i<=size(I[1]); i++) |
---|
2308 | { |
---|
2309 | for(j=1; j<=I[2][i]; j++){ |
---|
2310 | J[size(J)+1] = I[1][i]; |
---|
2311 | } |
---|
2312 | } |
---|
2313 | return (multarrAdd(A,J)); |
---|
2314 | } |
---|
2315 | } |
---|
2316 | } |
---|
2317 | return(A); |
---|
2318 | } |
---|
2319 | |
---|
2320 | |
---|
2321 | // adds defining polys to the arrangement |
---|
2322 | static proc multarrAddIdeal(multarr A, ideal I){ |
---|
2323 | for(int k=1; k<=size(I); k++){ |
---|
2324 | A = multarrAddPoly(A,I[k]); |
---|
2325 | } return (A); |
---|
2326 | } |
---|
2327 | |
---|
2328 | |
---|
2329 | // union of two multarrangements |
---|
2330 | static proc multarrAddArr(multarr A, multarr B){ |
---|
2331 | return (multarrAddIdeal(A, multIdeal(B))); |
---|
2332 | } |
---|
2333 | |
---|
2334 | |
---|
2335 | //============================================================================// |
---|
2336 | //------------------------------- TYPE CASTING -------------------------------// |
---|
2337 | //============================================================================// |
---|
2338 | |
---|
2339 | |
---|
2340 | // computes the defining polynomial |
---|
2341 | static proc multarr2poly(multarr A) |
---|
2342 | "USAGE: multArrQPoly(multarr A); |
---|
2343 | RETURN: [poly] q: defining polynomial of a multiarrangement with multiplicities |
---|
2344 | SEE ALSO: arr2poly, multarr2poly |
---|
2345 | KEYWORDS: multiarrangement; defining polynomial |
---|
2346 | EXAMPLE: example multArrQPolys;" |
---|
2347 | |
---|
2348 | { |
---|
2349 | poly q=1; |
---|
2350 | for(int i=1;i<=size(A.l);i=i+1) |
---|
2351 | { |
---|
2352 | q=q*A.l[i]^A.m[i]; |
---|
2353 | } |
---|
2354 | return(q); |
---|
2355 | } |
---|
2356 | |
---|
2357 | example |
---|
2358 | { |
---|
2359 | "EXAMPLE: Computing the Q-Poly for a multiarrangement"; echo = 2; |
---|
2360 | ring R = 0,(x,y,z),dp; |
---|
2361 | multarr A = ideal(x,y,z,x,y,x-y,x-z,x-z,y-z); A; |
---|
2362 | poly q=multarr2poly(A); q; |
---|
2363 | } |
---|
2364 | |
---|
2365 | |
---|
2366 | // converts simple arrangement to multiarrangement |
---|
2367 | proc arr2multarr(arr A, intvec v) |
---|
2368 | "USAGE: multArrFromIntvec(arr A, intvec v); |
---|
2369 | RETURN: [multarr] multiarrangement MA, which is the arrangement A with multiplicities v |
---|
2370 | NOTE: the size of v must match the number of hyperplanes of the arrangement A |
---|
2371 | SEE ALSO: arr2multarr, multarr2arr |
---|
2372 | KEYWORDS: multiarrangement; arrangement; constructor |
---|
2373 | EXAMPLE: example arr2multarr" |
---|
2374 | |
---|
2375 | { |
---|
2376 | if(ncols(A.l)!=size(v)){ |
---|
2377 | ERROR("Vector's size does not match the hyperplanes."); |
---|
2378 | } |
---|
2379 | multarr B; |
---|
2380 | B.l=A.l; |
---|
2381 | B.m=v; |
---|
2382 | return(B); |
---|
2383 | } |
---|
2384 | |
---|
2385 | example |
---|
2386 | { |
---|
2387 | "EXAMPLE:"; echo = 2; |
---|
2388 | ring R = 0,(x,y,z),dp; |
---|
2389 | arr A = arrTypeB(3); A; |
---|
2390 | intvec v=2:9; v; |
---|
2391 | multarr MA=arr2multarr(A,v); |
---|
2392 | MA; |
---|
2393 | } |
---|
2394 | |
---|
2395 | |
---|
2396 | // converts multiarrangement to simple arrangement |
---|
2397 | proc multarr2arr(multarr A) |
---|
2398 | "USAGE: multarr2arr(multarr A, intvec v); |
---|
2399 | RETURN: [arr] arrangement A, with all multiplicities removed |
---|
2400 | SEE ALSO: arr2multarr, multarr2arr |
---|
2401 | KEYWORDS: multiarrangement; arrangement; constructor |
---|
2402 | EXAMPLE: example multarr2arr" |
---|
2403 | |
---|
2404 | { |
---|
2405 | arr B; |
---|
2406 | B.l=A.l; |
---|
2407 | return(B); |
---|
2408 | } |
---|
2409 | |
---|
2410 | example |
---|
2411 | { |
---|
2412 | "EXAMPLE:"; echo = 2; |
---|
2413 | ring r = 0,(x,y,z),dp; |
---|
2414 | multarr A=x2y3z5; A; |
---|
2415 | arr AS = multarr2arr(A); AS; |
---|
2416 | } |
---|
2417 | |
---|
2418 | |
---|
2419 | //============================================================================// |
---|
2420 | //------------------------------- PRINTING -----------------------------------// |
---|
2421 | //============================================================================// |
---|
2422 | |
---|
2423 | |
---|
2424 | // prints arrangement in the console |
---|
2425 | static proc multarrPrint(multarr A) |
---|
2426 | "USAGE: A; A arr |
---|
2427 | RETURN: [] better readable output in the console as newstruct print. |
---|
2428 | SEE ALSO: arrPrint, multarrPrint |
---|
2429 | KEYWORDS: print |
---|
2430 | EXAMPLE: example multarrPrint;" |
---|
2431 | |
---|
2432 | { |
---|
2433 | for(int j=1;j<=ncols(A.l);j++){ |
---|
2434 | print("_["+string(j)+"]=("+string(A.l[j])+")^"+ string(A.m[j])); |
---|
2435 | } |
---|
2436 | } |
---|
2437 | |
---|
2438 | example |
---|
2439 | { |
---|
2440 | "EXAMPLE:"; echo = 2; |
---|
2441 | ring R = 0,(x,y,z),dp; |
---|
2442 | multarr A = ideal(x2,y3,z); |
---|
2443 | A; |
---|
2444 | } |
---|
2445 | |
---|
2446 | |
---|
2447 | // number of hyperplanes with multiplicities |
---|
2448 | static proc multarrSize(multarr A) |
---|
2449 | "USAGE: size(A); A multarr |
---|
2450 | RETURN: [int] Number of hyperplanes with multiplicities |
---|
2451 | SEE ALSO: arrSize |
---|
2452 | KEYWORDS: multiarrangement; size; number; hyperplanes |
---|
2453 | EXAMPLE: example multarrSize;" |
---|
2454 | |
---|
2455 | { |
---|
2456 | return (sum(A.m)); |
---|
2457 | } |
---|
2458 | |
---|
2459 | example |
---|
2460 | { |
---|
2461 | "EXAMPLE:"; echo = 2; |
---|
2462 | ring R = 0,(x,y,z),dp; |
---|
2463 | multarr A = ideal(x2,y3,z); |
---|
2464 | A; |
---|
2465 | } |
---|
2466 | |
---|
2467 | |
---|
2468 | //============================================================================// |
---|
2469 | //------------------------ RESTRICTION & DELETION ----------------------------// |
---|
2470 | //============================================================================// |
---|
2471 | |
---|
2472 | |
---|
2473 | // decrements the multiplicity of a hyperplane by one |
---|
2474 | static proc multarrDelete(multarr A, int k) |
---|
2475 | "USAGE: multarrDelete(A, k); arrangement A, integer k; |
---|
2476 | RETURN: [multarr] the hyperplane Multiarrangement A', i.e. the multiarrangement |
---|
2477 | with multiplicity of H_k decremented by one. If m(H_k)=1, then the hyperplane |
---|
2478 | H_k is deleted |
---|
2479 | SEE ALSO: arrDelete, multarrDelete, |
---|
2480 | KEYWORDS: multiarrangement; delete |
---|
2481 | EXAMPLE: example multarrDelete;" |
---|
2482 | |
---|
2483 | { |
---|
2484 | if(k>ncols(A.l)){ |
---|
2485 | ERROR("There is no k-th hyperplane"); |
---|
2486 | } |
---|
2487 | |
---|
2488 | poly q = multarr2poly(A); |
---|
2489 | q = q/(A.l[k]); |
---|
2490 | multarr MA = q; |
---|
2491 | return (MA); |
---|
2492 | } |
---|
2493 | |
---|
2494 | example |
---|
2495 | { |
---|
2496 | "EXAMPLE:"; echo = 2; |
---|
2497 | ring R = 0,(x,y,z),dp; |
---|
2498 | multarr A =ideal(x2,y3,z); |
---|
2499 | multarr AD = multarrDelete(A,2); AD; |
---|
2500 | } |
---|
2501 | |
---|
2502 | |
---|
2503 | // restriction of A (as arr) to a hyperplane with multiplicities |
---|
2504 | proc multarrRestrict(arr A,intvec v, list #) |
---|
2505 | "USAGE: multarrRestrict(A, v); arrangement A, int/intvec v, optional argument "CC"; |
---|
2506 | RETURN: [multarr] the restricted hyperplane Multi-Arrangement (A^X) with multiplicities |
---|
2507 | i.e. counting how often one element of the restricted arrangement occurs |
---|
2508 | as intersetion of hyperplane of the first arrangement. This definition |
---|
2509 | is due to Guenter M. Ziegler. |
---|
2510 | NOTE: A has to be non-empty. |
---|
2511 | REMARKS: We restrict A to the flat X, defined by the equations in A[v]. |
---|
2512 | The restriction will only be performed, if the ideal defining |
---|
2513 | the flat X is monomial (i.e. X is an intersection of coordinate planes). |
---|
2514 | If the optional argument CC is given, the arrangement is transformed |
---|
2515 | in such a way that X has the above form. |
---|
2516 | SEE ALSO: multarrRestrict, multarrMultRestrict, arrRestrict |
---|
2517 | KEYWORDS: multiarrangement; restriction |
---|
2518 | EXAMPLE: example multarrRestrict; " |
---|
2519 | |
---|
2520 | { |
---|
2521 | ideal I=A.l; |
---|
2522 | I=I[v]; |
---|
2523 | option(redSB); |
---|
2524 | I=std(I); // defining equations for flat X |
---|
2525 | //option(none); |
---|
2526 | ideal Al=arr2ideal(A); |
---|
2527 | int i; |
---|
2528 | |
---|
2529 | if(isMonomial(I)==1){ |
---|
2530 | for(i=1;i<= size(I);i++){ |
---|
2531 | Al=subst(Al,I[i],0); |
---|
2532 | } |
---|
2533 | multarr AR=simplify(Al,3); |
---|
2534 | return(AR); |
---|
2535 | } |
---|
2536 | |
---|
2537 | if(size(#) == 0){ |
---|
2538 | ERROR("The flat X is not defined by a monomial ideal. " + |
---|
2539 | "Please do a coordinate change first. You can use arrCoordNormalize to " + |
---|
2540 | "transform the arrangement accordingly by adding the argument CC."); |
---|
2541 | } |
---|
2542 | if(#[1]!="CC"){ |
---|
2543 | ERROR("The flat X is not defined by a monomial ideal. " + |
---|
2544 | "Please do a coordinate change first. You can use arrCoordNormalize to " + |
---|
2545 | "transform the arrangement accordingly by adding the argument CC."); |
---|
2546 | } |
---|
2547 | if(size(v)==0){return(A);} |
---|
2548 | intvec w=v[1]; |
---|
2549 | intvec tmp; |
---|
2550 | |
---|
2551 | for(i=2;i<=size(v);i++){ |
---|
2552 | tmp=w,v[i]; |
---|
2553 | if(rank(matrix(A[w]))==size(tmp)){ |
---|
2554 | w=tmp; |
---|
2555 | } |
---|
2556 | } |
---|
2557 | |
---|
2558 | arr B=arrCoordNormalize(A,w); |
---|
2559 | ideal Bl=arr2ideal(B); |
---|
2560 | |
---|
2561 | for(i=1;i<= size(w);i++){ |
---|
2562 | Bl=subst(Bl,Bl[w[i]],0); |
---|
2563 | } |
---|
2564 | multarr BR=simplify(Bl,3); |
---|
2565 | return(BR); |
---|
2566 | } |
---|
2567 | |
---|
2568 | example |
---|
2569 | { |
---|
2570 | "EXAMPLE:"; echo = 2; |
---|
2571 | ring R = 0,x(1..5),dp; |
---|
2572 | arr A = arrEdelmanReiner(); A; |
---|
2573 | multarr AR = multarrRestrict(A,6,"CC"); AR; |
---|
2574 | } |
---|
2575 | |
---|
2576 | |
---|
2577 | // restriction of A (as multarr) to a hyperplane with multiplicities |
---|
2578 | proc multarrMultRestrict(multarr A,int k) |
---|
2579 | "USAGE: multarrMultRestrict(A, k); multiarrangement A, integer k; |
---|
2580 | RETURN: [multarr] the restricted hyperplane Multi-Arrangement (A^H_k) with multiplicities, i.e. |
---|
2581 | counting with multiplicities how often one element of the restricted arrangement occurs |
---|
2582 | as intersetion of hyperplane of the first multiarrangement. |
---|
2583 | This definition is due to Guenter M. Ziegler. |
---|
2584 | NOTE: A has to be non-empty. |
---|
2585 | REMARKS: The restriction will only be performed, if H_k = ker(x_i) for some i. |
---|
2586 | One can also restrict an arrangement with respect to any hyper- |
---|
2587 | plane k, but than a coordinate change is necessary first to make |
---|
2588 | H_k = ker(x_k). Since such a coordinate change is not unique, please |
---|
2589 | use arrCoordchange to do so. |
---|
2590 | SEE ALSO: multarrRestrict, multarrMultRestrict, arrRestrict |
---|
2591 | KEYWORDS: multiarrangement; restriction |
---|
2592 | EXAMPLE: example multarrMultRestrict; " |
---|
2593 | |
---|
2594 | { |
---|
2595 | ideal I = variables(A.l); |
---|
2596 | ideal J; |
---|
2597 | int j,i; |
---|
2598 | for(i=1;i<=size(A.l);i=i+1) |
---|
2599 | { |
---|
2600 | for(j=1;j<=A.m[i];j++){ |
---|
2601 | J[ncols(J)+1]=A.l[i]; |
---|
2602 | } |
---|
2603 | } |
---|
2604 | poly p = simplify(A.l[k],1); |
---|
2605 | int n = rvar(p); |
---|
2606 | if( n == 0 ){ |
---|
2607 | ERROR("H_" + string(k) + " = " + string(p) + |
---|
2608 | " is not of the form ker(x_i). Please do a coordinate change first." + |
---|
2609 | "You can use arrCoordinateChange to transform the arrangement accordingly."); |
---|
2610 | } |
---|
2611 | J = subst(J, p, 0); |
---|
2612 | multarr MA = simplify(J,3); |
---|
2613 | return(MA); |
---|
2614 | } |
---|
2615 | |
---|
2616 | example |
---|
2617 | { |
---|
2618 | "EXAMPLE:"; echo = 2; |
---|
2619 | ring R = 0,(x,y,z),dp; |
---|
2620 | multarr A =ideal(x2,y2,z2,(x-y)^3,(x-z)^2,(y-z)); A; |
---|
2621 | //The restriction of the multiarrangement is: |
---|
2622 | multarr AR = multarrMultRestrict(A,1); AR; |
---|
2623 | } |
---|
2624 | |
---|
2625 | |
---|
2626 | //============================================================================// |
---|
2627 | //----------------------------- #16 COMBINATORICS ---------------------------// |
---|
2628 | //============================================================================// |
---|
2629 | |
---|
2630 | /* |
---|
2631 | RELATIONS can have 3 values: |
---|
2632 | -1 : flat&hyperplane are parallel |
---|
2633 | 0 : hyperplane intersects the flat, but only further on in the lattice. |
---|
2634 | +1 : hyperplane is part of the flat. |
---|
2635 | |
---|
2636 | */ |
---|
2637 | // intvec: multiplicities of hyperplanes |
---|
2638 | |
---|
2639 | |
---|
2640 | proc arrLattice(arr ARR) |
---|
2641 | "USAGE: arrLattice(arr ARR) |
---|
2642 | RETURN: [arrposet] intersection poset of the arrangement |
---|
2643 | NOTE: The algorithm works by a bottom up approach, i.e. it calculates the |
---|
2644 | SEE ALSO: arrLattice, arrFlats |
---|
2645 | KEYWORDS: intersection lattice; poset; lattice |
---|
2646 | EXAMPLE: example arrFlats;" |
---|
2647 | |
---|
2648 | { |
---|
2649 | |
---|
2650 | //Performance test |
---|
2651 | system("--ticks-per-sec",1); |
---|
2652 | timer = 0; |
---|
2653 | int start = timer; |
---|
2654 | // start |
---|
2655 | dbprint(3-voice,newline); |
---|
2656 | dbprint(3-voice,"=== Computing poset ==="); |
---|
2657 | dbprint(3-voice,newline); |
---|
2658 | |
---|
2659 | // Initialization |
---|
2660 | list L,L2; |
---|
2661 | intvec u,v, NFLAT; |
---|
2662 | int i,j,l,m,n,rk, rk2,rk3, tic, toc, numberOfFlats,numberOfHplanes, counter, ntests; |
---|
2663 | matrix NewRel, M; |
---|
2664 | |
---|
2665 | // Initialization of matrix |
---|
2666 | ntests = 0; |
---|
2667 | M = matrix(ARR); // m x n+1 matrix |
---|
2668 | m = size(ARR); // m = # hplanes |
---|
2669 | n = nvars(basering); // n = # variables |
---|
2670 | numberOfHplanes = m; |
---|
2671 | |
---|
2672 | // Initialization of poset and flat |
---|
2673 | flat F; |
---|
2674 | arrposet P; P.A = ARR; |
---|
2675 | F.REL = 0:m; F.moebius = -1; |
---|
2676 | for(i=1; i<=n; i++){ P.r = insert(P.r,list());} |
---|
2677 | for(i=1; i<=m; i++){ F.REL[i] = 1; L[i] = F; F.REL[i] = 0;} |
---|
2678 | F.moebius = 0; |
---|
2679 | P.r[1] = L; |
---|
2680 | |
---|
2681 | // calculating r[i] step by step |
---|
2682 | for(rk=2; rk<=rank(ARR); rk++){ // maximal flat possible has rank A.v ofc |
---|
2683 | counter = 0; |
---|
2684 | tic = timer; |
---|
2685 | |
---|
2686 | // Initialization of current round |
---|
2687 | L = P.r[rk-1]; // flats from the last iteration give rise to current ones. |
---|
2688 | numberOfFlats = size(L); |
---|
2689 | L2 = list(); // new list for elts of rk= |
---|
2690 | NewRel = getOldRel(L,m); // contains all the information of the old arrangement |
---|
2691 | |
---|
2692 | // find the next leftmost child by looking for intersetion of F_j with hyperplanes |
---|
2693 | for(j=1; j<=numberOfFlats; j++){ // goes over the elts of L.r[i-1] |
---|
2694 | |
---|
2695 | NFLAT = L[j].REL; |
---|
2696 | |
---|
2697 | for(i=1; i<=numberOfHplanes; i++){ |
---|
2698 | if (NewRel[i,j] == 0){ |
---|
2699 | //test this hyperplane |
---|
2700 | NFLAT[i] = 1; |
---|
2701 | v = collectHplanes(NFLAT); |
---|
2702 | rk2 = rank(submat(M,v,1..n)); |
---|
2703 | ntests = ntests+2; |
---|
2704 | |
---|
2705 | // CHILD FOUND BY INTERSECTING WITH H_i |
---|
2706 | if(rk2 == rank(submat(M,v,1..n+1))){ |
---|
2707 | NewRel[i,j] = 1; counter++; |
---|
2708 | |
---|
2709 | // potentially there exist more hyperplanes which intersect in the same space |
---|
2710 | for(l=i+1; l<=m; l++){ |
---|
2711 | if(NewRel[l,j] == 0){ |
---|
2712 | NFLAT[l] = 1; |
---|
2713 | u = collectHplanes(NFLAT); |
---|
2714 | rk3 = rank(submat(M, u, 1..n)); |
---|
2715 | ntests = ntests +2; |
---|
2716 | |
---|
2717 | // FLAT EXISTS |
---|
2718 | if(rk3 == rank(submat(M,u,1..n+1)) ){ |
---|
2719 | if(rk2 == rk3){NewRel[l,j] = 1;} |
---|
2720 | else{NFLAT[l] = 0;} //itersecting with both yields higher rank flat |
---|
2721 | } |
---|
2722 | |
---|
2723 | // FLAT DOES NOT EXIST |
---|
2724 | else{NFLAT[l] = -1;} |
---|
2725 | } |
---|
2726 | } |
---|
2727 | |
---|
2728 | // Add new flat to list, reset NFLAT. |
---|
2729 | F.REL = NFLAT; |
---|
2730 | NFLAT = L[j].REL; // refresh. |
---|
2731 | |
---|
2732 | // Find Parents: (most expensive part it seems) |
---|
2733 | F.parents = j:1; //parent in any case |
---|
2734 | for(l=j+1; l<=numberOfFlats; l++){ |
---|
2735 | if(isParent(L[l].REL, F.REL)){ |
---|
2736 | F.parents = F.parents,l; |
---|
2737 | NewRel[1..m,l] = updateRel(F.REL, submat(NewRel,1..m,l)); |
---|
2738 | } |
---|
2739 | } |
---|
2740 | |
---|
2741 | L2 = insert(L2,F,size(L2)); |
---|
2742 | |
---|
2743 | // means that the flat is parallel |
---|
2744 | } else{ |
---|
2745 | L[j].REL[i] = -1; |
---|
2746 | NewRel[i,j] = -1; |
---|
2747 | NFLAT[i] = -1; |
---|
2748 | } |
---|
2749 | |
---|
2750 | } |
---|
2751 | } |
---|
2752 | } |
---|
2753 | |
---|
2754 | toc = timer; |
---|
2755 | dbprint(3-voice,"rank "+string(rk)+": found "+ string(counter) + |
---|
2756 | " flats in "+string(toc - tic)+"s"); |
---|
2757 | counter = 0; |
---|
2758 | P.r[rk] = L2; |
---|
2759 | } |
---|
2760 | |
---|
2761 | dbprint(3-voice,newline); |
---|
2762 | //dbprint(3-voice,"Computation time: "+string(timer - start)+"s"); |
---|
2763 | dbprint(3-voice,"Matrix tests: "+string(ntests)); |
---|
2764 | return (P); |
---|
2765 | } |
---|
2766 | |
---|
2767 | example |
---|
2768 | { |
---|
2769 | "EXAMPLE: Intersection lattice of the braid arrangement in 3 dimensions "; echo = 2; |
---|
2770 | ring r; |
---|
2771 | arrLattice(arrTypeB(3)); |
---|
2772 | } |
---|
2773 | |
---|
2774 | |
---|
2775 | static proc arrRank(arr A){ |
---|
2776 | int n = rank(matrix(A)); |
---|
2777 | int m = nvars(A); |
---|
2778 | if(n < m) {return (n);} |
---|
2779 | return (m); |
---|
2780 | } |
---|
2781 | |
---|
2782 | static proc getOldRel(list L, int m){ |
---|
2783 | int n = size(L); |
---|
2784 | matrix NewRel[m][n]; |
---|
2785 | |
---|
2786 | for(int i=1; i<=n; i++){ |
---|
2787 | NewRel[1..m,i] = L[i].REL; |
---|
2788 | } |
---|
2789 | |
---|
2790 | return (NewRel); |
---|
2791 | } |
---|
2792 | |
---|
2793 | static proc updateRel(intvec REL, matrix NewRel){ |
---|
2794 | for(int i=1; i<=size(REL); i++){ |
---|
2795 | if(NewRel[i,1] == 0 && REL[i] == 1){NewRel[i,1] = 1;} |
---|
2796 | } |
---|
2797 | |
---|
2798 | return (NewRel); |
---|
2799 | } |
---|
2800 | |
---|
2801 | // Flat F is a parent of Flat G if F[i] == 1 => G[i] == 1 for all i |
---|
2802 | static proc isParent(intvec parent, intvec child){ |
---|
2803 | int counter = 0; |
---|
2804 | for(int i=1; i<= size(parent); i++){ |
---|
2805 | if(parent[i]==1){ |
---|
2806 | if(child[i] == 1){counter++;} |
---|
2807 | else{ return (0); } |
---|
2808 | } |
---|
2809 | } |
---|
2810 | if(counter == 0){ return (0); } |
---|
2811 | return (1); |
---|
2812 | } |
---|
2813 | |
---|
2814 | // transforms the "rel" intvec into an intvec which contains the indices of the hyperplanes. |
---|
2815 | static proc collectHplanes(intvec RELATIONS){ |
---|
2816 | intvec result; |
---|
2817 | |
---|
2818 | for(int i=1; i<=size(RELATIONS); i++){ |
---|
2819 | if(RELATIONS[i] > 0){result = result,i;} |
---|
2820 | } |
---|
2821 | |
---|
2822 | result = result[2..size(result)]; |
---|
2823 | |
---|
2824 | return (result); |
---|
2825 | |
---|
2826 | } |
---|
2827 | |
---|
2828 | static proc getFlat(arrposet P, int i, int j){ |
---|
2829 | list L = P.r[i]; |
---|
2830 | return (L[j]); |
---|
2831 | } |
---|
2832 | |
---|
2833 | static proc setFlat(arrposet P, int i, int j, flat F){ |
---|
2834 | list L = P.r[i]; |
---|
2835 | L[j] = F; |
---|
2836 | P.r[i] = L; |
---|
2837 | return (P); |
---|
2838 | } |
---|
2839 | |
---|
2840 | static proc getFlag(arrposet P, int i, int j){ |
---|
2841 | flat F = getFlat(P,i,j); |
---|
2842 | return (F.flag); |
---|
2843 | } |
---|
2844 | |
---|
2845 | static proc setFlag(arrposet P, int i, int j, int flag){ |
---|
2846 | flat F = getFlat(P,i,j); |
---|
2847 | F.flag = flag; |
---|
2848 | return (setFlat(P,i,j,F)); |
---|
2849 | } |
---|
2850 | |
---|
2851 | //============================================================================// |
---|
2852 | //------------------------- #16 INTERSECTION-LATTICE -------------------------// |
---|
2853 | //============================================================================// |
---|
2854 | //============================================================================// |
---|
2855 | //-------------------------- MOEBIUS -------------------------------------// |
---|
2856 | //============================================================================// |
---|
2857 | |
---|
2858 | |
---|
2859 | // calculates the moebius values of the poset |
---|
2860 | proc moebius(arrposet P) |
---|
2861 | "USAGE: moebius(arrposet P) |
---|
2862 | RETURN: [arrposet] fills in the moebius values of the flats in the poset |
---|
2863 | SEE ALSO: moebius |
---|
2864 | KEYWORDS: moebius function |
---|
2865 | EXAMPLE: example arrCharPoly; shows an example" |
---|
2866 | |
---|
2867 | { |
---|
2868 | list L; |
---|
2869 | int i,j; |
---|
2870 | |
---|
2871 | for(i=2; i<=rank(P.A); i++){ |
---|
2872 | L = P.r[i]; |
---|
2873 | |
---|
2874 | for(j=1; j<=size(L); j++){ |
---|
2875 | arrposet Q = P; |
---|
2876 | export Q; |
---|
2877 | L[j].moebius = -(1 + moebiusRecursion(i, j)); |
---|
2878 | kill Q; |
---|
2879 | } |
---|
2880 | P.r[i] = L; |
---|
2881 | |
---|
2882 | } |
---|
2883 | |
---|
2884 | return (P); |
---|
2885 | } |
---|
2886 | |
---|
2887 | example |
---|
2888 | { |
---|
2889 | "EXAMPLE: We look at the moebius values of the braid arrangement in 4 dimensions: "; echo = 2; |
---|
2890 | ring R = 0,(x,y,z,t),dp; |
---|
2891 | arr A = arrBraid(4); |
---|
2892 | arrposet P = arrLattice(A); |
---|
2893 | P; |
---|
2894 | //As you can see the values are not calculated yet: |
---|
2895 | |
---|
2896 | printMoebius(P); |
---|
2897 | P = moebius(P); |
---|
2898 | |
---|
2899 | //Now all entries are initialized: |
---|
2900 | |
---|
2901 | printMoebius(P); |
---|
2902 | } |
---|
2903 | |
---|
2904 | |
---|
2905 | // moebiusrecursion(i,j) is the sum moebius(X) over {X | V>X>Flat_ij} |
---|
2906 | // Vss: all moebius values of rank < i are known and correctly saved in P |
---|
2907 | static proc moebiusRecursion(int i, int j){ |
---|
2908 | flat F = getFlat(Q, i, j); |
---|
2909 | int m = F.moebius; |
---|
2910 | |
---|
2911 | if(getFlag(Q,i,j) == 1){return (0);} |
---|
2912 | else { |
---|
2913 | Q = setFlag(Q,i,j,1); |
---|
2914 | if(i > 1){ |
---|
2915 | for(int k=1; k<=size(F.parents); k++){ |
---|
2916 | m = m + moebiusRecursion(i-1, F.parents[k]); |
---|
2917 | } |
---|
2918 | } |
---|
2919 | } |
---|
2920 | |
---|
2921 | return (m); |
---|
2922 | } |
---|
2923 | |
---|
2924 | |
---|
2925 | // pi(A,t) = sum[X in L(A); moebius(X)*(-t)^r(X)] |
---|
2926 | proc arrPoincare(def input) |
---|
2927 | "USAGE: arrPoincare(A); arr A |
---|
2928 | RETURN: [intvec] The Poincare polynomial as integer vector of the arrangement, which |
---|
2929 | is equal to the second kind Poincare-Series of the Orlik-Solomon Algebra. |
---|
2930 | SEE ALSO: arrPoincare, arrChambers, arrBoundedChambers |
---|
2931 | KEYWORDS: Poincare polynomial |
---|
2932 | EXAMPLE: example arrPoincare;" |
---|
2933 | |
---|
2934 | { |
---|
2935 | while(1){ |
---|
2936 | if(typeof(input) == "arr"){ |
---|
2937 | arr A = input; |
---|
2938 | arrposet P = arrLattice(A); |
---|
2939 | break; |
---|
2940 | } |
---|
2941 | if(typeof(input) == "arrposet"){ |
---|
2942 | arr A = input.A; |
---|
2943 | arrposet P = input; |
---|
2944 | break; |
---|
2945 | } |
---|
2946 | ERROR("Bad input type"); |
---|
2947 | } |
---|
2948 | |
---|
2949 | P = moebius(P); |
---|
2950 | int i, j, coeff, sign; |
---|
2951 | list L; |
---|
2952 | |
---|
2953 | intvec v = 1; |
---|
2954 | sign = 1; |
---|
2955 | for(i=1; i<= rank(A); i++){ |
---|
2956 | L = P.r[i]; |
---|
2957 | sign = (-1)*sign; |
---|
2958 | coeff = 0; |
---|
2959 | for(j=1;j<=size(L);j++){ |
---|
2960 | coeff = coeff + L[j].moebius; |
---|
2961 | } |
---|
2962 | coeff = sign*coeff; |
---|
2963 | v = v,coeff; |
---|
2964 | } |
---|
2965 | return (v); |
---|
2966 | } |
---|
2967 | |
---|
2968 | example |
---|
2969 | { |
---|
2970 | "EXAMPLE: The Poincare polynomial of the braid arrangement in k dimensions is given as: " + |
---|
2971 | "pi(A,t) = (1 + t)*...*(1 + (k-1)*t)"; echo = 2; |
---|
2972 | ring R = 0,(x,y,z,u,v),dp; |
---|
2973 | arr A = arrBraid(5); |
---|
2974 | intvec v = arrPoincare(A); |
---|
2975 | (1+x)*(1+2x)*(1+3x)*(1+4x); |
---|
2976 | v; |
---|
2977 | } |
---|
2978 | |
---|
2979 | |
---|
2980 | // X(A,t) = t^l * pi(A,-t^-1) |
---|
2981 | proc arrCharPoly(def input) |
---|
2982 | "USAGE: arrCharPoly(arr A) |
---|
2983 | RETURN: [intvec] coefficients of the characteristic polynomial of A in increasing order |
---|
2984 | REMARKS: The algorithm only returns the coefficients of the characteristic polynomial since they |
---|
2985 | are whole numbers but the basering could be something different. |
---|
2986 | SEE ALSO: arrCharPoly, arrPoincare |
---|
2987 | KEYWORDS: characteristic polynomial |
---|
2988 | EXAMPLE: example arrCharPoly; shows an example" |
---|
2989 | |
---|
2990 | { |
---|
2991 | intvec v = arrPoincare(input); |
---|
2992 | int l = size(v); |
---|
2993 | v = v[l..1]; //reverses order |
---|
2994 | for(int i=2; i<=l; i=i+2){ |
---|
2995 | v[i] = -v[i]; |
---|
2996 | } |
---|
2997 | |
---|
2998 | return (v); |
---|
2999 | } |
---|
3000 | |
---|
3001 | example |
---|
3002 | { |
---|
3003 | "EXAMPLE: The characteristic polynomial of the Braid arrangement in k dimensions is given as: " + |
---|
3004 | "X(A,t) = t*(t-1)*...*(t-(k-1)))"; echo = 2; |
---|
3005 | ring R = 0,(x,y,z,u,v),dp; |
---|
3006 | arr A = arrBraid(5); |
---|
3007 | intvec v = arrCharPoly(A); |
---|
3008 | x*(x-1)*(x-2)*(x-3)*(x-4); |
---|
3009 | v; |
---|
3010 | } |
---|
3011 | |
---|
3012 | |
---|
3013 | // number of chambers of the arrangement |
---|
3014 | proc arrChambers(arr A) |
---|
3015 | "USAGE: arrChambers(A); arr A |
---|
3016 | RETURN: [int] The number of chambers of an arrangement, which is equal to the |
---|
3017 | evaluation of the Poincare polynomial at 1. |
---|
3018 | SEE ALSO: arrPoincare, arrChambers, arrBoundedChambers |
---|
3019 | KEYWORDS: chambers |
---|
3020 | EXAMPLE: example arrChambers;" |
---|
3021 | |
---|
3022 | { |
---|
3023 | intvec v = arrPoincare(A); |
---|
3024 | return(sum(v)); |
---|
3025 | } |
---|
3026 | |
---|
3027 | example |
---|
3028 | { |
---|
3029 | "EXAMPLE: Computing the number of number of chambers in a simple arrangement"; echo = 2; |
---|
3030 | ring R = 0,(x,y),dp; |
---|
3031 | arr A = ideal(x,y,x+y-1); |
---|
3032 | arrChambers(A); |
---|
3033 | } |
---|
3034 | |
---|
3035 | |
---|
3036 | // number of bounded chambers of the arrangement |
---|
3037 | proc arrBoundedChambers(arr A) |
---|
3038 | "USAGE: arrBoundedChambers(A); arr A |
---|
3039 | RETURN: [int] The number of bounded chambers of an arrangement, which is equal to |
---|
3040 | the evaluation of the Poincare polynomial at -1. |
---|
3041 | SEE ALSO: arrPoincare, arrChambers, arrBoundedChambers |
---|
3042 | KEYWORDS: chambers |
---|
3043 | EXAMPLE: example arrBoundedChambers;" |
---|
3044 | |
---|
3045 | { |
---|
3046 | intvec v = arrPoincare(A); |
---|
3047 | for(int i=2; i<=size(v); i=i+2){ |
---|
3048 | v[i] = -v[i]; |
---|
3049 | } |
---|
3050 | return (sum(v)); |
---|
3051 | } |
---|
3052 | |
---|
3053 | example |
---|
3054 | { |
---|
3055 | "EXAMPLE: Computing the number of bounded chambers in a simple arrangement"; echo = 2; |
---|
3056 | ring R = 0,(x,y),dp; |
---|
3057 | arr A = ideal(x,y,x+y-1); |
---|
3058 | arrBoundedChambers(A); |
---|
3059 | } |
---|
3060 | |
---|
3061 | //============================================================================// |
---|
3062 | //------------------------------- OUTPUT ------------------------------------// |
---|
3063 | //============================================================================// |
---|
3064 | |
---|
3065 | static proc arrPrintPoset(arrposet P){ |
---|
3066 | print("Given Arrangement:"); |
---|
3067 | P.A; |
---|
3068 | print("Corresponding poset:"); |
---|
3069 | |
---|
3070 | string s; |
---|
3071 | int i,j; |
---|
3072 | list L; |
---|
3073 | for(i=1; i<= size(P.r); i++){ |
---|
3074 | print("====== rank "+string(i)+": "+string(size(P.r[i]))+" flats ======"); |
---|
3075 | s = ""; |
---|
3076 | L = P.r[i]; |
---|
3077 | for(j=1; j<=size(L); j++){ |
---|
3078 | s = s + " (" + string(collectHplanes(L[j].REL)) + "), "; |
---|
3079 | } |
---|
3080 | |
---|
3081 | print(s); |
---|
3082 | if(size(P.r[i]) == 0){break;} |
---|
3083 | } |
---|
3084 | } |
---|
3085 | |
---|
3086 | |
---|
3087 | proc printMoebius(arrposet P) |
---|
3088 | "USAGE: printMoebius(A); arr A |
---|
3089 | RETURN: [] displays the moebius values of all the flats in the poset |
---|
3090 | REMARKS: Mainly used for debugging. |
---|
3091 | EXAMPLE: example printMoebius;" |
---|
3092 | |
---|
3093 | { |
---|
3094 | print("Moebius values: "); |
---|
3095 | |
---|
3096 | string s; |
---|
3097 | int i,j; |
---|
3098 | list L; |
---|
3099 | for(i=1; i<= size(P.r); i++){ |
---|
3100 | print("====== rank "+string(i)+": "+string(size(P.r[i]))+" flats ======"); |
---|
3101 | s = ""; |
---|
3102 | L = P.r[i]; |
---|
3103 | for(j=1; j<=size(L); j++){ |
---|
3104 | s = s + " (" + string(L[j].moebius) + "), "; |
---|
3105 | } |
---|
3106 | |
---|
3107 | print(s); |
---|
3108 | if(size(P.r[i]) == 0){break;} |
---|
3109 | } |
---|
3110 | } |
---|
3111 | |
---|
3112 | example |
---|
3113 | { |
---|
3114 | "EXAMPLE: We look at the moebius values of the braid arrangement in 4 dimensions: "; echo = 2; |
---|
3115 | ring R = 0,(x,y,z,t),dp; |
---|
3116 | arr A = arrBraid(4); |
---|
3117 | arrposet P = arrLattice(A); |
---|
3118 | P; |
---|
3119 | //As you can see the values are not calculated yet: |
---|
3120 | |
---|
3121 | printMoebius(P); |
---|
3122 | P = moebius(P); |
---|
3123 | |
---|
3124 | //Now all entries are initialized: |
---|
3125 | |
---|
3126 | printMoebius(P); |
---|
3127 | } |
---|
3128 | |
---|
3129 | //============================================================================// |
---|
3130 | //-------------------------------- arrFlats Stuff --------------------------------// |
---|
3131 | //============================================================================// |
---|
3132 | // test via hyperplanes. |
---|
3133 | proc arrFlats(arr ARR) |
---|
3134 | "USAGE: size(A); A arr |
---|
3135 | RETURN: [arrposet] Intersection lattice |
---|
3136 | SEE ALSO: arrFlats |
---|
3137 | KEYWORDS: intersection lattice |
---|
3138 | EXAMPLE: example arrFlats;" |
---|
3139 | |
---|
3140 | { |
---|
3141 | print(newline); |
---|
3142 | print("=== Computing poset ==="); |
---|
3143 | print(newline); |
---|
3144 | |
---|
3145 | //Performance test |
---|
3146 | system("--ticks-per-sec",1000); |
---|
3147 | timer = 0; |
---|
3148 | int time = timer; |
---|
3149 | |
---|
3150 | // Initialization |
---|
3151 | list L,L2; |
---|
3152 | intvec u,v,w,src; |
---|
3153 | int i,j,k,l,m,n,d, rk,rk2; |
---|
3154 | int counter, ntests; ntests = 0; |
---|
3155 | matrix S,T; |
---|
3156 | intmat REL; |
---|
3157 | |
---|
3158 | // Initialization of matrix |
---|
3159 | matrix M = matrix(ARR); // m x n+1 matrix |
---|
3160 | m = size(ARR); |
---|
3161 | n = nvars(basering); |
---|
3162 | |
---|
3163 | // Initialization of P |
---|
3164 | arrflats P; P.A = ARR; |
---|
3165 | for(i=1; i<=n; i++){ P.r = insert(P.r,list());} |
---|
3166 | for(i=1; i<=m; i++){ L[i] = i;} |
---|
3167 | P.r[1] = L; |
---|
3168 | |
---|
3169 | // calculating r[i] step by step |
---|
3170 | for(i=2; i<=n; i++){ // maximal flat possible has rank A.v ofc |
---|
3171 | |
---|
3172 | counter = 0; |
---|
3173 | |
---|
3174 | L = P.r[i-1]; // flats from the last iteration give rise to current ones. |
---|
3175 | if(size(L) <= 1){break;} // finish early if no more intersections possible. |
---|
3176 | L2 = list(); // new list for elts of rk=i |
---|
3177 | |
---|
3178 | |
---|
3179 | // REL is an intmat that contains information about the relations between |
---|
3180 | // the i-1 flats and the Hyperplanes. |
---|
3181 | // REL[i,j] = 1 means that H_i intersects F_j |
---|
3182 | // REL[i,j] = 0 means that it has not been tested yet |
---|
3183 | // REL[i,j] = -1 means that H_i does NOT intersect F_j, i.e. they are parallel |
---|
3184 | // resetting the REL-matrix |
---|
3185 | REL = intmat(intvec(0),m,size(L)); |
---|
3186 | // fills in the trivial entries, i.e. each flat intersects all hplanes it is composed of. |
---|
3187 | for(j=1; j<=size(L); j++){ |
---|
3188 | u = L[j]; |
---|
3189 | for(k=1; k<= size(u); k++){ |
---|
3190 | REL[u[k],j] = 1; |
---|
3191 | } |
---|
3192 | } |
---|
3193 | |
---|
3194 | //find new flat of rank i |
---|
3195 | for(j=1; j<=size(L); j++){ // goes over the elts of L.r[i-1] |
---|
3196 | |
---|
3197 | // Looking for intersetion of F_j with hyperplanes |
---|
3198 | u = L[j]; |
---|
3199 | |
---|
3200 | for(k=1; k<=m; k++){ |
---|
3201 | if (REL[k,j] == 0 ){ |
---|
3202 | v = insertVal(u,k); |
---|
3203 | counter++; |
---|
3204 | ntests = ntests +2; |
---|
3205 | rk = rank(submat(M,v,1..n)); |
---|
3206 | if(rk == rank(submat(M,v,1..n+1))){ |
---|
3207 | //INTERSECTION FOUND |
---|
3208 | // potentially there exist more hyperplanes which intersect in the same |
---|
3209 | // space. For example (x,y,x+y) all intersect in (0,0) |
---|
3210 | |
---|
3211 | REL[k,j] = 1; |
---|
3212 | |
---|
3213 | for(l=k+1; l<=m; l++){ |
---|
3214 | if(REL[l,j] == 0){ |
---|
3215 | u = insertVal(v,l); |
---|
3216 | ntests = ntests +2; |
---|
3217 | rk2 = rank(submat(M, u, 1..n)); |
---|
3218 | |
---|
3219 | if( rk == rk2 && rk == rank(submat(M,u,1..n+1)) ){ |
---|
3220 | v = u; |
---|
3221 | REL[l,j] = 1; |
---|
3222 | } |
---|
3223 | }} |
---|
3224 | |
---|
3225 | // Add new flat to list, reset u. |
---|
3226 | L2 = insert(L2,v,size(L2)); |
---|
3227 | u = L[j]; //if k = m this is not necessary.... |
---|
3228 | } |
---|
3229 | }} |
---|
3230 | } |
---|
3231 | |
---|
3232 | print("rank: "+string(i)+", expected OPS: "+string(binomial(size(L),2)) |
---|
3233 | +", executed OPS: " + string(counter)); |
---|
3234 | counter = 0; |
---|
3235 | |
---|
3236 | //cleanup => doing this during computation increases speed!! |
---|
3237 | for(j=1; j<size(L2); j++){ |
---|
3238 | u = L2[j]; |
---|
3239 | for(k=j+1; k<= size(L2); k++){ |
---|
3240 | if(u == L2[k]){ |
---|
3241 | L2 = delete(L2,k); k--; counter++; |
---|
3242 | |
---|
3243 | } |
---|
3244 | } |
---|
3245 | } |
---|
3246 | |
---|
3247 | P.r[i] = L2; |
---|
3248 | print("Cleaned up: "+string(counter)+" hyperplanes"); |
---|
3249 | print(newline); |
---|
3250 | } |
---|
3251 | |
---|
3252 | print(newline); |
---|
3253 | //print("Computation time: "+string(timer - time)); |
---|
3254 | print("Matrix tests: "+string(ntests)); |
---|
3255 | |
---|
3256 | return (P); |
---|
3257 | } |
---|
3258 | |
---|
3259 | example |
---|
3260 | { |
---|
3261 | "EXAMPLE: "; echo = 2; |
---|
3262 | ring R = 0,(x(1..5)),dp; |
---|
3263 | arrFlats(arrBraid(5)); |
---|
3264 | } |
---|
3265 | |
---|
3266 | // inserst k in u in the correct lexicographical place |
---|
3267 | static proc insertVal(intvec u, int k){ |
---|
3268 | |
---|
3269 | if(k<u[1]){ |
---|
3270 | u = k,u; |
---|
3271 | return (u); |
---|
3272 | } |
---|
3273 | if(k>u[size(u)]){ |
---|
3274 | u = u,k; |
---|
3275 | return (u); |
---|
3276 | } |
---|
3277 | |
---|
3278 | int i = 2; |
---|
3279 | while(u[i]<k){i++;} |
---|
3280 | |
---|
3281 | u = u[1..(i-1)],k,u[i..size(u)]; |
---|
3282 | |
---|
3283 | |
---|
3284 | return(u); |
---|
3285 | |
---|
3286 | } |
---|
3287 | |
---|
3288 | static proc arrPrintFlats(arrflats P){ |
---|
3289 | print("Given Arrangement:"); |
---|
3290 | P.A; |
---|
3291 | print("Corresponding poset:"); |
---|
3292 | |
---|
3293 | string s; |
---|
3294 | int i,j; |
---|
3295 | list L; |
---|
3296 | for(i=1; i<= size(P.r); i++){ |
---|
3297 | print("====== rank "+string(i)+": "+string(size(P.r[i]))+" flats ======"); |
---|
3298 | s = ""; |
---|
3299 | L = P.r[i]; |
---|
3300 | for(j=1; j<=size(L); j++){ |
---|
3301 | s = s + " (" + string(L[j]) + "), "; |
---|
3302 | } |
---|
3303 | |
---|
3304 | print(s); |
---|
3305 | } |
---|
3306 | } |
---|
3307 | |
---|
3308 | |
---|
3309 | // Compares two posets |
---|
3310 | static proc cposet(arrflats P, arrflats Q){ |
---|
3311 | |
---|
3312 | list L = P.r; |
---|
3313 | list K = Q.r; |
---|
3314 | list L2, K2; |
---|
3315 | int i,j, isequal; |
---|
3316 | |
---|
3317 | for(i =1; i<=size(L); i++){ |
---|
3318 | L2 = L[i]; |
---|
3319 | K2 = K[i]; |
---|
3320 | isequal = 1; |
---|
3321 | for(j=1; j<= size(L2); j++){ |
---|
3322 | if(L2[j] != K2[j]){isequal = 0; break;} |
---|
3323 | } |
---|
3324 | print("@rank "+string(i)+" |P.r[i]| = "+string(size(L[i]))+", |Q.r[i]| = "+string(size(K[i]))); |
---|
3325 | print("Elements are the same: "+string(isequal)); |
---|
3326 | } |
---|
3327 | |
---|
3328 | } |
---|
3329 | |
---|
3330 | //============================================================================// |
---|
3331 | //-------------------------------- END OF CODE -------------------------------// |
---|
3332 | //============================================================================// |
---|
3333 | |
---|
3334 | //============================================================================// |
---|
3335 | //------------------------------- UNUSED PROCEDURES --------------------------// |
---|
3336 | //============================================================================// |
---|
3337 | |
---|
3338 | |
---|
3339 | static proc printParents(arrposet P) |
---|
3340 | "USAGE: printParents(A); arr A |
---|
3341 | RETURN: [] displays the parent-lists of all the flats in the poset |
---|
3342 | REMARKS: Mainly used for debugging. |
---|
3343 | SEE ALSO: printParents, printMoebius, printFlags |
---|
3344 | EXAMPLE: example printParents;" |
---|
3345 | |
---|
3346 | { |
---|
3347 | print("Given Arrangement:"); |
---|
3348 | P.A; |
---|
3349 | print("Corresponding poset:"); |
---|
3350 | |
---|
3351 | string s; |
---|
3352 | int i,j; |
---|
3353 | list L; |
---|
3354 | for(i=1; i<= size(P.r); i++){ |
---|
3355 | print("====== rank "+string(i)+": "+string(size(P.r[i]))+" flats ======"); |
---|
3356 | s = ""; |
---|
3357 | L = P.r[i]; |
---|
3358 | for(j=1; j<=size(L); j++){ |
---|
3359 | s = s + " (" + string(L[j].parents) + "), "; |
---|
3360 | } |
---|
3361 | |
---|
3362 | print(s); |
---|
3363 | if(size(P.r[i]) == 0){break;} |
---|
3364 | } |
---|
3365 | } |
---|
3366 | |
---|
3367 | static proc printFlags(arrposet P) |
---|
3368 | "USAGE: printFlags(A); arr A |
---|
3369 | RETURN: [] displays the flag bits of all the flats in the poset |
---|
3370 | REMARKS: Mainly used for debugging. |
---|
3371 | SEE ALSO: printParents, printMoebius, printFlags |
---|
3372 | EXAMPLE: example printParents;" |
---|
3373 | { |
---|
3374 | print("Flag values:"); |
---|
3375 | |
---|
3376 | string s; |
---|
3377 | int i,j; |
---|
3378 | list L; |
---|
3379 | for(i=1; i<= size(P.r); i++){ |
---|
3380 | print("====== rank "+string(i)+": "+string(size(P.r[i]))+" flats ======"); |
---|
3381 | s = ""; |
---|
3382 | L = P.r[i]; |
---|
3383 | for(j=1; j<=size(L); j++){ |
---|
3384 | s = s + " (" + string(L[j].flag) + "), "; |
---|
3385 | } |
---|
3386 | |
---|
3387 | print(s); |
---|
3388 | if(size(P.r[i]) == 0){break;} |
---|
3389 | } |
---|
3390 | } |
---|
3391 | |
---|
3392 | // expects 2 intvecs of stric ascending order |
---|
3393 | // returns merged intvec |
---|
3394 | static proc mergeIV(intvec u, intvec v){ |
---|
3395 | intvec res; |
---|
3396 | int k; |
---|
3397 | int m = 1; |
---|
3398 | int n = 1; |
---|
3399 | for(k=1; m<=size(u) && n<=size(v); k++){ |
---|
3400 | while(1){ |
---|
3401 | if(u[m] < v[n]){res[k] = u[m]; m++; break;} |
---|
3402 | if(u[m] > v[n]){res[k] = v[n]; n++; break;} |
---|
3403 | if(u[m] == v[n]){res[k] = u[m]; m++; n++; break;} |
---|
3404 | ERROR("Something went wrong in proc mergeIV"); |
---|
3405 | } |
---|
3406 | } |
---|
3407 | while(m <= size(u)){res[k] = u[m]; m++; k++;} |
---|
3408 | while(n <= size(v)){res[k] = v[n]; n++; k++;} |
---|
3409 | return (res); |
---|
3410 | } |
---|
3411 | |
---|
3412 | // checks if v[i] < u[i] for all i |
---|
3413 | static proc isSmaller(intvec u, intvec v){ |
---|
3414 | for(int i= 1; i<=size(u); i++){ |
---|
3415 | if(u[i] > v[i]){ |
---|
3416 | return (0); |
---|
3417 | } |
---|
3418 | } |
---|
3419 | |
---|
3420 | return (1); |
---|
3421 | } |
---|
3422 | |
---|
3423 | static proc new2old(arrposet P){ |
---|
3424 | arrflats Q; |
---|
3425 | Q.A = P.A; |
---|
3426 | int i,j; |
---|
3427 | list L = P.r; |
---|
3428 | for(i=1; i<=size(L); i++){ |
---|
3429 | list L2 = L[i]; |
---|
3430 | for(j=1; j<=size(L2); j++){ |
---|
3431 | L2[j] = collectHplanes(L2[j].REL); |
---|
3432 | } |
---|
3433 | L[i] = L2; |
---|
3434 | } |
---|
3435 | |
---|
3436 | Q.r = L; |
---|
3437 | |
---|
3438 | return (Q); |
---|
3439 | } |
---|
3440 | /********************************************************************* |
---|
3441 | newstruct("arr","ideal l"); |
---|
3442 | |
---|
3443 | Defines a hyperplane arrangement by a list of linear polynomials such that the hyperplanes are the |
---|
3444 | varieties of those polynomials. |
---|
3445 | |
---|
3446 | supported operators: |
---|
3447 | A = input defines an arrangement by the input which may consist of the types arr/ideal/poly/matrix/list |
---|
3448 | A + input adds arrangement defined by the right hand side to the left hand side |
---|
3449 | A - input deletes hyperplanes from the arrangement |
---|
3450 | <, >, <=, >=, ==, != set theoretical comparison of two arrangements |
---|
3451 | A[int] access to a single hyperplane |
---|
3452 | A[intvec] subarrangement defined by the indexed hyperplanes |
---|
3453 | A; prints the arrangement on the screen |
---|
3454 | |
---|
3455 | supported type conversions: |
---|
3456 | matrix(A) returns coefficient matrix of the defining polynomials |
---|
3457 | poly(A) returns the defining polynomial which is the product of all the single polynomials |
---|
3458 | list(A) returns the defining polynomials as a list |
---|
3459 | ideal(A) returns the defining polynomials as an ideal |
---|
3460 | type2arr(input) returns arrangement defined by the input which may consist of the types arr/ideal/poly/matrix/list |
---|
3461 | |
---|
3462 | supported inherited functions: |
---|
3463 | delete(A, intvec) deletes indexed hyperplanes from the arrangement |
---|
3464 | homog(A) checks whether the defining polys are all homogeneous <=> arr is central |
---|
3465 | homog(A, rvar) homogenizes the def. polys with respect to the given ring variable |
---|
3466 | size(A) number of hyperplanes in the arrangement |
---|
3467 | subst(A,rvar,poly,...) substitutes ringvariables with given polynomials, though they need to remain linear |
---|
3468 | variables(A) ideal generated by the ring variables that A depends on |
---|
3469 | nvars(A) number of ring variables that A depends on |
---|
3470 | |
---|
3471 | |
---|
3472 | |
---|
3473 | newstruct("multarr","ideal l, intvec m"); |
---|
3474 | |
---|
3475 | Defines a hyperplane arrangement with multiplicities by a list of linear polynomials such that the hyperplanes are the |
---|
3476 | varieties of those polynomials and an intvec in which the multiplicities are saved. |
---|
3477 | |
---|
3478 | supported operators: |
---|
3479 | M = input defines an multarr by the input which may consist of the types multarr/ideal/poly/list |
---|
3480 | M + input adds arrangement defined by the right hand side to the left hand side |
---|
3481 | M; prints the multarr on the screen |
---|
3482 | |
---|
3483 | supported type conversions: |
---|
3484 | poly(M) returns defining polynomial with multiplicities |
---|
3485 | arr2multarr(A,intvec) returns multarr with multiplicities defined by the intvec. |
---|
3486 | multarr2arr(M) returns the internal arrangement (all multiplicities set to 1) |
---|
3487 | delete(M,int) decrements the multiplicity of the hyperplane defined by the index by 1 |
---|
3488 | size(M) returns number of hyperplanes counting multiplicities. |
---|
3489 | *-----------------------------------------------------------------------*/ |
---|