Changeset 18bd9c in git for Singular/LIB/finvar.lib
- Timestamp:
- Jul 21, 1999, 6:52:04 PM (25 years ago)
- Branches:
- (u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
- Children:
- f92cf85a06c9548f084c06b5adcdb9808ec2c69b
- Parents:
- 1f75da5a9bbe84595567f0e42d2bca935d67bbce
- File:
-
- 1 edited
-
Singular/LIB/finvar.lib (modified) (127 diffs)
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/finvar.lib
r1f75da r18bd9c 1 // $Id: finvar.lib,v 1.1 7 1999-07-19 14:07:31 obachmanExp $1 // $Id: finvar.lib,v 1.18 1999-07-21 16:52:00 Singular Exp $ 2 2 // author: Agnes Eileen Heydtmann, email:agnes@math.uni-sb.de 3 3 // last change: 98/11/05 4 4 ////////////////////////////////////////////////////////////////////////////// 5 version="$Id: finvar.lib,v 1.1 7 1999-07-19 14:07:31 obachmanExp $"5 version="$Id: finvar.lib,v 1.18 1999-07-21 16:52:00 Singular Exp $" 6 6 info=" 7 7 LIBRARY: finvar.lib LIBRARY TO CALCULATE INVARIANT RINGS & MORE … … 72 72 /////////////////////////////////////////////////////////////////////////////// 73 73 proc unique (list #) 74 { for (int i=1;i<size(#);i =i+1)74 { for (int i=1;i<size(#);i++) 75 75 { if (#[i]==#[size(#)]) 76 76 { return(0); … … 85 85 RETURNS: the i-th cyclotomic polynomial (type <poly>) as one in the first ring 86 86 variable 87 THEORY: x^i-1 is divided by the j-th cyclotomic polynomial where j takes on 87 THEORY: x^i-1 is divided by the j-th cyclotomic polynomial where j takes on 88 88 the value of proper divisors of i 89 89 EXAMPLE: example cyclotomic; shows an example … … 119 119 { break; 120 120 } 121 j =j+1;121 j++; 122 122 } 123 123 min=min/leadcoef(min); // making sure that the leading … … 154 154 EXAMPLE: example group_reynolds; shows an example 155 155 " 156 { int ch=char(basering); // the existance of the Reynolds operator 156 { int ch=char(basering); // the existance of the Reynolds operator 157 157 // is dependent on the characteristic of 158 158 // the base field … … 198 198 while (TEST<>I) 199 199 { TEST=TEST*G(1); 200 o1 =o1+1;200 o1++; 201 201 } 202 202 } 203 203 int i=1; 204 204 // -------------- doubles among the generators should be avoided ------------- 205 for (int j=2;j<=gen_num;j =j+1)// this loop adds the parameters to the205 for (int j=2;j<=gen_num;j++) // this loop adds the parameters to the 206 206 { // group, leaving out doubles and 207 207 // checking whether the parameters are … … 217 217 } 218 218 if (unique(G(1..i),#[j])) 219 { i =i+1;219 { i++; 220 220 matrix G(i)=#[j]; 221 221 if (ch<>0 && minpoly==0) // finding out of which order the group … … 224 224 while (TEST<>I) 225 225 { TEST=TEST*G(i); 226 o2 =o2+1;226 o2++; 227 227 } 228 228 o1=o1*o2/gcd(o1,o2); // lowest common multiple of the element … … 246 246 { l=0; // l is the number of products we get in 247 247 // one going 248 for (m=g-j+1;m<=g;m =m+1)249 { for (k=1;k<=i;k =k+1)250 { l =l+1;248 for (m=g-j+1;m<=g;m++) 249 { for (k=1;k<=i;k++) 250 { l++; 251 251 matrix P(l)=G(k)*G(m); // possible new element 252 252 } 253 253 } 254 254 j=0; 255 for (k=1;k<=l;k =k+1)255 for (k=1;k<=l;k++) 256 256 { if (unique(G(1..g),P(k))) 257 { j =j+1;// a new factor for next run258 g =g+1;257 { j++; // a new factor for next run 258 g++; 259 259 matrix G(g)=P(k); // a new group element - 260 260 if (ch<>0 && minpoly==0 && i<>1) //finding out of which order the group … … 263 263 while (TEST<>I) 264 264 { TEST=TEST*G(g); 265 o2 =o2+1;265 o2++; 266 266 } 267 267 o1=o1*o2/gcd(o1,o2); // lowest common multiple of the element … … 332 332 { int i=0; 333 333 while((n/root^i)<>1) 334 { i =i+1;334 { i++; 335 335 } 336 336 return(i); … … 342 342 G1,G2,...: nxn <matrices> generating a finite matrix group, ringname: 343 343 a <string> giving a name for a new ring of characteristic 0 for the 344 Molien series in case of prime characteristic, lcm: an <int> giving 344 Molien series in case of prime characteristic, lcm: an <int> giving 345 345 the lowest common multiple of the elements' orders in case of prime 346 346 characteristic, minpoly==0 and a non-cyclic group, flags: an optional … … 363 363 Molien series (named M) is stored 364 364 DISPLAY: information if the third component of flags does not equal 0 365 THEORY: In characteristic 0 the terms 1/det(1-xE) for all group elements of 365 THEORY: In characteristic 0 the terms 1/det(1-xE) for all group elements of 366 366 the Molien series are computed in a straight forward way. In prime 367 367 characteristic a Brauer lift is involved. The returned matrix gives 368 enumerator and denominator of the expanded version where common 368 enumerator and denominator of the expanded version where common 369 369 factors have been canceled. 370 370 EXAMPLE: example molien; shows an example … … 538 538 // new term of the Molien series 539 539 //------------ computing 1/det(1+xE) for all E in the group ------------------ 540 for (int j=1;j<=g;j =j+1)540 for (int j=1;j<=g;j++) 541 541 { if (not(typeof(#[j])=="matrix")) 542 542 { "ERROR: the parameters must be a list of matrices and maybe an <intvec>"; … … 634 634 int i, j, k; 635 635 // -------------- finding all the terms of the Molien series ----------------- 636 for (i=1;i<=g;i =i+1)636 for (i=1;i<=g;i++) 637 637 { setring br; 638 638 if (not(typeof(#[i])=="matrix")) … … 649 649 F=L[1]; 650 650 C=L[2]; 651 for (j=2;j<=ncols(F);j =j+1)651 for (j=2;j<=ncols(F);j++) 652 652 { F[j]=-1*(F[j]-x); // F[j] is now an eigenvalue of G(i), 653 653 // it is a power of a primitive r-th root … … 672 672 // { break; 673 673 // } 674 // k =k+1;674 // k++; 675 675 // } 676 676 setring `newring`; … … 680 680 { if (i<>g) 681 681 { Mstring=string(M); 682 for (j=1;j<=size(Mstring);j =j+1)682 for (j=1;j<=size(Mstring);j++) 683 683 { if (Mstring[j]=="e") 684 684 { interval=0; … … 758 758 //----------------------- simulating characteristic 0 ------------------------ 759 759 string chst=charstr(br); 760 for (int i=1;i<=size(chst);i =i+1)760 for (int i=1;i<=size(chst);i++) 761 761 { if (chst[i]==",") 762 762 { break; … … 794 794 string links, rechts; 795 795 //----------------- finding all terms of the Molien series ------------------- 796 for (i=1;i<=g;i =i+1)796 for (i=1;i<=g;i++) 797 797 { setring br; 798 798 if (not(typeof(#[i])=="matrix")) … … 805 805 } 806 806 string stM(i)=string(#[i]); 807 for (j=1;j<=size(stM(i));j =j+1)807 for (j=1;j<=size(stM(i));j++) 808 808 { if (stM(i)[j]==" 809 809 ") … … 1048 1048 // new term of the Molien series 1049 1049 int i=1; 1050 for (int j=2;j<=gen_num;j =j+1)// this loop adds the parameters to the1050 for (int j=2;j<=gen_num;j++) // this loop adds the parameters to the 1051 1051 { // group, leaving out doubles and 1052 1052 // checking whether the parameters are … … 1062 1062 } 1063 1063 if (unique(G(1..i),#[j])) 1064 { i =i+1;1064 { i++; 1065 1065 matrix G(i)=#[j]; 1066 1066 A(1)=concat(A(1),#[j]*vars); // adding ring homomorphisms to A(1) - … … 1098 1098 { l=0; // l is the number of products we get in 1099 1099 // one going 1100 for (m=g-j+1;m<=g;m =m+1)1101 { for (k=1;k<=i;k =k+1)1102 { l =l+1;1100 for (m=g-j+1;m<=g;m++) 1101 { for (k=1;k<=i;k++) 1102 { l++; 1103 1103 matrix P(l)=G(k)*G(m); // possible new element 1104 1104 } 1105 1105 } 1106 1106 j=0; 1107 for (k=1;k<=l;k =k+1)1107 for (k=1;k<=l;k++) 1108 1108 { if (unique(G(1..g),P(k))) 1109 { j =j+1;// a new factor for next run1110 g =g+1;1109 { j++; // a new factor for next run 1110 g++; 1111 1111 matrix G(g)=P(k); // a new group element - 1112 1112 A(1)=concat(A(1),P(k)*vars); // adding new mapping to A(1) … … 1196 1196 // operator 1197 1197 string chst=charstr(br); 1198 for (int i=1;i<=size(chst);i =i+1)1198 for (int i=1;i<=size(chst);i++) 1199 1199 { if (chst[i]==",") 1200 1200 { break; … … 1226 1226 string links, rechts; 1227 1227 string stM(1)=string(#[1]); 1228 for (o=1;o<=size(stM(1));o =o+1)1228 for (o=1;o<=size(stM(1));o++) 1229 1229 { if (stM(1)[o]==" 1230 1230 ") … … 1244 1244 // new term of the Molien series 1245 1245 i=1; 1246 for (int j=2;j<=gen_num;j =j+1) // this loop adds the parameters to the1246 for (int j=2;j<=gen_num;j++) // this loop adds the parameters to the 1247 1247 { // group, leaving out doubles and 1248 1248 // checking whether the parameters are … … 1259 1259 } 1260 1260 if (unique(G(1..i),#[j])) 1261 { i =i+1;1261 { i++; 1262 1262 matrix G(i)=#[j]; 1263 1263 A(1)=concat(A(1),G(i)*vars); // adding ring homomorphisms to A(1) 1264 1264 string stM(i)=string(G(i)); 1265 for (o=1;o<=size(stM(i));o =o+1)1265 for (o=1;o<=size(stM(i));o++) 1266 1266 { if (stM(i)[o]==" 1267 1267 ") … … 1308 1308 { l=0; // l is the number of products we get in 1309 1309 // one going 1310 for (m=g-j+1;m<=g;m =m+1)1311 { for (k=1;k<=i;k =k+1)1312 { l =l+1;1310 for (m=g-j+1;m<=g;m++) 1311 { for (k=1;k<=i;k++) 1312 { l++; 1313 1313 matrix P(l)=G(k)*G(m); // possible new element 1314 1314 } 1315 1315 } 1316 1316 j=0; 1317 for (k=1;k<=l;k =k+1)1317 for (k=1;k<=l;k++) 1318 1318 { if (unique(G(1..g),P(k))) 1319 { j =j+1;// a new factor for next run1320 g =g+1;1319 { j++; // a new factor for next run 1320 g++; 1321 1321 matrix G(g)=P(k); // a new group element - 1322 1322 A(1)=concat(A(1),G(g)*vars); // adding new mapping to A(1) 1323 1323 string stM(g)=string(G(g)); 1324 for (o=1;o<=size(stM(g));o =o+1)1324 for (o=1;o<=size(stM(g));o++) 1325 1325 { if (stM(g)[o]==" 1326 1326 ") … … 1433 1433 (first n if there is no third parameter given, otherwise the next n 1434 1434 terms depending on a previous calculation) and an intermediate result 1435 (type <poly>) of the calculation to be used as third parameter in a 1435 (type <poly>) of the calculation to be used as third parameter in a 1436 1436 next run of partial_molien 1437 1437 THEORY: The following calculation is implemented: … … 1480 1480 poly A(2)=fetch(br,A(2)); // fetching A(2) and M[1,2] into R 1481 1481 poly den=fetch(br,A(1)); 1482 for (int i=1; i<=n; i =i+1) // getting n terms and adding them up1482 for (int i=1; i<=n; i++) // getting n terms and adding them up 1483 1483 { min=lead(A(2)); 1484 1484 A(1)=A(1)+min; … … 1534 1534 map pREY; 1535 1535 matrix rowREY[1][n]; 1536 for (int i=1;i<=m;i =i+1)1536 for (int i=1;i<=m;i++) 1537 1537 { rowREY=REY[i,1..n]; 1538 1538 pREY=br,ideal(rowREY); // f is now the i-th ring homomorphism … … 1564 1564 shoud be, G1,G2,...: <matrices> generating a finite matrix group 1565 1565 RETURNS: the basis (type <ideal>) of the space of invariants of degree g 1566 THEORY: A general polynomial of degree g is generated and the generators of 1566 THEORY: A general polynomial of degree g is generated and the generators of 1567 1567 the matrix group applied. The difference ought to be 0 and this way a 1568 1568 system of linear equations is created. It is solved by computing … … 1582 1582 int n=nvars(br); 1583 1583 //---------------------- checking that the input is ok ----------------------- 1584 for (i=1;i<=a;i =i+1)1584 for (i=1;i<=a;i++) 1585 1585 { if (typeof(#[i])=="matrix") 1586 1586 { if (nrows(#[i])==n && ncols(#[i])==n) … … 1615 1615 int j,k; 1616 1616 //------------------- building the system of linear equations ---------------- 1617 for (i=1;i<=a;i =i+1)1617 for (i=1;i<=a;i++) 1618 1618 { I=ideal(matrix(vars)*transpose(imap(br,G(i)))); 1619 1619 I=I,p(1..m); 1620 1620 f=T,I; 1621 1621 Pnew=f(P); 1622 for (j=1;j<=m;j =j+1)1622 for (j=1;j<=m;j++) 1623 1623 { temp=P/mon[j]-Pnew/mon[j]; 1624 for (k=1;k<=m;k =k+1)1624 for (k=1;k<=m;k++) 1625 1625 { S[m*(i-1)+j,k]=temp/p(k); 1626 1626 } … … 1635 1635 ideal I=ideal(matrix(mon)*s); // I contains a basis of homogeneous 1636 1636 if (I[1]<>0) // invariants of degree d 1637 { for (i=1;i<=ncols(I);i =i+1)1637 { for (i=1;i<=ncols(I);i++) 1638 1638 { I[i]=I[i]/leadcoef(I[i]); // setting leading coefficients to 1 1639 1639 } … … 1758 1758 { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; echo=2; 1759 1759 ring R=0,(x,y,z),dp; 1760 matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; 1760 matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; 1761 1761 intvec flags=0,1,0; 1762 1762 matrix REY,M=reynolds_molien(A,flags); … … 1781 1781 if (step_fac>j) 1782 1782 { B=compress(evaluate_reynolds(REY,mon)); 1783 for (i=1;i<=ncols(B);i =i+1)// those are taken our that are o mod sP1783 for (i=1;i<=ncols(B);i++) // those are taken our that are o mod sP 1784 1784 { if (reduce(B[i],sP)==0) 1785 1785 { B[i]=0; … … 1801 1801 part_mon=mon[lower_bound..upper_bound]; 1802 1802 part_mon=compress(evaluate_reynolds(REY,part_mon)); 1803 for (i=1;i<=ncols(part_mon);i =i+1)1803 for (i=1;i<=ncols(part_mon);i++) 1804 1804 { if (reduce(part_mon[i],sP)==0) 1805 1805 { part_mon[i]=0; … … 1821 1821 proc next_vector(intmat vec) 1822 1822 { int n=ncols(vec); // p: >0, n: <0, p0: >=0, n0: <=0 1823 for (int i=1;i<=n;i =i+1)// finding out which is the first1823 for (int i=1;i<=n;i++) // finding out which is the first 1824 1824 { if (vec[1,i]<>0) // component <>0 1825 1825 { break; … … 1853 1853 return(new); 1854 1854 } // case left: 1,*,...,* 1855 for (i=2;i<=n;i =i+1)1855 for (i=2;i<=n;i++) 1856 1856 { if (vec[1,i]>0) // make first positive negative and 1857 1857 { vec[1,i]=-vec[1,i]; // return … … 1862 1862 } // positive 1863 1863 } 1864 for (i=2;i<=n-1;i =i+1)// case: 1,p,...,p after 1,n,...,n1864 for (i=2;i<=n-1;i++) // case: 1,p,...,p after 1,n,...,n 1865 1865 { if (vec[1,i]>0) 1866 1866 { vec[1,2]=vec[1,i]-1; // shuffleing things around... … … 1878 1878 /////////////////////////////////////////////////////////////////////////////// 1879 1879 // Maps integers to elements of the base field. It is only called if the base 1880 // field is of prime characteristic. If the base field has q elements 1880 // field is of prime characteristic. If the base field has q elements 1881 1881 // (depending on minpoly) 1..q is mapped to those q elements. 1882 1882 /////////////////////////////////////////////////////////////////////////////// … … 1899 1899 while (1) 1900 1900 { if (i<p^j) // finding an upper bound on i 1901 { for (k=0;k<j-1;k =k+1)1901 { for (k=0;k<j-1;k++) 1902 1902 { out=out+((i/p^k)%p)*a^k; // finding how often p^k is contained in 1903 1903 } // i … … 1908 1908 return(out); 1909 1909 } 1910 j =j+1;1910 j++; 1911 1911 } 1912 1912 } … … 1956 1956 } 1957 1957 P[j]=test_poly; // test_poly ist added to primary 1958 j =j+1;// invariants1958 j++; // invariants 1959 1959 } 1960 1960 return(P,sP,CI,dB); … … 2005 2005 { setring br; 2006 2006 new=next_vector(new); 2007 for (j=1;j<=cd;j =j+1)2007 for (j=1;j<=cd;j++) 2008 2008 { pnew[1,j]=int_number_map(new[1,j]); // mapping an integer into br 2009 2009 } 2010 2010 if (unique(vec(1..counter),pnew)) // checking whether we tried pnew before 2011 { counter =counter+1;2011 { counter++; 2012 2012 matrix vec(counter)=pnew; // keeping track of the ones we tried - 2013 2013 test_poly=(vec(counter)*tB)[1,1]; // linear combination - … … 2045 2045 return(P,sP,CI,dB-d+1); 2046 2046 } 2047 k =k+1;2047 k++; 2048 2048 } 2049 2049 return(P,sP,CI,dB); … … 2158 2158 } // i.e. we can take all of B 2159 2159 else 2160 { for(j=i+1;j>i+dif;j =j+1)2160 { for(j=i+1;j>i+dif;j++) 2161 2161 { CI=CI+ideal(var(j)^d); 2162 2162 } … … 2166 2166 } 2167 2167 if (v) 2168 { for (j=1;j<=dif;j =j+1)2168 { for (j=1;j<=dif;j++) 2169 2169 { " We find: "+string(P[i+j]); 2170 2170 } … … 2307 2307 } 2308 2308 else // i.e. we can take all of B 2309 { for(j=i+1;j>i+dif;j =j+1)2309 { for(j=i+1;j>i+dif;j++) 2310 2310 { CI=CI+ideal(var(j)^d); 2311 2311 } … … 2315 2315 } 2316 2316 if (v) 2317 { for (j=1;j<=size(P)-i;j =j+1)2317 { for (j=1;j<=size(P)-i;j++) 2318 2318 { " We find: "+string(P[i+j]); 2319 2319 } … … 2422 2422 while(1) // repeat until n primary invariants are 2423 2423 { // found - 2424 d =d+1;// degree where we'll search2424 d++; // degree where we'll search 2425 2425 if (v) 2426 2426 { " Computing primary invariants in degree "+string(d)+":"; … … 2444 2444 } 2445 2445 else // i.e. we can take all of B 2446 { for(j=i+1;j<=i+dif;j =j+1)2446 { for(j=i+1;j<=i+dif;j++) 2447 2447 { CI=CI+ideal(var(j)^d); 2448 2448 } … … 2452 2452 } 2453 2453 if (v) 2454 { for (j=1;j<=dif;j =j+1)2454 { for (j=1;j<=dif;j++) 2455 2455 { " We find: "+string(P[i+j]); 2456 2456 } … … 2565 2565 while(1) // repeat until n primary invariants are 2566 2566 { // found - 2567 d =d+1;// degree where we'll search2567 d++; // degree where we'll search 2568 2568 if (v) 2569 2569 { " Computing primary invariants in degree "+string(d)+":"; … … 2588 2588 } 2589 2589 else // i.e. we can take all of B 2590 { for(j=i+1;j<=i+dif;j =j+1)2590 { for(j=i+1;j<=i+dif;j++) 2591 2591 { CI=CI+ideal(var(j)^d); 2592 2592 } … … 2596 2596 } 2597 2597 if (v) 2598 { for (j=1;j<=size(P)-i;j =j+1)2598 { for (j=1;j<=size(P)-i;j++) 2599 2599 { " We find: "+string(P[i+j]); 2600 2600 } … … 2674 2674 // as the number of primary invariants, 2675 2675 // we should get 2676 for (int i=1;i<=gen_num;i =i+1)2676 for (int i=1;i<=gen_num;i++) 2677 2677 { if (typeof(#[i])=="matrix") 2678 2678 { if (nrows(#[i])<>n or ncols(#[i])<>n) … … 2712 2712 while(1) // repeat until n primary invariants are 2713 2713 { // found - 2714 d =d+1;// degree where we'll search2714 d++; // degree where we'll search 2715 2715 if (v) 2716 2716 { " Computing primary invariants in degree "+string(d)+":"; … … 2734 2734 } 2735 2735 else // i.e. we can take all of B 2736 { for(j=i+1;j<=i+dif;j =j+1)2736 { for(j=i+1;j<=i+dif;j++) 2737 2737 { CI=CI+ideal(var(j)^d); 2738 2738 } … … 2742 2742 } 2743 2743 if (v) 2744 { for (j=1;j<=size(P)-i;j =j+1)2744 { for (j=1;j<=size(P)-i;j++) 2745 2745 { " We find: "+string(P[i+j]); 2746 2746 } … … 2841 2841 v=0; 2842 2842 } 2843 for (int i=1;i<=gen_num;i =i+1)2843 for (int i=1;i<=gen_num;i++) 2844 2844 { if (typeof(#[i])=="matrix") 2845 2845 { if (nrows(#[i])<>n or ncols(#[i])<>n) … … 2981 2981 matrix test_matrix[1][dif]; // the linear combination to test 2982 2982 intvec h; // Hilbert series 2983 for (j=i+1;j<=i+dif;j =j+1)2983 for (j=i+1;j<=i+dif;j++) 2984 2984 { CI=CI+ideal(var(j)^d); // homogeneous polynomial of the same 2985 2985 // degree as the one we're looking for … … 3035 3035 { " Give a new <int> > "+string(max)+" that bounds the range of coefficients:"; 3036 3036 answer=read(""); 3037 for (j=1;j<=size(answer)-1;j =j+1)3038 { for (k=0;k<=9;k =k+1)3037 for (j=1;j<=size(answer)-1;j++) 3038 { for (k=0;k<=9;k++) 3039 3039 { if (answer[j]==string(k)) 3040 3040 { break; … … 3079 3079 ideal TEST; 3080 3080 while (dif>0) 3081 { for (j=i+1;j<=i+dif;j =j+1)3081 { for (j=i+1;j<=i+dif;j++) 3082 3082 { CI=CI+ideal(var(j)^d); // homogeneous polynomial of the same 3083 3083 // degree as the one we're looking for … … 3144 3144 { " Give a new <int> > "+string(max)+" that bounds the range of coefficients:"; 3145 3145 answer=read(""); 3146 for (j=1;j<=size(answer)-1;j =j+1)3147 { for (k=0;k<=9;k =k+1)3146 for (j=1;j<=size(answer)-1;j++) 3147 { for (k=0;k<=9;k++) 3148 3148 { if (answer[j]==string(k)) 3149 3149 { break; … … 3291 3291 } // dif invariants - 3292 3292 else // i.e. we can take all of B 3293 { for(j=i+1;j>i+dif;j =j+1)3293 { for(j=i+1;j>i+dif;j++) 3294 3294 { CI=CI+ideal(var(j)^d); 3295 3295 } … … 3303 3303 } 3304 3304 if (v) 3305 { for (j=1;j<=dif;j =j+1)3305 { for (j=1;j<=dif;j++) 3306 3306 { " We find: "+string(P[i+j]); 3307 3307 } … … 3445 3445 } 3446 3446 else // i.e. we can take all of B 3447 { for(j=i+1;j>i+dif;j =j+1)3447 { for(j=i+1;j>i+dif;j++) 3448 3448 { CI=CI+ideal(var(j)^d); 3449 3449 } … … 3457 3457 } 3458 3458 if (v) 3459 { for (j=1;j<=size(P)-i;j =j+1)3459 { for (j=1;j<=size(P)-i;j++) 3460 3460 { " We find: "+string(P[i+j]); 3461 3461 } … … 3565 3565 while(1) // repeat until n primary invariants are 3566 3566 { // found - 3567 d =d+1;// degree where we'll search3567 d++; // degree where we'll search 3568 3568 if (v) 3569 3569 { " Computing primary invariants in degree "+string(d)+":"; … … 3587 3587 } 3588 3588 else // i.e. we can take all of B 3589 { for(j=i+1;j<=i+dif;j =j+1)3589 { for(j=i+1;j<=i+dif;j++) 3590 3590 { CI=CI+ideal(var(j)^d); 3591 3591 } … … 3599 3599 } 3600 3600 if (v) 3601 { for (j=1;j<=dif;j =j+1)3601 { for (j=1;j<=dif;j++) 3602 3602 { " We find: "+string(P[i+j]); 3603 3603 } … … 3713 3713 while(1) // repeat until n primary invariants are 3714 3714 { // found - 3715 d =d+1;// degree where we'll search3715 d++; // degree where we'll search 3716 3716 if (v) 3717 3717 { " Computing primary invariants in degree "+string(d)+":"; … … 3735 3735 } 3736 3736 else // i.e. we can take all of B 3737 { for(j=i+1;j<=i+dif;j =j+1)3737 { for(j=i+1;j<=i+dif;j++) 3738 3738 { CI=CI+ideal(var(j)^d); 3739 3739 } … … 3747 3747 } 3748 3748 if (v) 3749 { for (j=1;j<=size(P)-i;j =j+1)3749 { for (j=1;j<=size(P)-i;j++) 3750 3750 { " We find: "+string(P[i+j]); 3751 3751 } … … 3834 3834 // as the number of primary invariants, 3835 3835 // we should get 3836 for (int i=1;i<=gen_num;i =i+1)3836 for (int i=1;i<=gen_num;i++) 3837 3837 { if (typeof(#[i])=="matrix") 3838 3838 { if (nrows(#[i])<>n or ncols(#[i])<>n) … … 3871 3871 while(1) // repeat until n primary invariants are 3872 3872 { // found - 3873 d =d+1;// degree where we'll search3873 d++; // degree where we'll search 3874 3874 if (v) 3875 3875 { " Computing primary invariants in degree "+string(d)+":"; … … 3892 3892 } 3893 3893 else // i.e. we can take all of B 3894 { for(j=i+1;j<=i+dif;j =j+1)3894 { for(j=i+1;j<=i+dif;j++) 3895 3895 { CI=CI+ideal(var(j)^d); 3896 3896 } … … 3904 3904 } 3905 3905 if (v) 3906 { for (j=1;j<=size(P)-i;j =j+1)3906 { for (j=1;j<=size(P)-i;j++) 3907 3907 { " We find: "+string(P[i+j]); 3908 3908 } … … 4015 4015 } 4016 4016 } 4017 for (int i=1;i<=gen_num;i =i+1)4017 for (int i=1;i<=gen_num;i++) 4018 4018 { if (typeof(#[i])=="matrix") 4019 4019 { if (nrows(#[i])<>n or ncols(#[i])<>n) … … 4168 4168 intmat PP[s][1]; 4169 4169 intmat TEST[s][1]; 4170 for (int i=1;i<=s;i =i+1)4170 for (int i=1;i<=s;i++) 4171 4171 { if (i<0) 4172 4172 { "ERROR: The entries of <intvec> may not be <= 0"; … … 4180 4180 intmat PPd_neu_gross[s][nc]; 4181 4181 PPd_neu_gross[i..s,1..nc]=PPd_neu[1..s-i+1,1..nc]; 4182 for (j=1;j<=nc;j =j+1)4182 for (j=1;j<=nc;j++) 4183 4183 { PPd_neu_gross[i,j]=PPd_neu_gross[i,j]+1; 4184 4184 } … … 4221 4221 irreducible secondary invariants (type <matrix>) 4222 4222 DISPLAY: information if v does not equal 0 4223 THEORY: The secondary invariants are calculated by finding a basis (in terms 4224 of monomials) of the basering modulo the primary invariants, mapping 4223 THEORY: The secondary invariants are calculated by finding a basis (in terms 4224 of monomials) of the basering modulo the primary invariants, mapping 4225 4225 those to invariants with the Reynolds operator and using these images 4226 4226 or their power products s. t. they are linearly independent modulo … … 4271 4271 //- finding the polynomial giving number and degrees of secondary invariants - 4272 4272 poly p=1; 4273 for (j=1;j<=n;j =j+1)// calculating the denominator of the4273 for (j=1;j<=n;j++) // calculating the denominator of the 4274 4274 { p=p*(1-var(1)^deg(P[j])); // Hilbert series of the ring generated 4275 4275 } // by the primary invariants - … … 4302 4302 ideal S=1; // 1 is the first secondary invariant - 4303 4303 //--------------------- generating secondary invariants ---------------------- 4304 for (i=2;i<=m;i =i+1)// going through dimmat -4304 for (i=2;i<=m;i++) // going through dimmat - 4305 4305 { if (int(dimmat[i,1])<>0) // when it is == 0 we need to find 0 4306 4306 { // elements in the current degree (i-1) … … 4314 4314 } // secondary invariants 4315 4315 if (size(ideal(PP))<>0) 4316 { for (j=1;j<=ncols(PP);j =j+1)// going through all the power products4316 { for (j=1;j<=ncols(PP);j++) // going through all the power products 4317 4317 { pp=1; 4318 for (k=1;k<=nrows(PP);k =k+1)4318 for (k=1;k<=nrows(PP);k++) 4319 4319 { pp=pp*IS[1,k]^PP[k,j]; 4320 4320 } 4321 4321 if (reduce(pp,TEST)<>0) 4322 4322 { S=S,pp; 4323 counter =counter+1;4323 counter++; 4324 4324 if (v) 4325 4325 { " We find: "+string(pp); … … 4350 4350 { " We find: "+string(B[1]); 4351 4351 } 4352 for (j=2;j<=int(dimmat[i,1]);j =j+1)4352 for (j=2;j<=int(dimmat[i,1]);j++) 4353 4353 { deg_vec=deg_vec,i-1; 4354 4354 if (v) … … 4358 4358 } 4359 4359 else 4360 { for (j=1;j<=int(dimmat[i,1]);j =j+1)4360 { for (j=1;j<=int(dimmat[i,1]);j++) 4361 4361 { deg_vec=deg_vec,i-1; 4362 4362 if (v) … … 4371 4371 { // invariants that are linearly 4372 4372 // independent modulo TEST 4373 j =j+1;4373 j++; 4374 4374 if (reduce(B[j],TEST)<>0) // B[j] should be added 4375 4375 { S=S,B[j]; … … 4381 4381 { deg_vec=deg_vec,i-1; 4382 4382 } 4383 counter =counter+1;4383 counter++; 4384 4384 if (v) 4385 4385 { " We find: "+string(B[j]); … … 4481 4481 intvec deg_dim_vec; 4482 4482 //- finding the polynomial giving number and degrees of secondary invariants - 4483 for (j=1;j<=n;j =j+1)4483 for (j=1;j<=n;j++) 4484 4484 { deg_dim_vec[j]=deg(P[j]); 4485 4485 } 4486 4486 setring `ring_name`; 4487 4487 poly p=1; 4488 for (j=1;j<=n;j =j+1)// calculating the denominator of the4488 for (j=1;j<=n;j++) // calculating the denominator of the 4489 4489 { p=p*(1-var(1)^deg_dim_vec[j]); // Hilbert series of the ring generated 4490 4490 } // by the primary invariants - … … 4505 4505 m=nrows(dimmat); // m-1 is the highest degree 4506 4506 deg_dim_vec=1; 4507 for (j=2;j<=m;j =j+1)4507 for (j=2;j<=m;j++) 4508 4508 { deg_dim_vec=deg_dim_vec,int(dimmat[j,1]); 4509 4509 } … … 4522 4522 ideal S=1; // 1 is the first secondary invariant 4523 4523 //------------------- generating secondary invariants ------------------------ 4524 for (i=2;i<=m;i =i+1)// going through deg_dim_vec -4524 for (i=2;i<=m;i++) // going through deg_dim_vec - 4525 4525 { if (deg_dim_vec[i]<>0) // when it is == 0 we need to find 0 4526 4526 { // elements in the current degree (i-1) … … 4534 4534 } // irreducible secondary invariants 4535 4535 if (size(ideal(PP))<>0) 4536 { for (j=1;j<=ncols(PP);j =j+1)// going through all of those4536 { for (j=1;j<=ncols(PP);j++) // going through all of those 4537 4537 { pp=1; 4538 for (k=1;k<=nrows(PP);k =k+1)4538 for (k=1;k<=nrows(PP);k++) 4539 4539 { pp=pp*IS[1,k]^PP[k,j]; 4540 4540 } 4541 4541 if (reduce(pp,TEST)<>0) 4542 4542 { S=S,pp; 4543 counter =counter+1;4543 counter++; 4544 4544 if (v) 4545 4545 { " We find: "+string(pp); … … 4570 4570 { " We find: "+string(B[1]); 4571 4571 } 4572 for (j=2;j<=deg_dim_vec[i];j =j+1)4572 for (j=2;j<=deg_dim_vec[i];j++) 4573 4573 { deg_vec=deg_vec,i-1; 4574 4574 if (v) … … 4578 4578 } 4579 4579 else 4580 { for (j=1;j<=deg_dim_vec[i];j =j+1)4580 { for (j=1;j<=deg_dim_vec[i];j++) 4581 4581 { deg_vec=deg_vec,i-1; 4582 4582 if (v) … … 4591 4591 { // invariants that are linearly 4592 4592 // independent modulo TEST 4593 j =j+1;4593 j++; 4594 4594 if (reduce(B[j],TEST)<>0) // B[j] should be added 4595 4595 { S=S,B[j]; … … 4601 4601 { deg_vec=deg_vec,i-1; 4602 4602 } 4603 counter =counter+1;4603 counter++; 4604 4604 if (v) 4605 4605 { " We find: "+string(B[j]); … … 4643 4643 P: a 1xn <matrix> with primary invariants, REY: a gxn <matrix> 4644 4644 representing the Reynolds operator, deg_vec: an optional <intvec> 4645 listing some degrees where no non-trivial homogeneous invariants can 4645 listing some degrees where no non-trivial homogeneous invariants can 4646 4646 be found, v: an optional <int> 4647 4647 ASSUME: n is the number of variables of the basering, g the size of the group, … … 4716 4716 int j, m, d; 4717 4717 int max=1; 4718 for (j=1;j<=n;j =j+1)4718 for (j=1;j<=n;j++) 4719 4719 { max=max*deg(P[j]); 4720 4720 } … … 4738 4738 //--------------------- generating secondary invariants ---------------------- 4739 4739 while (counter<>max) 4740 { i =i+1;4740 { i++; 4741 4741 if (deg_vec[k]<>i) 4742 4742 { if (v) … … 4748 4748 // and that don't reduce to 0 modulo sP 4749 4749 TEST=sP; 4750 for (j=1;j<=ncols(B);j =j+1)4750 for (j=1;j<=ncols(B);j++) 4751 4751 { // that are linearly independent modulo 4752 4752 // TEST 4753 4753 if (reduce(B[j],TEST)<>0) // B[j] should be added 4754 4754 { S=S,B[j]; 4755 counter =counter+1;4755 counter++; 4756 4756 if (v) 4757 4757 { " We find: "+string(B[j]); … … 4775 4775 } 4776 4776 else 4777 { k =k+1;4777 { k++; 4778 4778 } 4779 4779 } … … 4874 4874 int j, m, d; 4875 4875 int max=1; 4876 for (j=1;j<=n;j =j+1)4876 for (j=1;j<=n;j++) 4877 4877 { max=max*deg(P[j]); 4878 4878 } … … 4900 4900 //------------------- generating secondary invariants ------------------------ 4901 4901 while (counter<>max) 4902 { i =i+1;4902 { i++; 4903 4903 if (deg_vec[l]<>i) 4904 4904 { if (v) … … 4911 4911 // invariants 4912 4912 if (size(ideal(PP))<>0) 4913 { for (j=1;j<=ncols(PP);j =j+1)// going through all those power products4913 { for (j=1;j<=ncols(PP);j++) // going through all those power products 4914 4914 { pp=1; 4915 for (k=1;k<=nrows(PP);k =k+1)4915 for (k=1;k<=nrows(PP);k++) 4916 4916 { pp=pp*IS[1,k]^PP[k,j]; 4917 4917 } 4918 4918 if (reduce(pp,TEST)<>0) 4919 4919 { S=S,pp; 4920 counter =counter+1;4920 counter++; 4921 4921 if (v) 4922 4922 { " We find: "+string(pp); … … 4938 4938 // operator that are linearly independent 4939 4939 // and that don't reduce to 0 modulo sP 4940 for (j=1;j<=ncols(B);j =j+1)4940 for (j=1;j<=ncols(B);j++) 4941 4941 { if (reduce(B[j],TEST)<>0) // B[j] should be added 4942 4942 { S=S,B[j]; … … 4948 4948 { irreducible_deg_vec=irreducible_deg_vec,i; 4949 4949 } 4950 counter =counter+1;4950 counter++; 4951 4951 if (v) 4952 4952 { " We find: "+string(B[j]); … … 4971 4971 } 4972 4972 else 4973 { l =l+1;4973 { l++; 4974 4974 } 4975 4975 } … … 5022 5022 } 5023 5023 int v=#[size(#)]; 5024 for (i=1;i<=gen_num;i =i+1)5024 for (i=1;i<=gen_num;i++) 5025 5025 { if (typeof(#[i])<>"matrix") 5026 5026 { "ERROR: These parameters should be generators of the finite matrix group."; … … 5036 5036 { int v=0; 5037 5037 int gen_num=size(#); 5038 for (i=1;i<=gen_num;i =i+1)5038 for (i=1;i<=gen_num;i++) 5039 5039 { if (typeof(#[i])<>"matrix") 5040 5040 { "ERROR: These parameters should be generators of the finite matrix group."; … … 5095 5095 matrix M(1)[gen_num][k]; // M(1) will contain a module 5096 5096 ideal B; 5097 for (i=1;i<=gen_num;i =i+1)5097 for (i=1;i<=gen_num;i++) 5098 5098 { B=ideal(matrix(maxideal(1))*transpose(#[i])); // image of the various 5099 5099 // variables under the i-th generator - … … 5116 5116 M[1..k,1..m]=matrix(f(M(2))); // will contain a module - 5117 5117 matrix P=f(P); // primary invariants in the new ring 5118 for (i=1;i<=n;i =i+1)5119 { for (j=1;j<=k;j =j+1)5118 for (i=1;i<=n;i++) 5119 { for (j=1;j<=k;j++) 5120 5120 { M[j,m+(i-1)*k+j]=y(i)-P[1,i]; 5121 5121 } … … 5134 5134 S=matrix(trivialS)*matrix(M(2)); // S now contains the secondary 5135 5135 // invariants 5136 for (i=1; i<=m;i =i+1)5136 for (i=1; i<=m;i++) 5137 5137 { S[1,i]=S[1,i]/leadcoef(S[1,i]); // making elements nice 5138 5138 } … … 5140 5140 if (v) 5141 5141 { " These are the secondary invariants: "; 5142 for (i=1;i<=m;i =i+1)5142 for (i=1;i<=m;i++) 5143 5143 { " "+string(S[1,i]); 5144 5144 } … … 5386 5386 (i.e. we will not even attempt to compute the Reynolds operator or 5387 5387 Molien series), the second component should give the size of intervals 5388 between canceling common factors in the expansion of the Molien 5388 between canceling common factors in the expansion of the Molien 5389 5389 series, 5390 5390 0 (the default) means only once after generating all terms, in prime … … 5450 5450 } 5451 5451 } 5452 for (int i=1;i<=gen_num;i =i+1)5452 for (int i=1;i<=gen_num;i++) 5453 5453 { if (typeof(#[i])=="matrix") 5454 5454 { if (nrows(#[i])<>n or ncols(#[i])<>n) … … 5619 5619 execute "minpoly=number("+mp+");"; 5620 5620 ideal I=ideal(imap(br,F)); 5621 for (int i=1;i<=m;i =i+1)5621 for (int i=1;i<=m;i++) 5622 5622 { I[i]=I[i]-y(i); 5623 5623 } … … 5626 5626 execute "minpoly=number("+mp+");"; 5627 5627 ideal vars; 5628 for (i=2;i<=n;i =i+1)5628 for (i=2;i<=n;i++) 5629 5629 { vars[i]=0; 5630 5630 } … … 5662 5662 THEORY: A Groebner basis of the ideal of algebraic relations of the invariant 5663 5663 ring generators is calculated, then one of the basis elements plus the 5664 ideal generators. The variables of the original ring are eliminated 5664 ideal generators. The variables of the original ring are eliminated 5665 5665 and the polynomials that are left define the relative orbit variety 5666 5666 with respect to I. … … 5681 5681 ideal J=ideal(imap(br,F)); 5682 5682 ideal I=imap(br,I); 5683 for (int i=1;i<=m;i =i+1)5683 for (int i=1;i<=m;i++) 5684 5684 { J[i]=J[i]-y(i); 5685 5685 } … … 5689 5689 J=std(J); 5690 5690 ideal vars; 5691 //for (i=1;i<=n;i =i+1)5691 //for (i=1;i<=n;i++) 5692 5692 //{ vars[i]=0; 5693 5693 //} … … 5697 5697 ideal G=emb(J); 5698 5698 J=J-G; 5699 for (i=1;i<=ncols(G);i =i+1)5699 for (i=1;i<=ncols(G);i++) 5700 5700 { if (J[i]<>0) 5701 5701 { G[i]=0; … … 5706 5706 execute "minpoly=number("+mp+");"; 5707 5707 ideal vars; 5708 for (i=2;i<=n;i =i+1)5708 for (i=2;i<=n;i++) 5709 5709 { vars[i]=0; 5710 5710 } … … 5753 5753 execute "ring R=("+charstr(br)+"),("+varstr(br)+","+varstr(E)+"),lp;"; 5754 5754 ideal F=imap(br,F); 5755 for (int i=1;i<=n;i =i+1)5755 for (int i=1;i<=n;i++) 5756 5756 { F=0,F; 5757 5757 }
Note: See TracChangeset
for help on using the changeset viewer.
