Changeset 808a9f3 in git for Singular/LIB/primdec.lib
- Timestamp:
- Jun 20, 2006, 1:55:10 PM (18 years ago)
- Branches:
- (u'spielwiese', '4a9821a93ffdc22a6696668bd4f6b8c9de3e6c5f')
- Children:
- 88b2dd92d922dd2e8a9f9759e7edde2f1e1d2829
- Parents:
- 951a6327b56dcd11021810ab04eba401d1432b98
- File:
-
- 1 edited
-
Singular/LIB/primdec.lib (modified) (26 diffs)
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/primdec.lib
r951a63 r808a9f3 1 1 /////////////////////////////////////////////////////////////////////////////// 2 version="$Id: primdec.lib,v 1.11 8 2006-06-19 15:08:32Singular Exp $";2 version="$Id: primdec.lib,v 1.119 2006-06-20 11:55:10 Singular Exp $"; 3 3 category="Commutative Algebra"; 4 4 info=" … … 7 7 @* Wolfram Decker, decker@math.uni-sb.de (SY) 8 8 @* Hans Schoenemann, hannes@mathematik.uni-kl.de (SY) 9 @* Santiago Laplagne, slaplagn@dm.uba.ar (GTZ) 9 10 10 11 OVERVIEW: … … 1021 1022 } 1022 1023 1023 if(voice>=8){ primary=extF(primary)};1024 if(voice>=8){ primary=extF(primary) } 1024 1025 1025 1026 //test whether all ideals in the decomposition are primary and … … 2043 2044 static proc minAssPrimes(ideal i, list #) 2044 2045 "USAGE: minAssPrimes(i); i ideal 2045 minAssPrimes(i,1); i ideal (to use also the factorizing Groebner) 2046 Optional parameters in list #: (can be entered in any order) 2047 0, "facstd" -> uses facstd to first decompose the ideal 2048 1, "noFacstd" -> does not use facstd (default) 2049 "SL" -> the new algorithm is used (default) 2050 "GTZ" -> the old algorithm is used 2046 2051 RETURN: list = the minimal associated prime ideals of i 2047 2052 EXAMPLE: example minAssPrimes; shows an example 2048 2053 " 2049 2054 { 2050 def P=basering; 2051 if(size(i)==0){return(list(ideal(0)));} 2052 list q=simplifyIdeal(i); 2053 list re=maxideal(1); 2054 int j,a,k; 2055 intvec op=option(get); 2056 map phi=P,q[2]; 2057 2058 if(npars(P)==0){option(redSB);} 2055 string algorithm; // Algorithm to be used 2056 string facstdOption; // To uses proc facstd 2057 int j; // Counter 2058 def P = basering; 2059 if(size(i) == 0){return(list(ideal(0)));} 2060 2061 // Set input parameters 2062 algorithm = "SL"; // Default: SL algorithm 2063 facstdOption = "noFacstd"; // Default: facstd is not used 2064 if(size(#) > 0){ 2065 int valid; 2066 for(j = 1; j <= size(#); j++) 2067 { 2068 valid = 0; 2069 if((typeof(#[j]) == "int") or (typeof(#[j]) == "number")) 2070 { 2071 if (#[j] == 0) {facstdOption = "noFacstd"; valid = 1;} // facstd is not used. 2072 if (#[j] == 1) {facstdOption = "facstd"; valid = 1;} // facstd is used. 2073 } 2074 if(typeof(#[j]) == "string") 2075 { 2076 if(#[j] == "GTZ" || #[j] == "SL") 2077 { 2078 algorithm = #[j]; 2079 valid = 1; 2080 } 2081 if(#[j] == "noFacstd" || #[j] == "facstd") 2082 { 2083 facstdOption = #[j]; 2084 valid = 1; 2085 } 2086 } 2087 if(valid == 0) 2088 { 2089 dbprint(1, "Warning! The following input parameter was not recognized:", #[j]); 2090 } 2091 } 2092 } 2093 2094 list q = simplifyIdeal(i); 2095 list re = maxideal(1); 2096 int a, k; 2097 intvec op = option(get); 2098 map phi = P,q[2]; 2099 2100 list result; 2101 2102 if(npars(P) == 0){option(redSB);} 2059 2103 2060 2104 if(attrib(i,"isSB")!=1) … … 2073 2117 } 2074 2118 } 2075 if(dim(i)==-1){return(ideal(1));} 2076 if((dim(i)==0)&&(npars(P)==0)) 2077 { 2078 int di=vdim(i); 2119 2120 if(dim(i) == -1){re = ideal(1); return(re);} 2121 if((dim(i) == 0) && (npars(P) == 0)) 2122 { 2123 int di = vdim(i); 2079 2124 execute ("ring gnir=("+charstr(P)+"),("+varstr(P)+"),lp;"); 2080 ideal J =std(imap(P,i));2081 attrib(J, "isSB",1);2082 if(vdim(J) !=di)2083 { 2084 J =fglm(P,i);2085 } 2086 list pr =triangMH(J,2);2087 list qr, re;2088 2089 for(k =1;k<=size(pr);k++)2125 ideal J = std(imap(P,i)); 2126 attrib(J, "isSB", 1); 2127 if(vdim(J) != di) 2128 { 2129 J = fglm(P, i); 2130 } 2131 list pr = triangMH(J, 2); 2132 list qr, re; 2133 2134 for(k = 1; k <= size(pr); k++) 2090 2135 { 2091 2136 if(primT(pr[k])) 2092 2137 { 2093 re[size(re) +1]=pr[k];2138 re[size(re) + 1] = pr[k]; 2094 2139 } 2095 2140 else 2096 2141 { 2097 attrib(pr[k], "isSB",1);2098 qr=decomp(pr[k],2);2099 for(j=1;j<=size(qr)/2;j++)2142 attrib(pr[k], "isSB", 1); 2143 // Lines changed 2144 if (algorithm == "GTZ") 2100 2145 { 2101 re[size(re)+1]=qr[2*j]; 2146 qr = decomp(pr[k], 2); 2147 } 2148 else 2149 { 2150 qr = newMinAssSL(pr[k]); 2151 } 2152 for(j = 1; j <= size(qr) / 2; j++) 2153 { 2154 re[size(re) + 1] = qr[2 * j]; 2102 2155 } 2103 2156 } 2104 2157 } 2105 2158 setring P; 2106 re =imap(gnir,re);2107 option(set, op);2159 re = imap(gnir, re); 2160 option(set, op); 2108 2161 return(phi(re)); 2109 2162 } 2110 2163 2111 if((size(#)==0)||(dim(i)==0)) 2112 { 2113 re[1]=decomp(i,2); 2114 re=union(re); 2115 option(set,op); 2116 return(phi(re)); 2117 } 2118 2119 q=facstd(i); 2164 // Lines changed 2165 if ((facstdOption == "noFacstd") || (dim(i) == 0)) 2166 { 2167 if (algorithm == "GTZ") 2168 { 2169 re[1] = decomp(i, 2); 2170 } 2171 else 2172 { 2173 re[1] = newMinAssSL(i); 2174 } 2175 re = union(re); 2176 option(set, op); 2177 return(phi(re)); 2178 } 2179 2180 q = facstd(i); 2120 2181 2121 2182 /* 2122 if((size(q) ==1)&&(dim(i)>1))2183 if((size(q) == 1) && (dim(i) > 1)) 2123 2184 { 2124 2185 execute ("ring gnir=("+charstr(P)+"),("+varstr(P)+"),lp;"); 2125 2126 list p=facstd(fetch(P,i)); 2127 if(size(p)>1) 2128 { 2129 a=1; 2186 list p = facstd(fetch(P, i)); 2187 if(size(p) > 1) 2188 { 2189 a = 1; 2130 2190 setring P; 2131 q =fetch(gnir,p);2191 q = fetch(gnir,p); 2132 2192 } 2133 2193 else … … 2140 2200 2141 2201 option(set,op); 2142 for(j=1;j<=size(q);j++) 2143 { 2144 if(a==0){attrib(q[j],"isSB",1);} 2145 re[j]=decomp(q[j],2); 2146 } 2147 re=union(re); 2202 // Debug 2203 dbprint(printlevel - voice, "Components returned by facstd", size(q), q); 2204 for(j = 1; j <= size(q); j++) 2205 { 2206 if(a == 0){attrib(q[j], "isSB", 1);} 2207 // Debug 2208 dbprint(printlevel - voice, "We compute the decomp of component", j); 2209 // Lines changed 2210 if (algorithm == "GTZ") 2211 { 2212 re[j] = decomp(q[j], 2); 2213 } 2214 else 2215 { 2216 re[j] = newMinAssSL(q[j]); 2217 } 2218 // Debug 2219 dbprint(printlevel - voice, "Number of components obtained for this component:", size(re[j]) / 2); 2220 dbprint(printlevel - voice, "re[j]:", re[j]); 2221 } 2222 re = union(re); 2223 2148 2224 return(phi(re)); 2149 2225 } … … 5420 5496 /////////////////////////////////////////////////////////////////////////////// 5421 5497 proc minAssGTZ(ideal i,list #) 5422 "USAGE: minAssGTZ(i); i ideal 5423 minAssGTZ(i,1); i ideal does not use the factorizing Groebner 5498 "USAGE: minAssGTZ(i); i ideal 5499 Optional parameters in list #: (can be entered in any order) 5500 0, "facstd" -> uses facstd to first decompose the ideal (default) 5501 1, "noFacstd" -> does not use facstd 5502 "SL" -> the new algorithm is used (default) 5503 "GTZ" -> the old algorithm is used 5424 5504 RETURN: a list, the minimal associated prime ideals of i. 5425 5505 NOTE: Designed for characteristic 0, works also in char k > 0 based … … 5428 5508 " 5429 5509 { 5510 int j; 5511 string algorithm; 5512 string facstdOption; 5513 int useFac; 5514 5515 // Set input parameters 5516 algorithm = "SL"; // Default: SL algorithm 5517 facstdOption = "facstd"; 5518 if(size(#) > 0){ 5519 int valid; 5520 for(j = 1; j <= size(#); j++) 5521 { 5522 valid = 0; 5523 if((typeof(#[j]) == "int") or (typeof(#[j]) == "number")) 5524 { 5525 if (#[j] == 1) {facstdOption = "noFacstd"; valid = 1;} // facstd is not used. 5526 if (#[j] == 0) {facstdOption = "facstd"; valid = 1;} // facstd is used. 5527 } 5528 if(typeof(#[j]) == "string") 5529 { 5530 if((#[j] == "GTZ") || (#[j] == "SL")) 5531 { 5532 algorithm = #[j]; 5533 valid = 1; 5534 } 5535 if((#[j] == "noFacstd") || (#[j] == "facstd")) 5536 { 5537 facstdOption = #[j]; 5538 valid = 1; 5539 } 5540 } 5541 if(valid == 0) 5542 { 5543 dbprint(1, "Warning! The following input parameter was not recognized:", #[j]); 5544 } 5545 } 5546 } 5547 5430 5548 if(ord_test(basering)!=1) 5431 5549 { … … 5438 5556 return(algeDeco(i,2)); 5439 5557 } 5440 if(size(#)==0) 5441 { 5442 return(minAssPrimes(i,1)); 5443 } 5444 return(minAssPrimes(i)); 5558 5559 list result; 5560 result = minAssPrimes(i, facstdOption, algorithm); 5561 return(result); 5445 5562 } 5446 5563 example … … 5544 5661 5545 5662 // Set input parameters 5546 algorithm = "SL"; // Default: SL algorithm 5547 il = 0; // Default: Full radical (not only equidim part) 5548 if (homog(i) == 1) { // Default: facStd is used, except if the ideal is homogeneous. 5549 useFac = 0; 5550 } else { 5551 useFac = 1; 5552 }; 5553 if(size(#) > 0){ 5554 int valid; 5555 for(j = 1; j <= size(#); j++){ 5556 valid = 0; 5557 if((typeof(#[j]) == "int") or (typeof(#[j]) == "number")) { 5558 il = #[j]; // If il == 1, equiRadical is computed 5559 valid = 1; 5560 }; 5561 if(typeof(#[j]) == "string"){ 5562 if(#[j] == "KL") { 5563 algorithm = "KL"; 5564 valid = 1}; 5565 if(#[j] == "KLdp") { 5566 algorithm = "KLdp"; 5567 valid = 1}; 5568 if(#[j] == "SL") { 5569 algorithm = "SL"; 5570 valid = 1}; 5571 if(#[j] == "noFacstd") { 5572 useFac = 0; 5573 valid = 1}; 5574 if(#[j] == "facstd") { 5575 useFac = 1; 5576 valid = 1}; 5577 if(#[j] == "equiRad") { 5578 il = 1; 5579 valid = 1}; 5580 if(#[j] == "fullRad") { 5581 il = 0; 5582 valid = 1}; 5583 }; 5584 if(valid == 0) { 5585 dbprint(1, "Warning! The following input parameter was not recognized:", #[j]); 5586 }; 5587 }; 5588 }; 5663 algorithm = "SL"; // Default: SL algorithm 5664 il = 0; // Default: Full radical (not only equidim part) 5665 if (homog(i) == 1) 5666 { // Default: facStd is used, except if the ideal is homogeneous. 5667 useFac = 0; 5668 } 5669 else 5670 { 5671 useFac = 1; 5672 } 5673 if(size(#) > 0) 5674 { 5675 int valid; 5676 for(j = 1; j <= size(#); j++) 5677 { 5678 valid = 0; 5679 if((typeof(#[j]) == "int") or (typeof(#[j]) == "number")) 5680 { 5681 il = #[j]; // If il == 1, equiRadical is computed 5682 valid = 1; 5683 } 5684 if(typeof(#[j]) == "string") 5685 { 5686 if(#[j] == "KL") 5687 { 5688 algorithm = "KL"; 5689 valid = 1 5690 } 5691 if(#[j] == "KLdp") 5692 { 5693 algorithm = "KLdp"; 5694 valid = 1 5695 } 5696 if(#[j] == "SL") 5697 { 5698 algorithm = "SL"; 5699 valid = 1 5700 } 5701 if(#[j] == "noFacstd") 5702 { 5703 useFac = 0; 5704 valid = 1 5705 } 5706 if(#[j] == "facstd") 5707 { 5708 useFac = 1; 5709 valid = 1 5710 } 5711 if(#[j] == "equiRad") 5712 { 5713 il = 1; 5714 valid = 1 5715 } 5716 if(#[j] == "fullRad") 5717 { 5718 il = 0; 5719 valid = 1 5720 } 5721 } 5722 if(valid == 0) 5723 { 5724 dbprint(1, "Warning! The following input parameter was not recognized:", #[j]); 5725 } 5726 } 5727 } 5589 5728 5590 5729 if(size(i) == 0){return(ideal(0));} … … 5612 5751 { 5613 5752 pr = facstd(i); 5614 } else { 5615 pr = i 5616 }; 5753 } 5754 else 5755 { 5756 pr = i 5757 } 5617 5758 option(set, op); 5618 5759 int s = size(pr); 5619 if(useFac == 1) { 5760 if(useFac == 1) 5761 { 5620 5762 dbprint(printlevel - voice, "Number of components returned by facstd: ", s); 5621 } ;5763 } 5622 5764 for(j = 1; j <= s; j++) 5623 5765 { 5624 attrib(pr[s + 1 - j], "isSB", 1); 5625 if((size(reduce(rad, pr[s + 1 - j], 1)) != 0) && ((dim(pr[s + 1 - j]) == di) || !il)) 5626 { 5627 // SL Debug messages 5628 dbprint(printlevel-voice, "We shall compute the radical of ", pr[s + 1 - j]); 5629 dbprint(printlevel-voice, "The dimension is: ", dim(pr[s+1-j])); 5630 5631 if(algorithm == "KL") { 5632 rad = intersect(rad, radicalKL(pr[s + 1 - j], rad, il)); 5633 } 5634 if(algorithm == "KLdp") { 5635 rad = intersect(rad, newRadicalKL(pr[s + 1 - j], rad, il)); 5636 } 5637 if(algorithm == "SL") { 5638 rad = intersect(rad, newRadicalSL(pr[s + 1 - j], il)); 5639 }; 5640 } else 5641 { 5642 // SL Debug 5643 dbprint(printlevel-voice, "The radical of this component is not needed."); 5644 dbprint(printlevel-voice, "size(reduce(rad, pr[s + 1 - j], 1))", size(reduce(rad, pr[s + 1 - j], 1))); 5645 dbprint(printlevel-voice, "dim(pr[s + 1 - j])", dim(pr[s + 1 - j])); 5646 dbprint(printlevel-voice, "il", il); 5647 }; 5766 attrib(pr[s + 1 - j], "isSB", 1); 5767 if((size(reduce(rad, pr[s + 1 - j], 1)) != 0) && ((dim(pr[s + 1 - j]) == di) || !il)) 5768 { 5769 // SL Debug messages 5770 dbprint(printlevel-voice, "We shall compute the radical of ", pr[s + 1 - j]); 5771 dbprint(printlevel-voice, "The dimension is: ", dim(pr[s+1-j])); 5772 5773 if(algorithm == "KL") 5774 { 5775 rad = intersect(rad, radicalKL(pr[s + 1 - j], rad, il)); 5776 } 5777 if(algorithm == "KLdp") 5778 { 5779 rad = intersect(rad, newRadicalKL(pr[s + 1 - j], rad, il)); 5780 } 5781 if(algorithm == "SL") 5782 { 5783 rad = intersect(rad, newRadicalSL(pr[s + 1 - j], il)); 5784 } 5785 } 5786 else 5787 { 5788 // SL Debug 5789 dbprint(printlevel-voice, "The radical of this component is not needed."); 5790 dbprint(printlevel-voice, "size(reduce(rad, pr[s + 1 - j], 1))", size(reduce(rad, pr[s + 1 - j], 1))); 5791 dbprint(printlevel-voice, "dim(pr[s + 1 - j])", dim(pr[s + 1 - j])); 5792 dbprint(printlevel-voice, "il", il); 5793 } 5648 5794 } 5649 5795 return(interred(phi(rad))); … … 5670 5816 // and radicalKL removed. 5671 5817 // 5672 static proc newRadicalKL(ideal I, ideal ser, list #) { 5673 // I The ideal for which the radical is computed 5674 // ideal ser Same as in radicalKL (used to reduce components already obtained) 5675 // list # If #[1] = 1, equiradical is computed. 5676 5677 // I needs to be a Groebner basis. 5678 if (attrib(I, "isSB") != 1) { 5679 I = groebner(I); 5680 }; 5681 5682 ideal rad; // The radical 5683 int allIndep = 1; // All max independent sets are used 5684 5685 list result = newRadicalReduction(I, ser, allIndep, #); 5686 int done = result[3]; 5687 rad = result[1]; 5688 if (done == 0) { 5689 rad = intersect(rad, newRadicalKL(result[2], ideal(1), #)); 5690 }; 5691 return(rad); 5692 }; 5693 5818 static proc newRadicalKL(ideal I, ideal ser, list #) 5819 { 5820 // I The ideal for which the radical is computed 5821 // ideal ser Same as in radicalKL (used to reduce components already obtained) 5822 // list # If #[1] = 1, equiradical is computed. 5823 5824 // I needs to be a Groebner basis. 5825 if (attrib(I, "isSB") != 1) 5826 { 5827 I = groebner(I); 5828 } 5829 5830 ideal rad; // The radical 5831 int allIndep = 1; // All max independent sets are used 5832 5833 list result = newRadicalReduction(I, ser, allIndep, #); 5834 int done = result[3]; 5835 rad = result[1]; 5836 if (done == 0) 5837 { 5838 rad = intersect(rad, newRadicalKL(result[2], ideal(1), #)); 5839 } 5840 return(rad); 5841 } 5694 5842 5695 5843 /////////////////////////////////////////////////////////////////////////////// … … 5706 5854 // obtained in intermediate steps. 5707 5855 { 5708 ideal rad = 1; 5709 ideal equiRad = 1; 5710 list primes; 5711 int k; // Counter 5712 int il; // If il = 1, only the equiradical is required. 5713 int iDim; // The dimension of I 5714 int stop = 0; // Checks if the radical has been obtained 5715 5716 if (attrib(I, "isSB") != 1) { 5717 I = groebner(I); 5718 }; 5719 iDim = dim(I); 5720 5721 // Checks if only equiradical is required 5722 if (size(#) > 0) { 5723 il = #[1]; 5724 }; 5725 5726 while(stop == 0) { 5727 dbprint (printlevel-voice, "// We call radLoopR to find new prime ideals."); 5728 primes = newRadicalSLIteration(I, rad); // A list of primes or intersections of primes, not included in P 5729 dbprint (printlevel - voice, "// Output of Iteration Step:"); 5730 dbprint (printlevel - voice, primes); 5731 if (size(primes) > 0) { 5732 dbprint (printlevel - voice, "// We intersect P with the ideal just obtained."); 5733 for(k = 1; k <= size(primes); k++) { 5734 rad = intersect(rad, primes[k]); 5735 if (il == 1){ 5736 if (attrib(primes[k], "isSB") != 1) { 5737 primes[k] = groebner(primes[k]); 5738 }; 5739 if (iDim == dim(primes[k])) { 5740 equiRad = intersect(equiRad, primes[k]); 5741 }; 5742 }; 5743 }; 5744 } else { 5745 stop = 1 5746 }; 5747 }; 5748 if (il == 0) { 5749 return(rad); 5750 } else { 5751 return(equiRad); 5752 }; 5753 }; 5856 ideal rad = 1; 5857 ideal equiRad = 1; 5858 list primes; 5859 int k; // Counter 5860 int il; // If il = 1, only the equiradical is required. 5861 int iDim; // The dimension of I 5862 int stop = 0; // Checks if the radical has been obtained 5863 5864 if (attrib(I, "isSB") != 1) 5865 { 5866 I = groebner(I); 5867 } 5868 iDim = dim(I); 5869 5870 // Checks if only equiradical is required 5871 if (size(#) > 0) 5872 { 5873 il = #[1]; 5874 } 5875 5876 while(stop == 0) 5877 { 5878 dbprint (printlevel-voice, "// We call radLoopR to find new prime ideals."); 5879 primes = newRadicalSLIteration(I, rad); // A list of primes or intersections of primes, not included in P 5880 dbprint (printlevel - voice, "// Output of Iteration Step:"); 5881 dbprint (printlevel - voice, primes); 5882 if (size(primes) > 0) 5883 { 5884 dbprint (printlevel - voice, "// We intersect P with the ideal just obtained."); 5885 for(k = 1; k <= size(primes); k++) 5886 { 5887 rad = intersect(rad, primes[k]); 5888 if (il == 1) 5889 { 5890 if (attrib(primes[k], "isSB") != 1) 5891 { 5892 primes[k] = groebner(primes[k]); 5893 } 5894 if (iDim == dim(primes[k])) 5895 { 5896 equiRad = intersect(equiRad, primes[k]); 5897 } 5898 } 5899 } 5900 } 5901 else 5902 { 5903 stop = 1; 5904 } 5905 } 5906 if (il == 0) 5907 { 5908 return(rad); 5909 } 5910 else 5911 { 5912 return(equiRad); 5913 } 5914 } 5754 5915 5755 5916 ////////////////////////////////////////////////////////////////////////// … … 5766 5927 // This is not used in the new algorithm. It is part of KL algorithm 5767 5928 5768 static proc newRadicalReduction(ideal I, ideal ser, int allIndep, list #) { 5929 static proc newRadicalReduction(ideal I, ideal ser, int allIndep, list #) 5930 { 5769 5931 // allMaximal 1 -> Indicates that the reduction to the zerodim case 5770 5932 // must be done for all indep set of the leading terms ideal … … 5774 5936 // It is used to set the value of done.) 5775 5937 5776 5777 5778 5779 5780 5781 5782 5783 5784 5785 5786 5787 5788 5789 5790 5791 5792 5793 5794 5795 5796 5797 5798 5799 5800 5801 5802 5803 5804 5938 attrib(I, "isSB", 1); // I needs to be a reduced standard basis 5939 list indep, fett; 5940 intvec @w, @hilb, op; 5941 int @wr, @n, @m, lauf, di; 5942 ideal fac, @h, collectrad, lsau; 5943 poly @q; 5944 string @va, quotring; 5945 5946 def @P = basering; 5947 int jdim = dim(I); // Computes the dimension of I 5948 int homo = homog(I); // Finds out if I is homogeneous 5949 ideal rad = ideal(1); // The unit ideal 5950 ideal te = ser; 5951 if(size(#) > 0) 5952 { 5953 @wr = #[1]; 5954 } 5955 if(homo == 1) 5956 { 5957 for(@n = 1; @n <= nvars(basering); @n++) 5958 { 5959 @w[@n] = ord(var(@n)); 5960 } 5961 @hilb = hilb(I, 1, @w); 5962 } 5963 5964 // SL 2006.04.11 1 Debug messages 5965 dbprint(printlevel-voice, "//Computes the radical of the ideal:", I); 5966 // SL 2006.04.11 2 Debug messages 5805 5967 5806 5968 //--------------------------------------------------------------------------- … … 5808 5970 //--------------------------------------------------------------------------- 5809 5971 5810 5811 5812 5813 5972 if (jdim==-1) 5973 { 5974 return(ideal(1), ideal(1), 1); 5975 } 5814 5976 5815 5977 //--------------------------------------------------------------------------- … … 5817 5979 //--------------------------------------------------------------------------- 5818 5980 5819 5820 5821 5822 5823 5824 5825 5826 5827 5828 5829 5830 5831 5832 5833 5834 5835 5836 5837 5838 5839 5840 5841 5842 5843 5844 5845 5981 if (jdim==0) 5982 { 5983 return(zeroRad(I), ideal(1), 1); 5984 } 5985 5986 //------------------------------------------------------------------------- 5987 //search for a maximal independent set indep,i.e. 5988 //look for subring such that the intersection with the ideal is zero 5989 //j intersected with K[var(indep[3]+1),...,var(nvar)] is zero, 5990 //indep[1] is the new varstring, indep[2] the string for the block-ordering 5991 //------------------------------------------------------------------------- 5992 5993 // SL 2006-04-24 1 If allIndep = 0, then it only computes one maximal 5994 // independent set. 5995 // This looks better for the new algorithm but not for KL 5996 // algorithm 5997 list parameters = allIndep; 5998 indep = newMaxIndependSet(I, parameters); 5999 // SL 2006-04-24 2 6000 6001 for(@m = 1; @m <= size(indep); @m++) 6002 { 6003 if((indep[@m][1] == varstr(basering)) && (@m == 1)) 6004 //this is the good case, nothing to do, just to have the same notations 6005 //change the ring 6006 { 6007 execute("ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),(" 5846 6008 +ordstr(basering)+");"); 5847 ideal @j = fetch(@P, I); 5848 attrib(@j, "isSB", 1); 6009 ideal @j = fetch(@P, I); 6010 attrib(@j, "isSB", 1); 6011 } 6012 else 6013 { 6014 @va = string(maxideal(1)); 6015 6016 execute("ring gnir1 = (" + charstr(basering) + "), (" + indep[@m][1] + "),(" 6017 + indep[@m][2] + ");"); 6018 execute("map phi = @P," + @va + ";"); 6019 if(homo == 1) 6020 { 6021 ideal @j = std(phi(I), @hilb, @w); 5849 6022 } 5850 6023 else 5851 6024 { 5852 @va = string(maxideal(1)); 5853 5854 execute("ring gnir1 = (" + charstr(basering) + "), (" + indep[@m][1] + "),(" 5855 + indep[@m][2] + ");"); 5856 execute("map phi = @P," + @va + ";"); 5857 if(homo == 1) 5858 { 5859 ideal @j = std(phi(I), @hilb, @w); 5860 } 5861 else 5862 { 5863 ideal @j = groebner(phi(I)); 5864 } 5865 } 5866 if((deg(@j[1]) == 0) || (dim(@j) < jdim)) 5867 { 5868 setring @P; 5869 break; 5870 } 5871 for (lauf = 1; lauf <= size(@j); lauf++) 5872 { 5873 fett[lauf] = size(@j[lauf]); 5874 } 5875 //------------------------------------------------------------------------ 5876 // We have now the following situation: 5877 // j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass 5878 // to this quotientring, j is there still a standardbasis, the 5879 // leading coefficients of the polynomials there (polynomials in 5880 // K[var(nnp+1),..,var(nva)]) are collected in the list h, 5881 // we need their LCM, gh, because of the following: 5882 // let (j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..rest..] 5883 // intersected with K[var(1),...,var(nva)] is (j:gh^n) 5884 // on the other hand j = ((j, gh^n) intersected with (j : gh^n)) 5885 5886 //------------------------------------------------------------------------ 5887 // The arrangement for the quotientring K(var(nnp+1),..,var(nva))[..rest..] 5888 // and the map phi:K[var(1),...,var(nva)] -----> 5889 // K(var(nnpr+1),..,var(nva))[..the rest..] 5890 //------------------------------------------------------------------------ 5891 quotring = newPrepareQuotientRing(nvars(basering) - indep[@m][3]); 5892 //------------------------------------------------------------------------ 5893 // We pass to the quotientring K(var(nnp+1),..,var(nva))[..the rest..] 5894 //------------------------------------------------------------------------ 5895 5896 execute(quotring); 5897 5898 // @j considered in the quotientring 5899 ideal @j = imap(gnir1, @j); 5900 5901 kill gnir1; 5902 5903 // j is a standardbasis in the quotientring but usually not minimal 5904 // here it becomes minimal 5905 5906 @j = clearSB(@j, fett); 5907 5908 // We need later LCM(h[1],...) = gh for saturation 5909 ideal @h; 5910 if(deg(@j[1]) > 0) 5911 { 5912 for(@n = 1; @n <= size(@j); @n++) 5913 { 5914 @h[@n] = leadcoef(@j[@n]); 5915 } 5916 op = option(get); 5917 option(redSB); 5918 @j = interred(@j); //to obtain a reduced standardbasis 5919 attrib(@j, "isSB", 1); 5920 option(set, op); 5921 5922 // SL 1 Debug messages 5923 dbprint(printlevel - voice, "zero_rad", basering, @j, dim(groebner(@j))); 5924 ideal zero_rad = zeroRad(@j); 5925 dbprint(printlevel - voice, "zero_rad passed"); 5926 // SL 2 5927 } 5928 else 5929 { 5930 ideal zero_rad = ideal(1); 5931 } 5932 5933 // We need the intersection of the ideals in the list quprimary with the 5934 // polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal 5935 // but fi polynomials, then the intersection of q with the polynomialring 5936 // is the saturation of the ideal generated by f1,...,fr with respect to 5937 // h which is the lcm of the leading coefficients of the fi considered in 6025 ideal @j = groebner(phi(I)); 6026 } 6027 } 6028 if((deg(@j[1]) == 0) || (dim(@j) < jdim)) 6029 { 6030 setring @P; 6031 break; 6032 } 6033 for (lauf = 1; lauf <= size(@j); lauf++) 6034 { 6035 fett[lauf] = size(@j[lauf]); 6036 } 6037 //------------------------------------------------------------------------ 6038 // We have now the following situation: 6039 // j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass 6040 // to this quotientring, j is there still a standardbasis, the 6041 // leading coefficients of the polynomials there (polynomials in 6042 // K[var(nnp+1),..,var(nva)]) are collected in the list h, 6043 // we need their LCM, gh, because of the following: 6044 // let (j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..rest..] 6045 // intersected with K[var(1),...,var(nva)] is (j:gh^n) 6046 // on the other hand j = ((j, gh^n) intersected with (j : gh^n)) 6047 6048 //------------------------------------------------------------------------ 6049 // The arrangement for the quotientring K(var(nnp+1),..,var(nva))[..rest..] 6050 // and the map phi:K[var(1),...,var(nva)] -----> 6051 // K(var(nnpr+1),..,var(nva))[..the rest..] 6052 //------------------------------------------------------------------------ 6053 quotring = newPrepareQuotientRing(nvars(basering) - indep[@m][3]); 6054 //------------------------------------------------------------------------ 6055 // We pass to the quotientring K(var(nnp+1),..,var(nva))[..the rest..] 6056 //------------------------------------------------------------------------ 6057 6058 execute(quotring); 6059 6060 // @j considered in the quotientring 6061 ideal @j = imap(gnir1, @j); 6062 6063 // j is a standardbasis in the quotientring but usually not minimal 6064 // here it becomes minimal 6065 6066 @j = clearSB(@j, fett); 6067 6068 // We need later LCM(h[1],...) = gh for saturation 6069 ideal @h; 6070 if(deg(@j[1]) > 0) 6071 { 6072 for(@n = 1; @n <= size(@j); @n++) 6073 { 6074 @h[@n] = leadcoef(@j[@n]); 6075 } 6076 op = option(get); 6077 option(redSB); 6078 @j = interred(@j); //to obtain a reduced standardbasis 6079 attrib(@j, "isSB", 1); 6080 option(set, op); 6081 6082 // SL 1 Debug messages 6083 dbprint(printlevel - voice, "zero_rad", basering, @j, dim(groebner(@j))); 6084 ideal zero_rad = zeroRad(@j); 6085 dbprint(printlevel - voice, "zero_rad passed"); 6086 // SL 2 6087 } 6088 else 6089 { 6090 ideal zero_rad = ideal(1); 6091 } 6092 6093 // We need the intersection of the ideals in the list quprimary with the 6094 // polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal 6095 // but fi polynomials, then the intersection of q with the polynomialring 6096 // is the saturation of the ideal generated by f1,...,fr with respect to 6097 // h which is the lcm of the leading coefficients of the fi considered in 5938 6098 // the quotientring: this is coded in saturn 5939 6099 … … 6026 6186 } else { 6027 6187 done = 0; 6028 } ;6188 } 6029 6189 // SL 2006.04.21 2 6030 6190 … … 6033 6193 return(result); 6034 6194 // SL 2006.04.21 2 6035 } ;6195 } 6036 6196 6037 6197 … … 6055 6215 good = 1 - rad_con(P[k], I); 6056 6216 k++; 6057 } ;6217 } 6058 6218 k--; 6059 6219 if (good == 0) { … … 6061 6221 list emptyList = list(); 6062 6222 return (emptyList); 6063 } ;6223 } 6064 6224 dbprint(printlevel - voice, "// That one was good!"); 6065 6225 dbprint(printlevel - voice, "// We saturate I with respect to this element."); … … 6069 6229 dbprint(printlevel - voice, "// The polynomial is 1, the saturation in not actually computed."); 6070 6230 ideal J = I; 6071 } ;6231 } 6072 6232 6073 6233 // We now call proc radicalNew; … … 6082 6242 6083 6243 return(result[1]); 6084 } ;6244 } 6085 6245 6086 6246 … … 6118 6278 } else { 6119 6279 allMaximal = 1; 6120 } ;6280 } 6121 6281 6122 6282 int nMax; … … 6125 6285 } else { 6126 6286 nMax = 1; 6127 } ;6287 } 6128 6288 6129 6289 for(n = 1; n <= nMax; n++) … … 6455 6615 pr; 6456 6616 } 6617 /////////////////////////////////////////////////////////////////////////////// 6618 static proc newDecompStep(ideal i, list #) 6619 "USAGE: newDecompStep(i); i ideal (for primary decomposition) 6620 newDecompStep(i,1); (for the associated primes of dimension of i) 6621 newDecompStep(i,2); (for the minimal associated primes) 6622 newDecompStep(i,3); (for the absolute primary decomposition) 6623 "oneIndep"; (for using only one max indep set) 6624 "intersect"; (returns alse the intersection of the components founded) 6625 6626 RETURN: list = list of primary ideals and their associated primes 6627 (at even positions in the list) 6628 (resp. a list of the minimal associated primes) 6629 NOTE: Algorithm of Gianni/Trager/Zacharias 6630 EXAMPLE: example newDecompStep; shows an example 6631 " 6632 { 6633 intvec op,@vv; 6634 def @P = basering; 6635 list primary,indep,ltras; 6636 intvec @vh,isat,@w; 6637 int @wr,@k,@n,@m,@n1,@n2,@n3,homo,seri,keepdi,abspri,ab,nn; 6638 ideal peek=i; 6639 ideal ser,tras; 6640 list data; 6641 list result; 6642 intvec @hilb; 6643 int isS=(attrib(i,"isSB")==1); 6644 6645 // Debug 6646 dbprint(printlevel - voice, "newDecompStep, v2.0"); 6647 6648 string indepOption = "allIndep"; 6649 string intersectOption = "noIntersect"; 6650 6651 if(size(#)>0) 6652 { 6653 int count = 1; 6654 if(typeof(#[count]) == "string") { 6655 if ((#[count] == "oneIndep") or (#[count] == "allIndep")){ 6656 indepOption = #[count]; 6657 count++; 6658 } 6659 } 6660 if(typeof(#[count]) == "string") { 6661 if ((#[count] == "intersect") or (#[count] == "noIntersect")){ 6662 intersectOption = #[count]; 6663 count++; 6664 } 6665 } 6666 if((typeof(#[count]) == "int") or (typeof(#[count]) == "number")) 6667 { 6668 if ((#[count]==1)||(#[count]==2)||(#[count]==3)) 6669 { 6670 @wr=#[count]; 6671 if(@wr==3){abspri = 1; @wr = 0;} 6672 count++; 6673 } 6674 } 6675 if(size(#)>count) 6676 { 6677 seri=1; 6678 peek=#[count + 1]; 6679 ser=#[count + 2]; 6680 } 6681 } 6682 if(abspri) 6683 { 6684 list absprimary,abskeep,absprimarytmp,abskeeptmp; 6685 } 6686 homo=homog(i); 6687 if(homo==1) 6688 { 6689 if(attrib(i,"isSB")!=1) 6690 { 6691 //ltras=mstd(i); 6692 tras=groebner(i); 6693 ltras=tras,tras; 6694 attrib(ltras[1],"isSB",1); 6695 } 6696 else 6697 { 6698 ltras=i,i; 6699 attrib(ltras[1],"isSB",1); 6700 } 6701 tras = ltras[1]; 6702 attrib(tras,"isSB",1); 6703 if(dim(tras)==0) 6704 { 6705 primary[1]=ltras[2]; 6706 primary[2]=maxideal(1); 6707 if(@wr>0) 6708 { 6709 list l; 6710 l[1]=maxideal(1); 6711 l[2]=maxideal(1); 6712 if (intersectOption == "intersect") { 6713 return(list(l, maxideal(1))); 6714 } else { 6715 return(l); 6716 } 6717 } 6718 if (intersectOption == "intersect") { 6719 return(list(primary, primary[1])); 6720 } else { 6721 return(primary); 6722 } 6723 6724 } 6725 for(@n=1;@n<=nvars(basering);@n++) 6726 { 6727 @w[@n]=ord(var(@n)); 6728 } 6729 @hilb=hilb(tras,1,@w); 6730 intvec keephilb=@hilb; 6731 } 6732 6733 //---------------------------------------------------------------- 6734 //i is the zero-ideal 6735 //---------------------------------------------------------------- 6736 6737 if(size(i)==0) 6738 { 6739 primary=i,i; 6740 if (intersectOption == "intersect") { 6741 return(list(primary, i)); 6742 } else { 6743 return(primary); 6744 } 6745 } 6746 6747 //---------------------------------------------------------------- 6748 //pass to the lexicographical ordering and compute a standardbasis 6749 //---------------------------------------------------------------- 6750 6751 int lp=islp(); 6752 6753 execute("ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),(C,lp);"); 6754 op=option(get); 6755 option(redSB); 6756 6757 ideal ser=fetch(@P,ser); 6758 if(homo==1) 6759 { 6760 if(!lp) 6761 { 6762 ideal @j=std(fetch(@P,i),@hilb,@w); 6763 } 6764 else 6765 { 6766 ideal @j=fetch(@P,tras); 6767 attrib(@j,"isSB",1); 6768 } 6769 } 6770 else 6771 { 6772 if(lp&&isS) 6773 { 6774 ideal @j=fetch(@P,i); 6775 attrib(@j,"isSB",1); 6776 } 6777 else 6778 { 6779 ideal @j=groebner(fetch(@P,i)); 6780 } 6781 } 6782 option(set,op); 6783 if(seri==1) 6784 { 6785 ideal peek=fetch(@P,peek); 6786 attrib(peek,"isSB",1); 6787 } 6788 else 6789 { 6790 ideal peek=@j; 6791 } 6792 if((size(ser)==0)&&(!abspri)) 6793 { 6794 ideal fried; 6795 @n=size(@j); 6796 for(@k=1;@k<=@n;@k++) 6797 { 6798 if(deg(lead(@j[@k]))==1) 6799 { 6800 fried[size(fried)+1]=@j[@k]; 6801 @j[@k]=0; 6802 } 6803 } 6804 if(size(fried)==nvars(basering)) 6805 { 6806 setring @P; 6807 primary[1]=i; 6808 primary[2]=i; 6809 if (intersectOption == "intersect") { 6810 return(list(primary, i)); 6811 } else { 6812 return(primary); 6813 } 6814 } 6815 if(size(fried)>0) 6816 { 6817 string newva; 6818 string newma; 6819 for(@k=1;@k<=nvars(basering);@k++) 6820 { 6821 @n1=0; 6822 for(@n=1;@n<=size(fried);@n++) 6823 { 6824 if(leadmonom(fried[@n])==var(@k)) 6825 { 6826 @n1=1; 6827 break; 6828 } 6829 } 6830 if(@n1==0) 6831 { 6832 newva=newva+string(var(@k))+","; 6833 newma=newma+string(var(@k))+","; 6834 } 6835 else 6836 { 6837 newma=newma+string(0)+","; 6838 } 6839 } 6840 newva[size(newva)]=")"; 6841 newma[size(newma)]=";"; 6842 execute("ring @deirf=("+charstr(gnir)+"),("+newva+",lp;"); 6843 execute("map @kappa=gnir,"+newma); 6844 ideal @j= @kappa(@j); 6845 @j=simplify(@j, 2); 6846 attrib(@j,"isSB",1); 6847 result = newDecompStep(@j, indepOption, intersectOption, @wr); 6848 if (intersectOption == "intersect") { 6849 list pr = result[1]; 6850 ideal intersection = result[2]; 6851 } else { 6852 list pr = result; 6853 } 6854 6855 setring gnir; 6856 list pr=imap(@deirf,pr); 6857 for(@k=1;@k<=size(pr);@k++) 6858 { 6859 @j=pr[@k]+fried; 6860 pr[@k]=@j; 6861 } 6862 if (intersectOption == "intersect") { 6863 ideal intersection = imap(@deirf, intersection); 6864 @j = intersection + fried; 6865 intersection = @j; 6866 } 6867 setring @P; 6868 if (intersectOption == "intersect") { 6869 return(list(imap(gnir,pr), imap(gnir,intersection))); 6870 } else { 6871 return(imap(gnir,pr)); 6872 } 6873 } 6874 } 6875 //---------------------------------------------------------------- 6876 //j is the ring 6877 //---------------------------------------------------------------- 6878 6879 if (dim(@j)==-1) 6880 { 6881 setring @P; 6882 primary=ideal(1),ideal(1); 6883 if (intersectOption == "intersect") { 6884 return(list(primary, ideal(1))); 6885 } else { 6886 return(primary); 6887 } 6888 } 6889 6890 //---------------------------------------------------------------- 6891 // the case of one variable 6892 //---------------------------------------------------------------- 6893 6894 if(nvars(basering)==1) 6895 { 6896 6897 list fac=factor(@j[1]); 6898 list gprimary; 6899 poly generator; 6900 ideal gIntersection; 6901 for(@k=1;@k<=size(fac[1]);@k++) 6902 { 6903 if(@wr==0) 6904 { 6905 gprimary[2*@k-1]=ideal(fac[1][@k]^fac[2][@k]); 6906 gprimary[2*@k]=ideal(fac[1][@k]); 6907 } 6908 else 6909 { 6910 gprimary[2*@k-1]=ideal(fac[1][@k]); 6911 gprimary[2*@k]=ideal(fac[1][@k]); 6912 } 6913 if (intersectOption == "intersect") { 6914 generator = generator * fac[1][@k]; 6915 } 6916 } 6917 if (intersectOption == "intersect") { 6918 gIntersection = generator; 6919 } 6920 6921 setring @P; 6922 primary=fetch(gnir,gprimary); 6923 if (intersectOption == "intersect") { 6924 ideal intersection = fetch(gnir,gIntersection); 6925 } 6926 6927 //HIER 6928 if(abspri) 6929 { 6930 list resu,tempo; 6931 string absotto; 6932 for(ab=1;ab<=size(primary)/2;ab++) 6933 { 6934 absotto= absFactorize(primary[2*ab][1],77); 6935 tempo=primary[2*ab-1],primary[2*ab],absotto,string(var(nvars(basering))); 6936 resu[ab]=tempo; 6937 } 6938 primary=resu; 6939 intersection = 1; 6940 for(ab=1;ab<=size(primary);ab++) { 6941 intersection = intersect(intersection, primary[ab][2]); 6942 } 6943 } 6944 if (intersectOption == "intersect") { 6945 return(list(primary, intersection)); 6946 } else { 6947 return(primary); 6948 } 6949 } 6950 6951 //------------------------------------------------------------------ 6952 //the zero-dimensional case 6953 //------------------------------------------------------------------ 6954 if (dim(@j)==0) 6955 { 6956 op=option(get); 6957 option(redSB); 6958 list gprimary= newZero_decomp(@j,ser,@wr); 6959 6960 setring @P; 6961 primary=fetch(gnir,gprimary); 6962 6963 if(size(ser)>0) 6964 { 6965 primary=cleanPrimary(primary); 6966 } 6967 //HIER 6968 if(abspri) 6969 { 6970 list resu,tempo; 6971 string absotto; 6972 for(ab=1;ab<=size(primary)/2;ab++) 6973 { 6974 absotto= absFactorize(primary[2*ab][1],77); 6975 tempo=primary[2*ab-1],primary[2*ab],absotto,string(var(nvars(basering))); 6976 resu[ab]=tempo; 6977 } 6978 primary=resu; 6979 } 6980 if (intersectOption == "intersect") { 6981 return(list(primary, fetch(gnir,@j))); 6982 } else { 6983 return(primary); 6984 } 6985 6986 } 6987 6988 poly @gs,@gh,@p; 6989 string @va,quotring; 6990 list quprimary,htprimary,collectprimary,lsau,lnew,allindep,restindep; 6991 ideal @h; 6992 int jdim=dim(@j); 6993 list fett; 6994 int lauf,di,newtest; 6995 //------------------------------------------------------------------ 6996 //search for a maximal independent set indep,i.e. 6997 //look for subring such that the intersection with the ideal is zero 6998 //j intersected with K[var(indep[3]+1),...,var(nvar] is zero, 6999 //indep[1] is the new varstring and indep[2] the string for block-ordering 7000 //------------------------------------------------------------------ 7001 if(@wr!=1) 7002 { 7003 allindep = newMaxIndependSetLp(@j, indepOption); 7004 for(@m=1;@m<=size(allindep);@m++) 7005 { 7006 if(allindep[@m][3]==jdim) 7007 { 7008 di++; 7009 indep[di]=allindep[@m]; 7010 } 7011 else 7012 { 7013 lauf++; 7014 restindep[lauf]=allindep[@m]; 7015 } 7016 } 7017 } 7018 else 7019 { 7020 indep = newMaxIndependSetLp(@j, indepOption); 7021 } 7022 7023 ideal jkeep=@j; 7024 if(ordstr(@P)[1]=="w") 7025 { 7026 execute("ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),("+ordstr(@P)+");"); 7027 } 7028 else 7029 { 7030 execute( "ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),(C,dp);"); 7031 } 7032 7033 if(homo==1) 7034 { 7035 if((ordstr(@P)[3]=="d")||(ordstr(@P)[1]=="d")||(ordstr(@P)[1]=="w") 7036 ||(ordstr(@P)[3]=="w")) 7037 { 7038 ideal jwork=imap(@P,tras); 7039 attrib(jwork,"isSB",1); 7040 } 7041 else 7042 { 7043 ideal jwork=std(imap(gnir,@j),@hilb,@w); 7044 } 7045 7046 } 7047 else 7048 { 7049 ideal jwork=groebner(imap(gnir,@j)); 7050 } 7051 list hquprimary; 7052 poly @p,@q; 7053 ideal @h,fac,ser; 7054 //Aenderung================ 7055 ideal @Ptest=1; 7056 //========================= 7057 di=dim(jwork); 7058 keepdi=di; 7059 7060 ser = 1; 7061 7062 setring gnir; 7063 for(@m=1; @m<=size(indep); @m++) 7064 { 7065 data[1] = indep[@m]; 7066 result = newReduction(@j, ser, @hilb, @w, jdim, abspri, @wr, data); 7067 quprimary = quprimary + result[1]; 7068 if(abspri) { 7069 absprimary = absprimary + result[2]; 7070 abskeep = abskeep + result[3]; 7071 } 7072 @h = result[5]; 7073 ser = result[4]; 7074 if(size(@h)>0) 7075 { 7076 //--------------------------------------------------------------- 7077 //we change to @Phelp to have the ordering dp for saturation 7078 //--------------------------------------------------------------- 7079 7080 setring @Phelp; 7081 @h=imap(gnir,@h); 7082 //Aenderung================================== 7083 if(defined(@LL)){kill @LL;} 7084 list @LL=minSat(jwork,@h); 7085 @Ptest=intersect(@Ptest,@LL[1]); 7086 ser = intersect(ser, @LL[1]); 7087 //=========================================== 7088 7089 if(@wr!=1) 7090 { 7091 //Aenderung================================== 7092 @q=@LL[2]; 7093 //=========================================== 7094 //@q=minSat(jwork,@h)[2]; 7095 } 7096 else 7097 { 7098 fac=ideal(0); 7099 for(lauf=1;lauf<=ncols(@h);lauf++) 7100 { 7101 if(deg(@h[lauf])>0) 7102 { 7103 fac=fac+factorize(@h[lauf],1); 7104 } 7105 } 7106 fac=simplify(fac,4); 7107 @q=1; 7108 for(lauf=1;lauf<=size(fac);lauf++) 7109 { 7110 @q=@q*fac[lauf]; 7111 } 7112 } 7113 jwork = std(jwork,@q); 7114 keepdi = dim(jwork); 7115 if(keepdi < di) 7116 { 7117 setring gnir; 7118 @j = imap(@Phelp, jwork); 7119 ser = imap(@Phelp, ser); 7120 break; 7121 } 7122 if(homo == 1) 7123 { 7124 @hilb = hilb(jwork, 1, @w); 7125 } 7126 7127 setring gnir; 7128 ser = imap(@Phelp, ser); 7129 @j = imap(@Phelp, jwork); 7130 } 7131 } 7132 7133 if((size(quprimary)==0)&&(@wr==1)) 7134 { 7135 @j=ideal(1); 7136 quprimary[1]=ideal(1); 7137 quprimary[2]=ideal(1); 7138 } 7139 if((size(quprimary)==0)) 7140 { 7141 keepdi = di - 1; 7142 quprimary[1]=ideal(1); 7143 quprimary[2]=ideal(1); 7144 } 7145 //--------------------------------------------------------------- 7146 //notice that j=sat(j,gh) intersected with (j,gh^n) 7147 //we finished with sat(j,gh) and have to start with (j,gh^n) 7148 //--------------------------------------------------------------- 7149 if((deg(@j[1])!=0)&&(@wr!=1)) 7150 { 7151 if(size(quprimary)>0) 7152 { 7153 setring @Phelp; 7154 ser=imap(gnir,ser); 7155 7156 hquprimary=imap(gnir,quprimary); 7157 if(@wr==0) 7158 { 7159 //Aenderung==================================================== 7160 //HIER STATT DURCHSCHNITT SATURIEREN! 7161 ideal htest=@Ptest; 7162 /* 7163 ideal htest=hquprimary[1]; 7164 for (@n1=2;@n1<=size(hquprimary)/2;@n1++) 7165 { 7166 htest=intersect(htest,hquprimary[2*@n1-1]); 7167 } 7168 */ 7169 //============================================================= 7170 } 7171 else 7172 { 7173 ideal htest=hquprimary[2]; 7174 7175 for (@n1=2;@n1<=size(hquprimary)/2;@n1++) 7176 { 7177 htest=intersect(htest,hquprimary[2*@n1]); 7178 } 7179 } 7180 7181 if(size(ser)>0) 7182 { 7183 ser=intersect(htest,ser); 7184 } 7185 else 7186 { 7187 ser=htest; 7188 } 7189 setring gnir; 7190 ser=imap(@Phelp,ser); 7191 } 7192 if(size(reduce(ser,peek,1))!=0) 7193 { 7194 for(@m=1;@m<=size(restindep);@m++) 7195 { 7196 // if(restindep[@m][3]>=keepdi) 7197 // { 7198 isat=0; 7199 @n2=0; 7200 7201 if(restindep[@m][1]==varstr(basering)) 7202 //the good case, nothing to do, just to have the same notations 7203 //change the ring 7204 { 7205 execute("ring gnir1 = ("+charstr(basering)+"),("+ 7206 varstr(basering)+"),("+ordstr(basering)+");"); 7207 ideal @j=fetch(gnir,jkeep); 7208 attrib(@j,"isSB",1); 7209 } 7210 else 7211 { 7212 @va=string(maxideal(1)); 7213 execute("ring gnir1 = ("+charstr(basering)+"),("+ 7214 restindep[@m][1]+"),(" +restindep[@m][2]+");"); 7215 execute("map phi=gnir,"+@va+";"); 7216 op=option(get); 7217 option(redSB); 7218 if(homo==1) 7219 { 7220 ideal @j=std(phi(jkeep),keephilb,@w); 7221 } 7222 else 7223 { 7224 ideal @j=groebner(phi(jkeep)); 7225 } 7226 ideal ser=phi(ser); 7227 option(set,op); 7228 } 7229 7230 for (lauf=1;lauf<=size(@j);lauf++) 7231 { 7232 fett[lauf]=size(@j[lauf]); 7233 } 7234 //------------------------------------------------------------------ 7235 //we have now the following situation: 7236 //j intersected with K[var(nnp+1),..,var(nva)] is zero so we may 7237 //pass to this quotientring, j is their still a standardbasis, the 7238 //leading coefficients of the polynomials there (polynomials in 7239 //K[var(nnp+1),..,var(nva)]) are collected in the list h, 7240 //we need their ggt, gh, because of the following: 7241 //let (j:gh^n)=(j:gh^infinity) then 7242 //j*K(var(nnp+1),..,var(nva))[..the rest..] 7243 //intersected with K[var(1),...,var(nva)] is (j:gh^n) 7244 //on the other hand j=(j,gh^n) intersected with (j:gh^n) 7245 7246 //------------------------------------------------------------------ 7247 7248 //the arrangement for the quotientring 7249 // K(var(nnp+1),..,var(nva))[..the rest..] 7250 //and the map phi:K[var(1),...,var(nva)] ----> 7251 //--->K(var(nnpr+1),..,var(nva))[..the rest..] 7252 //------------------------------------------------------------------ 7253 7254 quotring=prepareQuotientring(nvars(basering)-restindep[@m][3]); 7255 7256 //------------------------------------------------------------------ 7257 //we pass to the quotientring K(var(nnp+1),..,var(nva))[..rest..] 7258 //------------------------------------------------------------------ 7259 7260 execute(quotring); 7261 7262 // @j considered in the quotientring 7263 ideal @j=imap(gnir1,@j); 7264 ideal ser=imap(gnir1,ser); 7265 7266 kill gnir1; 7267 7268 //j is a standardbasis in the quotientring but usually not minimal 7269 //here it becomes minimal 7270 @j=clearSB(@j,fett); 7271 attrib(@j,"isSB",1); 7272 7273 //we need later ggt(h[1],...)=gh for saturation 7274 ideal @h; 7275 7276 for(@n=1;@n<=size(@j);@n++) 7277 { 7278 @h[@n]=leadcoef(@j[@n]); 7279 } 7280 //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..rest..] 7281 7282 op=option(get); 7283 option(redSB); 7284 list uprimary= newZero_decomp(@j,ser,@wr); 7285 //HIER 7286 if(abspri) 7287 { 7288 ideal II; 7289 ideal jmap; 7290 map sigma; 7291 nn=nvars(basering); 7292 map invsigma=basering,maxideal(1); 7293 for(ab=1;ab<=size(uprimary)/2;ab++) 7294 { 7295 II=uprimary[2*ab]; 7296 attrib(II,"isSB",1); 7297 if(deg(II[1])!=vdim(II)) 7298 { 7299 jmap=randomLast(50); 7300 sigma=basering,jmap; 7301 jmap[nn]=2*var(nn)-jmap[nn]; 7302 invsigma=basering,jmap; 7303 II=groebner(sigma(II)); 7304 } 7305 absprimarytmp[ab]= absFactorize(II[1],77); 7306 II=var(nn); 7307 abskeeptmp[ab]=string(invsigma(II)); 7308 invsigma=basering,maxideal(1); 7309 } 7310 } 7311 option(set,op); 7312 7313 //we need the intersection of the ideals in the list quprimary with 7314 //the polynomialring, i.e. let q=(f1,...,fr) in the quotientring 7315 //such an ideal but fi polynomials, then the intersection of q with 7316 //the polynomialring is the saturation of the ideal generated by 7317 //f1,...,fr with respect toh which is the lcm of the leading 7318 //coefficients of the fi considered in the quotientring: 7319 //this is coded in saturn 7320 7321 list saturn; 7322 ideal hpl; 7323 7324 for(@n=1;@n<=size(uprimary);@n++) 7325 { 7326 hpl=0; 7327 for(@n1=1;@n1<=size(uprimary[@n]);@n1++) 7328 { 7329 hpl=hpl,leadcoef(uprimary[@n][@n1]); 7330 } 7331 saturn[@n]=hpl; 7332 } 7333 //------------------------------------------------------------------ 7334 //we leave the quotientring K(var(nnp+1),..,var(nva))[..rest..] 7335 //back to the polynomialring 7336 //------------------------------------------------------------------ 7337 setring gnir; 7338 collectprimary=imap(quring,uprimary); 7339 lsau=imap(quring,saturn); 7340 @h=imap(quring,@h); 7341 7342 kill quring; 7343 7344 7345 @n2=size(quprimary); 7346 //================NEU========================================= 7347 if(deg(quprimary[1][1])<=0){ @n2=0; } 7348 //============================================================ 7349 7350 @n3=@n2; 7351 7352 for(@n1=1;@n1<=size(collectprimary)/2;@n1++) 7353 { 7354 if(deg(collectprimary[2*@n1][1])>0) 7355 { 7356 @n2++; 7357 quprimary[@n2]=collectprimary[2*@n1-1]; 7358 lnew[@n2]=lsau[2*@n1-1]; 7359 @n2++; 7360 lnew[@n2]=lsau[2*@n1]; 7361 quprimary[@n2]=collectprimary[2*@n1]; 7362 if(abspri) 7363 { 7364 absprimary[@n2/2]=absprimarytmp[@n1]; 7365 abskeep[@n2/2]=abskeeptmp[@n1]; 7366 } 7367 } 7368 } 7369 7370 7371 //here the intersection with the polynomialring 7372 //mentioned above is really computed 7373 7374 for(@n=@n3/2+1;@n<=@n2/2;@n++) 7375 { 7376 if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n])) 7377 { 7378 quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1]; 7379 quprimary[2*@n]=quprimary[2*@n-1]; 7380 } 7381 else 7382 { 7383 if(@wr==0) 7384 { 7385 quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1]; 7386 } 7387 quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1]; 7388 } 7389 } 7390 if(@n2>=@n3+2) 7391 { 7392 setring @Phelp; 7393 ser=imap(gnir,ser); 7394 hquprimary=imap(gnir,quprimary); 7395 for(@n=@n3/2+1;@n<=@n2/2;@n++) 7396 { 7397 if(@wr==0) 7398 { 7399 ser=intersect(ser,hquprimary[2*@n-1]); 7400 } 7401 else 7402 { 7403 ser=intersect(ser,hquprimary[2*@n]); 7404 } 7405 } 7406 setring gnir; 7407 ser=imap(@Phelp,ser); 7408 } 7409 7410 // } 7411 } 7412 //HIER 7413 if(abspri) 7414 { 7415 list resu,tempo; 7416 for(ab=1;ab<=size(quprimary)/2;ab++) 7417 { 7418 if (deg(quprimary[2*ab][1])!=0) 7419 { 7420 tempo=quprimary[2*ab-1],quprimary[2*ab], 7421 absprimary[ab],abskeep[ab]; 7422 resu[ab]=tempo; 7423 } 7424 } 7425 quprimary=resu; 7426 @wr=3; 7427 } 7428 if(size(reduce(ser,peek,1))!=0) 7429 { 7430 if(@wr>0) 7431 { 7432 // The following line was dropped to avoid the recursion step: 7433 //htprimary=newDecompStep(@j,@wr,peek,ser); 7434 htprimary = list(); 7435 } 7436 else 7437 { 7438 // The following line was dropped to avoid the recursion step: 7439 //htprimary=newDecompStep(@j,peek,ser); 7440 htprimary = list(); 7441 } 7442 // here we collect now both results primary(sat(j,gh)) 7443 // and primary(j,gh^n) 7444 @n=size(quprimary); 7445 if (deg(quprimary[1][1])<=0) { @n=0; } 7446 for (@k=1;@k<=size(htprimary);@k++) 7447 { 7448 quprimary[@n+@k]=htprimary[@k]; 7449 } 7450 } 7451 } 7452 } 7453 else 7454 { 7455 if(abspri) 7456 { 7457 list resu,tempo; 7458 for(ab=1;ab<=size(quprimary)/2;ab++) 7459 { 7460 tempo=quprimary[2*ab-1],quprimary[2*ab], 7461 absprimary[ab],abskeep[ab]; 7462 resu[ab]=tempo; 7463 } 7464 quprimary=resu; 7465 } 7466 } 7467 //--------------------------------------------------------------------------- 7468 //back to the ring we started with 7469 //the final result: primary 7470 //--------------------------------------------------------------------------- 7471 7472 setring @P; 7473 primary=imap(gnir,quprimary); 7474 7475 if (intersectOption == "intersect") { 7476 return(list(primary, imap(gnir, ser))); 7477 } else { 7478 return(primary); 7479 } 7480 } 7481 example 7482 { "EXAMPLE:"; echo = 2; 7483 ring r = 32003,(x,y,z),lp; 7484 poly p = z2+1; 7485 poly q = z4+2; 7486 ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4; 7487 list pr= newDecompStep(i); 7488 pr; 7489 testPrimary( pr, i); 7490 } 7491 7492 // 7493 proc newReduction(ideal @j, ideal ser, intvec @hilb, intvec @w, int jdim, int abspri, int @wr, list data) 7494 { 7495 string @va; 7496 string quotring; 7497 intvec op; 7498 intvec @vv; 7499 def gnir = basering; 7500 ideal isat=0; 7501 int @n; 7502 int @n1 = 0; 7503 int @n2 = 0; 7504 int @n3 = 0; 7505 int homo = homog(@j); 7506 int lauf; 7507 int @k; 7508 list fett; 7509 int keepdi; 7510 list collectprimary; 7511 list lsau; 7512 list lnew; 7513 ideal @h; 7514 7515 list indepInfo = data[1]; 7516 list quprimary = list(); 7517 7518 7519 //if(abspri) 7520 //{ 7521 int ab; 7522 list absprimarytmp,abskeeptmp; 7523 list absprimary, abskeep; 7524 //} 7525 // Debug 7526 dbprint(printlevel - voice, "newReduction, v2.0"); 7527 7528 if((indepInfo[1]==varstr(basering))) // &&(@m==1) 7529 //this is the good case, nothing to do, just to have the same notations 7530 //change the ring 7531 { 7532 execute("ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),(" 7533 +ordstr(basering)+");"); 7534 ideal @j = fetch(gnir, @j); 7535 attrib(@j,"isSB",1); 7536 ideal ser = fetch(gnir, ser); 7537 } 7538 else 7539 { 7540 @va=string(maxideal(1)); 7541 //Aenderung============== 7542 //if(@m==1) 7543 //{ 7544 //@j=fetch(@P,i); 7545 //} 7546 //======================= 7547 execute("ring gnir1 = ("+charstr(basering)+"),("+indepInfo[1]+"),(" 7548 +indepInfo[2]+");"); 7549 execute("map phi=gnir,"+@va+";"); 7550 op=option(get); 7551 option(redSB); 7552 if(homo==1) 7553 { 7554 ideal @j=std(phi(@j),@hilb,@w); 7555 } 7556 else 7557 { 7558 ideal @j=groebner(phi(@j)); 7559 } 7560 ideal ser=phi(ser); 7561 7562 option(set,op); 7563 } 7564 if((deg(@j[1])==0)||(dim(@j)<jdim)) 7565 { 7566 setring gnir; 7567 break; 7568 } 7569 for (lauf=1;lauf<=size(@j);lauf++) 7570 { 7571 fett[lauf]=size(@j[lauf]); 7572 } 7573 //------------------------------------------------------------------------ 7574 //we have now the following situation: 7575 //j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass 7576 //to this quotientring, j is their still a standardbasis, the 7577 //leading coefficients of the polynomials there (polynomials in 7578 //K[var(nnp+1),..,var(nva)]) are collected in the list h, 7579 //we need their ggt, gh, because of the following: let 7580 //(j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..the rest..] 7581 //intersected with K[var(1),...,var(nva)] is (j:gh^n) 7582 //on the other hand j=(j,gh^n) intersected with (j:gh^n) 7583 7584 //------------------------------------------------------------------------ 7585 7586 //arrangement for quotientring K(var(nnp+1),..,var(nva))[..the rest..] and 7587 //map phi:K[var(1),...,var(nva)] --->K(var(nnpr+1),..,var(nva))[..rest..] 7588 //------------------------------------------------------------------------ 7589 7590 quotring=prepareQuotientring(nvars(basering)-indepInfo[3]); 7591 7592 //--------------------------------------------------------------------- 7593 //we pass to the quotientring K(var(nnp+1),..,var(nva))[..the rest..] 7594 //--------------------------------------------------------------------- 7595 7596 ideal @jj=lead(@j); //!! vorn vereinbaren 7597 execute(quotring); 7598 7599 ideal @jj=imap(gnir1,@jj); 7600 @vv=clearSBNeu(@jj,fett); //!! vorn vereinbaren 7601 setring gnir1; 7602 @k=size(@j); 7603 for (lauf=1;lauf<=@k;lauf++) 7604 { 7605 if(@vv[lauf]==1) 7606 { 7607 @j[lauf]=0; 7608 } 7609 } 7610 @j=simplify(@j,2); 7611 setring quring; 7612 // @j considered in the quotientring 7613 ideal @j=imap(gnir1,@j); 7614 7615 ideal ser=imap(gnir1,ser); 7616 7617 kill gnir1; 7618 7619 //j is a standardbasis in the quotientring but usually not minimal 7620 //here it becomes minimal 7621 7622 attrib(@j,"isSB",1); 7623 7624 //we need later ggt(h[1],...)=gh for saturation 7625 ideal @h; 7626 if(deg(@j[1])>0) 7627 { 7628 for(@n=1;@n<=size(@j);@n++) 7629 { 7630 @h[@n]=leadcoef(@j[@n]); 7631 } 7632 //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..] 7633 op=option(get); 7634 option(redSB); 7635 7636 int zeroMinAss = @wr; 7637 if (@wr == 2) {zeroMinAss = 1;} 7638 list uprimary= newZero_decomp(@j, ser, zeroMinAss); 7639 7640 //HIER 7641 if(abspri) 7642 { 7643 ideal II; 7644 ideal jmap; 7645 map sigma; 7646 nn=nvars(basering); 7647 map invsigma=basering,maxideal(1); 7648 for(ab=1;ab<=size(uprimary)/2;ab++) 7649 { 7650 II=uprimary[2*ab]; 7651 attrib(II,"isSB",1); 7652 if(deg(II[1])!=vdim(II)) 7653 { 7654 jmap=randomLast(50); 7655 sigma=basering,jmap; 7656 jmap[nn]=2*var(nn)-jmap[nn]; 7657 invsigma=basering,jmap; 7658 II=groebner(sigma(II)); 7659 } 7660 absprimarytmp[ab]= absFactorize(II[1],77); 7661 II=var(nn); 7662 abskeeptmp[ab]=string(invsigma(II)); 7663 invsigma=basering,maxideal(1); 7664 } 7665 } 7666 option(set,op); 7667 } 7668 else 7669 { 7670 list uprimary; 7671 uprimary[1]=ideal(1); 7672 uprimary[2]=ideal(1); 7673 } 7674 //we need the intersection of the ideals in the list quprimary with the 7675 //polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal 7676 //but fi polynomials, then the intersection of q with the polynomialring 7677 //is the saturation of the ideal generated by f1,...,fr with respect to 7678 //h which is the lcm of the leading coefficients of the fi considered in 7679 //in the quotientring: this is coded in saturn 7680 7681 list saturn; 7682 ideal hpl; 7683 7684 for(@n=1;@n<=size(uprimary);@n++) 7685 { 7686 uprimary[@n]=interred(uprimary[@n]); // temporary fix 7687 hpl=0; 7688 for(@n1=1;@n1<=size(uprimary[@n]);@n1++) 7689 { 7690 hpl=hpl,leadcoef(uprimary[@n][@n1]); 7691 } 7692 saturn[@n]=hpl; 7693 } 7694 7695 //-------------------------------------------------------------------- 7696 //we leave the quotientring K(var(nnp+1),..,var(nva))[..the rest..] 7697 //back to the polynomialring 7698 //--------------------------------------------------------------------- 7699 setring gnir; 7700 7701 collectprimary=imap(quring,uprimary); 7702 lsau=imap(quring,saturn); 7703 @h=imap(quring,@h); 7704 7705 kill quring; 7706 7707 @n2=size(quprimary); 7708 @n3=@n2; 7709 7710 for(@n1=1;@n1<=size(collectprimary)/2;@n1++) 7711 { 7712 if(deg(collectprimary[2*@n1][1])>0) 7713 { 7714 @n2++; 7715 quprimary[@n2]=collectprimary[2*@n1-1]; 7716 lnew[@n2]=lsau[2*@n1-1]; 7717 @n2++; 7718 lnew[@n2]=lsau[2*@n1]; 7719 quprimary[@n2]=collectprimary[2*@n1]; 7720 if(abspri) 7721 { 7722 absprimary[@n2/2]=absprimarytmp[@n1]; 7723 abskeep[@n2/2]=abskeeptmp[@n1]; 7724 } 7725 } 7726 } 7727 7728 //here the intersection with the polynomialring 7729 //mentioned above is really computed 7730 for(@n=@n3/2+1;@n<=@n2/2;@n++) 7731 { 7732 if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n])) 7733 { 7734 quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1]; 7735 quprimary[2*@n]=quprimary[2*@n-1]; 7736 } 7737 else 7738 { 7739 if(@wr==0) 7740 { 7741 quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1]; 7742 } 7743 quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1]; 7744 } 7745 } 7746 7747 return(quprimary, absprimary, abskeep, ser, @h); 7748 } 7749 7750 7751 //////////////////////////////////////////////////////////////////////////// 7752 7753 7754 7755 7756 /////////////////////////////////////////////////////////////////////////////// 7757 // Based on minAssGTZ 7758 7759 proc newMinAss(ideal i,list #) 7760 "USAGE: newMinAss(i); i ideal 7761 newMinAss(i,1); i ideal does not use the factorizing Groebner 7762 RETURN: a list, the minimal associated prime ideals of i. 7763 NOTE: Designed for characteristic 0, works also in char k > 0 based 7764 on an algorithm of Yokoyama 7765 EXAMPLE: example newMinAss; shows an example 7766 " 7767 { 7768 string algorithm = "SL"; 7769 if(ord_test(basering) != 1) 7770 { 7771 ERROR( 7772 "// Not implemented for this ordering, please change to global ordering." 7773 ); 7774 } 7775 if(minpoly!=0) 7776 { 7777 return(algeDeco(i, 2)); 7778 } 7779 if(size(#)==0) 7780 { 7781 return(minAssPrimes(i, "facstd", algorithm)); 7782 } 7783 return(minAssPrimes(i, algorithm)); 7784 } 7785 example 7786 { "EXAMPLE:"; echo = 2; 7787 ring r = 0, (x, y, z), dp; 7788 poly p = z2 + 1; 7789 poly q = z3 + 2; 7790 ideal i = p * q^2, y - z2; 7791 list pr = newMinAss(i); 7792 pr; 7793 } 7794 7795 7796 7797 7798 7799 7800 7801 7802 /////////////////////////////////////////////////////////////////////////////// 7803 // 7804 // Computes the minimal associated primes of I via the new algorithm, 7805 // using primary decomposition in the zero dimensional case. 7806 // For reduction to the zerodimensional case, it uses the procedure 7807 // decomp of primdec.lib, with some modifications to avoid the recursion. 7808 // 7809 7810 proc newMinAssSL(ideal I) 7811 // Input = I, ideal 7812 // Output = primaryDec where primaryDec is the list of the minimal 7813 // associated primes and the primary components corresponding to these primes. 7814 7815 { 7816 ideal P = 1; 7817 list pd = list(); 7818 int k; 7819 int stop = 0; 7820 list primaryDec = list(); 7821 7822 while (stop == 0) { 7823 // Debug 7824 dbprint(printlevel - voice, "// We call newMinAssSLIteration to find new prime ideals!"); 7825 pd = newMinAssSLIteration(I, P); 7826 // Debug 7827 dbprint(printlevel - voice, "// Output of newMinAssSLIteration:"); 7828 dbprint(printlevel - voice, pd); 7829 if (size(pd[1]) > 0) { 7830 primaryDec = primaryDec + pd[1]; 7831 // Debug 7832 dbprint(printlevel - voice, "// We intersect the prime ideals obtained."); 7833 /* 7834 for (k = 1; k <= size(pd); k++) { 7835 P = intersect(P, pd[k]); 7836 } 7837 */ 7838 P = intersect(P, pd[2]); 7839 // Debug 7840 dbprint(printlevel - voice, "// Intersection finished."); 7841 } else { 7842 stop = 1; 7843 } 7844 } 7845 //return(insert(primaryDec, P)); 7846 // Returns only the primary components, not the radical. 7847 return(primaryDec); 7848 } 7849 7850 /////////////////////////////////////////////////////////////////////////////// 7851 // Given an ideal I and an ideal P (intersection of some minimal prime ideals 7852 // associated to I), it calculates new minimal prime ideals associated to I 7853 // which were not used to calculate P. 7854 // This version uses Primary Decomposition in the zerodimensional case. 7855 proc newMinAssSLIteration(ideal I, ideal P); 7856 { 7857 int k = 1; 7858 int good = 0; 7859 list primaryDec = list(); 7860 // Debug 7861 dbprint (printlevel-voice, "// We search for an element in P - sqrt(I)."); 7862 while ((k <= size(P)) and (good == 0)) { 7863 good = 1 - rad_con(P[k], I); 7864 k++; 7865 } 7866 k--; 7867 if (good == 0) { 7868 // Debug 7869 dbprint (printlevel - voice, "// No element was found, P = sqrt(I)."); 7870 return (list(primaryDec, ideal(0))); 7871 } 7872 // Debug 7873 dbprint (printlevel - voice, "// We found h = ", P[k]); 7874 dbprint (printlevel - voice, "// We calculate the saturation of I with respect to the element just founded."); 7875 ideal J = sat(I, P[k])[1]; 7876 7877 // Uses decomp from primdec, modified to avoid the recursion. 7878 // Debug 7879 dbprint(printlevel - voice, "// We do the reduction to the zerodimensional case, via decomp."); 7880 7881 primaryDec = newDecompStep(J, "oneIndep", "intersect", 2); 7882 // Debug 7883 dbprint(printlevel - voice, "// Proc decomp has found", size(primaryDec) / 2, "new primary components."); 7884 7885 return(primaryDec); 7886 } 7887 7888 7889 7890 /////////////////////////////////////////////////////////////////////////////////// 7891 // Based on maxIndependSet 7892 // Added list # as parameter 7893 // If the first element of # is 0, the output is only 1 max indep set. 7894 // If no list is specified or #[1] = 1, the output is all the max indep set of the 7895 // leading terms ideal. This is the original output of maxIndependSet 7896 7897 // The ordering given in the output has been changed to block dp instead of lp. 7898 7899 proc newMaxIndependSetLp(ideal j, list #) 7900 // #[1]=0 --> computes only one max indep set 7901 "USAGE: maxIndependentSet(i); i ideal 7902 RETURN: list = #1. new varstring with the maximal independent set at the end, 7903 #2. ordstring with the corresponding block ordering, 7904 #3. the number of independent variables 7905 // SL #3 explanation modified, 7906 // it was: integer where the independent set starts in the varstring 7907 NOTE: 7908 EXAMPLE: example maxIndependentSetLp; shows an example 7909 " 7910 { 7911 int n, k, di; 7912 list resu, hilf; 7913 string var1, var2; 7914 list v = indepSet(j, 0); 7915 7916 // SL 2006.04.21 1 Lines modified to use only one independent Set 7917 string indepOption; 7918 if (size(#) > 0) { 7919 indepOption = #[1]; 7920 } else { 7921 indepOption = "allIndep"; 7922 } 7923 7924 int nMax; 7925 if (indepOption == "allIndep") { 7926 nMax = size(v); 7927 } else { 7928 nMax = 1; 7929 } 7930 7931 for(n = 1; n <= nMax; n++) 7932 // SL 2006.04.21 2 7933 { 7934 di = 0; 7935 var1 = ""; 7936 var2 = ""; 7937 for(k = 1; k <= size(v[n]); k++) 7938 { 7939 if(v[n][k] != 0) 7940 { 7941 di++; 7942 var2 = var2 + "var(" + string(k) + "), "; 7943 } 7944 else 7945 { 7946 var1 = var1 + "var(" + string(k) + "), "; 7947 } 7948 } 7949 if(di > 0) 7950 { 7951 var1 = var1 + var2; 7952 var1 = var1[1..size(var1) - 2]; // The "- 2" removes the trailer comma 7953 hilf[1] = var1; 7954 // SL 2006.21.04 1 The order is now block dp instead of lp 7955 //hilf[2] = "dp(" + string(nvars(basering) - di) + "), dp(" + string(di) + ")"; 7956 // SL 2006.21.04 2 7957 // For decomp, lp ordering is needed. Nothing is changed. 7958 hilf[2] = "lp"; 7959 hilf[3] = di; 7960 resu[n] = hilf; 7961 } 7962 else 7963 { 7964 resu[n] = varstr(basering), ordstr(basering), 0; 7965 } 7966 } 7967 return(resu); 7968 } 7969 example 7970 { "EXAMPLE:"; echo = 2; 7971 ring s1 = (0, x, y), (a, b, c, d, e, f, g), lp; 7972 ideal i = ea - fbg, fa + be, ec - fdg, fc + de; 7973 i = std(i); 7974 list l = newMaxIndependSetLp(i); 7975 l; 7976 i = i, g; 7977 l = newMaxIndependSetLp(i); 7978 l; 7979 7980 ring s = 0, (x, y, z), lp; 7981 ideal i = z, yx; 7982 list l = newMaxIndependSetLp(i); 7983 l; 7984 } 7985 7986 7987 /////////////////////////////////////////////////////////////////////////////// 7988 7989 proc newZero_decomp (ideal j, ideal ser, int @wr, list #) 7990 "USAGE: newZero_decomp(j,ser,@wr); j,ser ideals, @wr=0 or 1 7991 (@wr=0 for primary decomposition, @wr=1 for computaion of associated 7992 primes) 7993 if #[1] = "nest", then #[2] indicates the nest level (number of recursive calls) 7994 When the nest level is high it indicates that the computation is difficult, 7995 and different methods are applied. 7996 RETURN: list = list of primary ideals and their radicals (at even positions 7997 in the list) if the input is zero-dimensional and a standardbases 7998 with respect to lex-ordering 7999 If ser!=(0) and ser is contained in j or if j is not zero-dimen- 8000 sional then ideal(1),ideal(1) is returned 8001 NOTE: Algorithm of Gianni/Trager/Zacharias 8002 EXAMPLE: example newZero_decomp; shows an example 8003 " 8004 { 8005 def @P = basering; 8006 int uytrewq; 8007 int nva = nvars(basering); 8008 int @k,@s,@n,@k1,zz; 8009 list primary,lres0,lres1,act,@lh,@wh; 8010 map phi,psi,phi1,psi1; 8011 ideal jmap,jmap1,jmap2,helpprim,@qh,@qht,ser1; 8012 intvec @vh,@hilb; 8013 string @ri; 8014 poly @f; 8015 8016 // Debug 8017 dbprint(printlevel - voice, "proc newZero_decomp"); 8018 8019 8020 if (dim(j)>0) 8021 { 8022 primary[1]=ideal(1); 8023 primary[2]=ideal(1); 8024 return(primary); 8025 } 8026 j=interred(j); 8027 8028 attrib(j,"isSB",1); 8029 8030 int nestLevel = 0; 8031 if (size(#) > 0){ 8032 if (typeof(#[1]) == "string") { 8033 if (#[1] == "nest") { 8034 nestLevel = #[2]; 8035 } 8036 # = list(); 8037 } 8038 } 8039 8040 if(vdim(j)==deg(j[1])) 8041 { 8042 act=factor(j[1]); 8043 for(@k=1;@k<=size(act[1]);@k++) 8044 { 8045 @qh=j; 8046 if(@wr==0) 8047 { 8048 @qh[1]=act[1][@k]^act[2][@k]; 8049 } 8050 else 8051 { 8052 @qh[1]=act[1][@k]; 8053 } 8054 primary[2*@k-1]=interred(@qh); 8055 @qh=j; 8056 @qh[1]=act[1][@k]; 8057 primary[2*@k]=interred(@qh); 8058 attrib( primary[2*@k-1],"isSB",1); 8059 8060 if((size(ser)>0)&&(size(reduce(ser,primary[2*@k-1],1))==0)) 8061 { 8062 primary[2*@k-1]=ideal(1); 8063 primary[2*@k]=ideal(1); 8064 } 8065 } 8066 return(primary); 8067 } 8068 8069 if(homog(j)==1) 8070 { 8071 primary[1]=j; 8072 if((size(ser)>0)&&(size(reduce(ser,j,1))==0)) 8073 { 8074 primary[1]=ideal(1); 8075 primary[2]=ideal(1); 8076 return(primary); 8077 } 8078 if(dim(j)==-1) 8079 { 8080 primary[1]=ideal(1); 8081 primary[2]=ideal(1); 8082 } 8083 else 8084 { 8085 primary[2]=maxideal(1); 8086 } 8087 return(primary); 8088 } 8089 8090 //the first element in the standardbase is factorized 8091 if(deg(j[1])>0) 8092 { 8093 act=factor(j[1]); 8094 testFactor(act,j[1]); 8095 } 8096 else 8097 { 8098 primary[1]=ideal(1); 8099 primary[2]=ideal(1); 8100 return(primary); 8101 } 8102 8103 //with the factors new ideals (hopefully the primary decomposition) 8104 //are created 8105 if(size(act[1])>1) 8106 { 8107 if(size(#)>1) 8108 { 8109 primary[1]=ideal(1); 8110 primary[2]=ideal(1); 8111 primary[3]=ideal(1); 8112 primary[4]=ideal(1); 8113 return(primary); 8114 } 8115 for(@k=1;@k<=size(act[1]);@k++) 8116 { 8117 if(@wr==0) 8118 { 8119 primary[2*@k-1]=std(j,act[1][@k]^act[2][@k]); 8120 8121 } 8122 else 8123 { 8124 primary[2*@k-1]=std(j,act[1][@k]); 8125 } 8126 if((act[2][@k]==1)&&(vdim(primary[2*@k-1])==deg(act[1][@k]))) 8127 { 8128 primary[2*@k] = primary[2*@k-1]; 8129 } 8130 else 8131 { 8132 8133 primary[2*@k] = primaryTest(primary[2*@k-1],act[1][@k]); 8134 8135 } 8136 } 8137 } 8138 else 8139 { 8140 primary[1]=j; 8141 if((size(#)>0)&&(act[2][1]>1)) 8142 { 8143 act[2]=1; 8144 primary[1]=std(primary[1],act[1][1]); 8145 } 8146 if(@wr!=0) 8147 { 8148 primary[1]=std(j,act[1][1]); 8149 } 8150 if((act[2][1]==1)&&(vdim(primary[1])==deg(act[1][1]))) 8151 { 8152 primary[2]=primary[1]; 8153 } 8154 else 8155 { 8156 primary[2]=primaryTest(primary[1],act[1][1]); 8157 } 8158 } 8159 8160 if(size(#)==0) 8161 { 8162 primary=splitPrimary(primary,ser,@wr,act); 8163 } 8164 8165 if((voice>=6)&&(char(basering)<=181)) 8166 { 8167 primary=splitCharp(primary); 8168 } 8169 8170 if((@wr==2)&&(npars(basering)>0)&&(voice>=6)&&(char(basering)>0)) 8171 { 8172 //the prime decomposition of Yokoyama in characteristic p 8173 list ke,ek; 8174 @k=0; 8175 while(@k<size(primary)/2) 8176 { 8177 @k++; 8178 if(size(primary[2*@k])==0) 8179 { 8180 ek=insepDecomp(primary[2*@k-1]); 8181 primary=delete(primary,2*@k); 8182 primary=delete(primary,2*@k-1); 8183 @k--; 8184 } 8185 ke=ke+ek; 8186 } 8187 for(@k=1;@k<=size(ke);@k++) 8188 { 8189 primary[size(primary)+1]=ke[@k]; 8190 primary[size(primary)+1]=ke[@k]; 8191 } 8192 } 8193 8194 if(nestLevel > 1){primary=extF(primary)} 8195 8196 //test whether all ideals in the decomposition are primary and 8197 //in general position 8198 //if not after a random coordinate transformation of the last 8199 //variable the corresponding ideal is decomposed again. 8200 if((npars(basering)>0)&&(nestLevel > 1)) 8201 { 8202 poly randp; 8203 for(zz=1;zz<nvars(basering);zz++) 8204 { 8205 randp=randp 8206 +(random(0,5)*par(1)^2+random(0,5)*par(1)+random(0,5))*var(zz); 8207 } 8208 randp=randp+var(nvars(basering)); 8209 } 8210 @k=0; 8211 while(@k<(size(primary)/2)) 8212 { 8213 @k++; 8214 if (size(primary[2*@k])==0) 8215 { 8216 for(zz=1;zz<size(primary[2*@k-1])-1;zz++) 8217 { 8218 attrib(primary[2*@k-1],"isSB",1); 8219 if(vdim(primary[2*@k-1])==deg(primary[2*@k-1][zz])) 8220 { 8221 primary[2*@k]=primary[2*@k-1]; 8222 } 8223 } 8224 } 8225 } 8226 8227 @k=0; 8228 ideal keep; 8229 while(@k<(size(primary)/2)) 8230 { 8231 @k++; 8232 if (size(primary[2*@k])==0) 8233 { 8234 8235 jmap=randomLast(100); 8236 jmap1=maxideal(1); 8237 jmap2=maxideal(1); 8238 @qht=primary[2*@k-1]; 8239 if((npars(basering)>0)&&(nestLevel > 1)) 8240 { 8241 jmap[size(jmap)]=randp; 8242 } 8243 8244 for(@n=2;@n<=size(primary[2*@k-1]);@n++) 8245 { 8246 if(deg(lead(primary[2*@k-1][@n]))==1) 8247 { 8248 for(zz=1;zz<=nva;zz++) 8249 { 8250 if(lead(primary[2*@k-1][@n])/var(zz)!=0) 8251 { 8252 jmap1[zz]=-1/leadcoef(primary[2*@k-1][@n])*primary[2*@k-1][@n] 8253 +2/leadcoef(primary[2*@k-1][@n])*lead(primary[2*@k-1][@n]); 8254 jmap2[zz]=primary[2*@k-1][@n]; 8255 @qht[@n]=var(zz); 8256 8257 } 8258 } 8259 jmap[nva]=subst(jmap[nva],lead(primary[2*@k-1][@n]),0); 8260 } 8261 } 8262 if(size(subst(jmap[nva],var(1),0)-var(nva))!=0) 8263 { 8264 // jmap[nva]=subst(jmap[nva],var(1),0); 8265 //hier geaendert +untersuchen!!!!!!!!!!!!!! 8266 } 8267 phi1=@P,jmap1; 8268 phi=@P,jmap; 8269 for(@n=1;@n<=nva;@n++) 8270 { 8271 jmap[@n]=-(jmap[@n]-2*var(@n)); 8272 } 8273 psi=@P,jmap; 8274 psi1=@P,jmap2; 8275 @qh=phi(@qht); 8276 8277 //=================== the new part ============================ 8278 8279 @qh=groebner(@qh); 8280 8281 //============================================================= 8282 // if(npars(@P)>0) 8283 // { 8284 // @ri= "ring @Phelp =" 8285 // +string(char(@P))+", 8286 // ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);"; 8287 // } 8288 // else 8289 // { 8290 // @ri= "ring @Phelp =" 8291 // +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);"; 8292 // } 8293 // execute(@ri); 8294 // ideal @qh=homog(imap(@P,@qht),@t); 8295 // 8296 // ideal @qh1=std(@qh); 8297 // @hilb=hilb(@qh1,1); 8298 // @ri= "ring @Phelp1 =" 8299 // +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);"; 8300 // execute(@ri); 8301 // ideal @qh=homog(imap(@P,@qh),@t); 8302 // kill @Phelp; 8303 // @qh=std(@qh,@hilb); 8304 // @qh=subst(@qh,@t,1); 8305 // setring @P; 8306 // @qh=imap(@Phelp1,@qh); 8307 // kill @Phelp1; 8308 // @qh=clearSB(@qh); 8309 // attrib(@qh,"isSB",1); 8310 //============================================================= 8311 8312 ser1=phi1(ser); 8313 @lh=newZero_decomp (@qh,phi(ser1),@wr, list("nest", nestLevel + 1)); 8314 8315 kill lres0; 8316 list lres0; 8317 if(size(@lh)==2) 8318 { 8319 helpprim=@lh[2]; 8320 lres0[1]=primary[2*@k-1]; 8321 ser1=psi(helpprim); 8322 lres0[2]=psi1(ser1); 8323 if(size(reduce(lres0[2],lres0[1],1))==0) 8324 { 8325 primary[2*@k]=primary[2*@k-1]; 8326 continue; 8327 } 8328 } 8329 else 8330 { 8331 8332 lres1=psi(@lh); 8333 lres0=psi1(lres1); 8334 } 8335 8336 //=================== the new part ============================ 8337 8338 primary=delete(primary,2*@k-1); 8339 primary=delete(primary,2*@k-1); 8340 @k--; 8341 if(size(lres0)==2) 8342 { 8343 lres0[2]=groebner(lres0[2]); 8344 } 8345 else 8346 { 8347 for(@n=1;@n<=size(lres0)/2;@n++) 8348 { 8349 if(specialIdealsEqual(lres0[2*@n-1],lres0[2*@n])==1) 8350 { 8351 lres0[2*@n-1]=groebner(lres0[2*@n-1]); 8352 lres0[2*@n]=lres0[2*@n-1]; 8353 attrib(lres0[2*@n],"isSB",1); 8354 } 8355 else 8356 { 8357 lres0[2*@n-1]=groebner(lres0[2*@n-1]); 8358 lres0[2*@n]=groebner(lres0[2*@n]); 8359 } 8360 } 8361 } 8362 primary=primary+lres0; 8363 8364 //============================================================= 8365 // if(npars(@P)>0) 8366 // { 8367 // @ri= "ring @Phelp =" 8368 // +string(char(@P))+", 8369 // ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);"; 8370 // } 8371 // else 8372 // { 8373 // @ri= "ring @Phelp =" 8374 // +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);"; 8375 // } 8376 // execute(@ri); 8377 // list @lvec; 8378 // list @lr=imap(@P,lres0); 8379 // ideal @lr1; 8380 // 8381 // if(size(@lr)==2) 8382 // { 8383 // @lr[2]=homog(@lr[2],@t); 8384 // @lr1=std(@lr[2]); 8385 // @lvec[2]=hilb(@lr1,1); 8386 // } 8387 // else 8388 // { 8389 // for(@n=1;@n<=size(@lr)/2;@n++) 8390 // { 8391 // if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1) 8392 // { 8393 // @lr[2*@n-1]=homog(@lr[2*@n-1],@t); 8394 // @lr1=std(@lr[2*@n-1]); 8395 // @lvec[2*@n-1]=hilb(@lr1,1); 8396 // @lvec[2*@n]=@lvec[2*@n-1]; 8397 // } 8398 // else 8399 // { 8400 // @lr[2*@n-1]=homog(@lr[2*@n-1],@t); 8401 // @lr1=std(@lr[2*@n-1]); 8402 // @lvec[2*@n-1]=hilb(@lr1,1); 8403 // @lr[2*@n]=homog(@lr[2*@n],@t); 8404 // @lr1=std(@lr[2*@n]); 8405 // @lvec[2*@n]=hilb(@lr1,1); 8406 // 8407 // } 8408 // } 8409 // } 8410 // @ri= "ring @Phelp1 =" 8411 // +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);"; 8412 // execute(@ri); 8413 // list @lr=imap(@Phelp,@lr); 8414 // 8415 // kill @Phelp; 8416 // if(size(@lr)==2) 8417 // { 8418 // @lr[2]=std(@lr[2],@lvec[2]); 8419 // @lr[2]=subst(@lr[2],@t,1); 8420 // 8421 // } 8422 // else 8423 // { 8424 // for(@n=1;@n<=size(@lr)/2;@n++) 8425 // { 8426 // if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1) 8427 // { 8428 // @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]); 8429 // @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1); 8430 // @lr[2*@n]=@lr[2*@n-1]; 8431 // attrib(@lr[2*@n],"isSB",1); 8432 // } 8433 // else 8434 // { 8435 // @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]); 8436 // @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1); 8437 // @lr[2*@n]=std(@lr[2*@n],@lvec[2*@n]); 8438 // @lr[2*@n]=subst(@lr[2*@n],@t,1); 8439 // } 8440 // } 8441 // } 8442 // kill @lvec; 8443 // setring @P; 8444 // lres0=imap(@Phelp1,@lr); 8445 // kill @Phelp1; 8446 // for(@n=1;@n<=size(lres0);@n++) 8447 // { 8448 // lres0[@n]=clearSB(lres0[@n]); 8449 // attrib(lres0[@n],"isSB",1); 8450 // } 8451 // 8452 // primary[2*@k-1]=lres0[1]; 8453 // primary[2*@k]=lres0[2]; 8454 // @s=size(primary)/2; 8455 // for(@n=1;@n<=size(lres0)/2-1;@n++) 8456 // { 8457 // primary[2*@s+2*@n-1]=lres0[2*@n+1]; 8458 // primary[2*@s+2*@n]=lres0[2*@n+2]; 8459 // } 8460 // @k--; 8461 //============================================================= 8462 } 8463 } 8464 return(primary); 8465 } 8466 example 8467 { "EXAMPLE:"; echo = 2; 8468 ring r = 0,(x,y,z),lp; 8469 poly p = z2+1; 8470 poly q = z4+2; 8471 ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4; 8472 i=std(i); 8473 list pr= newZero_decomp(i,ideal(0),0); 8474 pr; 8475 } 8476 /////////////////////////////////////////////////////////////////////////////// 8477 6457 8478 //////////////////////////////////////////////////////////////////////////// 6458 8479 /*
Note: See TracChangeset
for help on using the changeset viewer.
