Changeset b6ae8c in git
 Timestamp:
 Mar 15, 2011, 7:17:52 PM (13 years ago)
 Branches:
 (u'fiekerDuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', '3720ae8bfcff4a4649ee98a15552089151d2d59b')
 Children:
 199cd68d23f96331d3327f7a0b11a46accefa1db
 Parents:
 4f7d761d43192f9c176ef0f44b0f61e891545eb8
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/multigrading.lib
r4f7d76 rb6ae8c 1 // * Mode: C++; tabwidth: 2; indenttabsmode: nil; cbasicoffset: 2 * 2 // Vimodeline: vim: filetype=c:syntax:shiftwidth=2:tabstop=8:textwidth=0:expandtab 1 3 /////////////////////////////////////////////////////////////////// 2 4 version="$Id$"; … … 5 7 LIBRARY: multigrading.lib Multigraded Rings 6 8 7 AUTHORS: Rene Birkner, rbirkner@math.fuberlin.de 9 AUTHORS: Benjamin Bechtold, benjamin.bechtold@googlemail.com 10 @* Rene Birkner, rbirkner@math.fuberlin.de 8 11 @* Lars Kastner, lkastner@math.fuberlin.de 12 @* Simon Keicher, keicher@mail.mathematik.unituebingen.de 9 13 @* Oleksandr Motsak, U@D, where U={motsak}, D={mathematik.unikl.de} 10 14 … … 12 16 OVERVIEW: This library allows one to virtually add multigradings to Singular. 13 17 For more see http://code.google.com/p/convexsingular/wiki/Multigrading 14 For theoretical references see: 18 For theoretical references see: 15 19 E. Miller, B. Sturmfels: 'Combinatorial Commutative Algebra' and 16 20 M. Kreuzer, L. Robbiano: 'Computational Commutative Algebra'. 17 21 18 22 NOTE: 'mDegBasis' relies on 4ti2 for computing Hilbert Bases. 23 All groups are finitely generated Abelian 19 24 20 25 PROCEDURES: 21 setBaseMultigrading(M, T); attach multiweights/torsionmatrices to the basering26 setBaseMultigrading(M,L); attach multiweights/grading group matrices to the basering 22 27 getVariableWeights([R]); get matrix of multidegrees of vars attached to a ring 23 getTorsion([R]); get torsion matrix attached to a ring 28 29 getGradingGroup([R]); get grading group attached to a ring 30 getLattice([R[,choice]]); get grading group' lattice attached to a ring (or its NF) 31 32 createGroup(S,L); create a group generated by S, with relations L 33 createQuotientGroup(L); create a group generated by the unit matrix whith relations L 34 createTorsionFreeGroup(S); create a group generated by S which is torsionfree 35 printGroup(G); print a group 36 37 areIsomorphicGroups(G,H); test wheter G an H are isomorphic groups (TODO Tuebingen) 38 isGroup(G); test whether G is a valid group 39 isGroupHomomorphism(L1,L2,A); test wheter A defines a group homomrphism from L1 to L2 40 41 isGradedRingHomomorphism(R,f,A); test graded ring homomorph 42 createGradedRingHomomorphism(R,f,A); create a graded ring homomorph 24 43 25 44 setModuleGrading(M,v); attach multiweights of units to a module and return it 26 45 getModuleGrading(M); get multiweights of module units (attached to M) 27 46 47 isSublattice(A,B); test whether A is a sublattice of B 48 imageLattice(P,L); computes an integral basis for the image of the lattice L under the linear map P. 49 intRank(A); computes the rank of the intmat A 50 kernelLattice(P); computes an integral basis for the kernel of the linear map P. 51 latticeBasis(B); compute an integral basis of the lattice B 52 preimageLattice(P,L); computes an integral basis for the preimage of the lattice L under the linear map P. 53 projectLattices(B); compute a linear map of lattices having the primitive span of B as its kernel. 54 intersectLattices(A,B); compute an integral basis for the intersection of the lattices A and B. 55 isIntegralSurjective(P); test whether the map P of lattices is surjective. 56 isPrimitiveSublattice(A); test whether A generates a primitive sublattice. 57 intInverse(A); compute the integral inverse matrix of the intmat A 58 intAdjoint(A,i,j); delete row i and column j of the intmat A. 59 integralSection(P); for a given linear surjective map P of lattices this procedure returns an integral section of P. 60 primitiveSpan(A); compute a basis for the minimal primitive sublattice that contains the given vectors (by A). 61 62 factorgroup(G,H); create the group G mod H 63 productgroup(G,H); create the group G x H 64 28 65 mDeg(A); compute the multidegree of A 29 66 mDegBasis(d); compute all monomials of multidegree d 30 mDegPartition(p); compute the multigradedhomogenous components of p 31 32 isTorsionFree(); test whether the current multigrading is torsion free 33 isTorsionElement(p); test whether p has zero multidegree 34 isHomogenous(a); test whether 'a' is multigradedhomogenous 67 mDegPartition(p); compute the multigradedhomogeneous components of p 68 69 isTorsionFree(); test whether the current multigrading is free 70 isPositive(); test whether the current multigrading is positive 71 isZeroElement(p); test whether p has zero multidegree 72 areZeroElements(M); test whether an integer matrix M considered as a collection of columns has zero multidegree 73 isHomogeneous(a); test whether 'a' is multigradedhomogeneous 35 74 36 75 equalMDeg(e1,e2[,V]); test whether e1==e2 in the current multigrading … … 38 77 mDegGroebner(M); compute the multigraded GB/SB of M 39 78 mDegSyzygy(M); compute the multigraded syzygies of M 79 mDegModulo(I,J); compute the multigraded 'modulo' module of I and J 40 80 mDegResolution(M,l[,m]); compute the multigraded resolution of M 41 42 defineHomogenous(p); get a torsion matrix wrt which p becomes homogenous 81 mDegTensor(m,n); compute the tensor product of multigraded modules m,n 82 mDegTor(i,m,n); compute the Tor_i(m,n) for multigraded modules m,n 83 84 defineHomogeneous(p); get a grading group wrt which p becomes homogeneous 43 85 pushForward(f); find the finest grading on the image ring, homogenizing f 44 45 hermite(A); compute the Hermite Normal Form of a matrix 86 gradiator(h); coarsens grading of the ring until h becomes homogeneous 87 88 hermiteNormalForm(A); compute the Hermite Normal Form of a matrix 89 smithNormalForm(A,#); compute matrices D,P,Q with D=P*A*Q and D is the smith normal form of A 46 90 47 91 hilbertSeries(M); compute the multigraded Hilbert Series of M 92 evalHilbertSeries(h,v); evaluate hilberts series h by substituting v[i] for t_(i) 93 94 lll(A); applies LLL(.) of lll.lib which only works for lists on a matrix A 48 95 49 96 (parameters in square brackets [] are optional) 50 97 51 KEYWORDS: multigrad eding, multidegree, multiweights, multigradedhomogenous98 KEYWORDS: multigrading, multidegree, multiweights, multigradedhomogeneous, integral linear algebra 52 99 "; 53 100 … … 56 103 57 104 LIB "standard.lib"; // for groebner 105 LIB "lll.lib"; // for lll_matrix 106 LIB "matrix.lib"; // for mDegTor 107 108 /******************************************************/ 109 110 static proc concatintmat(intmat A, intmat B) 111 { 112 113 if ( nrows(A) != nrows(B) ) 114 { 115 ERROR("matrices A and B have different number of rows."); 116 } 117 118 intmat At = transpose(A); 119 intmat Bt = transpose(B); 120 121 intmat Ct[nrows(At) + nrows(Bt)][ncols(At)] = At, Bt; 122 123 return(transpose(Ct)); 124 } 125 126 127 /******************************************************/ 128 proc createGradedRingHomomorphism(def src, ideal Im, def A) 129 "USAGE: createGradedRingHomomorphism(R, f, A); ring R, ideal f, group homomorphism A 130 PURPOSE: create a multigraded group ring homomorphism defined by 131 a ring map from R to the current ring, given by generators images f 132 and a group homomorphism A between grading groups 133 RETURN: graded ring homorphism 134 EXAMPLE: example createGradedRingHomomorphism; shows an example 135 " 136 { 137 string isGRH = "isGRH"; 138 139 if( !isGradedRingHomomorphism(def src, ideal Im, def A) ) 140 { 141 ERROR("Input data is wrong"); 142 } 143 144 list h; 145 h[3] = A; 146 147 // map f = src, Im; 148 h[2] = Im; // f? 149 h[1] = src; 150 151 attrib(h, isGRH, (1==1)); // mark it "a graded ring homomorphism" 152 153 return(h); 154 } 155 example 156 { 157 "EXAMPLE:"; echo=2; 158 159 // TODO! 160 161 } 162 163 164 /******************************************************/ 165 proc isGradedRingHomomorphism(def src, ideal Im, def A) 166 "USAGE: isGradedRingHomomorphism(R, f, A); ring R, ideal f, group homomorphism A 167 PURPOSE: test a multigraded group ring homomorphism defined by 168 a ring map from R to the current ring, given by generators images f 169 and a group homomorphism A between grading groups 170 RETURN: int, 1 for TRUE, 0 otherwise 171 EXAMPLE: example isGradedRingHomomorphism; shows an example 172 " 173 { 174 def dst = basering; 175 176 intmat result_degs = mDeg(Im); 177 print(result_degs); 178 179 setring src; 180 181 intmat input_degs = mDeg(maxideal(1)); 182 print(input_degs); 183 184 def image_degs = A * input_degs; 185 print( image_degs ); 186 187 def df = image_degs  result_degs; 188 print(df); 189 190 setring dst; 191 192 return (areZeroElements( df )); 193 } 194 example 195 { 196 "EXAMPLE:"; echo=2; 197 198 ring r = 0, (x, y, z), dp; 199 intmat S1[3][3] = 200 1, 0, 0, 201 0, 1, 0, 202 0, 0, 1; 203 intmat L1[3][1] = 204 0, 205 0, 206 0; 207 208 def G1 = createGroup(S1, L1); // (S1 + L1)/L1 209 printGroup(G1); 210 211 setBaseMultigrading(S1, L1); // to change... 212 213 ring R = 0, (a, b, c), dp; 214 intmat S2[2][3] = 215 1, 0, 216 0, 1; 217 intmat L2[2][1] = 218 0, 219 2; 220 221 def G2 = createGroup(S2, L2); 222 printGroup(G2); 223 224 setBaseMultigrading(S2, L2); // to change... 225 226 227 map F = r, a, b, c; 228 intmat A[nrows(L2)][nrows(L1)] = 229 1, 0, 0, 230 3, 2, 6; 231 232 // graded ring homomorphism is given by (compatible): 233 print(F); 234 print(A); 235 236 isGradedRingHomomorphism(r, ideal(F), A); 237 def h = createGradedRingHomomorphism(r, ideal(F), A); 238 239 print(h); 240 241 // not a homo.. 242 intmat B[nrows(L2)][nrows(L1)] = 243 1, 1, 1, 244 0, 0, 0; 245 print(B); 246 247 isGradedRingHomomorphism(r, ideal(F), B); 248 def hh = createGradedRingHomomorphism(r, ideal(F), B); 249 250 if( defined(hh) ) { ERROR("That input was not valid"); } 251 } 252 253 254 proc createQuotientGroup(intmat L) 255 " 256 L  relations 257 TODO: bad name 258 " 259 { 260 int r = nrows(L); int i; 261 intmat S[r][r]; // SQUARE!!! 262 for(i = r; i > 0; i){ S[i, i] = 1; } 263 return (createGroup(S,L)); 264 } 265 266 proc createTorsionFreeGroup(intmat S) 267 " 268 S  generators 269 TODO: bad name 270 " 271 { 272 int r = nrows(S); int i; 273 intmat L[r][1] = 0; 274 return (createGroup(S,L)); 275 } 276 277 278 /******************************************************/ 279 proc createGroup(intmat S, intmat L) 280 "USAGE: createGroup(S, L); S, L are integer matrices 281 PURPOSE: create the group of the form (S+L)/L, i.e. 282 S specifies generators, L specifies relations. 283 RETURN: group 284 EXAMPLE: example createGroup; shows an example 285 " 286 { 287 string isGroup = "isGroup"; 288 string attrGroupHNF = "hermite"; 289 string attrGroupSNF = "smith"; 290 291 292 /* 293 if( size(#) > 0 ) 294 { 295 if( typeof(#[1]) == "intmat" ) 296 { 297 intmat S = #[1]; 298 } else { ERROR("Wrong optional argument: 1"); } 299 300 if( size(#) > 1 ) 301 { 302 if( typeof(#[2]) == "intmat" ) 303 { 304 intmat L = #[2]; 305 } else { ERROR("Wrong optional argument: 2"); } 306 } 307 } 308 if( !defined(S) ) 309 {} 310 */ 311 312 if( nrows(L) != nrows(S) ) 313 { 314 ERROR("Incompatible matrices!"); 315 } 316 317 def H = attrib(L, attrGroupHNF); 318 if( !defined(H)  typeof(H) != "intmat") 319 { 320 attrib(L, attrGroupHNF, hermiteNormalForm(L)); 321 } else { kill H; } 322 323 def HH = attrib(L, attrGroupSNF); 324 if( !defined(HH)  typeof(HH) != "intmat") 325 { 326 attrib(L, attrGroupSNF, smithNormalForm(L)); 327 } else { kill HH; } 328 329 list G; // Please, note the order: Generators + Relations: 330 G[1] = S; G[2] = L; 331 332 attrib(G, isGroup, (1==1)); // mark it "a group" 333 334 return (G); 335 } 336 example 337 { 338 "EXAMPLE:"; echo=2; 339 340 intmat S[3][3] = 341 1, 0, 0, 342 0, 1, 0, 343 0, 0, 1; 344 345 intmat L[3][2] = 346 1, 1, 347 1, 3, 348 1, 5; 349 350 def G = createGroup(S, L); // (S+L)/L 351 352 printGroup(G); 353 354 kill S, L, G; 355 356 ///////////////////////////////////////////////// 357 intmat S[2][3] = 358 1, 2, 1, 359 1, 1, 0; 360 361 intmat L[2][1] = 362 0, 363 2; 364 365 def G = createGroup(S, L); // (S+L)/L 366 367 printGroup(G); 368 369 kill S, L, G; 370 371 //  extreme case  // 372 intmat S[1][3] = 373 1, 1, 10; 374 375 // Torsion: 376 intmat L[1][1] = 377 0; 378 379 def G = createGroup(S, L); // (S+L)/L 380 381 printGroup(G); 382 } 383 384 385 /******************************************************/ 386 proc printGroup(def G) 387 "USAGE: printGroup(G); G is a group 388 PURPOSE: prints the group G 389 RETURN: nothing 390 EXAMPLE: example printGroup; shows an example 391 " 392 { 393 "Generators: "; 394 print(G[1]); 395 396 "Relations: "; 397 print(G[2]); 398 399 // attrib(G[2]); 400 } 401 example 402 { 403 "EXAMPLE:"; echo=2; 404 405 } 406 407 /******************************************************/ 408 proc areIsomorphicGroups(def G, def H) 409 "USAGE: areIsomorphicGroups(G, H); G and H are groups 410 PURPOSE: ? 411 RETURN: int, 1 for TRUE, 0 otherwise 412 EXAMPLE: example areIsomorphicGroups; shows an example 413 " 414 { 415 return (1); // TRUE 416 } 417 example 418 { 419 "EXAMPLE:"; echo=2; 420 421 } 422 423 424 proc isGroup(def G) 425 "test whether G is a valid group" 426 { 427 string isGroup = "isGroup"; 428 429 // valid? 430 if( typeof(G) != "list" ){ return(0); } 431 432 def a = attrib(G, isGroup); 433 434 ///// TODO for Hans: fix attr^2 bug in Singular! 435 436 // if( !defined(a) ) { return(0); } 437 // if( typeof(a) != "int" ) { return(0); } 438 // if( !a ){ return(0); } 439 440 441 if( size(G) != 2 ){ return(0); } 442 if( typeof(G[1]) != "intmat" ){ return(0); } 443 if( typeof(G[2]) != "intmat" ){ return(0); } 444 if( nrows(G[1]) != nrows(G[2]) ){ return(0); } 445 446 return(1==1); 447 } 448 449 58 450 59 451 /******************************************************/ 60 452 proc setBaseMultigrading(intmat M, list #) 61 "USAGE: setBaseMultigrading(M[, T]); M, T are integer matrices62 PURPOSE: attaches weights of variables and torsionto the basering.453 "USAGE: setBaseMultigrading(M[, G]); M is an intege matrix, G is a group (or lattice) 454 PURPOSE: attaches weights of variables and grading group to the basering. 63 455 NOTE: M encodes the weights of variables columnwise. 64 The torsion is given by the lattice spanned by the columns of the integer65 matrix T in Z^nrows(M) over Z.66 456 RETURN: nothing 67 457 EXAMPLE: example setBaseMultigrading; shows an example … … 69 459 { 70 460 string attrMgrad = "mgrad"; 71 string attrTorsion = "torsion"; 72 string attrTorsionHNF = "torsion_hermite"; 73 74 75 attrib(basering, attrMgrad, M); 76 77 if( size(#) > 0 and typeof(#[1]) == "intmat" ) 78 { 79 def T = #[1]; 80 } else 81 { 82 intmat T[nrows(M)][1]; 83 } 84 85 if( nrows(T) == nrows(M) ) 86 { 87 attrib(basering, attrTorsion, T); 88 // def H; 89 // attrib(basering, attrTorsionHNF, H); 90 } 461 string attrGradingGroup = "gradingGroup"; 462 463 if( size(#) > 0 ) 464 { 465 if( typeof(#[1]) == "intmat" ) 466 { 467 def L = createGroup(M, #[1]); 468 } 469 470 if( isGroup(#[1]) ) 471 { 472 def L = #[1]; 473 474 if( !isSublattice(M, L[1]) ) 475 { 476 ERROR("Multigrading is not contained in the grading group!"); 477 } 478 } 479 } 91 480 else 92 481 { 93 ERROR("Incompatible matrix sizes!"); 94 } 482 def L = createTorsionFreeGroup(M); 483 } 484 485 if( !defined(L) ){ ERROR("Wrong arguments: no group given?"); } 486 487 488 489 attrib(basering, attrMgrad, M); 490 attrib(basering, attrGradingGroup, L); 491 492 ideal Q = ideal(basering); 493 if( !isHomogeneous(Q) ) // easy now, but would be hard before setting ring attributes! 494 { 495 "Warning: your quotient ideal is not homogenous (multigrading was set anyway)!"; 496 } 497 498 499 500 95 501 } 96 502 example 97 503 { 98 504 "EXAMPLE:"; echo=2; 99 505 100 506 ring R = 0, (x, y, z), dp; 101 507 … … 106 512 0, 0, 1; 107 513 108 // Torsion:514 // GradingGroup: 109 515 intmat L[3][2] = 110 516 1, 1, 111 1, 3, 517 1, 3, 112 518 1, 5; 113 519 114 520 // attaches M & L to R (==basering): 115 521 setBaseMultigrading(M, L); // Grading: Z^3/L … … 117 523 // Weights are accessible via "getVariableWeights()": 118 524 getVariableWeights(); 119 120 // Test all possible usages: 525 526 // Test all possible usages: 121 527 (getVariableWeights() == M) && (getVariableWeights(R) == M) && (getVariableWeights(basering) == M); 122 528 123 // Torsion is accessible via "getTorsion()":124 get Torsion();125 126 // Test all possible usages: 127 (get Torsion() == L) && (getTorsion(R) == L) && (getTorsion(basering) == L);128 129 // And its hermite NF via get Torsion("hermite"):130 get Torsion("hermite");131 132 // Test all possible usages: 133 intmat H = hermite (L);134 (get Torsion("hermite") == H) && (getTorsion(R, "hermite") == H) && (getTorsion(basering, "hermite") == H);529 // Grading group is accessible via "getLattice()": 530 getLattice(); 531 532 // Test all possible usages: 533 (getLattice() == L) && (getLattice(R) == L) && (getLattice(basering) == L); 534 535 // And its hermite NF via getLattice("hermite"): 536 getLattice("hermite"); 537 538 // Test all possible usages: 539 intmat H = hermiteNormalForm(L); 540 (getLattice("hermite") == H) && (getLattice(R, "hermite") == H) && (getLattice(basering, "hermite") == H); 135 541 136 542 kill L, M; … … 147 553 0, 148 554 2; 149 555 150 556 // attaches M & L to R (==basering): 151 557 setBaseMultigrading(M, L); // Grading: Z + (Z/2Z) … … 154 560 getVariableWeights() == M; 155 561 156 // Torsion is accessible via "get Torsion()":157 get Torsion() == L;562 // Torsion is accessible via "getLattice()": 563 getLattice() == L; 158 564 159 565 kill L, M; … … 167 573 intmat L[1][1] = 168 574 0; 169 575 170 576 // attaches M & L to R (==basering): 171 577 setBaseMultigrading(M); // Grading: Z^3 … … 174 580 getVariableWeights() == M; 175 581 176 // Torsion is accessible via "get Torsion()":177 get Torsion() == L;582 // Torsion is accessible via "getLattice()": 583 getLattice() == L; 178 584 } 179 585 … … 208 614 def M = attrib(R, attrMgrad); 209 615 if( typeof(M) == "intmat"){ return (M); } 210 ERROR( "Sorry no multigrading matrix!" ); 616 ERROR( "Sorry no multigrading matrix!" ); 211 617 } 212 618 example 213 619 { 214 620 "EXAMPLE:"; echo=2; 215 621 216 622 ring R = 0, (x, y, z), dp; 217 623 … … 222 628 0, 0, 1; 223 629 224 // Torsion:630 // Grading group: 225 631 intmat L[3][2] = 226 632 1, 1, 227 1, 3, 633 1, 3, 228 634 1, 5; 229 635 230 636 // attaches M & L to R (==basering): 231 637 setBaseMultigrading(M, L); // Grading: Z^3/L … … 243 649 1, 1, 0; 244 650 245 // Torsion:651 // Grading group: 246 652 intmat L[2][1] = 247 653 0, 248 654 2; 249 655 250 656 // attaches M & L to R (==basering): 251 657 setBaseMultigrading(M, L); // Grading: Z + (Z/2Z) … … 262 668 1, 1, 10; 263 669 264 // Torsion:670 // Grading group: 265 671 intmat L[1][1] = 266 672 0; 267 673 268 674 // attaches M & L to R (==basering): 269 675 setBaseMultigrading(M); // Grading: Z^3 … … 273 679 } 274 680 275 /******************************************************/ 276 proc getTorsion(list #) 277 "USAGE: getTorsion([R[,opt]]) 278 PURPOSE: get associated torsion matrix, i.e. generators (cols) of the Torsion group 279 RETURN: intmat, the torsion matrix, or its hermite normal form 280 if an optional argument (\"hermite\") is given 281 EXAMPLE: example getTorsion; shows an example 282 " 283 { 284 string attrTorsion = "torsion"; 285 string attrTorsionHNF = "torsion_hermite"; 681 682 proc getGradingGroup(list #) 683 "USAGE: getGradingGroup([R]) 684 PURPOSE: get associated grading group 685 RETURN: group, the grading group 686 EXAMPLE: example getGradingGroup; shows an example 687 " 688 { 689 string attrGradingGroup = "gradingGroup"; 286 690 287 691 int i = 1; … … 293 697 def R = #[i]; 294 698 i++; 295 } 296 } 699 } 700 } 297 701 298 702 if( !defined(R) ) … … 301 705 } 302 706 303 if( size(#) >= i ) 304 { 305 if( #[i] == "hermite" ) 306 { 307 if( typeof(attrib(R, attrTorsionHNF)) != "intmat" ) 308 { 309 def M = getTorsion(R); 310 if( typeof(M) != "intmat") 311 { 312 ERROR( "Sorry no torsion matrix!" ); 313 } 314 M = hermite(M); 315 attrib(R, attrTorsionHNF, M); // this might not work with R != basering... 316 } 317 return (attrib(R, attrTorsionHNF)); 318 } 319 } 320 321 def M = attrib(R, attrTorsion); 322 if( typeof(M) != "intmat") 323 { 324 ERROR( "Sorry no torsion matrix!" ); 325 } 326 return (M); 707 def G = attrib(R, attrGradingGroup); 708 709 if( !isGroup(G) ) 710 { 711 ERROR("Sorry no grading group!"); 712 } 713 714 return(G); 327 715 } 328 716 example 329 717 { 330 718 "EXAMPLE:"; echo=2; 331 719 332 720 ring R = 0, (x, y, z), dp; 333 721 … … 341 729 intmat L[3][2] = 342 730 1, 1, 343 1, 3, 731 1, 3, 344 732 1, 5; 345 733 346 734 // attaches M & L to R (==basering): 347 735 setBaseMultigrading(M, L); // Grading: Z^3/L 348 736 349 // Torsion is accessible via "getTorsion()":350 getTorsion() == L; 351 352 // its hermite NF:353 print(getTorsion("hermite"));354 355 kill L, M ;737 def G = getGradingGroup(); 738 739 printGroup( G ); 740 741 G[1] == M; G[2] == L; 742 743 kill L, M, G; 356 744 357 745 //  isomorphic multigrading  // … … 366 754 0, 367 755 2; 368 756 369 757 // attaches M & L to R (==basering): 370 758 setBaseMultigrading(M, L); // Grading: Z + (Z/2Z) 371 759 372 // Torsion is accessible via "getTorsion()": 373 getTorsion() == L; 374 375 // its hermite NF: 376 print(getTorsion("hermite")); 377 378 kill L, M; 379 760 def G = getGradingGroup(); 761 762 printGroup( G ); 763 764 G[1] == M; G[2] == L; 765 766 kill L, M, G; 380 767 //  extreme case  // 381 768 … … 387 774 intmat L[1][1] = 388 775 0; 389 776 390 777 // attaches M & L to R (==basering): 391 778 setBaseMultigrading(M); // Grading: Z^3 392 779 393 // Torsion is accessible via "getTorsion()": 394 getTorsion() == L; 780 def G = getGradingGroup(); 781 782 printGroup( G ); 783 784 G[1] == M; G[2] == L; 785 786 kill L, M, G; 787 } 788 789 790 /******************************************************/ 791 proc getLattice(list #) 792 "USAGE: getLattice([R[,opt]]) 793 PURPOSE: get associated grading group matrix, i.e. generators (cols) of the grading group 794 RETURN: intmat, the grading group matrix, or 795 its hermite normal form if an optional argument (\"hermiteNormalForm\") is given or 796 smith normal form if an optional argument (\"smith\") is given 797 EXAMPLE: example getLattice; shows an example 798 " 799 { 800 int i = 1; 801 if( size(#) >= i ) 802 { 803 if( ( typeof(#[i]) == "ring" ) or ( typeof(#[i]) == "qring" ) ) 804 { 805 i++; 806 } 807 } 808 809 string attrGradingGroupHNF = "hermite"; 810 string attrGradingGroupSNF = "smith"; 811 812 def G = getGradingGroup(#); 813 814 // printGroup(G); 815 816 817 818 def T = G[2]; 819 820 if( size(#) >= i ) 821 { 822 if( #[i] == "hermite" ) 823 { 824 def M = attrib(T, attrGradingGroupHNF); 825 if( (!defined(M)) or (typeof(M) != "intmat") ) 826 { 827 M = hermiteNormalForm(T); 828 } 829 return (M); 830 } 831 832 if( #[i] == "smith" ) 833 { 834 def M = attrib(T, attrGradingGroupSNF); 835 if( (!defined(M)) or (typeof(M) != "intmat") ) 836 { 837 M = smithNormalForm(T); 838 } 839 return (M); 840 } 841 } 842 843 return(T); 844 } 845 example 846 { 847 "EXAMPLE:"; echo=2; 848 849 ring R = 0, (x, y, z), dp; 850 851 // Weights of variables 852 intmat M[3][3] = 853 1, 0, 0, 854 0, 1, 0, 855 0, 0, 1; 856 857 // Torsion: 858 intmat L[3][2] = 859 1, 1, 860 1, 3, 861 1, 5; 862 863 // attaches M & L to R (==basering): 864 setBaseMultigrading(M, L); // Grading: Z^3/L 865 866 // Torsion is accessible via "getLattice()": 867 getLattice() == L; 395 868 396 869 // its hermite NF: 397 print(getTorsion("hermite")); 398 } 870 print(getLattice("hermite")); 871 872 kill L, M; 873 874 //  isomorphic multigrading  // 875 876 // Weights of variables 877 intmat M[2][3] = 878 1, 2, 1, 879 1, 1, 0; 880 881 // Torsion: 882 intmat L[2][1] = 883 0, 884 2; 885 886 // attaches M & L to R (==basering): 887 setBaseMultigrading(M, L); // Grading: Z + (Z/2Z) 888 889 // Torsion is accessible via "getLattice()": 890 getLattice() == L; 891 892 // its hermite NF: 893 print(getLattice("hermite")); 894 895 kill L, M; 896 897 //  extreme case  // 898 899 // Weights of variables 900 intmat M[1][3] = 901 1, 1, 10; 902 903 // Torsion: 904 intmat L[1][1] = 905 0; 906 907 // attaches M & L to R (==basering): 908 setBaseMultigrading(M); // Grading: Z^3 909 910 // Torsion is accessible via "getLattice()": 911 getLattice() == L; 912 913 // its hermite NF: 914 print(getLattice("hermite")); 915 } 916 917 proc getGradedGenerator(def m, int i) 918 " 919 returns m[i], but with grading 920 " 921 { 922 if( typeof(m) == "ideal" ) 923 { 924 return (m[i]); 925 } 926 927 if( typeof(m) == "module" ) 928 { 929 def v = getModuleGrading(m); 930 931 return ( setModuleGrading(m[i],v) ); 932 } 933 934 ERROR("m is expected to be an ideal or a module"); 935 } 936 399 937 400 938 /******************************************************/ … … 410 948 411 949 def V = attrib(m, attrModuleGrading); 412 950 413 951 if( typeof(V) != "intmat" ) 414 952 { … … 419 957 return (VV); 420 958 } 421 959 422 960 ERROR("Sorry: vector or module need modulegradingmatrix! See 'getModuleGrading'."); 423 961 } … … 432 970 ERROR("Sorry wrong width of V: " + string(ncols(V))); 433 971 } 434 972 435 973 return (V); 436 974 } … … 438 976 { 439 977 "EXAMPLE:"; echo=2; 440 978 441 979 ring R = 0, (x,y), dp; 442 980 intmat M[2][2]= … … 446 984 1, 2, 3, 4, 0, 447 985 0, 10, 20, 30, 1; 448 986 449 987 setBaseMultigrading(M, T); 450 988 451 989 ideal I = x, y, xy^5; 452 isHomogen ous(I);453 990 isHomogeneous(I); 991 454 992 intmat V = mDeg(I); print(V); 455 993 456 994 module S = syz(I); print(S); 457 995 458 996 S = setModuleGrading(S, V); 459 997 460 998 getModuleGrading(S) == V; 461 462 vector v = setModuleGrading(S[1], V);999 1000 vector v = getGradedGenerator(S, 1); 463 1001 getModuleGrading(v) == V; 464 isHomogen ous(v);465 print( mDeg(v) ); 466 467 isHomogen ous(S);1002 isHomogeneous(v); 1003 print( mDeg(v) ); 1004 1005 isHomogeneous(S); 468 1006 print( mDeg(S) ); 469 1007 } … … 474 1012 PURPOSE: attaches the multiweights of free module generators to 'm' 475 1013 WARNING: The method does not verify whether the multigrading makes the 476 module/vector homogen ous. One can do that using isHomogenous(m).1014 module/vector homogeneous. One can do that using isHomogeneous(m). 477 1015 EXAMPLE: example setModuleGrading; shows an example 478 1016 " … … 484 1022 if(nrows(G) != nrows(R)){ ERROR("Incompatible gradings.");} 485 1023 if(ncols(G) < nrows(m)){ ERROR("Multigrading does not fit to module.");} 486 1024 487 1025 attrib(m, attrModuleGrading, G); 488 1026 return(m); … … 491 1029 { 492 1030 "EXAMPLE:"; echo=2; 493 1031 494 1032 ring R = 0, (x,y), dp; 495 1033 intmat M[2][2]= … … 499 1037 1, 2, 3, 4, 0, 500 1038 0, 10, 20, 30, 1; 501 1039 502 1040 setBaseMultigrading(M, T); 503 1041 504 1042 ideal I = x, y, xy^5; 505 1043 intmat V = mDeg(I); 506 1044 507 1045 // V == M; modulo T 508 1046 print(V); 509 1047 510 1048 module S = syz(I); 511 1049 512 1050 S = setModuleGrading(S, V); 513 1051 getModuleGrading(S) == V; 514 1052 515 1053 print(S); 516 517 vector v = S[1]; v = setModuleGrading(v, V);1054 1055 vector v = getGradedGenerator(S, 1); 518 1056 getModuleGrading(v) == V; 519 1057 520 print( mDeg(v) ); 521 522 isHomogen ous(S);1058 print( mDeg(v) ); 1059 1060 isHomogeneous(S); 523 1061 524 1062 print( mDeg(S) ); … … 526 1064 527 1065 528 1066 proc mDegTensor(module m, module n){ 1067 matrix M = m; 1068 matrix N = n; 1069 intmat gm = getModuleGrading(m); 1070 intmat gn = getModuleGrading(n); 1071 int grows = nrows(gm); 1072 int mr = nrows(M); 1073 int mc = ncols(M); 1074 if(rank(M) == 0){ mc = 0;} 1075 int nr = nrows(N); 1076 int nc = ncols(N); 1077 if(rank(N) == 0){ nc = 0;} 1078 intmat gresult[nrows(gm)][mr*nr]; 1079 matrix result[mr*nr][mr*nc+mc*nr]; 1080 int i, j; 1081 int column = 1; 1082 for(i = 1; i<=mr; i++){ 1083 for(j = 1; j<=nr; j++){ 1084 gresult[1..grows,(i1)*nr+j] = gm[1..grows,i]+gn[1..grows,j]; 1085 } 1086 } 1087 //gresult; 1088 if( nc!=0 ){ 1089 for(i = 1; i<=mr; i++) 1090 { 1091 result[((i1)*nr+1)..(i*nr),((i1)*nc+1)..(i*nc)] = N[1..nr,1..nc]; 1092 } 1093 } 1094 list rownumbers, colnumbers; 1095 //print(result); 1096 if( mc!=0 ){ 1097 for(j = 1; j<=nr; j++) 1098 { 1099 rownumbers = nr*(0..(mr1))+j*(1:mr); 1100 colnumbers = ((mr*nc+j):mc)+nr*(0..(mc1)); 1101 result[rownumbers[1..mr],colnumbers[1..mc] ] = M[1..mr,1..mc]; 1102 } 1103 } 1104 module res = result; 1105 res = setModuleGrading(res, gresult); 1106 //getModuleGrading(res); 1107 return(res); 1108 } 1109 example 1110 { 1111 "EXAMPLE: ";echo=2; 1112 ring r = 0,(x),dp; 1113 intmat g[2][1]=1,1; 1114 setBaseMultigrading(g); 1115 matrix m[5][3]=1,2,3,4,5,6,7,8,9,10,11,12,13,14,15; 1116 matrix n[3][2]=x,x2,x3,x4,x5,x6; 1117 module mm = m; 1118 module nn = n; 1119 intmat gm[2][5]=1,2,3,4,5,0,0,0,0,0; 1120 intmat gn[2][3]=0,0,0,1,2,3; 1121 mm = setModuleGrading(mm, gm); 1122 nn = setModuleGrading(nn, gn); 1123 module mmtnn = mDegTensor(mm, nn); 1124 print(mmtnn); 1125 getModuleGrading(mmtnn); 1126 LIB "homolog.lib"; 1127 module tt = tensorMod(mm,nn); 1128 print(tt); 1129 1130 kill m, mm, n, nn, gm, gn; 1131 1132 matrix m[7][3] = x, x1,x+2, 3x, 4x, x5, x6, x7, x8, 9, 10, 11x, 12 x, 13x, 14x, x15, (x4)^2, x17, 18x, 19x, 20x, 21x; 1133 matrix n[2][4] = 1, 2, 3, 4, x, x2, x3, x4; 1134 module mm = m; 1135 module nn = n; 1136 print(mm); 1137 print(nn); 1138 intmat gm[2][7] = 1, 2, 3, 4, 5, 6, 7, 0, 0, 0, 0, 0, 0, 0; 1139 intmat gn[2][2] = 0, 0, 1, 2; 1140 mm = setModuleGrading(mm, gm); 1141 nn = setModuleGrading(nn, gn); 1142 module mmtnn = mDegTensor(mm, nn); 1143 print(mmtnn); 1144 getModuleGrading(mmtnn); 1145 matrix a = mmtnn; 1146 matrix b = tensorMod(mm, nn); 1147 print(ab); 1148 1149 } 1150 1151 proc mDegTor(int i, module m, module n) 1152 { 1153 def res = mDegResolution(n, 0, 1); 1154 //print(res); 1155 list l = res; 1156 if(size(l)<i){ return(0);} 1157 else 1158 { 1159 1160 matrix fd[nrows(m)][0]; 1161 matrix fd2[nrows(l[i+1])][0]; 1162 matrix fd3[nrows(l[i])][0]; 1163 1164 module freedim = fd; 1165 module freedim2 = fd2; 1166 module freedim3 = fd3; 1167 1168 freedim = setModuleGrading(freedim,getModuleGrading(m)); 1169 freedim2 = setModuleGrading(freedim2,getModuleGrading(l[i+1])); 1170 freedim3 = setModuleGrading(freedim3, getModuleGrading(l[i])); 1171 1172 module mimag = mDegTensor(freedim3, m); 1173 //"mimag ok."; 1174 module mf = mDegTensor(l[i], freedim); 1175 //"mf ok."; 1176 module mim1 = mDegTensor(freedim2 ,m); 1177 module mim2 = mDegTensor(l[i+1],freedim); 1178 //"mim1+2 ok."; 1179 module mker = mDegModulo(mf,mimag); 1180 //"mker ok."; 1181 module mim = mim1,mim2; 1182 mim = setModuleGrading(mim, getModuleGrading(mim1)); 1183 //"mim: r: ",nrows(mim)," c: ",ncols(mim); 1184 //"mim1: r: ",nrows(mim1)," c: ",ncols(mim1); 1185 //"mim2: r: ",nrows(mim2)," c: ",ncols(mim2); 1186 //matrix mimmat = mim; 1187 //matrix mimmat1[16][4]=mimmat[1..16,25..28]; 1188 //print(mimmat1matrix(mim2)); 1189 return(mDegModulo(mker,mim)); 1190 //return(0); 1191 } 1192 return(0); 1193 } 1194 example 1195 { 1196 "EXAMPLE: ";echo=2; 1197 LIB "homolog.lib"; 1198 ring r = 0,(x_(1..4)),dp; 1199 intmat g[2][4]=1,1,0,0,0,1,1,1; 1200 setBaseMultigrading(g); 1201 ideal i = maxideal(1); 1202 module m = mDegSyzygy(i); 1203 module rt = Tor(2,m,m); 1204 module mDegT = mDegTor(2,m,m); 1205 print(matrix(rt)matrix(mDegT)); 1206 /* 1207 ring r = 0,(x),dp; 1208 intmat g[2][1]=1,1; 1209 setBaseMultigrading(g); 1210 matrix m[5][3]=1,2,3,4,5,6,7,8,9,10,11,12,13,14,15; 1211 matrix n[3][2]=x,x2,x3,x4,x5,x6; 1212 module mm = m; 1213 module nn = n; 1214 intmat gm[2][5]=1,1,1,1,1,1,1,1,1,1,1; 1215 intmat gn[2][3]=0,2,4,0,2,4; 1216 mm = setModuleGrading(mm, gm); 1217 nn = setModuleGrading(nn, gn); 1218 isHomogeneous(mm,"checkGens"); 1219 isHomogeneous(nn,"checkGens"); 1220 mDegTor(1,mm, nn); 1221 1222 kill m, mm, n, nn, gm, gn; 1223 1224 matrix m[7][3] = x, x1,x+2, 3x, 4x, x5, x6, x7, x8, 9, 10, 11x, 12 x, 13x, 14x, x15, (x4)^2, x17, 18x, 19x, 20x, 21x; 1225 matrix n[2][4] = 1, 2, 3, 4, x, x2, x3, x4; 1226 module mm = m; 1227 module nn = n; 1228 print(mm); 1229 print(nn); 1230 intmat gm[2][7] = 1, 2, 3, 4, 5, 6, 7, 0, 0, 0, 0, 0, 0, 0; 1231 intmat gn[2][2] = 0, 0, 1, 2; 1232 mm = setModuleGrading(mm, gm); 1233 nn = setModuleGrading(nn, gn); 1234 module mmtnn = mDegTensor(mm, nn); 1235 */ 1236 } 1237 1238 1239 /******************************************************/ 1240 proc isGroupHomomorphism(def L1, def L2, intmat A) 1241 "USAGE: gisGoupHomomorphism(L1,L2,A); L1 and L2 are groups, A is an integer matrix 1242 PURPOSE: checks whether A defines a group homomorphism phi: L1 > L2 1243 RETURN: int, 1 if A defines the homomorphism and 0 otherwise 1244 EXAMPLE: example isGroupHomomorphism; shows an example 1245 " 1246 { 1247 // TODO: L1, L2 1248 if( (ncols(A) != nrows(L1)) or (nrows(A) != nrows(L2)) ) 1249 { 1250 ERROR("Incompatible sizes!"); 1251 } 1252 1253 intmat im = A * L1; 1254 1255 return (areZeroElements(im, L2)); 1256 } 1257 example 1258 { 1259 "EXAMPLE:"; echo=2; 1260 1261 intmat L1[4][1]= 1262 0, 1263 0, 1264 0, 1265 2; 1266 1267 intmat L2[3][2]= 1268 0, 0, 1269 2, 0, 1270 0, 3; 1271 1272 intmat A[3][4] = 1273 1, 2, 3, 0, 1274 7, 0, 0, 0, 1275 1, 2, 0, 3; 1276 print( A ); 1277 1278 isGroupHomomorphism(L1, L2, A); 1279 1280 intmat B[3][4] = 1281 1, 2, 3, 0, 1282 7, 0, 0, 0, 1283 1, 2, 0, 2; 1284 print( B ); 1285 1286 isGroupHomomorphism(L1, L2, B); // Not a homomorphism! 1287 } 529 1288 530 1289 /******************************************************/ 531 1290 proc isTorsionFree() 532 1291 "USAGE: isTorsionFree() 533 PURPOSE: Determines whether the multigrading attached to the current ring is torsionfree.1292 PURPOSE: Determines whether the multigrading attached to the current ring is free. 534 1293 RETURN: boolean, the result of the test 535 1294 EXAMPLE: example isTorsionFree; shows an example 536 1295 " 537 1296 { 538 539 intmat H = hermite(transpose(getTorsion("hermite"))); // TODO: ?cache it? //****** 1297 intmat H = smithNormalForm(getLattice()); // TODO: ?cache it? //****** 540 1298 541 1299 int i, j; … … 550 1308 551 1309 if(H[j, i]!=0) 552 { 1310 { 553 1311 d=d*H[j, i]; 554 1312 } … … 556 1314 557 1315 if( (d*d)==1 ) 558 { 1316 { 559 1317 return(1==1); 560 1318 } … … 572 1330 1, 2, 3, 4, 0, 573 1331 0,10,20,30, 1; 574 1332 575 1333 setBaseMultigrading(M,T); 576 577 // Is the resulting group torsionfree?1334 1335 // Is the resulting group free? 578 1336 isTorsionFree(); 579 1337 … … 582 1340 583 1341 ring R=0,(x,y,z),dp; 584 intmat A[3][3] = 1342 intmat A[3][3] = 585 1343 1,0,0, 586 1344 0,1,0, … … 591 1349 1,2,0,3; 592 1350 setBaseMultigrading(A,B); 593 // Is the resulting group torsionfree?1351 // Is the resulting group free? 594 1352 isTorsionFree(); 595 1353 … … 598 1356 599 1357 1358 static proc gcdcomb(int a, int b) 1359 { 1360 // a; 1361 // b; 1362 intvec av = a,1,0; 1363 intvec bv = b,0,1; 1364 intvec save; 1365 while(av[1]*bv[1] != 0) 1366 { 1367 bv = bv  (bv[1]  bv[1]%av[1])/av[1] * av; 1368 save = bv; 1369 bv = av; 1370 av = save; 1371 } 1372 if(bv[1] < 0) 1373 { 1374 bv = bv; 1375 } 1376 return(bv); 1377 } 1378 1379 1380 proc lll(def A) 1381 " 1382 The lll algorithm of lll.lib only works for lists of vectors. 1383 Maybe one should rescript it for matrices. This method will 1384 convert a matrix to a list, plug it into lll and make the result 1385 a matrix and return it. 1386 " 1387 { 1388 if(typeof(A) == "list") 1389 { 1390 int sizeA= size (A); 1391 if (sizeA == 0) 1392 { 1393 return (A); 1394 } 1395 if (typeof (A [1]) != "intvec") 1396 { 1397 ERROR("Unrecognized type."); 1398 } 1399 int columns= size (A [1]); 1400 int i; 1401 for (i= 2; i <= sizeA; i++) 1402 { 1403 if (typeof (A[i]) != "intvec") 1404 { 1405 ERROR("Unrecognized type."); 1406 } 1407 if (size (A [i]) != columns) 1408 { 1409 ERROR ("expected equal dimension"); 1410 } 1411 } 1412 int j; 1413 intmat m [columns] [sizeA]; 1414 for (i= 1; i <= sizeA; i++) 1415 { 1416 for (j= 1; j <= columns; j++) 1417 { 1418 m[i,j]= A[i] [j]; 1419 } 1420 } 1421 m= system ("LLL", m); 1422 list result= list(); 1423 1424 for (i= 1; i <= sizeA; i++) 1425 { 1426 intvec buf= intvec (m[i , 1]); 1427 for (j= 2; j <= columns; j++) 1428 { 1429 buf= buf,intvec (m [i, j]); 1430 } 1431 result= result+ list (buf); 1432 1433 } 1434 return(result); 1435 } 1436 else 1437 { 1438 if(typeof(A) == "intmat") 1439 { 1440 A= system ("LLL", A); 1441 return(A); 1442 } 1443 else 1444 { 1445 ERROR("Unrecognized type."); 1446 } 1447 } 1448 } 1449 1450 example 1451 { 1452 1453 "EXAMPLE:"; 1454 1455 ring R = 0,x,dp; 1456 intmat m[5][5]=13,25,37,83,294,12,33,9,0,64,77,12,34,6,1,43,2,88,91,100,46,32,37,42,15; 1457 lll(m); 1458 list l=intvec(13,25,37, 83, 294),intvec(12, 33, 9,0,64), intvec (77,12,34,6,1), intvec (43,2,88,91,100), intvec (46,32,37,42,15); 1459 lll(l); 1460 } 1461 1462 1463 proc smithNormalForm(intmat A, list #) 1464 " 1465 This method returns 3 Matrices P, D and Q such that D = P*A*Q. 1466 WARNING: This might not be what you expect. 1467 " 1468 { 1469 list l1 = hermiteNormalForm(A, 5); 1470 // l1; 1471 intmat B = transpose(l1[1]); 1472 list l2 = hermiteNormalForm(B, 5); 1473 // l2; 1474 intmat P = transpose(l2[2]); 1475 intmat D = transpose(l2[1]); 1476 intmat Q = l1[2]; 1477 int cc = ncols(D); 1478 int rr = nrows(D); 1479 intmat transform; 1480 int k = 1; 1481 int a, b, c; 1482 // D; 1483 intvec v; 1484 if((cc==1)(rr==1)){ 1485 if(size(#)==0) 1486 { 1487 return(D); 1488 } else 1489 { 1490 return(list(P,D,Q)); 1491 } 1492 } 1493 while(D[k+1,k+1] !=0){ 1494 if(D[k+1,k+1]%D[k,k]!=0){ 1495 b = D[k, k]; c = D[k+1, k+1]; 1496 v = gcdcomb(D[k,k],D[k+1,k+1]); 1497 transform = unitMatrix(cc); 1498 transform[k+1,k] = 1; 1499 a = v[3]*D[k+1,k+1]/v[1]; 1500 transform[k, k+1] = a; 1501 transform[k+1, k+1] = a+1; 1502 //det(transform); 1503 D = D*transform; 1504 Q = Q*transform; 1505 //D; 1506 transform = unitMatrix(rr); 1507 transform[k,k] = v[2]; 1508 transform[k,k+1] = v[3]; 1509 transform[k+1,k] = c/v[1]; 1510 transform[k+1,k+1] = b/v[1]; 1511 D = transform * D; 1512 P = transform * P; 1513 //" "; 1514 //D; 1515 //"small transform: ", det(transform); 1516 //transform; 1517 k=0; 1518 } 1519 k++; 1520 if((k==rr)  (k==cc)){ 1521 break; 1522 } 1523 } 1524 //"here is the size ",size(#); 1525 if(size(#) == 0){ 1526 return(D); 1527 } else { 1528 return(list(P, D, Q)); 1529 } 1530 } 1531 example 1532 { 1533 "EXAMPLE: "; echo=2; 1534 1535 intmat A[5][7] = 1536 1,0,1,0,2,9,71, 1537 0,24,248,32,96,448,3496, 1538 0,4,42,4,8,30,260, 1539 0,0,0,18,90,408,3168, 1540 0,0,0,32,224,1008,7872; 1541 1542 list l = smithNormalForm(A, 5); 1543 1544 l; 1545 l[1]*A*l[3]; 1546 det(l[1]); 1547 det(l[3]); 1548 } 1549 600 1550 601 1551 /******************************************************/ 602 proc hermite (intmat A)603 "USAGE: hermite ( A );604 PURPOSE: Computes the (lower triangular) Hermite Normal Form 1552 proc hermiteNormalForm(intmat A, list #) 1553 "USAGE: hermiteNormalForm( A ); 1554 PURPOSE: Computes the (lower triangular) Hermite Normal Form 605 1555 of the matrix A by column operations. 606 1556 RETURN: intmat, the Hermite Normal Form of A 607 EXAMPLE: example hermite; shows an example 608 " 609 { 610 intmat g; 611 int i,j,k, save, d, a1, a2, c, l; 612 c=0; 613 l=0; 614 for(i=1; ((i+l)<=nrows(A))&&((i+l)<=ncols(A)); i++) 615 { 616 //i; 617 if(A[i+l, i]==0) 618 { 619 for(j=i;j<=ncols(A); j++) 620 { 621 if(A[i+l, j]!=0) 1557 EXAMPLE: example hermiteNormalForm; shows an example 1558 " 1559 { 1560 1561 int row, column, i, j; 1562 int rr = nrows(A); 1563 int cc = ncols(A); 1564 intvec savev, gcdvec, v1, v2; 1565 intmat q = unitMatrix(cc); 1566 intmat transform; 1567 column = 1; 1568 for(row = 1; (row<=rr)&&(column<=cc); row++) 1569 { 1570 if(A[row,column]==0) 622 1571 { 623 for(k=1;k<=nrows(A);k++) 624 { 625 save=A[k, i]; 626 A[k, i]=A[k, j]; 627 A[k, j]=save; 628 } 629 break; 1572 for(j = column; j<=cc; j++) 1573 { 1574 if(A[row, j]!=0) 1575 { 1576 transform = unitMatrix(cc); 1577 transform[j,j] = 0; 1578 transform[column, column] = 0; 1579 transform[column,j] = 1; 1580 transform[j,column] = 1; 1581 q = q*transform; 1582 A = A*transform; 1583 break; 1584 } 1585 } 630 1586 } 631 if(j==ncols(A)){ c=1; l=l+1; } 632 } 633 } 634 635 if((i+l)>nrows(A)){ break; } 636 637 if(A[i+l, i]<0) 638 { 639 for(k=1;k<=nrows(A);k++){ A[k, i]=A[k, i]*(1); } 640 } 641 642 if(c==0) 643 { 644 for(j=i+1;j<=ncols(A);j++) 645 { 646 //A; 647 if(A[i+l, j]<0) 1587 if(A[row,column] == 0) 648 1588 { 649 for(k=1;k<=nrows(A);k++){ A[k, j]=(1)*A[k, j];} 1589 row++; 1590 continue; 650 1591 } 651 652 if(A[i+l, i]==1) 1592 for(j = column+1; j<=cc; j++) 653 1593 { 654 d=A[i+l, j]; 655 for(k=1;k<=nrows(A);k++) 656 { 657 A[k, j]=A[k, j]d*A[k, i]; 658 } 1594 if(A[row, j]!=0) 1595 { 1596 gcdvec = gcdcomb(A[row,column],A[row,j]); 1597 // gcdvec; 1598 // typeof(A[1..rr,column]); 1599 v1 = A[1..rr,column]; 1600 v2 = A[1..rr,j]; 1601 transform = unitMatrix(cc); 1602 transform[j,j] = v1[row]/gcdvec[1]; 1603 transform[column, column] = gcdvec[2]; 1604 transform[column,j] = v2[row]/gcdvec[1]; 1605 transform[j,column] = gcdvec[3]; 1606 q = q*transform; 1607 A = A*transform; 1608 // A; 1609 } 659 1610 } 660 else1611 if(A[row,column]<0) 661 1612 { 662 while((A[i+l, i] * A[i+l, j])!=0) 663 { 664 if(A[i+l, i]> A[i+l, j]) 665 { 666 667 for(k=1;k<=nrows(A);k++) 668 { 669 save=A[k, i]; 670 A[k, i]=A[k, j]; 671 A[k, j]=save; 672 } 673 } 674 a1=A[i+l, j]%A[i+l,i]; 675 a2=A[i+l, j]a1; 676 d=a2/A[i+l, i]; 677 for(k=1;k<=nrows(A);k++) 678 { 679 A[k, j]=A[k, j] d*A[k, i]; 680 } 681 } 1613 transform = unitMatrix(cc); 1614 transform[column,column] = 1; 1615 q = q*transform; 1616 A = A*transform; 1617 } 1618 for( j=1; j<column; j++){ 1619 if(A[row, j]!=0){ 1620 transform = unitMatrix(cc); 1621 transform[column, j] = (A[row,j]+A[row, j]%A[row, column])/A[row, column]; 1622 if(A[row,j]<0){ 1623 transform[column,j]=transform[column,j]+1;} 1624 q = q*transform; 1625 A = A*transform; 682 1626 } 683 1627 } 684 for(j=1;j<i;j++) 685 { 686 a1=A[i+l, j]%A[i+l,i]; 687 d=(A[i+l, j]a1)/A[i+l, i]; 688 for(k=1;k<=nrows(A);k++){ A[k, j]=A[k, j]d*A[k, i];} 689 } 690 } 691 c=0; 692 } 693 return( A); 1628 column++; 1629 } 1630 if(size(#) > 0){ 1631 return(list(A, q)); 1632 } 1633 return(A); 694 1634 } 695 1635 example … … 697 1637 "EXAMPLE:"; echo=2; 698 1638 699 intmat M[2][5] = 1639 intmat M[2][5] = 700 1640 1, 2, 3, 4, 0, 701 1641 0,10,20,30, 1; 702 1642 703 1643 // Hermite Normal Form of M: 704 print(hermite (M));705 706 intmat T[3][4] = 1644 print(hermiteNormalForm(M)); 1645 1646 intmat T[3][4] = 707 1647 3,3,3,3, 708 1648 2,1,3,0, … … 710 1650 711 1651 // Hermite Normal Form of T: 712 print(hermite (T));713 714 intmat A[4][5] = 1652 print(hermiteNormalForm(T)); 1653 1654 intmat A[4][5] = 715 1655 1,2,3,2,2, 716 1656 1,2,3,4,0, … … 719 1659 720 1660 // Hermite Normal Form of A: 721 print(hermite(A)); 722 } 723 724 725 /******************************************************/ 726 proc isTorsionElement(intvec mdeg) 727 "USAGE: isTorsionElement(intvec mdeg); 728 PURPOSE: For a integer vector mdeg representing the multidegree of some polynomial 729 or vector this method computes if the multidegree is contained in the torsion 730 group, i.e. if it is zero in the multigrading. 731 EXAMPLE: example isTorsionElement; shows an example 732 " 733 { 734 intmat H = getTorsion("hermite"); 735 int x, k, i; 736 737 int r = nrows(H); 738 int c = ncols(H); 739 740 int rr = nrows(mdeg); 741 742 for( i = 1; (i<=r) && (i<=c); i++) 743 { 744 if(H[i, i]!=0) 745 { 746 x = mdeg[i]%H[i, i]; 747 748 if(x!=0) 749 { 750 return(1==0); 751 } 752 753 x = mdeg[i] / H[i,i]; 754 755 for( k=1; k <= rr; k++) 756 { 757 mdeg[k] = mdeg[k]  x*H[k,i]; 758 } 759 } 760 } 761 762 return( mdeg == 0 ); 763 764 } 1661 print(hermiteNormalForm(A)); 1662 } 1663 1664 proc areZeroElements(intmat m, list #) 1665 "same as isZeroElement but for an integer matrix considered as a collection of columns" 1666 { 1667 int r = nrows(m); 1668 int i = ncols(m); 1669 1670 intvec v; 1671 1672 for( ; i > 0; i ) 1673 { 1674 v = m[1..r, i]; 1675 if( !isZeroElement(v, #) ) 1676 { 1677 return (0); 1678 } 1679 } 1680 return(1); 1681 } 1682 765 1683 example 766 1684 { … … 776 1694 1; 777 1695 1696 intmat tt[2][1]= 1697 1, 1698 1; 1699 778 1700 setBaseMultigrading(g,t); 779 1701 … … 785 1707 poly f = z2; 786 1708 787 intvec v1 = mDeg(a)  mDeg(b); 788 v1; 789 isTorsionElement(v1); 790 791 intvec v2 = mDeg(a)  mDeg(c); 792 v2; 793 isTorsionElement(v2); 794 795 intvec v3 = mDeg(e)  mDeg(f); 796 v3; 797 isTorsionElement(v3); 798 799 intvec v4 = mDeg(c)  mDeg(d); 800 v4; 801 isTorsionElement(v4); 1709 intmat m[5][2]=mDeg(a)mDeg(b),mDeg(b)mDeg(c),mDeg(c)mDeg(d),mDeg(d)mDeg(e),mDeg(e)mDeg(f); 1710 m=transpose(m); 1711 areZeroElementes(m); 1712 areZeroElementes(m,tt); 802 1713 } 803 1714 804 1715 805 1716 /******************************************************/ 806 proc defineHomogenous(poly f, list #) 807 "USAGE: defineHomogenous(f[, G]); polynomial f, integer matrix G 808 PURPOSE: Yields a matrix which has to be appended to the torsion matrix to make the 809 polynomial f homogenous in the grading by grad. 810 EXAMPLE: example defineHomogenous; shows an example 1717 proc isZeroElement(intvec mdeg, list #) 1718 "USAGE: isZeroElement(d, [T]); intvec d, group T 1719 PURPOSE: For a integer vector mdeg representing the multidegree of some polynomial 1720 or vector this method computes if the multidegree is contained in the grading group 1721 group (either set globally or given as an optional argument), i.e. if it is zero in the multigrading. 1722 EXAMPLE: example isZeroElement; shows an example 811 1723 " 812 1724 { … … 815 1727 if( typeof(#[1]) == "intmat" ) 816 1728 { 817 intmat grad = #[1]; 818 } 819 } 820 821 if( !defined(grad) ) 822 { 823 intmat grad = getVariableWeights(); 824 } 825 826 intmat newtor[nrows(grad)][size(f)1]; 827 int i,j; 828 intvec l = grad*leadexp(f); 1729 intmat H = hermiteNormalForm(#[1]); 1730 } else 1731 { 1732 if( typeof(#[1]) == "list" ) 1733 { 1734 list L = #[1]; 1735 intmat H = attrib(L, "hermite"); // todo 1736 } 1737 } 1738 1739 } 1740 if( !defined(H) ) 1741 { 1742 intmat H = getLattice("hermite"); 1743 } 1744 1745 int x, k, i, row; 1746 1747 int r = nrows(H); 1748 int c = ncols(H); 1749 1750 int rr = nrows(mdeg); 1751 row = 1; 829 1752 intvec v; 830 for(i=2; i <= size(f); i++) 831 { 832 v = grad * leadexp(f[i])  l; 833 for( j=1; j<=size(v); j++) 834 { 835 newtor[j,i1] = v[j]; 836 } 837 } 838 return(newtor); 839 } 840 example 841 { 842 "EXAMPLE:"; echo=2; 843 844 ring r =0,(x,y,z),dp; 845 intmat grad[2][3] = 846 1,0,1, 847 0,1,1; 848 849 setBaseMultigrading(grad); 850 851 poly f = x2y3z5+x3zx; 852 853 intmat M = defineHomogenous(f); 854 M; 855 defineHomogenous(f, grad) == M; 856 857 isHomogenous(f); 858 setBaseMultigrading(grad, M); 859 isHomogenous(f); 860 } 861 862 /******************************************************/ 863 proc pushForward(map f) 864 "USAGE: pushForward(f); 865 PURPOSE: Computes the finest grading of the image ring which makes the map f 866 a map of graded rings. The group map between the two grading groups is given 867 by transpose( (Id, 0) ). Pay attention that the group spanned by the columns of 868 the torsion matrix may not be a subgroup of the grading group. Still all columns 869 are needed to find the correct image of the preimage gradings. 870 EXAMPLE: example pushForward; shows an example 871 " 872 { 873 874 int k,i,j; 875 // f; 876 877 // listvar(); 878 def pre = preimage(f); 879 880 // "pre: "; pre; 881 882 intmat oldgrad=getVariableWeights(pre); 883 intmat oldtor=getTorsion(pre); 884 885 int n=nvars(pre); 886 int np=nvars(basering); 887 int p=nrows(oldgrad); 888 int pp=p+np; 889 890 intmat newgrad[pp][np]; 891 892 for(i=1;i<=np;i++){ newgrad[p+i,i]=1;} 893 894 //newgrad; 895 896 897 898 list newtor; 899 intmat toadd; 900 int columns=0; 901 902 intmat toadd1[pp][n]; 903 intvec v; 904 poly im; 905 906 for(i=1;i<=p;i++){ 907 for(j=1;j<=n;j++){ toadd1[i,j]=oldgrad[i,j];} 908 } 909 910 for(i=1;i<=n;i++){ 911 im=f[i]; 912 //im; 913 toadd = defineHomogenous(im, newgrad); 914 newtor=insert(newtor,toadd); 915 columns=columns+ncols(toadd); 916 917 v=leadexp(f[i]); 918 for(j=p+1;j<=p+np;j++){ toadd1[j,i]=v[jp];} 919 } 920 921 newtor=insert(newtor,toadd1); 922 columns=columns+ncols(toadd1); 923 924 925 if(typeof(basering)=="qring"){ 926 //"Entering qring"; 927 ideal a=ideal(basering); 928 for(i=1;i<=size(a);i++){ 929 toadd = defineHomogenous(a[i], newgrad); 930 //toadd; 931 columns=columns+ncols(toadd); 932 newtor=insert(newtor,toadd); 933 } 934 } 935 936 //newtor; 937 intmat imofoldtor[pp][ncols(oldtor)]; 938 for(i=1; i<=nrows(oldtor);i++){ 939 for(j=1; j<=ncols(oldtor); j++){ 940 imofoldtor[i,j]=oldtor[i,j]; 941 } 942 } 943 944 columns=columns+ncols(oldtor); 945 newtor=insert(newtor, imofoldtor); 946 947 intmat torsion[pp][columns]; 948 columns=0; 949 for(k=1;k<=size(newtor);k++){ 950 for(i=1;i<=pp;i++){ 951 for(j=1;j<=ncols(newtor[k]);j++){torsion[i,j+columns]=newtor[k][i,j];} 952 } 953 columns=columns+ncols(newtor[k]); 954 } 955 956 torsion=hermite(torsion); 957 intmat result[pp][pp]; 958 for(i=1;i<=pp;i++){ 959 for(j=1;j<=pp;j++){result[i,j]=torsion[i,j];} 960 } 961 962 setBaseMultigrading(newgrad, result); 963 964 } 965 example 966 { 967 "EXAMPLE:"; echo=2; 968 969 ring r = 0,(x,y,z),dp; 970 971 972 973 // Setting degrees for preimage ring.; 974 intmat grad[3][3] = 975 1,0,0, 976 0,1,0, 977 0,0,1; 978 979 setBaseMultigrading(grad); 980 981 // grading on r: 982 getVariableWeights(); 983 getTorsion(); 984 985 // only for the purpose of this example 986 if( voice > 1 ){ /*keepring(r);*/ export(r); } 987 988 ring R = 0,(a,b),dp; 989 ideal i = a2b2+a6b5+ab3,a7b+b15ab6+a6b6; 990 991 // The quotient ring by this ideal will become our image ring.; 992 qring Q = std(i); 993 994 listvar(); 995 996 map f = r,a2b6+b5+a3b+a2+ab,a2b73a2b5+b4+a,a6b6b3+a2; f; 997 998 999 // TODO: Unfortunately this is not a very spectacular example...: 1000 // Pushing forward f: 1001 pushForward(f); 1002 1003 // due to pushForward we have got new grading on Q 1004 getVariableWeights(); 1005 getTorsion(); 1006 1007 1008 // only for the purpose of this example 1009 if( voice > 1 ){ kill r; } 1010 1011 } 1012 1013 1014 /******************************************************/ 1015 proc equalMDeg(intvec exp1, intvec exp2, list #) 1016 "USAGE: equalMDeg(exp1, exp2[, V]); intvec exp1, exp2, intmat V 1017 PURPOSE: Tests if the exponent vectors of two monomials (given by exp1 and exp2) 1018 represent the same multidegree. 1019 NOTE: the integer matrix V encodes multidegrees of module components, 1020 if module component is present in exp1 and exp2 1021 EXAMPLE: example equalMDeg; shows an example 1022 " 1023 { 1024 if( size(exp1) != size(exp2) ) 1025 { 1026 ERROR("Sorry: we cannot compare exponents comming from a polynomial and a vector yet!"); 1027 } 1028 1029 if( exp1 == exp2) 1030 { 1031 return (1==1); 1032 } 1033 1034 1035 1036 intmat M = getVariableWeights(); 1037 1038 if( nrows(exp1) > ncols(M) ) // vectors => last exponent is the module component! 1039 { 1040 if( (size(#) == 0) or (typeof(#[1])!="intmat") ) 1041 { 1042 ERROR("Sorry: wrong or missing moduleunitweightsmatrix V!"); 1043 } 1044 intmat V = #[1]; 1045 1046 // typeof(V); print(V); 1047 1048 int N = ncols(M); 1049 int r = nrows(M); 1050 1051 intvec d = intvec(exp1[1..N])  intvec(exp2[1..N]); 1052 intvec dm = intvec(V[1..r, exp1[N+1]])  intvec(V[1..r, exp2[N+1]]); 1053 1054 intvec difference = M * d + dm; 1055 } 1056 else 1057 { 1058 intvec d = (exp1  exp2); 1059 intvec difference = M * d; 1060 } 1061 1062 if (isFreeRepresented()) // no torsion!? 1063 { 1064 return ( difference == 0); 1065 } 1066 return ( isTorsionElement( difference ) ); 1067 } 1068 example 1069 { 1070 "EXAMPLE:"; echo=2; 1753 for(i=1; (i<=r)&&(row<=r)&&(i<=c); i++) 1754 { 1755 while((H[row,i]==0)&&(row<=r)) 1756 { 1757 row++; 1758 if(row == (r+1)){ 1759 break; 1760 } 1761 } 1762 if(row<=r){ 1763 if(H[row,i]!=0) 1764 { 1765 v = H[1..r,i]; 1766 mdeg = mdeg(mdeg[row]mdeg[row]%v[row])/v[row]*v; 1767 } 1768 } 1769 } 1770 return( mdeg == 0 ); 1771 1772 } 1773 example 1774 { 1775 "EXAMPLE:"; echo=2; 1071 1776 1072 1777 ring r = 0,(x,y,z),dp; … … 1075 1780 1,0,1, 1076 1781 0,1,1; 1077 1078 1782 intmat t[2][1]= 1079 1783 2, 1080 1784 1; 1785 1786 intmat tt[2][1]= 1787 1, 1788 1; 1081 1789 1082 1790 setBaseMultigrading(g,t); … … 1089 1797 poly f = z2; 1090 1798 1799 intvec v1 = mDeg(a)  mDeg(b); 1800 v1; 1801 isZeroElement(v1); 1802 isZeroElement(v1, tt); 1803 1804 intvec v2 = mDeg(a)  mDeg(c); 1805 v2; 1806 isZeroElement(v2); 1807 isZeroElement(v2, tt); 1808 1809 intvec v3 = mDeg(e)  mDeg(f); 1810 v3; 1811 isZeroElement(v3); 1812 isZeroElement(v3, tt); 1813 1814 intvec v4 = mDeg(c)  mDeg(d); 1815 v4; 1816 isZeroElement(v4); 1817 isZeroElement(v4, tt); 1818 } 1819 1820 1821 /******************************************************/ 1822 proc defineHomogeneous(poly f, list #) 1823 "USAGE: defineHomogeneous(f[, G]); polynomial f, integer matrix G 1824 PURPOSE: Yields a matrix which has to be appended to the grading group matrix to make the 1825 polynomial f homogeneous in the grading by grad. 1826 EXAMPLE: example defineHomogeneous; shows an example 1827 " 1828 { 1829 if( size(#) > 0 ) 1830 { 1831 if( typeof(#[1]) == "intmat" ) 1832 { 1833 intmat grad = #[1]; 1834 } 1835 } 1836 1837 if( !defined(grad) ) 1838 { 1839 intmat grad = getVariableWeights(); 1840 } 1841 1842 intmat newgg[nrows(grad)][size(f)1]; 1843 int i,j; 1844 intvec l = grad*leadexp(f); 1845 intvec v; 1846 for(i=2; i <= size(f); i++) 1847 { 1848 v = grad * leadexp(f[i])  l; 1849 for( j=1; j<=size(v); j++) 1850 { 1851 newgg[j,i1] = v[j]; 1852 } 1853 } 1854 return(newgg); 1855 } 1856 example 1857 { 1858 "EXAMPLE:"; echo=2; 1859 1860 ring r =0,(x,y,z),dp; 1861 intmat grad[2][3] = 1862 1,0,1, 1863 0,1,1; 1864 1865 setBaseMultigrading(grad); 1866 1867 poly f = x2y3z5+x3zx; 1868 1869 intmat M = defineHomogeneous(f); 1870 M; 1871 defineHomogeneous(f, grad) == M; 1872 1873 isHomogeneous(f); 1874 setBaseMultigrading(grad, M); 1875 isHomogeneous(f); 1876 } 1877 1878 1879 proc gradiator(def h) 1880 PURPOSE: coarsens the grading of the basering until the polynom or ideal h becomes homogeneous. 1881 1882 { 1883 if(typeof(h)=="poly"){ 1884 intmat W = getVariableWeights(); 1885 intmat L = getLattice(); 1886 intmat toadd = defineHomogeneous(h); 1887 //h; 1888 //toadd; 1889 if(ncols(toadd) == 0) 1890 { 1891 return(1==1); 1892 } 1893 int rr = nrows(W); 1894 intmat newL[rr][ncols(L)+ncols(toadd)]; 1895 newL[1..rr,1..ncols(L)] = L[1..rr,1..ncols(L)]; 1896 newL[1..rr,(ncols(L)+1)..(ncols(L)+ncols(toadd))] = toadd[1..rr,1..ncols(toadd)]; 1897 setBaseMultigrading(W,newL); 1898 return(1==1); 1899 } 1900 if(typeof(h)=="ideal"){ 1901 int i; 1902 def s = (1==1); 1903 for(i=1;i<=size(h);i++){ 1904 s = s && gradiator(h[i]); 1905 } 1906 return(s); 1907 } 1908 return(1==0); 1909 } 1910 example 1911 { 1912 "EXAMPLE:"; echo=2; 1913 ring r = 0,(x,y,z),dp; 1914 intmat g[2][3] = 1,0,1,0,1,1; 1915 intmat l[2][1] = 3,0; 1916 1917 setBaseMultigrading(g,l); 1918 1919 getLattice(); 1920 1921 ideal i = y5+x4, 1922 y6+xz, 1923 x2y; 1924 gradiator(i); 1925 getLattice(); 1926 isHomogeneous(i); 1927 } 1928 1929 1930 proc pushForward(map f) 1931 "USAGE: pushForward(f); 1932 PURPOSE: Computes the finest grading of the image ring which makes the map f 1933 a map of graded rings. The group map between the two grading groups is given 1934 by transpose( (Id, 0) ). Pay attention that the group spanned by the columns of 1935 the grading group matrix may not be a subgroup of the grading group. Still all columns 1936 are needed to find the correct image of the preimage gradings. 1937 EXAMPLE: example pushForward; shows an example 1938 " 1939 { 1940 1941 int k,i,j; 1942 // f; 1943 1944 // listvar(); 1945 def pre = preimage(f); 1946 1947 // "pre: "; pre; 1948 1949 intmat oldgrad=getVariableWeights(pre); 1950 intmat oldtor=getLattice(pre); 1951 1952 int n=nvars(pre); 1953 int np=nvars(basering); 1954 int p=nrows(oldgrad); 1955 int pp=p+np; 1956 1957 intmat newgrad[pp][np]; 1958 1959 for(i=1;i<=np;i++){ newgrad[p+i,i]=1;} 1960 1961 //newgrad; 1962 1963 1964 1965 list newtor; 1966 intmat toadd; 1967 int columns=0; 1968 1969 intmat toadd1[pp][n]; 1970 intvec v; 1971 poly im; 1972 1973 for(i=1;i<=p;i++){ 1974 for(j=1;j<=n;j++){ toadd1[i,j]=oldgrad[i,j];} 1975 } 1976 1977 for(i=1;i<=n;i++){ 1978 im=f[i]; 1979 //im; 1980 toadd = defineHomogeneous(im, newgrad); 1981 newtor=insert(newtor,toadd); 1982 columns=columns+ncols(toadd); 1983 1984 v=leadexp(f[i]); 1985 for(j=p+1;j<=p+np;j++){ toadd1[j,i]=v[jp];} 1986 } 1987 1988 newtor=insert(newtor,toadd1); 1989 columns=columns+ncols(toadd1); 1990 1991 1992 if(typeof(basering)=="qring"){ 1993 //"Entering qring"; 1994 ideal a=ideal(basering); 1995 for(i=1;i<=size(a);i++){ 1996 toadd = defineHomogeneous(a[i], newgrad); 1997 //toadd; 1998 columns=columns+ncols(toadd); 1999 newtor=insert(newtor,toadd); 2000 } 2001 } 2002 2003 //newtor; 2004 intmat imofoldtor[pp][ncols(oldtor)]; 2005 for(i=1; i<=nrows(oldtor);i++){ 2006 for(j=1; j<=ncols(oldtor); j++){ 2007 imofoldtor[i,j]=oldtor[i,j]; 2008 } 2009 } 2010 2011 columns=columns+ncols(oldtor); 2012 newtor=insert(newtor, imofoldtor); 2013 2014 intmat gragr[pp][columns]; 2015 columns=0; 2016 for(k=1;k<=size(newtor);k++){ 2017 for(i=1;i<=pp;i++){ 2018 for(j=1;j<=ncols(newtor[k]);j++){gragr[i,j+columns]=newtor[k][i,j];} 2019 } 2020 columns=columns+ncols(newtor[k]); 2021 } 2022 2023 gragr=hermiteNormalForm(gragr); 2024 intmat result[pp][pp]; 2025 for(i=1;i<=pp;i++){ 2026 for(j=1;j<=pp;j++){result[i,j]=gragr[i,j];} 2027 } 2028 2029 setBaseMultigrading(newgrad, result); 2030 2031 } 2032 example 2033 { 2034 "EXAMPLE:"; echo=2; 2035 2036 ring r = 0,(x,y,z),dp; 2037 2038 2039 2040 // Setting degrees for preimage ring.; 2041 intmat grad[3][3] = 2042 1,0,0, 2043 0,1,0, 2044 0,0,1; 2045 2046 setBaseMultigrading(grad); 2047 2048 // grading on r: 2049 getVariableWeights(); 2050 getLattice(); 2051 2052 // only for the purpose of this example 2053 if( voice > 1 ){ /*keepring(r);*/ export(r); } 2054 2055 ring R = 0,(a,b),dp; 2056 ideal i = a2b2+a6b5+ab3,a7b+b15ab6+a6b6; 2057 2058 // The quotient ring by this ideal will become our image ring.; 2059 qring Q = std(i); 2060 2061 listvar(); 2062 2063 map f = r,a2b6+b5+a3b+a2+ab,a2b73a2b5+b4+a,a6b6b3+a2; f; 2064 2065 2066 // TODO: Unfortunately this is not a very spectacular example...: 2067 // Pushing forward f: 2068 pushForward(f); 2069 2070 // due to pushForward we have got new grading on Q 2071 getVariableWeights(); 2072 getLattice(); 2073 2074 2075 // only for the purpose of this example 2076 if( voice > 1 ){ kill r; } 2077 2078 } 2079 2080 2081 /******************************************************/ 2082 proc equalMDeg(intvec exp1, intvec exp2, list #) 2083 "USAGE: equalMDeg(exp1, exp2[, V]); intvec exp1, exp2, intmat V 2084 PURPOSE: Tests if the exponent vectors of two monomials (given by exp1 and exp2) 2085 represent the same multidegree. 2086 NOTE: the integer matrix V encodes multidegrees of module components, 2087 if module component is present in exp1 and exp2 2088 EXAMPLE: example equalMDeg; shows an example 2089 " 2090 { 2091 if( size(exp1) != size(exp2) ) 2092 { 2093 ERROR("Sorry: we cannot compare exponents comming from a polynomial and a vector yet!"); 2094 } 2095 2096 if( exp1 == exp2) 2097 { 2098 return (1==1); 2099 } 2100 2101 2102 2103 intmat M = getVariableWeights(); 2104 2105 if( nrows(exp1) > ncols(M) ) // vectors => last exponent is the module component! 2106 { 2107 if( (size(#) == 0) or (typeof(#[1])!="intmat") ) 2108 { 2109 ERROR("Sorry: wrong or missing moduleunitweightsmatrix V!"); 2110 } 2111 intmat V = #[1]; 2112 2113 // typeof(V); print(V); 2114 2115 int N = ncols(M); 2116 int r = nrows(M); 2117 2118 intvec d = intvec(exp1[1..N])  intvec(exp2[1..N]); 2119 intvec dm = intvec(V[1..r, exp1[N+1]])  intvec(V[1..r, exp2[N+1]]); 2120 2121 intvec difference = M * d + dm; 2122 } 2123 else 2124 { 2125 intvec d = (exp1  exp2); 2126 intvec difference = M * d; 2127 } 2128 2129 if (isFreeRepresented()) // no grading group!? 2130 { 2131 return ( difference == 0); 2132 } 2133 return ( isZeroElement( difference ) ); 2134 } 2135 example 2136 { 2137 "EXAMPLE:"; echo=2;printlevel=3; 2138 2139 ring r = 0,(x,y,z),dp; 2140 2141 intmat g[2][3]= 2142 1,0,1, 2143 0,1,1; 2144 2145 intmat t[2][1]= 2146 2, 2147 1; 2148 2149 setBaseMultigrading(g,t); 2150 2151 poly a = x10yz; 2152 poly b = x8y2z; 2153 poly c = x4z2; 2154 poly d = y5; 2155 poly e = x2y2; 2156 poly f = z2; 2157 1091 2158 1092 2159 equalMDeg(leadexp(a), leadexp(b)); … … 1116 2183 /******************************************************/ 1117 2184 static proc isFreeRepresented() 1118 "check whether the base muligrading is torsionfree (it is zero).1119 " 1120 { 1121 intmat T = get Torsion();2185 "check whether the base muligrading is free (it is zero). 2186 " 2187 { 2188 intmat T = getLattice(); 1122 2189 1123 2190 intmat Z[nrows(T)][ncols(T)]; 1124 2191 1125 return (T == Z); // no torsion!2192 return (T == Z); // no grading group! 1126 2193 } 1127 2194 1128 2195 1129 2196 /******************************************************/ 1130 proc isHomogen ous(def a, list #)1131 "USAGE: isHomogen ous(a[, f]); a polynomial/vector/ideal/module1132 RETURN: boolean, TRUE if a is (multi)homogen ous, and FALSE otherwise1133 EXAMPLE: example isHomogen ous; shows an example2197 proc isHomogeneous(def a, list #) 2198 "USAGE: isHomogeneous(a[, f]); a polynomial/vector/ideal/module 2199 RETURN: boolean, TRUE if a is (multi)homogeneous, and FALSE otherwise 2200 EXAMPLE: example isHomogeneous; shows an example 1134 2201 " 1135 2202 { … … 1137 2204 { 1138 2205 return ( size(mDegPartition(a)) <= 1 ) 1139 }1140 1141 if( typeof(a) == "vector" or typeof(a) == "module" )1142 {1143 intmat V = getModuleGrading(a);1144 2206 } 1145 2207 … … 1153 2215 for( int i = ncols(a); i > 0; i ) 1154 2216 { 1155 aa = a[i]; 1156 1157 if(defined(V)) 1158 { 1159 aa = setModuleGrading(aa, V); 1160 } 1161 1162 if(!isHomogenous(aa)) 2217 aa = getGradedGenerator(a, i); 2218 2219 if(!isHomogeneous(aa)) 1163 2220 { 1164 2221 return(0==1); … … 1171 2228 def g = groebner(a); // !!!! 1172 2229 1173 def b, aa; int j; 2230 def b, aa; int j; 1174 2231 for( int i = ncols(a); i > 0; i ) 1175 2232 { 1176 aa = a[i]; 1177 1178 if(defined(V)) 1179 { 1180 aa = setModuleGrading(aa, V); 1181 } 2233 aa = getGradedGenerator(a, i); 1182 2234 1183 2235 b = mDegPartition(aa); … … 1191 2243 } 1192 2244 return(1==1); 1193 } 2245 } 1194 2246 } 1195 2247 example … … 1200 2252 1201 2253 //Grading and Torsion matrices: 1202 intmat M[3][3] = 2254 intmat M[3][3] = 1203 2255 1,0,0, 1204 2256 0,1,0, … … 1217 2269 print(mDeg(_)); 1218 2270 1219 isHomogen ous(f); // f: is not homogenous2271 isHomogeneous(f); // f: is not homogeneous 1220 2272 1221 2273 poly g = 1xy2z3; 1222 isHomogen ous(g); // g: is homogenous2274 isHomogeneous(g); // g: is homogeneous 1223 2275 mDegPartition(g); 1224 2276 … … 1226 2278 ///////////////////////////////////////////////////////// 1227 2279 // new Torsion matrix: 1228 intmat T[3][4] = 2280 intmat T[3][4] = 1229 2281 3,3,3,3, 1230 2282 2,1,3,0, 1231 2283 1,2,0,3; 1232 2284 1233 2285 setBaseMultigrading(M,T); 1234 2286 1235 2287 f; 1236 isHomogen ous(f);2288 isHomogeneous(f); 1237 2289 mDegPartition(f); 1238 2290 1239 //  2291 //  1240 2292 g; 1241 isHomogen ous(g);2293 isHomogeneous(g); 1242 2294 mDegPartition(g); 1243 2295 … … 1246 2298 ring R = 0, (x,y,z), dp; 1247 2299 1248 intmat A[2][3] = 2300 intmat A[2][3] = 1249 2301 0,0,1, 1250 2302 3,2,1; 1251 intmat T[2][1] = 1252 1, 2303 intmat T[2][1] = 2304 1, 1253 2305 4; 1254 2306 setBaseMultigrading(A, T); 1255 2307 1256 isHomogen ous(ideal(x2  y3 xy +z, x*yz, x^3  y^2*z + x^2 y^3)); // 11257 isHomogen ous(ideal(x2  y3 xy +z, x*yz, x^3  y^2*z + x^2 y^3), "checkGens");1258 isHomogen ous(ideal(x+y, x2  y2)); // 02308 isHomogeneous(ideal(x2  y3 xy +z, x*yz, x^3  y^2*z + x^2 y^3)); // 1 2309 isHomogeneous(ideal(x2  y3 xy +z, x*yz, x^3  y^2*z + x^2 y^3), "checkGens"); 2310 isHomogeneous(ideal(x+y, x2  y2)); // 0 1259 2311 1260 2312 // Degree partition: … … 1262 2314 mDegPartition(x3 y2z + x2 y3 + z + 1); 1263 2315 1264 2316 1265 2317 module N = gen(1) + (x+y) * gen(2), z*gen(3); 1266 2318 1267 2319 intmat V[2][3] = 0; // 1, 2, 3, 4, 5, 6; // columnwise weights of components!!?? 1268 2320 1269 2321 vector v1, v2; 1270 2322 1271 2323 v1 = setModuleGrading(N[1], V); v1; 1272 2324 mDegPartition(v1); … … 1278 2330 1279 2331 N = setModuleGrading(N, V); 1280 isHomogen ous(N);2332 isHomogeneous(N); 1281 2333 print( mDeg(N) ); 1282 2334 1283 /////////////////////////////////////// 1284 1285 V = 1286 1, 2, 3, 2335 /////////////////////////////////////// 2336 2337 V = 2338 1, 2, 3, 1287 2339 4, 5, 6; 1288 2340 … … 1296 2348 1297 2349 N = setModuleGrading(N, V); 1298 isHomogen ous(N);2350 isHomogeneous(N); 1299 2351 print( mDeg(N) ); 1300 2352 1301 /////////////////////////////////////// 1302 1303 V = 1304 0, 0, 0, 2353 /////////////////////////////////////// 2354 2355 V = 2356 0, 0, 0, 1305 2357 4, 1, 0; 1306 2358 1307 2359 N = gen(1) + x * gen(2), z*gen(3); 1308 2360 N = setModuleGrading(N, V); print(N); 1309 isHomogen ous(N);2361 isHomogeneous(N); 1310 2362 print( mDeg(N) ); 1311 v1 = setModuleGrading(N[1], V); print(v1);2363 v1 = getGradedGenerator(N,1); print(v1); 1312 2364 mDegPartition(v1); 1313 2365 print( mDeg(_) ); 1314 2366 N = setModuleGrading(N, V); print(N); 1315 isHomogen ous(N);2367 isHomogeneous(N); 1316 2368 print( mDeg(N) ); 1317 2369 } … … 1339 2391 int r = nrows(M); 1340 2392 2393 if( (typeof(A) == "vector") or (typeof(A) == "module") ) 2394 { 2395 intmat V = getModuleGrading(A); 2396 2397 if( nrows(V) != r ) 2398 { 2399 ERROR("Sorry wrong mgradsize of V: " + string(nrows(V))); 2400 } 2401 } 2402 1341 2403 if( A == 0 ) 1342 2404 { … … 1345 2407 } 1346 2408 1347 if( (typeof(A) == "vector") or (typeof(A) == "module") )1348 {1349 intmat V = getModuleGrading(A);1350 1351 if( nrows(V) != r )1352 {1353 ERROR("Sorry wrong mgradsize of V: " + string(nrows(V)));1354 }1355 }1356 1357 2409 intvec m; m[r] = 0; 1358 2410 … … 1367 2419 A = A  lead(A); 1368 2420 while( size(A) > 0 ) 1369 { 2421 { 1370 2422 v = leadexp(A); // v; 1371 2423 m = max( m, M * v, r ); // ???? … … 1388 2440 A = A  lead(A); 1389 2441 while( size(A) > 0 ) 1390 { 2442 { 1391 2443 v = leadexp(A); // v; 1392 2444 … … 1413 2465 { 1414 2466 G[j, i] = d[j]; 1415 } 2467 } 1416 2468 } 1417 2469 return(G); … … 1425 2477 for( i = ncols(A); i > 0; i ) 1426 2478 { 1427 v = setModuleGrading(A[i], V);1428 1429 // G[1..r, i] 2479 v = getGradedGenerator(A, i); 2480 2481 // G[1..r, i] 1430 2482 d = mDeg(v); 1431 2483 … … 1433 2485 { 1434 2486 G[j, i] = d[j]; 1435 } 2487 } 1436 2488 1437 2489 } … … 1439 2491 return(G); 1440 2492 } 1441 2493 1442 2494 } 1443 2495 example … … 1453 2505 print(Ta); 1454 2506 1455 // attrib(A, " torsion", Ta); // to think about2507 // attrib(A, "gradingGroup", Ta); // to think about 1456 2508 1457 2509 // "poly:"; … … 1463 2515 1464 2516 setBaseMultigrading(A, Ta); 1465 2517 1466 2518 mDeg( x*x*y ); 1467 2519 1468 2520 mDeg( y*y*y*x ); 1469 2521 1470 2522 mDeg( x*y + x + 1 ); 1471 2523 … … 1477 2529 1478 2530 // "ideal:"; 1479 2531 1480 2532 ideal I = y*x*x, x*y*y*y; 1481 2533 print( mDeg(I) ); … … 1485 2537 1486 2538 // "vectors:"; 1487 2539 1488 2540 intmat B[2][2] = 0, 1, 1, 0; 1489 2541 print(B); 1490 2542 1491 2543 mDeg( setModuleGrading(y*y*y*gen(2), B )); 1492 2544 mDeg( setModuleGrading(x*x*gen(1), B )); 1493 2545 1494 2546 1495 2547 vector V = x*gen(1) + y*gen(2); 1496 2548 V = setModuleGrading(V, B); … … 1499 2551 vector v1 = setModuleGrading([0, 0, 0], B); 1500 2552 print( mDeg( v1 ) ); 1501 2553 1502 2554 vector v2 = setModuleGrading([0], B); 1503 2555 print( mDeg( v2 ) ); 1504 2556 1505 2557 // "module:"; 1506 2558 1507 2559 module D = x*gen(1), y*gen(2); 1508 2560 D; 1509 2561 D = setModuleGrading(D, B); 1510 2562 print( mDeg( D ) ); 1511 2563 1512 2564 1513 2565 module DD = [0, 0],[0, 0, 0]; … … 1518 2570 DDD = setModuleGrading(DDD, B); 1519 2571 print( mDeg( DDD ) ); 1520 } 2572 2573 }; 2574 2575 2576 2577 1521 2578 1522 2579 /******************************************************/ … … 1526 2583 EXAMPLE: example mDegPartition; shows an example 1527 2584 " 1528 { 1529 if( typeof(p) == "poly" ) 1530 { 1531 ideal I; 2585 { // TODO: What about an ideal or module??? 2586 2587 if( typeof(p) == "poly" ) 2588 { 2589 ideal I; 1532 2590 poly mp, t, tt; 1533 2591 } … … 1536 2594 if( typeof(p) == "vector" ) 1537 2595 { 1538 module I; 2596 module I; 1539 2597 vector mp, t, tt; 1540 2598 } … … 1547 2605 if( typeof(p) == "vector" ) 1548 2606 { 1549 intmat V = getModuleGrading(p); 2607 intmat V = getModuleGrading(p); 1550 2608 } 1551 2609 else … … 1554 2612 } 1555 2613 1556 if( size(p) > 1) 2614 if( size(p) > 1) 1557 2615 { 1558 2616 intvec m; … … 1561 2619 { 1562 2620 m = leadexp(p); 1563 mp = lead(p); 2621 mp = lead(p); 1564 2622 p = p  lead(p); 1565 2623 tt = p; t = 0; 1566 2624 1567 2625 while( size(tt) > 0 ) 1568 { 2626 { 1569 2627 // TODO: we do not cache matrices (M,T,H,V), which remain the same :( 1570 2628 // TODO: we need some lowlevel procedure with all these arguments...! 1571 if( equalMDeg( leadexp(tt), m, V ) ) 2629 if( equalMDeg( leadexp(tt), m, V ) ) 1572 2630 { 1573 2631 mp = mp + lead(tt); // "mp", mp; … … 1616 2674 1617 2675 mDegPartition(f); 1618 2676 1619 2677 vector v = xy*gen(1)x3y2*gen(2)+x4y*gen(3); 1620 2678 intmat B[2][3]=1,1,2,0,0,1; 1621 2679 v = setModuleGrading(v,B); 1622 2680 getModuleGrading(v); 1623 2681 1624 2682 mDegPartition(v, B); 1625 2683 } … … 1631 2689 { 1632 2690 intmat A[n][n]; 1633 2691 1634 2692 for( int i = n; i > 0; i ) 1635 2693 { … … 1646 2704 " 1647 2705 USAGE: finestMDeg(r); ring r 1648 RETURN: ring, r endowed with the finest multigrading 2706 RETURN: ring, r endowed with the finest multigrading 1649 2707 TODO: not yet... 1650 2708 " … … 1675 2733 1676 2734 1677 if( n > 0) 1678 { 1679 1680 intmat L[N][n]; 2735 if( n > 0) 2736 { 2737 2738 intmat L[N][n]; 1681 2739 // list L; 1682 2740 int j = n; … … 1686 2744 p = I[i]; 1687 2745 1688 if( size(p) > 1 ) 2746 if( size(p) > 1 ) 1689 2747 { 1690 2748 intvec m0 = leadexp(p); … … 1701 2759 1702 2760 print(L); 1703 setBaseMultigrading(A, L); 1704 } 2761 setBaseMultigrading(A, L); 2762 } 1705 2763 else 1706 2764 { … … 1739 2797 1740 2798 if( size(#) > 0 and typeof(#[1]) == "intmat" ) 1741 { 2799 { 1742 2800 attrib(F, "P", #[1]); 1743 2801 } … … 1797 2855 // check: 1798 2856 print(L*M); 1799 } 2857 }; 2858 2859 1800 2860 1801 2861 /******************************************************/ … … 1819 2879 " 1820 2880 { 2881 if( system("sh","which hilbert 2> /dev/null 1> /dev/null") != 0 ) 2882 { 2883 ERROR("Sorry: cannot find 'hilbert' command from 4ti2. Please install 4ti2!"); 2884 } 2885 1821 2886 // 1822 2887 // Initialization and Sanity Checks … … 1870 2935 // using standard unix commands 1871 2936 // 2937 2938 1872 2939 j=system("sh","hilbert q n sing4ti2"); ////////// be quiet + no loggin!!! 1873 2940 1874 j=system("sh", "awk \'BEGIN{ORS=\",\";}{print $0;}\' sing4ti2.hil " + 1875 " sed s/[\\\ \\\t\\\v\\\f]/,/g " + 1876 " sed s/,+/,/gsed s/,,/,/g " + 1877 " sed s/,,/,/g " + 2941 j=system("sh", "awk \'BEGIN{ORS=\",\";}{print $0;}\' sing4ti2.hil " + 2942 " sed s/[\\\ \\\t\\\v\\\f]/,/g " + 2943 " sed s/,+/,/gsed s/,,/,/g " + 2944 " sed s/,,/,/g " + 1878 2945 "> sing4ti2.converted" ); 1879 if( defined(keepfiles) <= 0) 1880 { 1881 j=system("sh",("rm f sing4ti2.hil sing4ti2."+fileending)); 1882 } 2946 2947 1883 2948 // 1884 2949 // reading output of 4ti2 … … 1889 2954 1890 2955 string s = read(ausg); 2956 2957 if( defined(keepfiles) <= 0) 2958 { 2959 j=system("sh",("rm f sing4ti2.hil sing4ti2.converted sing4ti2."+fileending)); 2960 } 2961 1891 2962 string ergstr = "intvec erglist = " + s + "0;"; 1892 2963 execute(ergstr); 1893 2964 1894 2965 // print(erglist); 1895 2966 1896 2967 int Rnc = erglist[1]; 1897 2968 int Rnr = erglist[2]; 1898 2969 1899 2970 intmat R[Rnr][Rnc]; 1900 2971 … … 1910 2981 } 1911 2982 } 2983 2984 1912 2985 1913 2986 return (R); … … 1931 3004 1, 1, 0, 0, 1, 0,1, 0, 0; 1932 3005 hilbert4ti2intmat(M); 1933 hermite (M);3006 hermiteNormalForm(M); 1934 3007 } 1935 3008 … … 1962 3035 1963 3036 /******************************************************/ 1964 proc mDegBasis(intvec d) 3037 proc mDegBasis(intvec d) 1965 3038 " 1966 3039 USAGE: multidegree d … … 1991 3064 } 1992 3065 1993 intmat T = get Torsion(R);3066 intmat T = getLattice(R); 1994 3067 1995 3068 if( isFreeRepresented() ) … … 2006 3079 2007 3080 intmat AA[nr][nc + 2 * n]; 2008 AA[1..nr, 1.. nc] = A[1..nr, 1.. nc]; 2009 AA[1..nr, nc + (1.. n)] = T[1..nr, 1.. n]; 2010 AA[1..nr, nc + n + (1.. n)] = T[1..nr, 1.. n]; 3081 AA[1..nr, 1.. nc] = A[1..nr, 1.. nc]; 3082 AA[1..nr, nc + (1.. n)] = T[1..nr, 1.. n]; 3083 AA[1..nr, nc + n + (1.. n)] = T[1..nr, 1.. n]; 2011 3084 2012 3085 2013 3086 // print ( AA ); 2014 3087 2015 intmat K = leftKernelZ(( AA ) ); // 3088 intmat K = leftKernelZ(( AA ) ); // 2016 3089 2017 3090 // print(K); … … 2022 3095 // "!"; 2023 3096 2024 intmat B = hilbert4ti2intmat(transpose(KK), 1); 3097 intmat B = hilbert4ti2intmat(transpose(KK), 1); 2025 3098 2026 3099 // "!"; print(B); … … 2033 3106 2034 3107 2035 int i; 3108 int i; 2036 3109 int nnr = nrows(B); 2037 3110 int nnc = ncols(B); … … 2088 3161 2089 3162 intmat g1[2][2]=1,0,0,1; 2090 intmat t[2][1]=2,0;3163 intmat l[2][1]=2,0; 2091 3164 intmat g2[2][2]=1,1,1,1; 2092 3165 intvec v1=4,0; 2093 3166 intvec v2=4,4; 2094 3167 2095 3168 intmat g3[1][2]=1,1; 2096 3169 setBaseMultigrading(g3); … … 2098 3171 v3; 2099 3172 mDegBasis(v3); 2100 2101 setBaseMultigrading(g1, t);3173 3174 setBaseMultigrading(g1,l); 2102 3175 mDegBasis(v1); 2103 3176 setBaseMultigrading(g2); 2104 3177 mDegBasis(v2); 2105 3178 2106 3179 intmat M[2][2] = 1, 1, 1, 1; 2107 3180 intvec d = 2, 2; … … 2117 3190 intmat M[2][3] = 1, 2, 1, 1, 1, 0; 2118 3191 2119 intmat T[2][1] = 0, 2;3192 intmat L[2][1] = 0, 2; 2120 3193 2121 3194 intvec d = 4, 1; 2122 3195 2123 setBaseMultigrading(M, T);3196 setBaseMultigrading(M, L); 2124 3197 2125 3198 mDegBasis(d); … … 2162 3235 2163 3236 proc mDegSyzygy(def I) 2164 "USAGE: mDegSyzygy(I); I is a poly/vector/ideal/module3237 "USAGE: mDegSyzygy(I); I is a ideal or a module 2165 3238 PURPOSE: computes the multigraded syzygy module of I 2166 3239 RETURNS: module, the syzygy module of I … … 2169 3242 " 2170 3243 { 2171 if( isHomogen ous(I, "checkGens") == 0)2172 { 2173 ERROR ("Sorry: inhomogen ous input!");2174 } 3244 if( isHomogeneous(I, "checkGens") == 0) 3245 { 3246 ERROR ("Sorry: inhomogeneous input!"); 3247 } 2175 3248 module S = syz(I); 2176 3249 S = setModuleGrading(S, mDeg(I)); … … 2180 3253 { 2181 3254 "EXAMPLE:"; echo=2; 2182 3255 2183 3256 2184 3257 ring r = 0,(x,y,z,w),dp; … … 2188 3261 setBaseMultigrading(MM); 2189 3262 module M = ideal( xwyz, x2zy3, xz2y2w, yw2z3); 2190 2191 3263 3264 2192 3265 intmat v[2][nrows(M)]= 2193 3266 1, 2194 3267 0; 2195 3268 2196 3269 M = setModuleGrading(M, v); 2197 3270 2198 isHomogen ous(M);3271 isHomogeneous(M); 2199 3272 "Multidegrees: "; print(mDeg(M)); 2200 3273 // Let's compute syzygies! … … 2203 3276 "Multidegrees: "; print(mDeg(S)); 2204 3277 2205 isHomogenous(S); 3278 isHomogeneous(S); 3279 } 3280 3281 3282 3283 proc mDegModulo(def I, def J) 3284 "USAGE: mDegModulo(I); I, J are ideals or modules 3285 PURPOSE: computes the multigraded 'modulo' module of I and J 3286 RETURNS: module, see 'modulo' command 3287 NOTE: I and J should have the same multigrading, and their 3288 generators must be multigraded homogeneous 3289 EXAMPLE: example mDegModulo; shows an example 3290 " 3291 { 3292 if( (isHomogeneous(I, "checkGens") == 0) or (isHomogeneous(J, "checkGens") == 0) ) 3293 { 3294 ERROR ("Sorry: inhomogeneous input!"); 3295 } 3296 module K = modulo(I, J); 3297 K = setModuleGrading(K, mDeg(I)); 3298 return (K); 3299 } 3300 example 3301 { 3302 "EXAMPLE:"; echo=2; 3303 3304 ring r = 0,(x,y,z),dp; 3305 intmat MM[2][3]= 3306 1,1,1, 3307 0,1,3; 3308 setBaseMultigrading(MM); 3309 3310 ideal h1 = x, y, z; 3311 ideal h2 = x; 3312 3313 "Multidegrees: "; print(mDeg(h1)); 3314 3315 // Let's compute modulo(h1, h2): 3316 def K = mDegModulo(h1, h2); K; 3317 3318 "Module Units Multigrading: "; print( getModuleGrading(K) ); 3319 "Multidegrees: "; print(mDeg(K)); 3320 3321 isHomogeneous(K); 2206 3322 } 2207 3323 … … 2210 3326 "USAGE: mDegGroebner(I); I is a poly/vector/ideal/module 2211 3327 PURPOSE: computes the multigraded standard/groebner basis of I 2212 NOTE: I must be multigraded homogeneous 3328 NOTE: I must be multigraded homogeneous 2213 3329 RETURNS: ideal/module, the computed basis 2214 3330 EXAMPLE: example mDegGroebner; shows an example 2215 3331 " 2216 3332 { 2217 if( isHomogen ous(I) == 0)2218 { 2219 ERROR ("Sorry: inhomogen ous input!");2220 } 3333 if( isHomogeneous(I) == 0) 3334 { 3335 ERROR ("Sorry: inhomogeneous input!"); 3336 } 2221 3337 2222 3338 def S = groebner(I); 2223 3339 2224 3340 if( typeof(I) == "module" or typeof(I) == "vector" ) 2225 3341 { 2226 S = setModuleGrading(S, getModuleGrading(I)); 3342 S = setModuleGrading(S, getModuleGrading(I)); 2227 3343 } 2228 3344 … … 2241 3357 setBaseMultigrading(MM); 2242 3358 2243 3359 2244 3360 module M = ideal( xwyz, x2zy3, xz2y2w, yw2z3); 2245 2246 3361 3362 2247 3363 intmat v[2][nrows(M)]= 2248 3364 1, 2249 3365 0; 2250 3366 2251 3367 M = setModuleGrading(M, v); 2252 3368 … … 2258 3374 "Multidegrees: "; print(mDeg(M)); 2259 3375 2260 isHomogen ous(M);3376 isHomogeneous(M); 2261 3377 2262 3378 ///////////////////////////////////////////////////////////////////////////// … … 2266 3382 "Multidegrees: "; print(mDeg(S)); 2267 3383 2268 isHomogen ous(S);3384 isHomogeneous(S); 2269 3385 2270 3386 ///////////////////////////////////////////////////////////////////////////// … … 2274 3390 "Multidegrees: "; print(mDeg(S)); 2275 3391 2276 isHomogen ous(S);3392 isHomogeneous(S); 2277 3393 } 2278 3394 … … 2281 3397 proc mDegResolution(def I, int ll, list #) 2282 3398 "USAGE: mDegResolution(I,l,[f]); I is poly/vector/ideal/module; l,f are integers 2283 PURPOSE: computes the multigraded resolution of I of the length l, 2284 or the whole resolution if l is zero. Returns minimal resolution if an optional 3399 PURPOSE: computes the multigraded resolution of I of the length l, 3400 or the whole resolution if l is zero. Returns minimal resolution if an optional 2285 3401 argument 1 is supplied 2286 3402 NOTE: input must have multigradedhomogeneous generators. 2287 The returned list is truncated beginning with the first zero differential. 3403 The returned list is truncated beginning with the first zero differential. 2288 3404 RETURNS: list, the computed resolution 2289 3405 EXAMPLE: example mDegResolution; shows an example 2290 3406 " 2291 3407 { 2292 if( isHomogen ous(I, "checkGens") == 0)2293 { 2294 ERROR ("Sorry: inhomogen ous input!");2295 } 3408 if( isHomogeneous(I, "checkGens") == 0) 3409 { 3410 ERROR ("Sorry: inhomogeneous input!"); 3411 } 2296 3412 2297 3413 def R = res(I, ll, #); list L = R; int l = size(L); 2298 3414 def V = getModuleGrading(I); 2299 3415 if( (typeof(I) == "module") or (typeof(I) == "vector") ) 2300 3416 { 2301 L[1] = setModuleGrading(L[1], getModuleGrading(I));2302 } 2303 2304 int i; 3417 L[1] = setModuleGrading(L[1], V); 3418 } 3419 3420 int i; 2305 3421 for( i = 2; i <= l; i++ ) 2306 3422 { … … 2313 3429 } 2314 3430 } 2315 3431 2316 3432 return (L); 2317 3433 2318 3434 2319 3435 } 2320 3436 example 2321 3437 { 2322 3438 "EXAMPLE:"; echo=2; 2323 3439 2324 3440 ring r = 0,(x,y,z,w),dp; 2325 3441 … … 2330 3446 setBaseMultigrading(M); 2331 3447 2332 3448 2333 3449 module m= ideal( xwyz, x2zy3, xz2y2w, yw2z3); 2334 2335 isHomogen ous(ideal( xwyz, x2zy3, xz2y2w, yw2z3), "checkGens");2336 3450 3451 isHomogeneous(ideal( xwyz, x2zy3, xz2y2w, yw2z3), "checkGens"); 3452 2337 3453 ideal A = xwyz, x2zy3, xz2y2w, yw2z3; 2338 3454 2339 3455 int j; 2340 3456 2341 3457 for(j=1; j<=ncols(A); j++) 2342 3458 { 2343 3459 mDegPartition(A[j]); 2344 3460 } 2345 3461 2346 3462 intmat v[2][1]= 2347 3463 1, 2348 3464 0; 2349 3465 2350 3466 m = setModuleGrading(m, v); 2351 3467 … … 2374 3490 2375 3491 ///////////////////////////////////////////////////////////////////////////// 2376 3492 2377 3493 L = mDegResolution(maxideal(1), 0, 1); 2378 3494 … … 2384 3500 "Multigrading: "; print(mDeg(L[j])); 2385 3501 } 2386 3502 2387 3503 kill v; 2388 3504 2389 3505 2390 3506 def h = hilbertSeries(m); … … 2393 3509 numerator1; 2394 3510 factorize(numerator1); 2395 3511 2396 3512 denominator1; 2397 3513 factorize(denominator1); … … 2407 3523 proc hilbertSeries(def I) 2408 3524 "USAGE: hilbertSeries(I); I is poly/vector/ideal/module 2409 PURPOSE: computes the multigraded Hilbert Series of M 2410 NOTE: input must have multigradedhomogeneous generators. 3525 PURPOSE: computes the multigraded Hilbert Series of M 3526 NOTE: input must have multigradedhomogeneous generators. 2411 3527 Multigrading should be positive. 2412 RETURNS: a ring in variables t_(i), s_(i), with polynomials 2413 numerator1 and denominator1 and muturally prime numerator2 3528 RETURNS: a ring in variables t_(i), s_(i), with polynomials 3529 numerator1 and denominator1 and muturally prime numerator2 2414 3530 and denominator2, quotients of which give the series. 2415 3531 EXAMPLE: example hilbertSeries; shows an example 2416 3532 " 2417 3533 { 2418 3534 2419 3535 if( !isFreeRepresented() ) 2420 3536 { 2421 ERROR("SORRY: ONLY TORSIONFREE CASE (POSITIVE GRADING)"); 2422 } 2423 3537 "Things might happen, since we are not free."; 3538 //ERROR("SORRY: ONLY TORSIONFREE CASE (POSITIVE GRADING)"); 3539 } 3540 2424 3541 int i, j, k, v; 2425 3542 2426 3543 intmat M = getVariableWeights(); 2427 3544 2428 3545 int cc = ncols(M); 2429 3546 int n = nrows(M); … … 2437 3554 2438 3555 int l = size(RES); 2439 3556 2440 3557 list L; L[l + 1] = 0; 2441 3558 … … 2444 3561 intmat zeros[n][1]; 2445 3562 L[1] = zeros; 2446 } 3563 } 2447 3564 else 2448 3565 { … … 2454 3571 L[j + 1] = mDeg(RES[j]); 2455 3572 } 2456 3573 2457 3574 l++; 2458 3575 2459 3576 ring R = 0,(t_(1..n),s_(1..n)),dp; 2460 2461 ideal units; 3577 3578 ideal units; 2462 3579 for( i=n; i>=1; i) 2463 3580 { 2464 3581 units[i] = (var(i) * var(n + i)  1); 2465 3582 } 2466 3583 2467 3584 qring Q = std(units); 2468 3585 2469 3586 // TODO: should not it be a quotient ring depending on Torsion??? 2470 3587 // I am not sure about what to do in the torsion case, but since … … 2475 3592 poly monom, summand, numerator; 2476 3593 poly denominator = 1; 2477 3594 2478 3595 for( i = 1; i <= cc; i++) 2479 3596 { 2480 3597 monom = 1; 2481 2482 3598 for( k = 1; k <= n; k++) 2483 3599 { … … 2487 3603 { 2488 3604 monom = monom * (var(k)^(v)); 2489 } 3605 } 2490 3606 else 2491 3607 { … … 2493 3609 } 2494 3610 } 2495 3611 2496 3612 if( monom == 1) 2497 3613 { … … 2501 3617 denominator = denominator * (1  monom); 2502 3618 } 2503 2504 for( j = 1; j<= l; j++) 3619 3620 for( j = 1; j<= l; j++) 2505 3621 { 2506 3622 summand = 0; … … 2516 3632 { 2517 3633 monom = monom * (var(k)^v); 2518 } 3634 } 2519 3635 else 2520 3636 { … … 2526 3642 numerator = numerator  (1)^j * summand; 2527 3643 } 2528 3644 2529 3645 if( denominator == 0 ) 2530 3646 { 2531 3647 ERROR("Multigrading not positive."); 2532 } 2533 3648 } 3649 2534 3650 poly denominator1 = denominator; 2535 3651 poly numerator1 = numerator; … … 2565 3681 "The s_(i)variables are defined to be the inverse of the t_(i)variables."; 2566 3682 "  "; 2567 3683 2568 3684 return(Q); 2569 3685 } … … 2571 3687 { 2572 3688 "EXAMPLE:"; echo=2; 2573 3689 2574 3690 ring r = 0,(x,y,z,w),dp; 2575 3691 intmat g[2][4]= … … 2577 3693 0,1,3,4; 2578 3694 setBaseMultigrading(g); 2579 3695 2580 3696 module M = ideal(xwyz, x2zy3, xz2y2w, yw2z3); 2581 3697 intmat V[2][1]= … … 2589 3705 factorize(numerator2); 2590 3706 factorize(denominator2); 2591 3707 2592 3708 kill g, h; setring r; 2593 3709 … … 2595 3711 1,2,3,4, 2596 3712 0,0,5,8; 2597 3713 2598 3714 setBaseMultigrading(g); 2599 3715 2600 3716 ideal I = x^2, y, z^3; 2601 3717 I = std(I); … … 2612 3728 mDeg(I); 2613 3729 def h = hilbertSeries(I); setring h; 2614 3730 2615 3731 factorize(numerator2); 2616 3732 factorize(denominator2); … … 2619 3735 //////////////////////////////////////////////// 2620 3736 ring R = 0,(x,y,z),dp; 2621 intmat W[2][3] = 3737 intmat W[2][3] = 2622 3738 1,1, 1, 2623 3739 0,0,1; 2624 3740 setBaseMultigrading(W); 2625 3741 ideal I = x3y,yz2,y2z,z4; 2626 3742 2627 3743 def h = hilbertSeries(I); setring h; 2628 3744 2629 3745 factorize(numerator2); 2630 3746 factorize(denominator2); … … 2633 3749 //////////////////////////////////////////////// 2634 3750 ring R = 0,(x,y,z,a,b,c),dp; 2635 intmat W[2][6] = 3751 intmat W[2][6] = 2636 3752 1,1, 1,1,1,1, 2637 3753 0,0,1,0,0,0; 2638 3754 setBaseMultigrading(W); 2639 3755 ideal I = x3y,yz2,y2z,z4; 2640 3756 2641 3757 def h = hilbertSeries(I); setring h; 2642 3758 2643 3759 factorize(numerator2); 2644 3760 factorize(denominator2); 2645 3761 2646 3762 kill R, W, h; 2647 3763 //////////////////////////////////////////////// 2648 3764 // This is example 5.3.9. from Robbianos book. 2649 3765 2650 3766 ring R = 0,(x,y,z,w),dp; 2651 intmat W[1][4] = 3767 intmat W[1][4] = 2652 3768 1,1, 1,1; 2653 3769 setBaseMultigrading(W); … … 2655 3771 2656 3772 hilb(std(I)); 2657 3773 2658 3774 def h = hilbertSeries(I); setring h; 2659 3775 2660 3776 numerator1; 2661 3777 denominator1; … … 2663 3779 factorize(numerator2); 2664 3780 factorize(denominator2); 2665 3781 2666 3782 2667 3783 kill h; … … 2672 3788 2673 3789 hilb(std(I2)); 2674 3790 2675 3791 def h = hilbertSeries(I2); setring h; 2676 3792 … … 2682 3798 //////////////////////////////////////////////// 2683 3799 setring R; 2684 3800 2685 3801 W = 2,2,2,2; 2686 3802 2687 3803 setBaseMultigrading(W); 2688 3804 … … 2694 3810 2695 3811 kill w; 2696 3812 2697 3813 2698 3814 def h = hilbertSeries(I2); setring h; 2699 3815 2700 3816 2701 3817 numerator1; denominator1; 2702 3818 kill h; 2703 3819 2704 3820 2705 3821 kill R, W; 2706 3822 … … 2716 3832 2717 3833 hilb(std(I)); 2718 3834 2719 3835 def h = hilbertSeries(I); setring h; 2720 3836 … … 2732 3848 2733 3849 numerator1; denominator1; 3850 3851 kill h; 3852 //////////////////////////////////////////////// 3853 setring R; 3854 3855 I = x^5; I; 3856 3857 hilb(std(I)); 3858 hilb(std(I), 1); 3859 3860 def h = hilbertSeries(I); setring h; 3861 3862 numerator1; denominator1; 3863 3864 3865 kill h; 3866 //////////////////////////////////////////////// 3867 setring R; 3868 3869 I = x^10; I; 3870 3871 hilb(std(I)); 3872 3873 def h = hilbertSeries(I); setring h; 3874 3875 numerator1; denominator1; 2734 3876 2735 3877 kill h; … … 2737 3879 setring R; 2738 3880 2739 I = x^5; I; 2740 2741 hilb(std(I)); 2742 hilb(std(I), 1); 2743 2744 def h = hilbertSeries(I); setring h; 3881 module M = 1; 3882 3883 M = setModuleGrading(M, W); 3884 3885 3886 hilb(std(M)); 3887 3888 def h = hilbertSeries(M); setring h; 2745 3889 2746 3890 numerator1; denominator1; 2747 2748 3891 2749 3892 kill h; … … 2751 3894 setring R; 2752 3895 2753 I = x^10; I; 2754 2755 hilb(std(I)); 2756 2757 def h = hilbertSeries(I); setring h; 3896 kill M; module M = x^5*gen(1); 3897 // intmat V[1][3] = 0; // TODO: this would lead to a wrong result!!!? 3898 intmat V[1][1] = 0; // all gen(i) of degree 0! 3899 3900 M = setModuleGrading(M, V); 3901 3902 hilb(std(M)); 3903 3904 def h = hilbertSeries(M); setring h; 2758 3905 2759 3906 numerator1; denominator1; … … 2761 3908 kill h; 2762 3909 //////////////////////////////////////////////// 2763 setring R; 2764 2765 module M = 1; 2766 2767 M = setModuleGrading(M, W); 2768 2769 2770 hilb(std(M)); 2771 2772 def h = hilbertSeries(M); setring h; 3910 setring R; 3911 3912 module N = x^5*gen(3); 3913 3914 kill V; 3915 3916 intmat V[1][3] = 0; // all gen(i) of degree 0! 3917 3918 N = setModuleGrading(N, V); 3919 3920 hilb(std(N)); 3921 3922 def h = hilbertSeries(N); setring h; 2773 3923 2774 3924 numerator1; denominator1; … … 2776 3926 kill h; 2777 3927 //////////////////////////////////////////////// 2778 setring R; 2779 2780 kill M; module M = x^5*gen(1); 2781 // intmat V[1][3] = 0; // TODO: this would lead to a wrong result!!!? 2782 intmat V[1][1] = 0; // all gen(i) of degree 0! 2783 2784 M = setModuleGrading(M, V); 2785 2786 hilb(std(M)); 2787 2788 def h = hilbertSeries(M); setring h; 2789 2790 numerator1; denominator1; 2791 2792 kill h; 2793 //////////////////////////////////////////////// 2794 setring R; 2795 2796 module N = x^5*gen(3); 2797 2798 kill V; 2799 2800 intmat V[1][3] = 0; // all gen(i) of degree 0! 2801 2802 N = setModuleGrading(N, V); 2803 2804 hilb(std(N)); 2805 2806 def h = hilbertSeries(N); setring h; 2807 2808 numerator1; denominator1; 2809 2810 kill h; 2811 //////////////////////////////////////////////// 2812 setring R; 2813 2814 3928 setring R; 3929 3930 2815 3931 module S = M + N; 2816 3932 2817 3933 S = setModuleGrading(S, V); 2818 3934 … … 2830 3946 } 2831 3947 2832 3948 proc evalHilbertSeries(def h, intvec v) 3949 " 3950 evaluate hilbert series h by substibuting v[i] for t_(i) (1/v[i] for s_(i)) 3951 return: int (h(v)) 3952 " 3953 { 3954 if( 2*size(v) != nvars(h) ) 3955 { 3956 ERROR("Wrong input/size!"); 3957 } 3958 3959 setring h; 3960 3961 if( defined(numerator2) and defined(denominator2) ) 3962 { 3963 poly n = numerator2; poly d = denominator2; 3964 } else 3965 { 3966 poly n = numerator1; poly d = denominator1; 3967 } 3968 3969 int N = size(v); 3970 int i; number k; 3971 ideal V; 3972 3973 for( i = N; i > 0; i  ) 3974 { 3975 k = v[i]; 3976 V[i] = var(i)  k; 3977 } 3978 3979 V = groebner(V); 3980 3981 n = NF(n, V); 3982 d = NF(d, V); 3983 3984 n; 3985 d; 3986 3987 if( d == 0 ) 3988 { 3989 ERROR("Sorry: denominator is zero!"); 3990 } 3991 3992 if( n == 0 ) 3993 { 3994 return (0); 3995 } 3996 3997 poly g = gcd(n, d); 3998 3999 if( g != leadcoef(g) ) 4000 { 4001 n = n / g; 4002 d = d / g; 4003 } 4004 4005 n; 4006 d; 4007 4008 4009 for( i = N; i > 0; i  ) 4010 { 4011 "i: ", i; 4012 n; 4013 d; 4014 4015 k = v[i]; 4016 k; 4017 4018 n = subst(n, var(i), k); 4019 d = subst(d, var(i), k); 4020 4021 if( k != 0 ) 4022 { 4023 k = 1/k; 4024 n = subst(n, var(N+i), k); 4025 d = subst(d, var(N+i), k); 4026 } 4027 } 4028 4029 n; 4030 d; 4031 4032 if( d == 0 ) 4033 { 4034 ERROR("Sorry: denominator is zero!"); 4035 } 4036 4037 if( n == 0 ) 4038 { 4039 return (0); 4040 } 4041 4042 poly g = gcd(n, d); 4043 4044 if( g != leadcoef(g) ) 4045 { 4046 n = n / g; 4047 d = d / g; 4048 } 4049 4050 n; 4051 d; 4052 4053 if( n != leadcoef(n)  d != leadcoef(d) ) 4054 { 4055 ERROR("Sorry cannot completely evaluate. Partial result: (" + string(n) + ")/(" + string(d) + ")"); 4056 } 4057 4058 n; 4059 d; 4060 4061 return (leadcoef(n)/leadcoef(d)); 4062 } 4063 example 4064 { 4065 "EXAMPLE:"; echo=2; 4066 4067 // TODO! 4068 4069 } 4070 4071 4072 proc isPositive() 4073 "USAGE: isPositive() 4074 PURPOSE: Computes whether the multigrading of the ring is positive. 4075 For computation theorem 8.6 of the Miller/Sturmfels book is used. 4076 RETURNS: true if the multigrading is positive 4077 EXAMPLE: example isPositive; shows an example 4078 " 4079 { 4080 ideal I = mDegBasis(0); 4081 ideal J = attrib(I,"ZeroPart"); 4082 /* 4083 I am not quite sure what this ZeroPart is anymore. I thought it 4084 should contain all monomials of degree 0, but then apparently 1 should 4085 be contained. It makes sense to exclude 1, but was this also the intention? 4086 */ 4087 return(J==0); 4088 } 4089 example 4090 { 4091 echo = 2; printlevel = 3; 4092 ring r = 0,(x,y),dp; 4093 intmat A[1][2]=1,1; 4094 setBaseMultigrading(A); 4095 isPositive(); 4096 4097 intmat A[1][2]=1,1; 4098 setBaseMultigrading(A); 4099 isPositive(A); 4100 } 2833 4101 2834 4102 /////////////////////////////////////////////////////////////////////////////// … … 2836 4104 proc testMultigradingLib () 2837 4105 { 2838 echo = 2; printlevel = 3;2839 4106 example setBaseMultigrading; 2840 4107 example setModuleGrading; 2841 4108 2842 4109 example getVariableWeights; 2843 example getTorsion; 4110 example getLattice; 4111 example getGradingGroup; 2844 4112 example getModuleGrading; 2845 4113 … … 2849 4117 2850 4118 2851 example hermite ;2852 example isHomogen ous;4119 example hermiteNormalForm; 4120 example isHomogeneous; 2853 4121 example isTorsionFree; 2854 4122 example pushForward; 2855 example defineHomogen ous;4123 example defineHomogeneous; 2856 4124 2857 4125 example equalMDeg; 2858 example is TorsionElement;4126 example isZeroElement; 2859 4127 2860 4128 example mDegResolution; 2861 2862 "// ******************* example hilbertSeries ************************//"; 4129 4130 "// ******************* example hilbertSeries ************************//"; 2863 4131 example hilbertSeries; 2864 4132 … … 2867 4135 2868 4136 "The End!"; 2869 2870 } 4137 } 4138 4139 4140 proc mDegTruncate(def M, intvec md) 4141 { 4142 "d: "; 4143 print(md); 4144 4145 "M: "; 4146 module LL = M; // + L for d+1 4147 LL; 4148 print(mDeg(LL)); 4149 4150 4151 intmat V = getModuleGrading(M); 4152 intvec vi; 4153 int s = nrows(M); 4154 int r = nrows(V); 4155 int i; 4156 module L; def B; 4157 for (i=s; i>0; i) 4158 { 4159 "comp: ", i; 4160 vi = V[1..r, i]; 4161 "v[i]: "; vi; 4162 4163 B = mDegBasis(md  vi); // ZeroPart is always the same... 4164 "B: "; B; 4165 4166 L = L, B*gen(i); 4167 } 4168 L = simplify(L, 2); 4169 L = setModuleGrading(L,V); 4170 4171 "L: "; L; 4172 print(mDeg(L)); 4173 4174 L = mDegModulo(L, LL); 4175 L = mDegGroebner(L); 4176 // L = minbase(prune(L)); 4177 4178 "??????????"; 4179 print(L); 4180 print(mDeg(L)); 4181 4182 V = getModuleGrading(L); 4183 4184 // take out other degrees 4185 for(i = ncols(L); i > 0; i ) 4186 { 4187 if( !equalMDeg( mDeg(getGradedGenerator(L, i)), md ) ) 4188 { 4189 L[i] = 0; 4190 } 4191 } 4192 4193 L = simplify(L, 2); 4194 L = setModuleGrading(L, V); 4195 print(L); 4196 print(mDeg(L)); 4197 4198 return(L); 4199 } 4200 example 4201 { 4202 "EXAMPLE:"; echo=2; 4203 4204 // TODO! 4205 ring r = 32003, (x,y), dp; 4206 4207 intmat M[2][2] = 4208 1, 0, 4209 0, 1; 4210 4211 setBaseMultigrading(M); 4212 4213 intmat V[2][1] = 4214 0, 4215 0; 4216 4217 "X:"; 4218 module h1 = x; 4219 h1 = setModuleGrading(h1, V); 4220 mDegTruncate(h1, mDeg(x)); 4221 mDegTruncate(h1, mDeg(y)); 4222 4223 "XY:"; 4224 module h2 = ideal(x, y); 4225 h2 = setModuleGrading(h2, V); 4226 mDegTruncate(h2, mDeg(x)); 4227 mDegTruncate(h2, mDeg(y)); 4228 mDegTruncate(h2, mDeg(xy)); 4229 } 4230 4231 4232 /******************************************************/ 4233 /* Some functions on lattices. 4234 TODO Tuebingen:  add functionality (see wiki) and 4235  adjust them to work for groups as well.*/ 4236 /******************************************************/ 4237 4238 4239 4240 /******************************************************/ 4241 proc imageLattice(intmat Q, intmat L) 4242 "USAGE: imageLattice(Q,L); Q and L are of type intmat 4243 PURPOSE: compute an integral basis for the image of the 4244 lattice L under the homomorphism of lattices Q. 4245 RETURN: intmat 4246 EXAMPLE: example imageLattice; shows an example 4247 " 4248 { 4249 intmat Mul = Q*L; 4250 intmat L = latticeBasis(Mul); 4251 4252 return(L); 4253 } 4254 example 4255 { 4256 "EXAMPLE:"; echo=2; 4257 4258 intmat Q[2][3] = 4259 1,2,3, 4260 3,2,1; 4261 4262 intmat L[3][2] = 4263 1,4, 4264 2,5, 4265 3,6; 4266 4267 // should be a 2x2 matrix with columns 4268 // [2,14], [0,36] 4269 imageLattice(Q,L); 4270 4271 } 4272 4273 /******************************************************/ 4274 proc intRank(intmat A) 4275 " 4276 USAGE: intRank(A); intmat A 4277 PURPOSE: compute the rank of the integral matrix A 4278 by computing a hermite normalform. 4279 RETURNS: int 4280 EXAMPLE: example intRank; shows an example 4281 " 4282 { 4283 intmat B = hermiteNormalForm(A); 4284 4285 // get number of zero columns 4286 int nzerocols = 0; 4287 int j; 4288 int i; 4289 int iszero; 4290 for ( j = 1; j <= ncols(B); j++ ) 4291 { 4292 iszero = 1; 4293 4294 for ( i = 1; i <= nrows(B); i++ ) 4295 { 4296 if ( B[i,j] != 0 ) 4297 { 4298 iszero = 0; 4299 break; 4300 } 4301 } 4302 4303 if ( iszero == 1 ) 4304 { 4305 nzerocols++; 4306 } 4307 } 4308 4309 // get number of zero rows 4310 int nzerorows = 0; 4311 4312 for ( i = 1; i <= nrows(B); i++ ) 4313 { 4314 iszero = 1; 4315 4316 for ( j = 1; j <= ncols(B); j++ ) 4317 { 4318 if ( B[i,j] != 0 ) 4319 { 4320 iszero = 0; 4321 break; 4322 } 4323 } 4324 4325 if ( iszero == 1 ) 4326 { 4327 nzerorows++; 4328 } 4329 } 4330 4331 int r = nrows(B)  nzerorows; 4332 4333 if ( ncols(B)  nzerocols < r ) 4334 { 4335 r = ncols(B)  nzerocols; 4336 } 4337 4338 return(r); 4339 } 4340 example 4341 { 4342 4343 intmat A[3][4] = 4344 1,0,1,0, 4345 1,2,0,0, 4346 0,0,0,0; 4347 4348 int r = intRank(A); 4349 4350 print(A); 4351 print(r); // Should be 2 4352 4353 kill A; 4354 4355 } 4356 4357 /*****************************************************/ 4358 4359 proc isSublattice(intmat L, intmat S) 4360 "USAGE: isSublattice(L, S); L, S are of tpye intmat 4361 PURPOSE: checks whether the lattice created by L is a 4362 sublattice of the lattice created by S. 4363 The procedure checks whether each generator of L is 4364 contained in S. 4365 RETURN: 0 if false, 1 if true 4366 EXAMPLE: example isSublattice; shows an example 4367 " 4368 { 4369 int a,b,g,i,j,k; 4370 intmat Ker; 4371 4372 // check whether each column v of L is contained in 4373 // the lattice generated by S 4374 for ( i = 1; i <= ncols(L); i++ ) 4375 { 4376 4377 // v is the ith column of L 4378 intvec v; 4379 for ( j = 1; j <= nrows(L); j++ ) 4380 { 4381 v[j] = L[j,i]; 4382 } 4383 4384 // concatenate B = [S,v] 4385 intmat B[nrows(L)][ncols(S) + 1]; 4386 4387 for ( a = 1; a <= nrows(S); a++ ) 4388 { 4389 for ( b = 1; b <= ncols(S); b++ ) 4390 { 4391 B[a,b] = S[a,b]; 4392 } 4393 } 4394 4395 for ( a = 1; a <= size(v); a++ ) 4396 { 4397 B[a,ncols(B)] = v[a]; 4398 } 4399 4400 4401 // check gcd 4402 Ker = kernelLattice(B); 4403 k = nrows(Ker); 4404 list R; // R is the last row 4405 4406 for ( j = 1; j <= ncols(Ker); j++ ) 4407 { 4408 R[j] = Ker[k,j]; 4409 } 4410 4411 g = R[1]; 4412 4413 for ( j = 2; j <= size(R); j++ ) 4414 { 4415 g = gcd(g,R[j]); 4416 } 4417 4418 if ( g != 1 ) 4419 { 4420 return(0); 4421 } 4422 4423 kill B, v, R; 4424 4425 } 4426 4427 return(1); 4428 } 4429 example 4430 { 4431 "EXAMPLE:"; echo=2; 4432 4433 //ring R = 0,(x,y),dp; 4434 intmat S2[3][2]= 4435 0, 2, 3, 4436 0, 1, 1, 4437 3, 0, 2; 4438 4439 intmat S1[3][3]= 4440 0, 6, 4441 0, 2, 4442 3, 4; 4443 4444 isSublattice(S1,S2); // Yes! 4445 4446 intmat S3[3][1] = 4447 0, 4448 0, 4449 1; 4450 4451 not(isSublattice(S3,S2)); // Yes! 4452 4453 } 4454 4455 /******************************************************/ 4456 4457 proc latticeBasis(intmat B) 4458 "USAGE: latticeBasis(B); intmat B 4459 PURPOSE: compute an integral basis for the lattice defined by 4460 the columns of B. 4461 RETURNS: intmat 4462 EXAMPLE: example latticeBasis; shows an example 4463 " 4464 { 4465 int n = ncols(B); 4466 int r = intRank(B); 4467 4468 if ( r == 0 ) 4469 { 4470 intmat H[nrows(B)][1]; 4471 int j; 4472 4473 for ( j = 1; j <= nrows(B); j++ ) 4474 { 4475 H[j,1] = 0; 4476 } 4477 } 4478 else 4479 { 4480 intmat H = hermiteNormalForm(B);; 4481 4482 if (r < n) 4483 { 4484 // delete columns r+1 to n 4485 // should be identical with the function 4486 // H = submat(H,1..nrows(H),1..r); 4487 // for matrices 4488 intmat Hdel[nrows(H)][r]; 4489 int k; 4490 int m; 4491 4492 for ( k = 1; k <= nrows(H); k++ ) 4493 { 4494 for ( m = 1; m <= r; m++ ) 4495 { 4496 Hdel[k,m] = H[k,m]; 4497 } 4498 } 4499 4500 H = Hdel; 4501 } 4502 } 4503 4504 return(H); 4505 } 4506 example 4507 { 4508 "EXAMPLE:"; echo=2; 4509 4510 intmat L[3][3] = 4511 1,4,8, 4512 2,5,10, 4513 3,6,12; 4514 4515 intmat B = latticeBasis(B); // should be a matrix with columns [1,2,3] and [0,3,6] 4516 4517 kill B,L; 4518 } 4519 4520 /******************************************************/ 4521 4522 proc kernelLattice(def P) 4523 " 4524 USAGE: kernelLattice(P); intmat P 4525 PURPOSE: compute a integral basis for the kernel of the 4526 homomorphism of lattices defined by the intmat P. 4527 RETURNS: intmat 4528 EXAMPLE: example kernelLattice; shows an example 4529 " 4530 { 4531 int n = ncols(P); 4532 int r = intRank(P); 4533 4534 if ( r == 0 ) 4535 { 4536 intmat U = unitMatrix(n); 4537 } 4538 else 4539 { 4540 if ( r == n ) 4541 { 4542 intmat U[n][1]; // now all entries are zero. 4543 } 4544 else 4545 { 4546 list L = hermiteNormalForm(P, "transform"); //hermite(P, "transform"); // now, Hermite = L[1] = A*L[2] 4547 intmat U = L[2]; 4548 4549 // delete columns 1 to r 4550 // should be identical with the function 4551 // U = submat(U,1..nrows(U),r+1..); 4552 // for matrices 4553 intmat Udel[nrows(U)][ncols(U)  r]; 4554 int k; 4555 int m; 4556 4557 for ( k = 1; k <= nrows(U); k++ ) 4558 { 4559 for ( m = r + 1; m <= ncols(U); m++ ) 4560 { 4561 Udel[k,m  r] = U[k,m]; 4562 } 4563 } 4564 4565 U = Udel; 4566 4567 } 4568 } 4569 4570 return(U); 4571 } 4572 example 4573 { 4574 "EXAMPLE"; echo = 2; 4575 4576 intmat LL[3][4] = 4577 1,4,7,10, 4578 2,5,8,11, 4579 3,6,9,12; 4580 4581 // should be a 4x2 matrix with colums 4582 // [1,2,1,0],[2,3,0,1] 4583 intmat B = kernelLattice(LL); 4584 4585 print(B); 4586 4587 kill B; 4588 4589 } 4590 4591 /*****************************************************/ 4592 4593 proc preimageLattice(def P, def B) 4594 " 4595 USAGE: preimageLattice(P, B); intmat P, intmat B 4596 PURPOSE: compute an integral basis for the preimage of B under 4597 the homomorphism of lattices defined by the intmat P. 4598 RETURNS: intmat 4599 EXAMPLE: example preimageLattice; shows an example 4600 " 4601 { 4602 // concatenate matrices: Con = [P,B] 4603 intmat Con[nrows(P)][ncols(P) + ncols(B)]; 4604 int i; 4605 int j; 4606 4607 for ( i = 1; i <= nrows(Con); i++ ) 4608 { 4609 for ( j = 1; j <= ncols(P); j++ ) // P first 4610 { 4611 Con[i,j] = P[i,j]; 4612 } 4613 } 4614 4615 for ( i = 1; i <= nrows(Con); i++ ) 4616 { 4617 for ( j = 1; j <= ncols(B); j++ ) // now B 4618 { 4619 Con[i,ncols(P) + j] =  B[i,j]; 4620 } 4621 } 4622 4623 4624 print(Con); 4625 4626 intmat L = kernelLattice(Con); 4627 4628 print(L); 4629 print(ncols(P)); 4630 print(ncols(L)); 4631 4632 // delete rows ncols(P)+1 to nrows(L) out of L 4633 intmat Del[ncols(P)][ncols(L)]; 4634 int k; 4635 int m; 4636 4637 for ( k = 1; k <= nrows(Del); k++ ) 4638 { 4639 for ( m = 1; m <= ncols(Del); m++ ) 4640 { 4641 Del[k,m] = L[k,m]; 4642 } 4643 } 4644 4645 L = latticeBasis(Del); 4646 4647 return(L); 4648 4649 } 4650 example 4651 { 4652 "EXAMPLE"; echo = 2; 4653 4654 intmat P[2][3] = 4655 2,6,10, 4656 4,8,12; 4657 4658 intmat B[2][1] = 4659 1, 4660 0; 4661 4662 // should be a 2x2 matrix with columns 4663 // [1,1,1] and [2,1,0] 4664 intmat L = preimageLattice(P,B); 4665 4666 kill B, P, L; 4667 4668 } 4669 4670 /******************************************************/ 4671 proc isPrimitiveSublattice(intmat A); 4672 "USAGE: isPrimitiveSublattice(A); intmat A 4673 PURPOSE: check whether the given set of integral vectors in ZZ^m, 4674 i.e. the columns of A, generate a primitive sublattice in ZZ^m 4675 (a direct summand of ZZ^m). 4676 RETURNS: int, where 0 is false and 1 is true. 4677 EXAMPLE: example isPrimitiveSublattice; shows an example 4678 " 4679 { 4680 intmat B = smithNormalForm(A); 4681 int r = intRank(B); 4682 4683 if ( r == 0 ) 4684 { 4685 return(1); 4686 } 4687 4688 if ( 1 < B[r,r] ) 4689 { 4690 return(0); 4691 } 4692 4693 return(1); 4694 } 4695 example 4696 { 4697 "EXAMPLE"; echo = 2; 4698 4699 intmat A[3][2] = 4700 1,4, 4701 2,5, 4702 3,6; 4703 4704 // should be 0 4705 int b = isPrimitiveSublattice(A); 4706 4707 kill A,b; 4708 } 4709 4710 /******************************************************/ 4711 proc isIntegralSurjective(intmat P); 4712 "USAGE: isIntegralSurjective(P); intmat P 4713 PURPOSE: test whether the given linear map P of lattices is 4714 surjective. 4715 RETURNS: int, where 0 is false and 1 is true. 4716 EXAMPLE: example isIntegralSurjective; shows an example 4717 " 4718 { 4719 int r = intRank(P); 4720 4721 if ( r < nrows(P) ) 4722 { 4723 return(0); 4724 } 4725 4726 if ( isPrimitiveSublattice(A) == 1 ) 4727 { 4728 return(1); 4729 } 4730 4731 return(0); 4732 } 4733 example 4734 { 4735 "EXAMPLE"; echo = 2; 4736 4737 intmat A[3][2] = 4738 1,3,5, 4739 2,4,6; 4740 4741 // should be 0 4742 int b = isIntegralSurjective(A); 4743 print(b); 4744 4745 kill A,b; 4746 } 4747 4748 /******************************************************/ 4749 proc projectLattice(intmat B) 4750 "USAGE: projectLattice(B); intmat B 4751 PURPOSE: A set of vectors in ZZ^m is given as the columns of B. 4752 Then this function provides a linear map ZZ^m > ZZ^n 4753 having the primitive span of B its kernel. 4754 RETURNS: intmat 4755 EXAMPLE: example projectLattice; shows an example 4756 " 4757 { 4758 int n = nrows(B); 4759 int r = intRank(B); 4760 4761 if ( r == 0 ) 4762 { 4763 intmat U = unitMatrix(n); 4764 } 4765 else 4766 { 4767 if ( r == n ) 4768 { 4769 intmat U[n][1]; // U now is the ndim zerovector 4770 } 4771 else 4772 { 4773 // we want a matrix with column operations so we transpose 4774 list L = hermiteNormalForm(B, "transform"); //hermite(transpose(B), "transform"); 4775 intmat U = transpose(L[2]); 4776 4777 // delete rows 1 to r 4778 intmat Udel[nrows(U)  r][ncols(U)]; 4779 int k; 4780 int m; 4781 4782 for ( k = 1; k <= nrows(U)  r ; k++ ) 4783 { 4784 for ( m = 1; m <= ncols(U); m++ ) 4785 { 4786 Udel[k,m] = U[k + r,m]; 4787 } 4788 } 4789 4790 U = Udel; 4791 4792 } 4793 } 4794 4795 return(U); 4796 } 4797 example 4798 { 4799 "EXAMPLE"; echo = 2; 4800 4801 intmat B[4][2] = 4802 1,5, 4803 2,6, 4804 3,7, 4805 4,8; 4806 4807 // should result in a (2x4)matrix with columns 4808 // [1, 2], [2, â3], [1, 0] and [0, 1]. 4809 intmat U = projectLattice(B); 4810 4811 } 4812 4813 /******************************************************/ 4814 proc intersectLattices(intmat A, intmat B) 4815 "USAGE: intersectLattices(A, B); intmat A, intmat B 4816 PURPOSE: compute an integral basis for the intersection of the 4817 lattices A and B. 4818 RETURNS: intmat 4819 EXAMPLE: example intersectLattices; shows an example 4820 " 4821 { 4822 // concatenate matrices: Con = [A,B] 4823 intmat Con[nrows(A)][ncols(A) + ncols(B)]; 4824 int i; 4825 int j; 4826 4827 for ( i = 1; i <= nrows(Con); i++ ) 4828 { 4829 for ( j = 1; j <= ncols(A); j++ ) // A first 4830 { 4831 Con[i,j] = A[i,j]; 4832 } 4833 } 4834 4835 for ( i = 1; i <= nrows(Con); i++ ) 4836 { 4837 for ( j = 1; j <= ncols(B); j++ ) // now B 4838 { 4839 Con[i,ncols(A) + j] =  B[i,j]; 4840 } 4841 } 4842 4843 intmat K = kernelLattice(Con); 4844 4845 // delete all rows in K from ncols(A)+1 onwards 4846 intmat Bas[ncols(A)][ncols(K)]; 4847 4848 for ( i = 1; i <= nrows(Bas); i++ ) 4849 { 4850 for ( j = 1; j <= ncols(Bas); j++ ) 4851 { 4852 Bas[i,j] = K[i,j]; 4853 } 4854 } 4855 4856 // take product in order to obtain the intersection 4857 intmat S = A * Bas; 4858 intmat Cut = hermiteNormalForm(S); //hermite(S); 4859 int r = intRank(Cut); 4860 4861 if ( r == 0 ) 4862 { 4863 intmat Cutdel[nrows(Cut)][1]; // is now the zerovector 4864 4865 Cut = Cutdel; 4866 } 4867 else 4868 { 4869 // delte columns from r+1 onwards 4870 intmat Cutdel[nrows(Cut)][r]; 4871 4872 for ( i = 1; i <= nrows(Cutdel); i++ ) 4873 { 4874 for ( j = 1; j <= r; j++ ) 4875 { 4876 Cutdel[i,j] = Cut[i,j]; 4877 } 4878 } 4879 4880 Cut = Cutdel; 4881 } 4882 4883 return(Cut); 4884 } 4885 example 4886 { 4887 "EXAMPLE"; echo = 2; 4888 4889 intmat A[3][2] = 4890 1,4, 4891 2,5, 4892 3,6; 4893 4894 intmat B[3][2] = 4895 6,9, 4896 7,10, 4897 8,11; 4898 4899 // should result in a (2x4)matrix with columns 4900 // [3, 3, 3], [0, 3, 6] 4901 intmat U = intersectLattices(A,B); 4902 4903 } 4904 4905 proc intInverse(intmat A); 4906 "USAGE: intInverse(A); intmat A 4907 PURPOSE: compute the integral inverse of the intmat A. 4908 If det(A) is neither 1 nor 1 an error is returned. 4909 RETURNS: intmat 4910 EXAMPLE: example intInverse; shows an example 4911 " 4912 { 4913 int d = det(A); 4914 4915 if ( d * d != 1 ) // is d = 1 or 1? Else: error 4916 { 4917 ERROR("determinant of the given intmat has to be 1 or 1."); 4918 } 4919 4920 int c; 4921 int i,j; 4922 intmat C[nrows(A)][ncols(A)]; 4923 intmat Ad; 4924 int s; 4925 4926 for ( i = 1; i <= nrows(C); i++ ) 4927 { 4928 for ( j = 1; j <= ncols(C); j++ ) 4929 { 4930 Ad = intAdjoint(A,i,j); 4931 s = 1; 4932 4933 if ( ((i + j) % 2) > 0 ) 4934 { 4935 s = 1; 4936 } 4937 4938 C[i,j] = d * s * intDet(Ad); // mult by d is equal to div by det 4939 } 4940 } 4941 4942 C = transpose(C); 4943 4944 return(C); 4945 } 4946 example 4947 { 4948 "EXAMPLE"; echo = 2; 4949 4950 intmat A[3][3] = 4951 1,1,3, 4952 3,2,0, 4953 0,0,1; 4954 4955 intmat B = intInverse(A); 4956 4957 // should be the unit matrix 4958 print(A * B); 4959 4960 kill A,B; 4961 } 4962 4963 4964 /******************************************************/ 4965 proc intAdjoint(intmat A, int indrow, int indcol) 4966 "USAGE: intAdjoint(A); intmat A 4967 PURPOSE: return the matrix where the given row and column are deleted. 4968 RETURNS: intmat 4969 EXAMPLE: example intAdjoint; shows an example 4970 " 4971 { 4972 int n = nrows(A); 4973 int m = ncols(A); 4974 int i, j; 4975 intmat B[n  1][m  1]; 4976 int a, b; 4977 4978 for ( i = 1; i < indrow; i++ ) 4979 { 4980 for ( j = 1; j < indcol; j++ ) 4981 { 4982 B[i,j] = A[i,j]; 4983 } 4984 for ( j = indcol + 1; j <= ncols(A); j++ ) 4985 { 4986 B[i,j  1] = A[i,j]; 4987 } 4988 } 4989 4990 for ( i = indrow + 1; i <= nrows(A); i++ ) 4991 { 4992 for ( j = 1; j < indcol; j++ ) 4993 { 4994 B[i  1,j] = A[i,j]; 4995 } 4996 for ( j = indcol+1; j <= ncols(A); j++ ) 4997 { 4998 B[i  1,j  1] = A[i,j]; 4999 } 5000 } 5001 5002 return(B); 5003 } 5004 example 5005 { 5006 "EXAMPLE"; echo = 2; 5007 5008 intmat A[2][3] = 5009 1,3,5, 5010 2,4,6; 5011 5012 intmat B = intAdjoint(A,2,2); 5013 print(B); 5014 5015 kill A,B; 5016 } 5017 5018 /******************************************************/ 5019 proc integralSection(intmat P); 5020 "USAGE: integralSection(P); intmat P 5021 PURPOSE: for a given linear surjective map P of lattices 5022 this procedure returns an integral section of P. 5023 RETURNS: intmat 5024 EXAMPLE: example integralSection; shows an example 5025 " 5026 { 5027 int m = nrows(P); 5028 int n = ncols(P); 5029 5030 if ( m == n ) 5031 { 5032 intmat U = intInverse(P); 5033 } 5034 else 5035 { 5036 intmat U = (hermiteNormalForm(P, "transform"))[2]; 5037 5038 // delete columns m+1 to n 5039 intmat Udel[nrows(U)][ncols(U)  (n  m)]; 5040 int k; 5041 int z; 5042 5043 for ( k = 1; k <= nrows(U); k++ ) 5044 { 5045 for ( z = 1; z <= m; z++ ) 5046 { 5047 Udel[k,z] = U[k,z]; 5048 } 5049 } 5050 5051 U = Udel; 5052 } 5053 5054 return(U); 5055 } 5056 example 5057 { 5058 "EXAMPLE"; echo = 2; 5059 5060 intmat P[2][4] = 5061 1,3,4,6, 5062 2,4,5,7; 5063 5064 // should be a matrix with two columns 5065 // for example: [â2, 1, 0, 0], [3, â3, 0, 1] 5066 intmat U = integralSection(P); 5067 5068 print(U); 5069 print(P * U); 5070 5071 kill U; 5072 } 5073 5074 5075 5076 /******************************************************/ 5077 proc factorgroup(G,H) 5078 "USAGE: factorgroup(G,H); list G, list H 5079 PURPOSE: returns a representation of the factor group G mod H using the first isomorphism thm 5080 RETURNS: list 5081 EXAMPLE: example factorgroup(G,H); shows an example 5082 " 5083 { 5084 intmat S1 = G[1]; 5085 intmat L1 = G[2]; 5086 intmat S2 = H[1]; 5087 intmat L2 = H[2]; 5088 5089 // check whether G,H are subgroups of a common group, i.e. whether L1 and L2 span the same lattice 5090 if ( !isSublattice(L1,L2)  !isSublattice(L2,L1)) 5091 { 5092 ERROR("G and H are not subgroups of a common group."); 5093 } 5094 5095 // check whether H is a subgroup of G, i.e. whether S2 is a sublattice of S1+L1 5096 intmat B = concatintmat(S1,L1); // check whether this gives the concatinated matrix 5097 if ( !isSublattice(S2,B) ) 5098 { 5099 ERROR("H is not a subgroup of G"); 5100 } 5101 // use first isomorphism thm to get the factor group 5102 intmat L = concatintmat(L1,S2); // check whether this gives the concatinated matrix 5103 list GmodH; 5104 GmodH[1]=S1; 5105 GmodH[2]=L; 5106 return(GmodH); 5107 } 5108 example 5109 { 5110 "EXAMPLE"; echo = 2; 5111 5112 intmat S1[2][2] = 5113 1,0, 5114 0,1; 5115 intmat L1[2][1] = 5116 2, 5117 0; 5118 5119 intmat S2[2][1] = 5120 1, 5121 0; 5122 intmat L2[2][1] = 5123 2, 5124 0; 5125 5126 list G = createGroup(S1,L1); 5127 list H = createGroup(S2,L2); 5128 5129 list N = factorgroup(G,H); 5130 print(N); 5131 5132 kill G,H,N,S1,L1,S2,L2; 5133 5134 } 5135 5136 /******************************************************/ 5137 proc productgroup(G,H) 5138 "USAGE: productgroup(G,H); list G, list H 5139 PURPOSE: returns a representation of the G x H 5140 RETURNS: list 5141 EXAMPLE: example productgroup(G,H); shows an example 5142 " 5143 { 5144 intmat S1 = G[1]; 5145 intmat L1 = G[2]; 5146 intmat S2 = H[1]; 5147 intmat L2 = H[2]; 5148 intmat OS1[nrows(S1)][ncols(S2)]; 5149 intmat OS2[nrows(S2)][ncols(S1)]; 5150 intmat OL1[nrows(L1)][ncols(L2)]; 5151 intmat OL2[nrows(L2)][ncols(L1)]; 5152 5153 // concatinate matrices to get S 5154 intmat A = concatintmat(S1,OS1); 5155 intmat B = concatintmat(OS2,S2); 5156 intmat At = transpose(A); 5157 intmat Bt = transpose(B); 5158 intmat St = concatintmat(At,Bt); 5159 intmat S = transpose(St); 5160 5161 // concatinate matrices to get L 5162 intmat C = concatintmat(L1,OL1); 5163 intmat D = concatintmat(OL2,L2); 5164 intmat Ct = transpose(C); 5165 intmat Dt = transpose(D); 5166 intmat Lt = concatintmat(Ct,Dt); 5167 intmat L = transpose(Lt); 5168 5169 list GxH; 5170 GxH[1]=S; 5171 GxH[2]=L; 5172 return(GxH); 5173 } 5174 example 5175 { 5176 "EXAMPLE"; echo = 2; 5177 5178 intmat S1[2][2] = 5179 1,0, 5180 0,1; 5181 intmat L1[2][1] = 5182 2, 5183 0; 5184 5185 intmat S2[2][2] = 5186 1,0, 5187 0,2; 5188 intmat L2[2][1] = 5189 0, 5190 3; 5191 5192 list G = createGroup(S1,L1); 5193 list H = createGroup(S2,L2); 5194 5195 list N = productgroup(G,H); 5196 print(N); 5197 5198 kill G,H,N,S1,L1,S2,L2; 5199 5200 } 5201 5202 /******************************************************/ 5203 proc primitiveSpan(intmat V); 5204 "USAGE: isIntegralSurjective(V); intmat V 5205 PURPOSE: compute an integral basis for the minimal primitive 5206 sublattice that contains the given vectors, i.e. the columns of V. 5207 RETURNS: int, where 0 is false and 1 is true. 5208 EXAMPLE: example isIntegralSurjective; shows an example 5209 " 5210 { 5211 int n = ncols(V); 5212 int m = nrows(V); 5213 int r = intRank(V); 5214 5215 5216 if ( r == 0 ) 5217 { 5218 intmat P[m][1]; // this is the mzerovector now 5219 } 5220 else 5221 { 5222 list L = smithNormalForm(V, "transform"); // L = [A,S,B] where S is the smithNF and S = A*S*B 5223 intmat P = intInverse(L[1]); 5224 5225 print(L); 5226 5227 if ( r < m ) 5228 { 5229 // delete columns r+1 to m in P: 5230 intmat Pdel[nrows(P)][r]; 5231 int i,j; 5232 5233 for ( i = 1; i <= nrows(Pdel); i++ ) 5234 { 5235 for ( j = 1; j <= ncols(Pdel); j++ ) 5236 { 5237 Pdel[i,j] = P[i,j]; 5238 } 5239 } 5240 5241 P = Pdel; 5242 } 5243 } 5244 5245 return(P); 5246 } 5247 example 5248 { 5249 "EXAMPLE"; echo = 2; 5250 5251 intmat V[2][4] = 5252 1,4, 5253 2,5, 5254 3,6; 5255 5256 // should return a (4x2)matrix with columns 5257 // [1, 2, 3] and [1, 1, 1] (or similar) 5258 intmat R = primitiveSpan(V); 5259 print(R); 5260 5261 kill V,R; 5262 } 5263 5264 /***********************************************************/
Note: See TracChangeset
for help on using the changeset viewer.