Changeset 67bd4c in git for Singular/LIB/primdec.lib
 Timestamp:
 Nov 5, 1997, 7:16:57 PM (26 years ago)
 Branches:
 (u'spielwiese', '2fa36c576e6a4ddbb1093b43c7f8e9835e17e52a')
 Children:
 8e8aca339d1600a8eb489af59da08ac141f0477d
 Parents:
 67a4495135101f4651e8ccc9b7618eb02713c97d
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/primdec.lib
r67a449 r67bd4c 1 // $Id: primdec.lib,v 1. 5 19970918 09:58:26Singular Exp $1 // $Id: primdec.lib,v 1.6 19971105 18:16:57 Singular Exp $ 2 2 /////////////////////////////////////////////////////// 3 3 // primdec.lib … … 7 7 ////////////////////////////////////////////////////// 8 8 9 LIBRARY: primdec.lib: PROCEDURE FOR PRIMARY DECOMPOSITION (G/T/Z) 10 11 minAssPrimes (ideal I, list choose) 12 // minimal associated primes 13 // The list choose must be either emty (minAssPrimes(I)) or 1 14 // (minAssPrimes(I,1)) 15 // In the second case the factorizing Buchberger Algorithm is used 16 // which in most cases may considerably speed up the algorithm 17 18 primdec (ideal I) 19 20 // Computes a complete primary decomposition via 9 LIBRARY: primdec.lib: PROCEDURE FOR PRIMARY DECOMPOSITION (I) 10 11 minAssPrimes (ideal I) 12 //computes the minimal associated primes of the ideal I 13 14 decomp (ideal I) 15 // Computes a complete primary decomposition via Gianni,Trager,Zacharias 21 16 22 17 radical(ideal I) 23 18 //computes the radical of the ideal I 24 19 20 equiRadical(ideal I) 21 //computes the radical of the equidimensional part of the ideal I 22 23 prepareAss(ideal I) 24 //computes the radicals of the equidimensional parts of I 25 25 26 LIB "random.lib"; 26 LIB " poly.lib";27 LIB "elim.lib"; 27 28 /////////////////////////////////////////////////////////////////////////////// 28 29 … … 399 400 int i; 400 401 ideal j=pr[1]; 401 for (i=2;i<=size(pr) div2;i++)402 for (i=2;i<=size(pr)/2;i++) 402 403 { 403 404 j=intersect(j,pr[2*i1]); … … 438 439 } 439 440 int k; 440 for (k=1;k<=size(l) div2;k=k+1)441 for (k=1;k<=size(l)/2;k=k+1) 441 442 { 442 443 " "; … … 485 486 486 487 proc primaryTest (ideal i, poly p) 487 USAGE: primaryTest(i,p); i ideal p poly488 RETURN: ideal = radical of i, if i is primary in general position,489 zerodimensional and the radical of i intersected with K[z]490 is (p), z the smallest variable with respect to the lexico491 graphical ordering, and 0 else492 NOTE: It is necessary that i is a standardbasis with respect to493 the lexicographical ordering and the first element in i is494 a power of p.495 EXAMPLE: example primaryTest; shows an example496 488 { 497 489 int m=1; … … 501 493 ideal h; 502 494 503 //the first generator of the prim ideal for the result504 495 ideal prm=p; 505 496 attrib(prm,"isSB",1); … … 527 518 //if not (0) is returned, else var(n)+h is added to prm 528 519 529 e=deg(lead(i[m])); 530 t=leadcoef(i[m])*e*var(n)+(i[m]lead(i[m]))/var(n)^(e1); 531 i[m]=poly(e)^e*leadcoef(i[m])^(e1)*i[m]; 532 533 if (reduce(i[m]t^e,prm) !=0) 534 { 535 return(ideal(0)); 536 } 537 h=interred(t); 538 t=h[1]; 520 e=deg(lead(i[m])); 521 // t=hilfe1(i,prm,m,n); 522 t=leadcoef(i[m])*e*var(n)+(i[m]lead(i[m]))/var(n)^(e1); 523 524 i[m]=poly(e)^e*leadcoef(i[m])^(e1)*i[m]; 525 526 if (reduce(i[m]t^e,prm) !=0) 527 { 528 return(ideal(0)); 529 } 530 h=interred(t); 531 t=h[1]; 539 532 540 533 prm = prm,t; … … 543 536 return(prm); 544 537 } 545 example 546 { "EXAMPLE:"; echo=2; 547 ring r = (0,a,b),(x,y,z),lp; 548 poly p = z^2+1; 549 ideal i = p^2,(a*yz^3)^3,(b*xyz+z4)^4; 550 ideal pr= primaryTest(i,p); 551 pr; 552 i = p^2,(yz3)^3,(xyz+z4)^4+1; 553 pr= primaryTest(i,p); 554 pr; 555 ring s=(0,e),(d,c,b,a,y,x,g,f),lp; 556 ideal i=f,g,x4,y,a,b3,c,d; 557 poly p=f; 558 ideal pr= primaryTest(i,p); 559 pr; 560 } 561 538 539 proc hilfe(ideal i,ideal prm,int m) 540 { 541 poly t; 542 int e; 543 544 if(size(i[m])==1) 545 { 546 t=var(n); 547 } 548 else 549 { 550 e=deg(lead(i[m])); 551 552 if(deg(poly(e))>=0) 553 { 554 t=leadcoef(i[m])*e*var(n)+(i[m]lead(i[m]))/var(n)^(e1); 555 i[m]=poly(e)^e*leadcoef(i[m])^(e1)*i[m]; 556 { 557 else 558 { 559 i[m]=i[m]/leadcoef(i[m]); 560 t=reduce(coef(i[m],var(n))[2,e+1],prm); 561 t=var(n)+factorize(t,1)[1]; 562 } 563 } 564 return(t); 565 } 566 proc hilfe1(ideal i,ideal prm,int m,int n) 567 { 568 poly t; 569 int e; 570 if(size(i[m])==1) 571 { 572 t=var(n); 573 } 574 else 575 { 576 e=deg(lead(i[m])); 577 i[m]=i[m]/leadcoef(i[m]); 578 t=reduce(coeffs(i[m],var(n))[1,1],prm); 579 if(size(t)==0){return(var(n));} 580 t=var(n)+factorize(t,1)[1]; 581 } 582 return(t); 583 } 562 584 563 585 /////////////////////////////////////////////////////////////////////////////// … … 569 591 int sl=size(l); 570 592 571 for(i=1;i<=sl div2;i++)593 for(i=1;i<=sl/2;i++) 572 594 { 573 595 if(sact[2][i]>1) … … 579 601 keepprime[i]=l[2*i1]; 580 602 } 581 603 } 582 604 i=0; 583 while(i<size(l) div2)605 while(i<size(l)/2) 584 606 { 585 607 i++; … … 600 622 } 601 623 j=0; 602 if(i<=sl div2)624 if(i<=sl/2) 603 625 { 604 626 j=1; … … 628 650 for(k=2;k<=r;k++) 629 651 { 630 keepprime[size(l) div2+k1]=interred(keepprime[i]+ideal(act[1][k]));652 keepprime[size(l)/2+k1]=interred(keepprime[i]+ideal(act[1][k])); 631 653 } 632 654 keepprime[i]=interred(keepprime[i]+ideal(act[1][1])); … … 655 677 { 656 678 l[s+2*k1]=keepresult[k]; 657 keepprime[s div2+k]=interred(keepresult[k]+ideal(act[1][k]));679 keepprime[s/2+k]=interred(keepresult[k]+ideal(act[1][k])); 658 680 if(vdim(keepresult[k])==deg(act[1][k])) 659 681 { … … 664 686 l[s+2*k]=ideal(0); 665 687 } 666 if((homog(keepresult[k])==1)(homog(keepprime[s div2+k])==1))688 if((homog(keepresult[k])==1)(homog(keepprime[s/2+k])==1)) 667 689 { 668 690 l[s+2*k]=maxideal(1); … … 688 710 l[s+2]=ideal(0); 689 711 } 690 keepprime[s div2+1]=interred(keepprime[i]+ideal(@f));691 if(homog(keepprime[s div2+1])==1)712 keepprime[s/2+1]=interred(keepprime[i]+ideal(@f)); 713 if(homog(keepprime[s/2+1])==1) 692 714 { 693 715 l[s+2]=maxideal(1); … … 712 734 return(l); 713 735 } 714 for(i=1;i<=size(l) div2;i++)736 for(i=1;i<=size(l)/2;i++) 715 737 { 716 738 if((size(l[2*i])==0)&&(specialIdealsEqual(keepprime[i],l[2*i1])!=1)) … … 900 922 901 923 @k=0; 902 while(@k<(size(primary) div 2)) 924 int zz; 925 while(@k<(size(primary)/2)) 903 926 { 904 927 @k++; 905 928 if (size(primary[2*@k])==0) 906 929 { 930 for(zz=1;zz<size(primary[2*@k1])1;zz++) 931 { 932 if(vdim(primary[2*@k1])==deg(primary[2*@k1][zz])) 933 { 934 primary[2*@k]=primary[2*@k1]; 935 } 936 } 937 } 938 } 939 @k=0; 940 while(@k<(size(primary)/2)) 941 { 942 @k++; 943 if (size(primary[2*@k])==0) 944 { 945 // "the univariat polynomials"; 946 // int qwe=timer; 947 // system("finduni",primary[2*@k1]); 948 907 949 jmap=randomLast(100); 950 // timerqwe; 951 // primary[2*@k1]; 952 // pause; 908 953 for(@n=2;@n<=size(primary[2*@k1]);@n++) 909 954 { … … 1002 1047 else 1003 1048 { 1004 for(@n=1;@n<=size(@lr) div2;@n++)1049 for(@n=1;@n<=size(@lr)/2;@n++) 1005 1050 { 1006 1051 if(specialIdealsEqual(@lr[2*@n1],@lr[2*@n])==1) … … 1037 1082 else 1038 1083 { 1039 for(@n=1;@n<=size(@lr) div2;@n++)1084 for(@n=1;@n<=size(@lr)/2;@n++) 1040 1085 { 1041 1086 if(specialIdealsEqual(@lr[2*@n1],@lr[2*@n])==1) … … 1067 1112 primary[2*@k1]=lres[1]; 1068 1113 primary[2*@k]=lres[2]; 1069 @s=size(primary) div2;1070 for(@n=1;@n<=size(lres) div21;@n++)1114 @s=size(primary)/2; 1115 for(@n=1;@n<=size(lres)/21;@n++) 1071 1116 { 1072 1117 primary[2*@s+2*@n1]=lres[2*@n+1]; … … 1111 1156 return(1) 1112 1157 } 1113 p= gcd(p,i[k]);1158 p=GCD(p,i[k]); 1114 1159 if(deg(p)==0) 1115 1160 { … … 1155 1200 if(size(@pa)==0) 1156 1201 { 1157 return(lcm (h));1202 return(lcmP(h)); 1158 1203 } 1159 1204 def bsr= basering; … … 1161 1206 execute "ring @r=0,("+@pa+","+varstr(bsr)+"),dp;"; 1162 1207 execute "ideal @i="+@id+";"; 1163 poly @p=lcm (@i);1208 poly @p=lcmP(@i); 1164 1209 string @ps=string(@p); 1165 1210 setring bsr; … … 1175 1220 ideal l=p,q; 1176 1221 poly pr= coeffLcm(l); 1222 pr; 1223 } 1224 1225 /////////////////////////////////////////////////////////////////////////////// 1226 1227 proc lcmP(ideal i) 1228 USAGE: lcm(i); i list of polynomials 1229 RETURN: poly = lcm(i[1],...,i[size(i)]) 1230 NOTE: 1231 EXAMPLE: example lcm; shows an example 1232 { 1233 int k,j; 1234 poly p,q; 1235 i=simplify(i,10); 1236 for(j=1;j<=size(i);j++) 1237 { 1238 if(deg(i[j])>0) 1239 { 1240 p=i[j]; 1241 break; 1242 } 1243 } 1244 if(deg(p)==1) 1245 { 1246 return(1); 1247 } 1248 for (k=j+1;k<=size(i);k++) 1249 { 1250 if(deg(i[k])!=0) 1251 { 1252 q=GCD(p,i[k]); 1253 if(deg(q)==0) 1254 { 1255 p=p*i[k]; 1256 } 1257 else 1258 { 1259 p=p/q; 1260 p=p*i[k]; 1261 } 1262 } 1263 } 1264 return(p); 1265 } 1266 example 1267 { "EXAMPLE:"; echo = 2; 1268 ring r = 0,(x,y,z),lp; 1269 poly p = (x+y)*(y+z); 1270 poly q = (z4+2)*(y+z); 1271 ideal l=p,q; 1272 poly pr= lcmP(l); 1273 pr; 1274 l=1,1,p,1,1,q,1; 1275 pr=lcmP(l); 1177 1276 pr; 1178 1277 } … … 1434 1533 } 1435 1534 1535 /////////////////////////////////////////////////////////////////////// 1536 1537 proc projdim(list l) 1538 { 1539 int i=size(l)+1; 1540 1541 while(i>2) 1542 { 1543 i; 1544 if((size(l[i])>0)&&(deg(l[i][1])>0)) 1545 { 1546 return(i); 1547 } 1548 } 1549 } 1550 1436 1551 /////////////////////////////////////////////////////////////////////////////// 1437 1552 proc cleanPrimary(list l) … … 1439 1554 int i,j; 1440 1555 list lh; 1441 for(i=1;i<=size(l) div2;i++)1556 for(i=1;i<=size(l)/2;i++) 1442 1557 { 1443 1558 if(deg(l[2*i1][1])>0) … … 1460 1575 EXAMPLE: example minAssPrimes; shows an example 1461 1576 { 1577 #[1]=1; 1462 1578 def @P=basering; 1463 execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),lp;"; 1464 // execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),dp;"; 1579 execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),(" 1580 +ordstr(basering)+");"; 1581 1582 1465 1583 ideal i=fetch(@P,i); 1466 1584 if(size(#)==0) … … 1478 1596 } 1479 1597 list @res,empty; 1598 ideal ser; 1480 1599 option(redSB); 1481 1600 list @pr=facstd(i); 1601 if(size(@pr)==1) 1602 { 1603 attrib(@pr[1],"isSB",1); 1604 if((dim(@pr[1])==0)&&(homog(@pr[1])==1)) 1605 { 1606 setring @P; 1607 list @res=maxideal(1); 1608 return(@res); 1609 } 1610 if(dim(@pr[1])>1) 1611 { 1612 setring @P; 1613 kill gnir; 1614 execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),lp;"; 1615 ideal i=fetch(@P,i); 1616 list @pr=facstd(i); 1617 ideal ser; 1618 } 1619 } 1482 1620 option(noredSB); 1483 1621 int j,k,odim,ndim,count; … … 1521 1659 else 1522 1660 { 1661 ser=ideal(1); 1523 1662 for(j=1;j<=size(@pr);j++) 1524 1663 { 1525 @res[j]=decomp(@pr[j],2); 1664 @res[j]=decomp(@pr[j],2); 1665 // @res[j]=decomp(@pr[j],2,@pr[j],ser); 1666 // for(k=1;k<=size(@res[j]);k++) 1667 // { 1668 // ser=intersect1(ser,@res[j][k]); 1669 // } 1526 1670 } 1527 1671 } … … 1546 1690 /////////////////////////////////////////////////////////////////////////////// 1547 1691 1548 proc union(list l )1692 proc union(list li) 1549 1693 { 1550 1694 int i,j,k; 1695 1696 def P=basering; 1697 1698 execute "ring ir = ("+charstr(basering)+"),("+varstr(basering)+"),lp;"; 1699 list l=fetch(P,li); 1551 1700 list @erg; 1552 i=0;1553 1701 1554 1702 for(k=1;k<=size(l);k++) 1555 1703 { 1556 for(j=1;j<=size(l[k]) div2;j++)1704 for(j=1;j<=size(l[k])/2;j++) 1557 1705 { 1558 1706 if(deg(l[k][2*j][1])!=0) … … 1611 1759 @wos[size(@wos)+1]=@erg[size(@erg)]; 1612 1760 } 1613 return(@wos); 1761 setring P; 1762 list @ser=fetch(ir,@wos); 1763 return(@ser); 1614 1764 } 1615 1765 /////////////////////////////////////////////////////////////////////////////// 1616 proc radical (ideal i)1766 proc radicalOld(ideal i) 1617 1767 { 1618 1768 list pr=minAssPrimes(i,1); … … 1624 1774 } 1625 1775 return(k); 1776 } 1777 /////////////////////////////////////////////////////////////////////////////// 1778 proc equiRadical(ideal i) 1779 { 1780 return(radical(i,1)); 1626 1781 } 1627 1782 /////////////////////////////////////////////////////////////////////////////// … … 1636 1791 { 1637 1792 def @P = basering; 1638 list primary,indep ;1793 list primary,indep,ltras; 1639 1794 intvec @vh,isat; 1640 1795 int @wr,@k,@n,@m,@n1,@n2,@n3,homo; … … 1665 1820 if(homo==1) 1666 1821 { 1667 tras=std(i); 1822 ltras=mstd(i); 1823 attrib(ltras[1],"isSB",1); 1824 tras=ltras[1]; 1668 1825 if(dim(tras)==0) 1669 1826 { 1670 primary[1]= tras;1827 primary[1]=ltras[2]; 1671 1828 primary[2]=maxideal(1); 1672 1829 if(@wr>0) … … 1679 1836 return(primary); 1680 1837 } 1681 intvec @hilb=hilb(tras,1);1838 intvec @hilb=hilb(tras,1); 1682 1839 } 1683 1840 … … 1737 1894 if(nvars(basering)==1) 1738 1895 { 1896 1739 1897 list fac=factor(@j[1]); 1740 1898 list gprimary; … … 1776 1934 return(primary); 1777 1935 } 1778 1779 1780 //1781 //search for a maximal independent set indep,i.e.1782 //look for subring such that the intersection with the ideal is zero1783 //j intersected with K[var(indep[3]+1),...,var(nvar] is zero,1784 //indep[1] is the new varstring and indep[2] the string for the blockordering1785 //1786 1936 1787 1937 poly @gs,@gh,@p; … … 1792 1942 list fett; 1793 1943 int lauf,di; 1944 // 1945 //search for a maximal independent set indep,i.e. 1946 //look for subring such that the intersection with the ideal is zero 1947 //j intersected with K[var(indep[3]+1),...,var(nvar] is zero, 1948 //indep[1] is the new varstring and indep[2] the string for the blockordering 1949 // 1794 1950 1795 1951 if(@wr!=1) … … 1959 2115 @n3=@n2; 1960 2116 1961 for(@n1=1;@n1<=size(collectprimary) div2;@n1++)2117 for(@n1=1;@n1<=size(collectprimary)/2;@n1++) 1962 2118 { 1963 2119 if(deg(collectprimary[2*@n1][1])>0) … … 1975 2131 //mentioned above is really computed 1976 2132 1977 for(@n=@n3 div 2+1;@n<=@n2 div2;@n++)2133 for(@n=@n3/2+1;@n<=@n2/2;@n++) 1978 2134 { 1979 2135 if(specialIdealsEqual(quprimary[2*@n1],quprimary[2*@n])) … … 2042 2198 quprimary[2]=ideal(1); 2043 2199 } 2200 2044 2201 // 2045 2202 //notice that j=sat(j,gh) intersected with (j,gh^n) … … 2055 2212 ideal htest=quprimary[1]; 2056 2213 2057 for (@n1=2;@n1<=size(quprimary) div2;@n1++)2214 for (@n1=2;@n1<=size(quprimary)/2;@n1++) 2058 2215 { 2059 2216 htest=intersect(htest,quprimary[2*@n11]); … … 2064 2221 ideal htest=quprimary[2]; 2065 2222 2066 for (@n1=2;@n1<=size(quprimary) div2;@n1++)2223 for (@n1=2;@n1<=size(quprimary)/2;@n1++) 2067 2224 { 2068 2225 htest=intersect(htest,quprimary[2*@n1]); … … 2193 2350 @n3=@n2; 2194 2351 2195 for(@n1=1;@n1<=size(collectprimary) div2;@n1++)2352 for(@n1=1;@n1<=size(collectprimary)/2;@n1++) 2196 2353 { 2197 2354 if(deg(collectprimary[2*@n1][1])>0) … … 2209 2366 //mentioned above is really computed 2210 2367 2211 for(@n=@n3 div 2+1;@n<=@n2 div2;@n++)2368 for(@n=@n3/2+1;@n<=@n2/2;@n++) 2212 2369 { 2213 2370 if(specialIdealsEqual(quprimary[2*@n1],quprimary[2*@n])) … … 2279 2436 testPrimary( pr, i); 2280 2437 } 2438 2439 /////////////////////////////////////////////////////////////////////////////// 2440 proc radical(ideal i,list #) 2441 { 2442 def @P=basering; 2443 int j,il; 2444 if(size(#)>0) 2445 { 2446 il=#[1]; 2447 } 2448 ideal re=1; 2449 option(redSB); 2450 list pr=facstd(i); 2451 2452 if(size(pr)==1) 2453 { 2454 attrib(pr[1],"isSB",1); 2455 if((dim(pr[1])==0)&&(homog(pr[1])==1)) 2456 { 2457 ideal @res=maxideal(1); 2458 return(@res); 2459 } 2460 if(dim(pr[1])>1) 2461 { 2462 execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),lp;"; 2463 ideal i=fetch(@P,i); 2464 list @pr=facstd(i); 2465 setring @P; 2466 pr=fetch(gnir,@pr); 2467 } 2468 } 2469 option(noredSB); 2470 int s=size(pr); 2471 if(s==1) 2472 { 2473 return(radicalEHV(i,ideal(1),il)); 2474 } 2475 intvec pos; 2476 pos[s]=0; 2477 2478 if(il==1) 2479 { 2480 int ndim,k; 2481 attrib(pr[1],"isSB",1); 2482 int odim=dim(pr[1]); 2483 int count=1; 2484 2485 for(j=2;j<=s;j++) 2486 { 2487 attrib(pr[j],"isSB",1); 2488 ndim=dim(pr[j]); 2489 if(ndim>odim) 2490 { 2491 for(k=count;k<=j1;k++) 2492 { 2493 pos[k]=1; 2494 } 2495 count=j; 2496 odim=ndim; 2497 } 2498 if(ndim<odim) 2499 { 2500 pos[j]=1; 2501 } 2502 } 2503 } 2504 2505 for(j=1;j<=s;j++) 2506 { 2507 if(pos[j]==0) 2508 { 2509 re=intersect1(re,radicalEHV(pr[s+1j],re,il)); 2510 } 2511 } 2512 return(re); 2513 } 2514 2515 proc intersect1(ideal i,ideal j) 2516 { 2517 def R=basering; 2518 execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+",@t),dp;"; 2519 ideal i=var(nvars(basering))*imap(R,i)+(var(nvars(basering))1)*imap(R,j); 2520 ideal j=eliminate(i,var(nvars(basering))); 2521 setring R; 2522 map phi=gnir,maxideal(1); 2523 return(phi(j)); 2524 } 2525 2526 2527 /////////////////////////////////////////////////////////////////////////////// 2528 proc radicalKL (list m,ideal ser,list #) 2529 USAGE: decomp(i); i ideal (for primary decomposition) (resp. 2530 decomp(i,1); (for the minimal associated primes) ) 2531 RETURN: list = list of primary ideals and their associated primes 2532 (at even positions in the list) 2533 (resp. a list of the minimal associated primes) 2534 NOTE: Algorithm of Gianni, Traeger, Zacharias 2535 EXAMPLE: example decomp; shows an example 2536 { 2537 ideal i=m[2]; 2538 // 2539 //i is the zeroideal 2540 // 2541 2542 if(size(i)==0) 2543 { 2544 return(ideal(0)); 2545 } 2546 2547 def @P = basering; 2548 list indep,allindep,restindep,fett,@mu; 2549 intvec @vh,isat; 2550 int @wr,@k,@n,@m,@n1,@n2,@n3,lauf,di; 2551 ideal @j,@j1,fac,@h,collectrad,htrad,lsau; 2552 ideal rad=ideal(1); 2553 ideal te=ser; 2554 poly @p,@q; 2555 string @va,quotring; 2556 int homo=homog(i); 2557 2558 if(size(#)>0) 2559 { 2560 @wr=#[1]; 2561 } 2562 @j=m[1]; 2563 @j1=m[2]; 2564 int jdim=dim(@j); 2565 if(size(reduce(ser,@j))==0) 2566 { 2567 return(ser); 2568 } 2569 if(homo==1) 2570 { 2571 if(jdim==0) 2572 { 2573 option(noredSB); 2574 return(maxideal(1)); 2575 } 2576 intvec @hilb=hilb(@j,1); 2577 } 2578 2579 2580 // 2581 //j is the ring 2582 // 2583 2584 if (jdim==1) 2585 { 2586 option(noredSB); 2587 return(ideal(0)); 2588 } 2589 2590 // 2591 // the case of one variable 2592 // 2593 2594 if(nvars(basering)==1) 2595 { 2596 fac=factorize(@j[1],1); 2597 poly @p=1; 2598 for(@k=1;@k<=size(fac);@k++) 2599 { 2600 @p=@p*fac[@k]; 2601 } 2602 option(noredSB); 2603 2604 return(ideal(@p)); 2605 } 2606 // 2607 //the case of a complete intersection 2608 // 2609 if(jdim+size(@j1)==nvars(basering)) 2610 { 2611 // ideal jac=minor(jacob(@j1),size(@j1)); 2612 // return(quotient(@j1,jac)); 2613 } 2614 2615 // 2616 //the zerodimensional case 2617 // 2618 2619 if (jdim==0) 2620 { 2621 @j1=system("finduni",@j); 2622 for(@k=1;@k<=size(@j1);@k++) 2623 { 2624 fac=factorize(cleardenom(@j1[@k]),1); 2625 @p=fac[1]; 2626 for(@n=2;@n<=size(fac);@n++) 2627 { 2628 @p=@p*fac[@n]; 2629 } 2630 @j=@j,@p; 2631 } 2632 @j=std(@j); 2633 option(noredSB); 2634 return(@j); 2635 } 2636 2637 // 2638 //search for a maximal independent set indep,i.e. 2639 //look for subring such that the intersection with the ideal is zero 2640 //j intersected with K[var(indep[3]+1),...,var(nvar] is zero, 2641 //indep[1] is the new varstring and indep[2] the string for the blockordering 2642 // 2643 2644 indep=maxIndependSet(@j); 2645 2646 di=dim(@j); 2647 2648 for(@m=1;@m<=size(indep);@m++) 2649 { 2650 if((indep[@m][1]==varstr(basering))&&(@m==1)) 2651 //this is the good case, nothing to do, just to have the same notations 2652 //change the ring 2653 { 2654 execute "ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),(" 2655 +ordstr(basering)+");"; 2656 ideal @j=fetch(@P,@j); 2657 attrib(@j,"isSB",1); 2658 } 2659 else 2660 { 2661 @va=string(maxideal(1)); 2662 execute "ring gnir1 = ("+charstr(basering)+"),("+indep[@m][1]+"),(" 2663 +indep[@m][2]+");"; 2664 execute "map phi=@P,"+@va+";"; 2665 if(homo==1) 2666 { 2667 ideal @j=std(phi(@j),@hilb); 2668 } 2669 else 2670 { 2671 ideal @j=std(phi(@j)); 2672 } 2673 } 2674 if((deg(@j[1])==0)(dim(@j)<jdim)) 2675 { 2676 setring @P; 2677 break; 2678 } 2679 for (lauf=1;lauf<=size(@j);lauf++) 2680 { 2681 fett[lauf]=size(@j[lauf]); 2682 } 2683 // 2684 //we have now the following situation: 2685 //j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass 2686 //to this quotientring, j is their still a standardbasis, the 2687 //leading coefficients of the polynomials there (polynomials in 2688 //K[var(nnp+1),..,var(nva)]) are collected in the list h, 2689 //we need their ggt, gh, because of the following: 2690 //let (j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..the rest..] 2691 //intersected with K[var(1),...,var(nva)] is (j:gh^n) 2692 //on the other hand j=(j,gh^n) intersected with (j:gh^n) 2693 2694 // 2695 2696 //the arrangement for the quotientring K(var(nnp+1),..,var(nva))[..the rest..] 2697 //and the map phi:K[var(1),...,var(nva)] >K(var(nnpr+1),..,var(nva))[..the rest..] 2698 // 2699 2700 quotring=prepareQuotientring(nvars(basering)indep[@m][3]); 2701 2702 // 2703 //we pass to the quotientring K(var(nnp+1),..,var(nva))[..the rest..] 2704 // 2705 2706 execute quotring; 2707 2708 // @j considered in the quotientring 2709 ideal @j=imap(gnir1,@j); 2710 2711 kill gnir1; 2712 2713 //j is a standardbasis in the quotientring but usually not minimal 2714 //here it becomes minimal 2715 2716 @j=clearSB(@j,fett); 2717 attrib(@j,"isSB",1); 2718 2719 //we need later ggt(h[1],...)=gh for saturation 2720 ideal @h; 2721 if(deg(@j[1])>0) 2722 { 2723 for(@n=1;@n<=size(@j);@n++) 2724 { 2725 @h[@n]=leadcoef(@j[@n]); 2726 } 2727 //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..] 2728 option(redSB); 2729 @j=interred(@j); 2730 attrib(@j,"isSB",1); 2731 list @mo=@j,@j; 2732 ideal zero_rad= radicalKL(@mo,ideal(1)); 2733 } 2734 else 2735 { 2736 ideal zero_rad=ideal(1); 2737 } 2738 2739 //we need the intersection of the ideals in the list quprimary with the 2740 //polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal 2741 //but fi polynomials, then the intersection of q with the polynomialring 2742 //is the saturation of the ideal generated by f1,...,fr with respect to 2743 //h which is the lcm of the leading coefficients of the fi considered in the 2744 //quotientring: this is coded in saturn 2745 2746 ideal hpl; 2747 2748 for(@n=1;@n<=size(zero_rad);@n++) 2749 { 2750 hpl=hpl,leadcoef(zero_rad[@n]); 2751 } 2752 2753 // 2754 //we leave the quotientring K(var(nnp+1),..,var(nva))[..the rest..] 2755 //back to the polynomialring 2756 // 2757 setring @P; 2758 2759 collectrad=imap(quring,zero_rad); 2760 lsau=simplify(imap(quring,hpl),2); 2761 @h=imap(quring,@h); 2762 2763 kill quring; 2764 2765 2766 //here the intersection with the polynomialring 2767 //mentioned above is really computed 2768 2769 collectrad=sat2(collectrad,lsau)[1]; 2770 2771 if(deg(@h[1])>=0) 2772 { 2773 fac=ideal(0); 2774 for(lauf=1;lauf<=ncols(@h);lauf++) 2775 { 2776 if(deg(@h[lauf])>0) 2777 { 2778 fac=fac+factorize(@h[lauf],1); 2779 } 2780 } 2781 fac=simplify(fac,4); 2782 @q=1; 2783 for(lauf=1;lauf<=size(fac);lauf++) 2784 { 2785 @q=@q*fac[lauf]; 2786 } 2787 2788 2789 @mu=mstd(quotient(@j+ideal(@q),rad)); 2790 @j=@mu[1]; 2791 attrib(@j,"isSB",1); 2792 2793 } 2794 if((deg(rad[1])>0)&&(deg(collectrad[1])>0)) 2795 { 2796 int xyz=timer; 2797 "bei collecterad"; 2798 rad=intersect(rad,collectrad); 2799 timerxyz; 2800 } 2801 else 2802 { 2803 if(deg(collectrad[1])>0) 2804 { 2805 rad=collectrad; 2806 } 2807 } 2808 2809 te=simplify(reduce(te*rad,@j),2); 2810 2811 if((dim(@j)<di)(size(te)==0)) 2812 { 2813 break; 2814 } 2815 if(homo==1) 2816 { 2817 @hilb=hilb(@j,1); 2818 } 2819 } 2820 2821 if(((@wr==1)&&(dim(@j)<di))(deg(@j[1])==0)(size(te)==0)) 2822 { 2823 return(rad); 2824 } 2825 ideal tes=radicalKL(@mu,rad,@wr); 2826 int sml=timer; 2827 "bei rad"; 2828 rad=intersect(rad,tes); 2829 timersml; 2830 // rad=intersect(rad,radicalKL(@mu,ideal(1),@wr)); 2831 2832 2833 option(noredSB); 2834 return(rad); 2835 } 2836 2837 2838 example 2839 { "EXAMPLE:"; echo = 2; 2840 } 2841 2842 2843 proc radicalEHV(ideal i,ideal re,list #) 2844 { 2845 ideal J,I,I0,radI0,L,radI1,I2,radI2; 2846 int l,il; 2847 if(size(#)>0) 2848 { 2849 il=#[1]; 2850 } 2851 option(redSB); 2852 list m=mstd(i); 2853 I=m[2]; 2854 option(noredSB); 2855 if(size(reduce(re,m[1]))==0) 2856 { 2857 return(re); 2858 } 2859 int cod=nvars(basering)dim(m[1]); 2860 if(nvars(basering)<9) 2861 { 2862 if(cod==size(m[2])) 2863 { 2864 J=minor(jacob(I),cod); 2865 return(quotient(I,J)); 2866 } 2867 2868 for(l=1;l<=cod;l++) 2869 { 2870 I0[l]=I[l]; 2871 } 2872 if(dim(std(I0))+cod==nvars(basering)) 2873 { 2874 J=minor(jacob(I0),cod); 2875 radI0=quotient(I0,J); 2876 L=quotient(radI0,I); 2877 radI1=quotient(radI0,L); 2878 2879 if(size(reduce(radI1,m[1]))==0) 2880 { 2881 return(I); 2882 } 2883 if(il==1) 2884 { 2885 return(radI1); 2886 } 2887 2888 I2=sat(I,radI1)[1]; 2889 2890 if(deg(I2[1])<=0) 2891 { 2892 return(radI1); 2893 } 2894 return(intersect(radI1,radicalEHV(I2,re,il))); 2895 } 2896 } 2897 return(radicalKL(m,re,il)); 2898 } 2899 2900 proc Ann(module M) 2901 { 2902 M=prune(M); //to obtain a small embedding 2903 return(quotient(M,freemodule(nrows(M)))); 2904 } 2905 2906 //computes the equidimensional part of the ideal i of codimension e 2907 proc int_ass_primary_e(ideal i, int e) 2908 { 2909 if(homog(i)!=1) 2910 { 2911 i=std(i); 2912 } 2913 list re=sres(i,0); //the resolution 2914 re=minres(re); //minimized resolution 2915 ideal ann=AnnExt_R(e,re); 2916 if(nvars(basering)dim(std(ann))!=e) 2917 { 2918 return(ideal(1)); 2919 } 2920 return(ann); 2921 } 2922 2923 //computes all equidimensional parts of the ideal i 2924 proc prepareAss(ideal i) 2925 { 2926 ideal j=std(i); 2927 int cod=nvars(basering)dim(j); 2928 int e; 2929 list er; 2930 ideal ann; 2931 if(homog(i)==1) 2932 { 2933 list re=sres(i,0); //the resolution 2934 re=minres(re); //minimized resolution 2935 } 2936 else 2937 { 2938 list re=mres(i,0); //fehler in sres 2939 } 2940 for(e=cod;e<=nvars(basering);e++) 2941 { 2942 ann=AnnExt_R(e,re); 2943 if(nvars(basering)dim(std(ann))==e) 2944 { 2945 er[size(er)+1]=equiRadical(ann); 2946 } 2947 } 2948 return(er); 2949 } 2950 2951 //computes the annihilator of Ext^n(R/i,R) with given resolution re 2952 //n is not necessarily the number of variables 2953 proc AnnExt_R(int n,list re) 2954 { 2955 if(n<nvars(basering)) 2956 { 2957 matrix f=transpose(re[n+1]); //Hom(_,R) 2958 module k=res(f,2)[2]; //the kernel 2959 matrix g=transpose(re[n]); //the image of Hom(_,R) 2960 ideal ann=quotient(g,k); //the anihilator 2961 } 2962 else 2963 { 2964 ideal ann=Ann(transpose(re[n])); 2965 } 2966 return(ann); 2967 } 2968
Note: See TracChangeset
for help on using the changeset viewer.