Changeset df39c9 in git
 Timestamp:
 Feb 8, 2007, 4:20:54 PM (16 years ago)
 Branches:
 (u'jengelhdatetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', '0604212ebb110535022efecad887940825b97c3f')
 Children:
 8b2345a5e557bd889ff1e0f5110102e0b4ae3c0c
 Parents:
 302a5d7a57af9072b9b2e8d50715cd8780fae77e
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/finvar.lib
r302a5d rdf39c9 1 1 /////////////////////////////////////////////////////////////////////////////// 2 version="$Id: finvar.lib,v 1.5 8 20070202 13:25:09 SingularExp $"2 version="$Id: finvar.lib,v 1.59 20070208 15:20:54 king Exp $" 3 3 category="Invariant theory"; 4 4 info=" … … 38 38 power_products() exponents for power products 39 39 secondary_char0() secondary invariants (s.i.) in char 0 40 irred_secondary_char0() irreducible secondary invariants in char 0 41 secondary_charp() secondary invariants in char p 42 secondary_no_molien() secondary invariants, without Molien series 40 irred_secondary_char0() irreducible s.i. in char 0 41 secondary_charp() s.i. in char p, with Molien series and 42 Reynolds operator 43 secondary_no_molien() s.i., without Molien series but with 44 Reynolds operator 45 irred_secondary_no_molien() irreducible s.i., without Molien series 43 46 but with Reynolds operator 44 47 secondary_and_irreducibles_no_molien() s.i. & irreducible s.i., without M.s. 45 secondary_not_cohen_macaulay() s.i. when invariant ringnot CohenMacaulay48 secondary_not_cohen_macaulay() s.i. when the invariant ring is not CohenMacaulay 46 49 orbit_variety() ideal of the orbit variety 47 50 rel_orbit_variety() ideal of a relative orbit variety (new version) … … 4305 4308 m=nrows(dimmat); // m1 is the highest degree 4306 4309 if (v) 4307 { " 4310 { "In degree 0 we have: 1"; 4308 4311 ""; 4309 4312 } … … 4321 4324 { // elements in the current degree (i1) 4322 4325 if (v) 4323 { " 4326 { "Searching in degree "+string(i1)+", we need to find "+string(int(dimmat[i,1]))+" invariant(s)..."; 4324 4327 } 4325 4328 TEST=sP; … … 4468 4471 { 4469 4472 // Internal parameters, whose choice might improve the performance  4470 int pieces = 15; // For generating reducible secondaries, blocks of #pieces# secondaries4473 int pieces = 3; // For generating reducible secondaries, blocks of #pieces# secondaries 4471 4474 // are formed. 4472 4475 // If this parameter is small, the memory consumption will decrease. … … 4566 4569 } 4567 4570 if (v) 4568 { " 4571 { "In degree 0 we have: 1"; 4569 4572 ""; 4570 4573 } … … 4580 4583 // This is the only explicit Groebner basis computation! 4581 4584 ideal Reductor,SaveRed; // sP union Reductor is a Groebner basis up to degree i1 4585 ideal sP_Reductor; // will contain sP union Reductor. 4586 int sPOffset = ncols(sP); 4582 4587 int SizeSave; 4583 4588 … … 4585 4590 // minimal degree of a nonconstant invariant factor. 4586 4591 list ISSort; // irr. sec. inv. sorted by degree 4592 int NrIS; 4587 4593 4588 4594 poly helpP; 4589 4595 ideal helpI; 4590 4596 ideal Indicator; // will tell us which candidates for sec. inv. we can choose 4591 ideal ReducedCandidates;4597 // ideal ReducedCandidates; 4592 4598 int helpint; 4593 4599 int k,k2,k3,minD; … … 4608 4614 { // elements in the current degree (i1) 4609 4615 if (v) 4610 { " 4616 { "Searching in degree "+string(i1)+", we need to find "+string(int(dimmat[i,1]))+ 4611 4617 " invariant(s)..."; 4612 4618 " Looking for Power Products..."; 4613 4619 } 4620 sP_Reductor = sP; 4614 4621 counter = 0; // we'll count up to dimmat[i,1] 4615 4622 Reductor = ideal(0); … … 4660 4667 ProdCand = simplify(Mul1*Mul2,4); 4661 4668 // else { ProdCand = Mul1*Mul2;} 4662 ReducedCandidates = reduce(ProdCand,sP);4669 // ReducedCandidates = reduce(ProdCand,sP); 4663 4670 // sP union SaveRed union Reductor is a homogeneous Groebner basis 4664 4671 // up to degree i1. 4665 // We first reduce by sP (which is fixed, so we can do it once for all), 4666 // then by SaveRed resp. by Reductor (which is modified during 4667 // the computations). 4668 Indicator = reduce(ReducedCandidates,SaveRed); 4669 // If Indicator[ii]==0 then ReducedCandidates it the reduction 4670 // of an invariant that is in the algebra generated by primary 4671 // invariants and previously computed secondary invariants. 4672 // Otherwise ProdCand[ii] can be taken as secondary invariant. 4672 // Indicator = reduce(ReducedCandidates,SaveRed); 4673 Indicator = reduce(ProdCand,sP_Reductor); 4674 // If Indicator[ii]<>0 then ProdCand[ii] can be taken as secondary invariant. 4673 4675 if (size(Indicator)<>0) 4674 4676 { for (ii=1;ii<=ncols(ProdCand);ii++) // going through all the power products … … 4681 4683 attrib(SSort[i1][k],"size",saveAttr); 4682 4684 if (v) 4683 { " We found",counter, " of", int(dimmat[i,1]),"secondaries in degree",i1;4685 { " We found",counter, "of", int(dimmat[i,1]),"secondaries in degree",i1; 4684 4686 } 4685 4687 if (int(dimmat[i,1])<>counter) 4686 4688 { Reductor = ideal(helpP); 4687 // Lemma: If G is a homogeneous Groebner basis up to degree i1 and p is a4688 // homogeneous polynomial of degree i1 then G union NF(p,G) is4689 // a homogeneous Groebner basis up to degree i1.4690 4689 attrib(Reductor, "isSB",1); 4691 // if Reductor becomes too large, we reduce the whole of Indicator by4692 // it, save it in SaveRed, and work with a smaller Reductor. This turns4693 // out to save a little time.4694 4690 Indicator=reduce(Indicator,Reductor); 4695 4691 SizeSave++; 4696 4692 SaveRed[SizeSave] = helpP; 4693 sP_Reductor[sPOffset+SizeSave] = helpP; 4694 attrib(sP_Reductor, "isSB",1); 4695 // Lemma: If G is a homogeneous Groebner basis up to degree i1 and p is a 4696 // homogeneous polynomial of degree i1 then G union NF(p,G) is 4697 // a homogeneous Groebner basis up to degree i1. Hence, sP_Reductor 4698 // is a homog. GB up to degree i1. 4697 4699 attrib(SaveRed, "isSB",1); 4700 attrib(Reductor, "isSB",1); 4698 4701 } 4699 4702 else … … 4709 4712 } 4710 4713 4714 NrIS = int(dimmat[i,1])counter; 4711 4715 // The remaining sec. inv. are irreducible! 4712 if ( int(dimmat[i,1])<>counter) // need more than all the power products4716 if (NrIS>0) // need more than all the power products 4713 4717 { if (v) 4714 { " There are irreducible secondary invariants in degree ", i1," !";4715 } 4716 mon = kbase(sP ,i1);4718 { " There are ",NrIS,"irreducible secondary invariants in degree ", i1; 4719 } 4720 mon = kbase(sP_Reductor,i1); 4717 4721 // mon = ideal(mon[ncols(mon)..1]); 4718 4722 ii = ncols(mon); 4719 4723 if ((i>IrrSwitch) or (ii>200)) // use sparse algorithm 4720 { if ( counter==0 && ii==int(dimmat[i,1])) // then we can use all of mon4724 { if (ii+counter==int(dimmat[i,1])) // then we can use all of mon 4721 4725 { B=normalize(evaluate_reynolds(REY,mon)); 4722 4726 IS=IS+B; … … 4730 4734 { ISSort[i1] = B; 4731 4735 } 4732 if (v) {" We found: "; print(B);}4736 if (v) {" We found all",NrIS,"irreducibles in degree ",i1;} 4733 4737 } 4734 4738 else … … 4741 4745 { ideal tmp = normalize(evaluate_reynolds(REY,ideal(mon[j+1..j+MonStep]))); 4742 4746 B[j+1..j+MonStep] = tmp[1..MonStep]; 4743 tmp = reduce(tmp,sP);4744 ReducedCandidates[j+1..j+MonStep] = tmp[1..MonStep];4745 tmp = reduce(tmp, SaveRed);4747 // tmp = reduce(tmp,sP); 4748 // ReducedCandidates[j+1..j+MonStep] = tmp[1..MonStep]; 4749 tmp = reduce(tmp,sP_Reductor); 4746 4750 Indicator[j+1..j+MonStep] = tmp[1..MonStep]; 4747 4751 kill tmp; … … 4750 4754 { ideal tmp = normalize(evaluate_reynolds(REY,ideal(mon[j+1..ii]))); 4751 4755 B[j+1..ii] = tmp[1..iij]; 4752 tmp = reduce(tmp,sP);4753 ReducedCandidates[j+1..ii] = tmp[1..iij];4754 tmp = reduce(tmp, SaveRed);4756 // tmp = reduce(tmp,sP); 4757 // ReducedCandidates[j+1..ii] = tmp[1..iij]; 4758 tmp = reduce(tmp,sP_Reductor); 4755 4759 Indicator[j+1..ii] = tmp[1..iij]; 4756 4760 kill tmp; … … 4771 4775 { ISSort[i1] = ideal(B[j]); 4772 4776 } 4773 if (v) 4774 { " We found the irreducible sec. inv. "+string(B[j]);4777 if (v) 4778 { " We found", irrcounter,"of", NrIS ,"irr. sec. inv. in degree ",i1; 4775 4779 } 4776 4780 Reductor = ideal(helpP); … … 4780 4784 SaveRed[SizeSave] = helpP; 4781 4785 attrib(SaveRed, "isSB",1); 4786 sP_Reductor[sPOffset+SizeSave] = helpP; 4787 attrib(sP_Reductor, "isSB",1); 4782 4788 } 4783 4789 B[j]=0; 4784 ReducedCandidates[j]=0;4790 // ReducedCandidates[j]=0; 4785 4791 Indicator[j]=0; 4786 4792 } … … 4788 4794 } // if i>IrrSwitch 4789 4795 else // use fast algorithm 4790 { B=sort_of_invariant_basis(sP ,REY,i1,int(dimmat[i,1])*6);4796 { B=sort_of_invariant_basis(sP_Reductor,REY,i1,int(dimmat[i,1])*6); 4791 4797 // B contains 4792 4798 // images of kbase(sP,i1) under the … … 4804 4810 { ISSort[i1] = B; 4805 4811 } 4806 if (v) {" We found: "; print(B);}4812 if (v) {" We found all",NrIS,"irreducibles in degree ",i1;} 4807 4813 } 4808 4814 else … … 4810 4816 j=0; // j goes through all of B  4811 4817 // Compare the comments on the computation of reducible sec. inv.! 4812 ReducedCandidates = reduce(B,sP);4813 Indicator = reduce( ReducedCandidates,SaveRed);4818 // ReducedCandidates = reduce(B,sP); 4819 Indicator = reduce(B,sP_Reductor); 4814 4820 while (int(dimmat[i,1])<>counter) 4815 4821 { j++; … … 4827 4833 { ISSort[i1] = ideal(B[j]); 4828 4834 } 4829 if (v) 4830 { " We found the irreducible sec. inv. "+string(B[j]);4835 if (v) 4836 { " We found", irrcounter,"of", NrIS,"irr. sec. inv. in degree ",i1; 4831 4837 } 4832 4838 Reductor = ideal(helpP); … … 4836 4842 SaveRed[SizeSave] = helpP; 4837 4843 attrib(SaveRed, "isSB",1); 4844 sP_Reductor[sPOffset+SizeSave] = helpP; 4845 attrib(sP_Reductor, "isSB",1); 4838 4846 } 4839 4847 B[j]=0; 4840 ReducedCandidates[j]=0;4848 // ReducedCandidates[j]=0; 4841 4849 Indicator[j]=0; 4842 4850 } … … 4892 4900 proc irred_secondary_char0 (matrix P, matrix REY, matrix M, list #) 4893 4901 " 4894 USAGE: irred_secondary_char0(P,REY,M[,v] );4902 USAGE: irred_secondary_char0(P,REY,M[,v][,\"PP\"]); 4895 4903 @* P: a 1xn <matrix> with homogeneous primary invariants, where 4896 4904 n is the number of variables of the basering; … … 4900 4908 Molien series; 4901 4909 @* v: an optional <int>; 4910 @* \"PP\": if this string occurs as (optional) parameter, then in 4911 all degrees power products of irr. sec. inv. will be computed. 4902 4912 RETURN: Irreducible homogeneous secondary invariants of the invariant ring 4903 4913 (type <matrix>) … … 4915 4925 independent modulo the primary invariants using Groebner basis techniques 4916 4926 (see paper \"Fast Computation of Secondary Invariants\" by S. King). 4917 The size of this set 4918 can be read off from the Molien series. Here, only irreducible secondary 4919 invariants are explicitly computed, which saves time and memory. 4927 The size of this set can be read off from the Molien series. Here, only 4928 irreducible secondary invariants are explicitly computed, which saves time and 4929 memory. 4930 @* Moreover, if no irr. sec. inv. in degree d1 have been found and unless the last 4931 optional paramter \"PP\" is used, a Groebner basis of primary invariants and 4932 irreducible secondary invariants up to degree d2 is computed, which allows to 4933 detect irr. sec. inv. in degree d without computing power products. 4920 4934 @* There are three internal parameters \"pieces\", \"MonStep\" and \"IrrSwitch\". 4921 4935 The default values of the parameters should be fine in most cases. However, … … 4927 4941 " 4928 4942 { // Internal parameters, whose choice might improve the performance  4929 int pieces = 15;// For generating reducible secondaries, blocks of #pieces# secondaries4943 int pieces = 3; // For generating reducible secondaries, blocks of #pieces# secondaries 4930 4944 // are formed. 4931 4945 // If this parameter is small, the memory consumption will decrease. 4932 4946 int MonStep = 15; // The Reynolds operator is applied to blocks of #MonStep# monomials. 4933 4947 // If this parameter is small, the memory consumption will decrease. 4934 int IrrSwitch = 7; // Up to degree #IrrSwitch#1, we use a method for generating4935 // irreducible secondaries that tends to produce a smaller output4936 // (which is good for subsequent computations). However, this method4937 // needs to much memory in high degrees, and so we use a sparser method4938 // from degree #IrrSwitch# on.4948 int IrrSwitch = 7; // Up to degree #IrrSwitch#1, or if there are few monomials in kbase(sP,i1), 4949 // we use a method for generating irreducible secondaries that tends to 4950 // produce a smaller output (which is good for subsequent computations). 4951 // However, this method needs to much memory in high degrees, and so we 4952 // use a sparser method from degree #IrrSwitch# upwards. 4939 4953 4940 4954 def br=basering; … … 4947 4961 int i; 4948 4962 if (size(#)>0) 4949 { if (typeof(#[ size(#)])=="int")4950 { int v=#[ size(#)];4963 { if (typeof(#[1])=="int") 4964 { int v=#[1]; 4951 4965 } 4952 4966 else … … 4956 4970 else 4957 4971 { int v=0; 4972 } 4973 4974 int UsePP; 4975 if (size(#)>0) 4976 { if (typeof(#[size(#)])=="string") 4977 { if (#[size(#)]=="PP") 4978 { UsePP = 1; 4979 } 4980 else 4981 { "ERROR: If the last optional parameter is a string, it should be \"PP\"."; 4982 return(matrix(ideal()),matrix(ideal())); 4983 } 4984 } 4958 4985 } 4959 4986 int n=nvars(br); // n is the number of variables, as well … … 4995 5022 m=nrows(dimmat); // m1 is the highest degree 4996 5023 if (v) 4997 { " We need to find";5024 { "There are"; 4998 5025 for (j=1;j<=m;j++) 4999 5026 { if (int(dimmat[j,1])<>1) … … 5006 5033 } 5007 5034 if (v) 5008 { " 5035 { "In degree 0 we have: 1"; 5009 5036 ""; 5010 5037 } … … 5019 5046 option(redSB); 5020 5047 ideal sP = groebner(ideal(P)); 5048 ideal TotalReductor = sP; // will eventually be groebner(ideal(P)+IS) 5049 ideal sP_Reductor; // will contain sP union Reductor. 5050 int sPOffset = ncols(sP); 5051 int LastNewIS=10; // if the Last New Irr. Sec. was found in degree i2 then 5052 // we compute TotalReductor again 5053 int TotalSet = 1; // It TotalSet is 0, we eventually need to recompute TotalReductor. 5021 5054 ideal Reductor,SaveRed; // sP union Reductor is a Groebner basis up to degree i1 5022 5055 int SizeSave; … … 5045 5078 // generating secondary invariants  5046 5079 for (i=2;i<=m;i++) // going through dimmat  5047 { if (int(dimmat[i,1])<>0) // when it is == 0 we need to find no 5080 { 5081 if ((LastNewIS == i3) and (UsePP==0)) 5082 { if (v) {"Computing Groebner basis for primaries and previously found irreducibles...";} 5083 TotalReductor = groebner(sP+IS); 5084 TotalSet = 1; 5085 } 5086 if (int(dimmat[i,1])<>0) // when it is == 0 we need to find no 5048 5087 { // elements in the current degree (i1) 5049 5088 if (v) 5050 { " Searching in degree "+string(i1)+", we need to find "+ 5051 string(int(dimmat[i,1]))+" invariant(s)..."; 5052 " Looking for Power Products..."; 5089 { "Searching in degree "+string(i1)+". There are "+ 5090 string(int(dimmat[i,1]))+" secondary invariant(s)..."; 5053 5091 } 5054 5092 counter = 0; // we'll count up to dimmat[i,1] … … 5060 5098 attrib(SaveRed,"isSB",1); 5061 5099 5100 if ( TotalSet==0 ) 5101 { 5102 if (v) 5103 { 5104 " Looking for Power Products..."; 5105 } 5106 sP_Reductor = sP; 5062 5107 // We start searching for reducible secondary invariants in degree i1, i.e., those 5063 5108 // that are power products of irreducible secondary invariants. … … 5126 5171 attrib(RedSSort[i1][k],"size",saveAttr); 5127 5172 if (v) 5128 { "We found",counter, " of", int(dimmat[i,1])," secondaries in degree",i1; 5173 { 5174 " We found reducible sec. inv. number "+string(counter); 5129 5175 } 5130 5176 if (int(dimmat[i,1])<>counter) 5131 5177 { Reductor = ideal(helpP); 5132 // Lemma: If G is a homogeneous Groebner basis up to degree i1 and p is a5133 // homogeneous polynomial of degree i1 then G union NF(p,G) is5134 // a homogeneous Groebner basis up to degree i1.5135 5178 attrib(Reductor, "isSB",1); 5136 // if Reductor becomes too large, we reduce the whole of Indicator by5137 // it, save it in SaveRed, and work with a smaller Reductor. This turns5138 // out to save a little time.5139 5179 Indicator=reduce(Indicator,Reductor); 5140 5180 SizeSave++; 5141 5181 SaveRed[SizeSave] = helpP; 5182 sP_Reductor[sPOffset+SizeSave] = helpP; 5183 attrib(sP_Reductor, "isSB",1); 5184 // Lemma: If G is a homogeneous Groebner basis up to degree i1 and p is a 5185 // homogeneous polynomial of degree i1 then G union NF(p,G) is 5186 // a homogeneous Groebner basis up to degree i1. Hence, sP_Reductor 5187 // is a homog. GB up to degree i1. 5142 5188 attrib(SaveRed, "isSB",1); 5143 5189 attrib(Reductor, "isSB",1); … … 5146 5192 { break; 5147 5193 } 5148 } 5149 } 5150 } 5194 } // new reducible sec. inv. found 5195 } // loop through ProdCand 5196 } // if there is some reducible sec. inv. 5151 5197 attrib(SaveRed, "isSB",1); 5152 5198 attrib(Reductor, "isSB",1); … … 5155 5201 } 5156 5202 } 5203 } // if TotalSet==0 5204 else 5205 { if (v) {" We don't compute Power Products!"; } 5206 // Instead of computing Power Products, we can compute a Groebner basis of the ideal 5207 // generated by the primary and previously found irr. sec. invariants up to degree i2. An invariant 5208 // polynomial of degree i1 belongs to this ideal if and only if it belongs to the subalgebra 5209 // generated by primary and irreducible secondary invariants up to degree i2. 5210 // Hence, if we find an invariant outside this ideal, it is an irreducible secondary 5211 // invariant of degree i1. 5212 } 5157 5213 5158 5214 // The remaining sec. inv. are irreducible! 5159 5215 if (int(dimmat[i,1])<>counter) // need more than all the power products 5216 // or work without power products 5160 5217 { if (v) 5161 { "There are irreducible secondary invariants in degree ", i1," !"; 5162 } 5163 mon = kbase(sP,i1); 5218 { " Looking for irreducible secondary invariants in degree ", i1; 5219 } 5220 // mon = kbase(sP,i1); 5221 if (TotalSet==0) 5222 { mon = kbase(sP_Reductor,i1); 5223 } 5224 else 5225 { mon = kbase(TotalReductor,i1); 5226 } 5164 5227 ii = ncols(mon); 5165 if ((i>IrrSwitch) or (ii>200)) // use sparse algorithm 5166 { if (counter==0 && ii==int(dimmat[i,1])) // then we can use all of mon 5228 if (mon[1]<>0) 5229 { if ((i>IrrSwitch) or (ii>200)) // use sparse algorithm 5230 { if (TotalSet==0 && ii+counter==int(dimmat[i,1])) // then we can use all of mon 5167 5231 { B=normalize(evaluate_reynolds(REY,mon)); 5168 5232 IS=IS+B; … … 5176 5240 { RedISSort[i1] = NF(B,sP); 5177 5241 } 5178 if (v) {" We found: "; print(B);} 5242 TotalSet=0; // i.e., TotalReductor needs to be recomputed 5243 // if (v) {" We found: "; print(B);} 5244 if (v) {" We found "+string(size(B))+" irred. sec. inv.";} 5179 5245 } 5180 5246 else … … 5182 5248 j=0; // j goes through all of mon  5183 5249 // Compare the comments on the computation of reducible sec. inv.! 5184 while (int(dimmat[i,1])<>counter) 5250 while (j < ii) // we can break the loop if counter==int(dimmat[i,1]) 5251 // while (int(dimmat[i,1])<>counter) 5185 5252 { if ((j mod MonStep) == 0) 5186 5253 { if ((j+MonStep) <= ii) … … 5189 5256 tmp = reduce(tmp,sP); 5190 5257 ReducedCandidates[j+1..j+MonStep] = tmp[1..MonStep]; 5258 if (TotalSet == 0) 5259 { 5191 5260 tmp = reduce(tmp,SaveRed); 5261 } 5262 else 5263 { 5264 tmp = reduce(tmp,TotalReductor); 5265 } 5192 5266 Indicator[j+1..j+MonStep] = tmp[1..MonStep]; 5193 5267 kill tmp; … … 5198 5272 tmp = reduce(tmp,sP); 5199 5273 ReducedCandidates[j+1..ii] = tmp[1..iij]; 5274 if (TotalSet == 0) 5275 { 5200 5276 tmp = reduce(tmp,SaveRed); 5277 } 5278 else 5279 { 5280 tmp = reduce(tmp,TotalReductor); 5281 } 5201 5282 Indicator[j+1..ii] = tmp[1..iij]; 5202 5283 kill tmp; … … 5206 5287 helpP = Indicator[j]; 5207 5288 if (helpP <>0) // B[j] should be added 5208 { counter++; irrcounter++; 5289 { TotalSet=0; // i.e., TotalReductor needs to be recomputed 5290 counter++; irrcounter++; 5209 5291 IS=IS,B[j]; 5210 5292 saveAttr = attrib(RedSSort[i1][i1],"size")+1; … … 5218 5300 } 5219 5301 if (v) 5220 { " We found the irreducible sec. inv. "+string(B[j]); 5302 { //" We found the irreducible sec. inv. "+string(B[j]); 5303 " We found irred. sec. inv. number "+string(irrcounter); 5221 5304 } 5222 5305 Reductor = ideal(helpP); … … 5227 5310 attrib(SaveRed, "isSB",1); 5228 5311 attrib(Reductor, "isSB",1); 5312 LastNewIS = i1; 5313 if (int(dimmat[i,1])==counter) { break; } 5229 5314 } 5230 5315 mon[j]=0; … … 5236 5321 } // if i>IrrSwitch 5237 5322 else // use fast algorithm 5238 { B=sort_of_invariant_basis(sP,REY,i1,int(dimmat[i,1])*6); 5323 { if (TotalSet == 0) 5324 { B=sort_of_invariant_basis(sP_Reductor,REY,i1,int(dimmat[i,1])*6); 5325 } 5326 else 5327 { B=sort_of_invariant_basis(TotalReductor,REY,i1,int(dimmat[i,1])*6); 5328 } 5329 // B=sort_of_invariant_basis(sP,REY,i1,int(dimmat[i,1])*6); 5239 5330 // B is a set of 5240 5331 // images of kbase(sP,i1) under the … … 5242 5333 // independent and that do not reduce to 5243 5334 // 0 modulo sP  5244 if ( counter==0 && ncols(B)==int(dimmat[i,1])) // then we can take all of B5335 if (TotalSet==0 && ncols(B)+counter==int(dimmat[i,1])) // then we can take all of B 5245 5336 { IS=IS+B; 5246 5337 saveAttr = attrib(RedSSort[i1][i1],"size")+int(dimmat[i,1]); … … 5253 5344 { RedISSort[i1] = B; 5254 5345 } 5255 if (v) {" We found: ";print(B);} 5346 // if (v) {" We found: ";print(B);} 5347 if (v) {" We found "+string(size(B))+" irred. sec. inv.";} 5256 5348 } 5257 5349 else … … 5266 5358 helpP = Indicator[j]; 5267 5359 if (helpP <>0) // B[j] should be added 5268 { counter++; irrcounter++; 5360 { TotalSet=0; // i.e., TotalReductor needs to be recomputed 5361 counter++; irrcounter++; 5269 5362 IS=IS,B[j]; 5270 5363 saveAttr = attrib(RedSSort[i1][i1],"size")+1; … … 5278 5371 } 5279 5372 if (v) 5280 { " We found the irreducible sec. inv. "+string(B[j]); 5373 { //" We found the irreducible sec. inv. "+string(B[j]); 5374 " We found irred. sec. inv. number "+string(irrcounter); 5281 5375 } 5282 5376 Reductor = ideal(helpP); … … 5287 5381 attrib(SaveRed, "isSB",1); 5288 5382 attrib(Reductor, "isSB",1); 5383 LastNewIS = i1; 5384 if (int(dimmat[i,1])==counter) { break; } 5289 5385 } 5290 5386 B[j]=0; … … 5294 5390 } 5295 5391 } // i<=IrrSwitch 5392 } // if mon[1]<>0 5296 5393 } // Computation of irreducible secondaries 5297 5394 if (v) … … 5424 5521 } 5425 5522 if (v) 5426 { " 5523 { "In degree 0 we have: 1"; 5427 5524 ""; 5428 5525 } … … 5441 5538 { // elements in the current degree (i1) 5442 5539 if (v) 5443 { " 5540 { "Searching in degree "+string(i1)+", we need to find "+string(deg_dim_vec[i])+" invariant(s)..."; 5444 5541 } 5445 5542 TEST=sP; … … 5587 5684 { 5588 5685 // Internal parameters, whose choice might improve the performance  5589 int pieces = 15; // For generating reducible secondaries, blocks of #pieces# secondaries5686 int pieces = 3; // For generating reducible secondaries, blocks of #pieces# secondaries 5590 5687 // are formed. 5591 5688 // If this parameter is small, the memory consumption will decrease. … … 5695 5792 } 5696 5793 if (v) 5697 { " 5794 { "In degree 0 we have: 1"; 5698 5795 ""; 5699 5796 } … … 5710 5807 // This is the only explicit Groebner basis computation! 5711 5808 ideal Reductor,SaveRed; // sP union Reductor is a Groebner basis up to degree i1 5809 ideal sP_Reductor; // will contain sP union Reductor. 5810 int sPOffset = ncols(sP); 5712 5811 int SizeSave; 5713 5812 … … 5715 5814 // minimal degree of a nonconstant invariant factor. 5716 5815 list ISSort; // irr. sec. inv. sorted by degree 5816 int NrIS; 5717 5817 5718 5818 poly helpP; … … 5738 5838 { // elements in the current degree (i1) 5739 5839 if (v) 5740 { " 5840 { "Searching in degree "+string(i1)+", we need to find "+string(deg_dim_vec[i])+ 5741 5841 " invariant(s)..."; 5742 5842 " Looking for Power Products..."; … … 5744 5844 counter = 0; // we'll count up to deg_dim_vec[i] 5745 5845 Reductor = ideal(0); 5846 sP_Reductor = sP; 5746 5847 helpint = 0; 5747 5848 SaveRed = Reductor; … … 5811 5912 attrib(SSort[i1][k],"size",saveAttr); 5812 5913 if (v) 5813 { " We found",counter, " of", deg_dim_vec[i],"secondaries in degree",i1;5914 { " We found", counter, "of", deg_dim_vec[i], "secondaries in degree",i1; 5814 5915 } 5815 5916 if (deg_dim_vec[i]<>counter) 5816 5917 { Reductor = ideal(helpP); 5918 attrib(Reductor, "isSB",1); 5919 Indicator=reduce(Indicator,Reductor); 5920 SizeSave++; 5921 sP_Reductor[sPOffset+SizeSave] = helpP; 5922 attrib(sP_Reductor, "isSB",1); 5817 5923 // Lemma: If G is a homogeneous Groebner basis up to degree i1 and p is a 5818 5924 // homogeneous polynomial of degree i1 then G union NF(p,G) is 5819 // a homogeneous Groebner basis up to degree i1. 5820 attrib(Reductor, "isSB",1); 5821 // if Reductor becomes too large, we reduce the whole of Indicator by 5822 // it, save it in SaveRed, and work with a smaller Reductor. This turns 5823 // out to save a little time. 5824 Indicator=reduce(Indicator,Reductor); 5825 SizeSave++; 5925 // a homogeneous Groebner basis up to degree i1. Hence, sP_Reductor 5926 // is a homog. GB up to degree i1. 5826 5927 SaveRed[SizeSave] = helpP; 5827 5928 attrib(SaveRed, "isSB",1); … … 5839 5940 } 5840 5941 5942 NrIS = deg_dim_vec[i]counter; 5841 5943 // The remaining sec. inv. are irreducible! 5842 if ( deg_dim_vec[i]<>counter) // need more than all the power products5944 if (NrIS>0) // need more than all the power products 5843 5945 { if (v) 5844 { " There are irreducible secondary invariants in degree ", i1," !";5845 } 5846 mon = kbase(sP ,i1);5946 { " There are ",NrIS,"irreducible secondary invariants in degree ", i1; 5947 } 5948 mon = kbase(sP_Reductor,i1); 5847 5949 ii = ncols(mon); 5848 5950 // mon = ideal(mon[ncols(mon)..1]); … … 5860 5962 { ISSort[i1] = B; 5861 5963 } 5862 if (v) {" We found: "; print(B);}5964 if (v) {" We found all",NrIS,"irreducibles in degree ",i1;} 5863 5965 } 5864 5966 else … … 5901 6003 { ISSort[i1] = ideal(B[j]); 5902 6004 } 5903 if (v) 5904 { " We found the irreducible sec. inv. "+string(B[j]);6005 if (v) 6006 { " We found", irrcounter,"of", NrIS ,"irr. sec. inv. in degree ",i1; 5905 6007 } 5906 6008 Reductor = ideal(helpP); … … 5934 6036 { ISSort[i1] = B; 5935 6037 } 5936 if (v) {" We found: "; print(B);}6038 if (v) {" We found all",NrIS,"irreducibles in degree ",i1;} 5937 6039 } 5938 6040 else … … 5957 6059 { ISSort[i1] = ideal(B[j]); 5958 6060 } 5959 if (v) 5960 { " We found the irreducible sec. inv. "+string(B[j]);6061 if (v) 6062 { " We found", irrcounter,"of", NrIS ,"irr. sec. inv. in degree ",i1; 5961 6063 } 5962 6064 Reductor = ideal(helpP); … … 6040 6142 hence the number of secondary invariants is the product of the degrees of 6041 6143 primary invariants divided by the group order. 6144 <secondary_and_irreducible_no_molien> should usually be faster and of more useful 6145 functionality. 6146 SEE ALSO: secondary_and_irreducible_no_molien 6042 6147 EXAMPLE: example secondary_no_molien; shows an example 6043 6148 " … … 6123 6228 { " We need to find "+string(max)+" secondary invariants."; 6124 6229 ""; 6125 " 6230 "In degree 0 we have: 1"; 6126 6231 ""; 6127 6232 } … … 6160 6265 if (deg_vec[k]<>i) 6161 6266 { if (v) 6162 { " 6267 { "Searching in degree "+string(i)+"..."; 6163 6268 } 6164 6269 mon = kbase(sP,i); … … 6194 6299 S[counter] = B[j]; 6195 6300 if (v) 6196 { " We found sec. inv. number "+string(counter)+" in degree "+string(i);6301 { " We found sec. inv. number "+string(counter)+" in degree "+string(i); 6197 6302 } 6198 6303 if (counter == max) {break;} … … 6226 6331 S[counter] = B[j]; 6227 6332 if (v) 6228 { " We found sec. inv. number "+string(counter)+" in degree "+string(i);6333 { " We found sec. inv. number "+string(counter)+" in degree "+string(i); 6229 6334 } 6230 6335 if (counter == max) {break;} … … 6292 6397 EXAMPLE: example secondary_and_irreducibles_no_molien; shows an example 6293 6398 " 6294 { int pieces = 15;// For generating reducible secondaries, blocks of #pieces# secondaries6399 { int pieces = 3; // For generating reducible secondaries, blocks of #pieces# secondaries 6295 6400 // are formed. 6296 6401 // If this parameter is small, the memory consumption will decrease. … … 6376 6481 { " We need to find "+string(max)+" secondary invariants."; 6377 6482 ""; 6378 " 6483 "In degree 0 we have: 1"; 6379 6484 ""; 6380 6485 } … … 6429 6534 if (deg_vec[l]<>i1) 6430 6535 { if (v) 6431 { " 6536 { "Searching in degree "+string(i1); 6432 6537 " Looking for Power Products..."; 6433 6538 } … … 6499 6604 attrib(SSort[i1][k],"size",saveAttr); 6500 6605 if (v) 6501 { " We found sec. inv. number",counter, "in degree",i1;6606 { " We found sec. inv. number",counter, "in degree",i1; 6502 6607 } 6503 6608 if (max<>counter) … … 6530 6635 if (max<>counter) 6531 6636 { if (v) 6532 { " Looking for irreducible secondary invariants in degree ", i1;6637 { " Looking for irreducible secondary invariants in degree ", i1; 6533 6638 } 6534 6639 mon = kbase(sP,i1); … … 6575 6680 } 6576 6681 if (v) 6577 { " We found the irreducible sec. inv. "+string(B[j]);6682 { " We found irreducible sec. inv. number", irrcounter, "in degree", i1; 6578 6683 } 6579 6684 if (counter == max) { break; } … … 6618 6723 } 6619 6724 if (v) 6620 { " We found the irreducible sec. inv. "+string(B[j]);6725 { " We found irreducible sec. inv. number", irrcounter, "in degree", i1; 6621 6726 } 6622 6727 Reductor = ideal(helpP); … … 6684 6789 matrix S,IS=secondary_and_irreducibles_no_molien(L[1..2],intvec(1,2),1); 6685 6790 print(S); 6791 print(IS); 6792 } 6793 /////////////////////////////////////////////////////////////////////////////// 6794 6795 proc irred_secondary_no_molien (matrix P, matrix REY, list #) 6796 "USAGE: irred_secondary_no_molien(P,REY[,deg_vec,v]); 6797 P: a 1xn <matrix> with primary invariants, REY: a gxn <matrix> 6798 representing the Reynolds operator, deg_vec: an optional <intvec> 6799 listing some degrees where no irreducible secondary invariants can 6800 be found, v: an optional <int> 6801 ASSUME: n is the number of variables of the basering, g the size of the group, 6802 REY is the 1st return value of group_reynolds(), reynolds_molien() or 6803 the second one of primary_invariants() 6804 RETURN: Irreducible secondary invariants of the invariant ring (type <matrix>) 6805 DISPLAY: information if v does not equal 0 6806 THEORY: Irred. secondary invariants are calculated by finding a basis (in terms of 6807 monomials) of the basering modulo primary and previously found secondary 6808 invariants, mapping those to invariants with the Reynolds operator. Among 6809 these images we pick a maximal subset that is linearly independent modulo 6810 the primary and previously found secondary invariants, using Groebner basis 6811 techniques. 6812 EXAMPLE: example irred_secondary_no_molien; shows an example 6813 " 6814 { int MonStep = 29; 6815 6816 def br=basering; 6817 // checking input and setting verbose  6818 if (size(#)==1 or size(#)==2) 6819 { if (typeof(#[size(#)])=="int") 6820 { if (size(#)==2) 6821 { if (typeof(#[size(#)1])=="intvec") 6822 { intvec deg_vec=#[size(#)1]; 6823 } 6824 else 6825 { "ERROR: the third parameter should be an <intvec>"; 6826 return(); 6827 } 6828 } 6829 int v=#[size(#)]; 6830 } 6831 else 6832 { if (size(#)==1) 6833 { if (typeof(#[size(#)])=="intvec") 6834 { intvec deg_vec=#[size(#)]; 6835 int v=0; 6836 } 6837 else 6838 { "ERROR: the third parameter should be an <intvec>"; 6839 return(); 6840 } 6841 } 6842 else 6843 { "ERROR: wrong list of parameters"; 6844 return(); 6845 } 6846 } 6847 } 6848 else 6849 { if (size(#)>2) 6850 { "ERROR: there are too many parameters"; 6851 return(); 6852 } 6853 int v=0; 6854 } 6855 int n=nvars(basering); // n is the number of variables, as well 6856 // as the size of the matrices, as well 6857 // as the number of primary invariants, 6858 // we should get 6859 if (ncols(P)<>n) 6860 { "ERROR: The first parameter ought to be the matrix of the primary"; 6861 " invariants." 6862 return(); 6863 } 6864 if (ncols(REY)<>n) 6865 { "ERROR: The second parameter ought to be the Reynolds operator."; 6866 return(); 6867 } 6868 if (char(br)<>0) 6869 { if (nrows(REY) mod char(br) == 0) 6870 { "ERROR: We need to be in the nonmodular case."; 6871 return(); 6872 } 6873 } 6874 if (v && voice==2) 6875 { ""; 6876 } 6877 if (defined(deg_vec)<>voice) 6878 { intvec deg_vec; 6879 } 6880 int l=1; 6881 6882 // initializing variables  6883 6884 int dgb = degBound; 6885 degBound = 0; 6886 int d,block, counter, totalcount, monsize,i,j, fin,hlp; 6887 poly helpP,lmp; 6888 block = 1; 6889 ideal mon, MON, Inv, IS,RedInv, NewIS, NewReductor; 6890 int NewIScounter; 6891 ideal sP = slimgb(ideal(P)); 6892 int Max = vdim(sP); 6893 if (Max<0) 6894 { "ERROR: The second parameter ought to be the Reynolds operator."; 6895 return(); 6896 } 6897 ideal sIS = sP; // sIS will be a Groebner basis up to degree d of P+irred. sec. inv. 6898 while (totalcount < Max) 6899 { d++; 6900 if (deg_vec[l]<>d) 6901 { if (v) { "Searching irred. sec. inv. in degree "+string(d); } 6902 NewIS = ideal(); 6903 NewIScounter=0; 6904 MON = kbase(sP,d); 6905 mon = kbase(sIS,d); 6906 monsize=ncols(mon); 6907 totalcount = totalcount + ncols(MON); 6908 if (mon[1]<>0) 6909 { 6910 if (v) { " We have "+string(monsize)+" candidates for irred. secondaries";} 6911 block=0; // Loops through the monomials 6912 while (block<monsize) 6913 { if ((block mod MonStep) == 0) 6914 { if ((block+MonStep) <= monsize) 6915 { ideal tmp = normalize(evaluate_reynolds(REY,ideal(mon[block+1..block+MonStep]))); 6916 Inv[block+1..block+MonStep] = tmp[1..MonStep]; 6917 tmp = reduce(tmp,sIS); 6918 RedInv[block+1..block+MonStep] = tmp[1..MonStep]; 6919 kill tmp; 6920 } 6921 else 6922 { ideal tmp = normalize(evaluate_reynolds(REY,ideal(mon[block+1..monsize]))); 6923 Inv[block+1..monsize] = tmp[1..monsizeblock]; 6924 tmp = reduce(tmp,sIS); 6925 RedInv[block+1..monsize] = tmp[1..monsizeblock]; 6926 kill tmp; 6927 } 6928 } 6929 block++; 6930 helpP = RedInv[block]; 6931 if (helpP<>0) 6932 { counter++; // found new irr. sec. inv.! 6933 if (v) { " We found irred. sec. inv. number "+string(counter)+" in degree "+string(d); } 6934 IS[counter] = Inv[block]; 6935 sIS = sIS,helpP; 6936 attrib(sIS,"isSB",1); 6937 NewIScounter++; 6938 NewIS[NewIScounter] = helpP; 6939 // for permutation of variables, the following might yield a speedup; 6940 // however, for general group actions it does not work: 6941 // for (i=block;i<=monsize;i++) 6942 // { if (Inv[block]<>reduce(Inv[block],std(mon[i]))) 6943 // { mon[i] = 0; 6944 // } 6945 // } 6946 6947 mon[block]=0; 6948 Inv[block]=0; 6949 RedInv[block]=0; 6950 NewReductor[1]=helpP; 6951 attrib(NewReductor,"isSB",1); 6952 RedInv = reduce(RedInv,NewReductor); 6953 } 6954 } // loop through monomials 6955 } // if mon[1]<>0 6956 if (totalcount < Max) 6957 { if (NewIScounter==0) 6958 { if (fin==0) 6959 { degBound = 0; 6960 sIS = groebner(sIS); 6961 fin=1; 6962 } 6963 } 6964 else 6965 { degBound = d+1; 6966 sIS = groebner(sIS+NewIS); 6967 degBound = 0; 6968 fin = 0; 6969 } 6970 } 6971 attrib(sIS,"isSB",1); 6972 } // if degree d is not excluded by deg_vec 6973 else 6974 { if (size(deg_vec)==l) 6975 { l=1; 6976 } 6977 else 6978 { l++; 6979 } 6980 } 6981 } // loop through degree 6982 degBound = dgb; 6983 return(matrix(IS)); 6984 } 6985 6986 example 6987 { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; echo=2; 6988 ring R=3,(x,y,z),dp; 6989 matrix A[3][3]=0,1,0,1,0,0,0,0,1; 6990 list L=primary_invariants(A,intvec(1,1,0)); 6991 // In that example, there are no secondary invariants 6992 // in degree 1 or 2. 6993 matrix IS=irred_secondary_no_molien(L[1..2],intvec(1,2),1); 6686 6994 print(IS); 6687 6995 } … … 6881 7189 RETURN: primary and secondary invariants (both of type <matrix>) generating 6882 7190 the invariant ring with respect to the matrix group generated by the 6883 matrices in the input and irreducible secondary invariants (type6884 <matrix>) if the Molien series was available7191 matrices in the input, and irreducible secondary invariants if we are 7192 in the nonmodular case. 6885 7193 DISPLAY: information about the various stages of the program if the third flag 6886 7194 does not equal 0 … … 6977 7285 list l=primary_charp_no_molien(L[1],v); 6978 7286 if (size(l)==2) 6979 { matrix S =secondary_no_molien(l[1],L[1],l[2],v);7287 { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],l[2],v); 6980 7288 } 6981 7289 else 6982 { matrix S =secondary_no_molien(l[1],L[1],v);7290 { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],v); 6983 7291 } 6984 return(l[1],S );7292 return(l[1],S,IS); 6985 7293 } 6986 7294 } … … 7003 7311 { list l=primary_char0_no_molien(L[1],v); 7004 7312 if (size(l)==2) 7005 { matrix S =secondary_no_molien(l[1],L[1],l[2],v);7313 { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],l[2],v); 7006 7314 } 7007 7315 else 7008 { matrix S =secondary_no_molien(l[1],L[1],v);7009 } 7010 return(l[1],S );7316 { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],v); 7317 } 7318 return(l[1],S,IS); 7011 7319 } 7012 7320 else … … 7014 7322 { list l=primary_charp_no_molien(L[1],v); // case 7015 7323 if (size(l)==2) 7016 { matrix S =secondary_no_molien(l[1],L[1],l[2],v);7324 { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],l[2],v); 7017 7325 } 7018 7326 else 7019 { matrix S =secondary_no_molien(l[1],L[1],v);7020 } 7021 return(l[1],S );7327 { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],v); 7328 } 7329 return(l[1],S,IS); 7022 7330 } 7023 7331 else // the modular case … … 7098 7406 among the coefficients) 7099 7407 RETURN: primary and secondary invariants (both of type <matrix>) generating 7100 invariant ring with respect to the matrix group generated by the7101 matrices in the input and irreducible secondary invariants (type7102 <matrix>) if the Molien series was available7408 the invariant ring with respect to the matrix group generated by the 7409 matrices in the input, and irreducible secondary invariants if we are 7410 in the nonmodular case. 7103 7411 DISPLAY: information about the various stages of the program if the third flag 7104 7412 does not equal 0 … … 7212 7520 list l=primary_charp_no_molien_random(L[1],max,v); 7213 7521 if (size(l)==2) 7214 { matrix S =secondary_no_molien(l[1],L[1],l[2],v);7522 { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],l[2],v); 7215 7523 } 7216 7524 else 7217 { matrix S =secondary_no_molien(l[1],L[1],v);7525 { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],v); 7218 7526 } 7219 return(l[1],S );7527 return(l[1],S,IS); 7220 7528 } 7221 7529 } … … 7238 7546 { list l=primary_char0_no_molien_random(L[1],max,v); 7239 7547 if (size(l)==2) 7240 { matrix S =secondary_no_molien(l[1],L[1],l[2],v);7548 { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],l[2],v); 7241 7549 } 7242 7550 else 7243 { matrix S =secondary_no_molien(l[1],L[1],v);7244 } 7245 return(l[1],S );7551 { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],v); 7552 } 7553 return(l[1],S,IS); 7246 7554 } 7247 7555 else … … 7249 7557 { list l=primary_charp_no_molien_random(L[1],max,v); // case 7250 7558 if (size(l)==2) 7251 { matrix S =secondary_no_molien(l[1],L[1],l[2],v);7559 { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],l[2],v); 7252 7560 } 7253 7561 else 7254 { matrix S =secondary_no_molien(l[1],L[1],v);7255 } 7256 return(l[1],S );7562 { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],v); 7563 } 7564 return(l[1],S,IS); 7257 7565 } 7258 7566 else // the modular case
Note: See TracChangeset
for help on using the changeset viewer.