1 | LIB "tst.lib"; |
---|
2 | tst_init(); |
---|
3 | |
---|
4 | LIB"ring.lib"; |
---|
5 | proc PrimaryTest(ideal i, poly p) |
---|
6 | "USAGE: PrimaryTest(i,p); i standard basis with respect to |
---|
7 | lp, p irreducible polynomial in K[var(n)], |
---|
8 | p^a=i[1] for some a; |
---|
9 | RETURN: an ideal, the radical of i if i is primary and in |
---|
10 | general position with respect to lp, |
---|
11 | the zero ideal else. |
---|
12 | " |
---|
13 | { |
---|
14 | int m,e; |
---|
15 | int n=nvars(basering); |
---|
16 | poly t; |
---|
17 | ideal prm=p; |
---|
18 | |
---|
19 | for(m=2;m<=size(i);m++) |
---|
20 | { |
---|
21 | if(size(ideal(leadexp(i[m])))==1) |
---|
22 | { |
---|
23 | n--; |
---|
24 | //----------------i[m] has a power of var(n) as leading term |
---|
25 | attrib(prm,"isSB",1); |
---|
26 | //--- ?? i[m]=(c*var(n)+h)^e modulo prm for h |
---|
27 | // in K[var(n+1),...], c in K ?? |
---|
28 | e=deg(lead(i[m])); |
---|
29 | t=leadcoef(i[m])*e*var(n)+(i[m]-lead(i[m])) |
---|
30 | /var(n)^(e-1); |
---|
31 | i[m]=poly(e)^e*leadcoef(i[m])^(e-1)*i[m]; |
---|
32 | //---if not (0) is returned, else c*var(n)+h is added to prm |
---|
33 | if (reduce(i[m]-t^e,prm,1) !=0) |
---|
34 | { |
---|
35 | return(ideal(0)); |
---|
36 | } |
---|
37 | prm = prm,cleardenom(simplify(t,1)); |
---|
38 | } |
---|
39 | } |
---|
40 | return(prm); |
---|
41 | } |
---|
42 | |
---|
43 | ring s=(0,x),(d,e,f,g),lp; |
---|
44 | ideal i=g^5,(x*f-g)^3,5*e-g^2,x*d^3; |
---|
45 | PrimaryTest(i,g); |
---|
46 | kill s; |
---|
47 | |
---|
48 | proc zeroDecomp(ideal i) |
---|
49 | "USAGE: zeroDecomp(i); i zerodimensional ideal |
---|
50 | RETURN: list l of lists of two ideals such that the |
---|
51 | intersection(l[j][1], j=1..)=i, the l[i][1] are |
---|
52 | primary and the l[i][2] their radicals |
---|
53 | NOTE: algorithm of Gianni/Trager/Zacharias |
---|
54 | " |
---|
55 | { |
---|
56 | def BAS = basering; |
---|
57 | //----the ordering is changed to the lexicographical one |
---|
58 | def R=changeord("lp"); setring R; |
---|
59 | ideal i=fetch(BAS,i); |
---|
60 | int n=nvars(R); |
---|
61 | int k; |
---|
62 | list result,rest; |
---|
63 | ideal primary,prim; |
---|
64 | option(redSB); |
---|
65 | |
---|
66 | //------the random coordinatechange and its inverse |
---|
67 | ideal m=maxideal(1); |
---|
68 | m[n]=0; |
---|
69 | poly p=(random(100,1,n)*transpose(m))[1,1]+var(n); |
---|
70 | m[n]=p; |
---|
71 | map phi=R,m; |
---|
72 | m[n]=2*var(n)-p; |
---|
73 | map invphi=R,m; |
---|
74 | ideal j=groebner(phi(i)); |
---|
75 | //-------------factorization of the first element in i |
---|
76 | list fac=factorize(j[1],2); |
---|
77 | //-------------computation of the primaries and primes |
---|
78 | for(k=1;k<=size(fac[1]);k++) |
---|
79 | { |
---|
80 | p=fac[1][k]^fac[2][k]; |
---|
81 | primary=groebner(j+p); |
---|
82 | prim=PrimaryTest(primary,fac[1][k]); |
---|
83 | //---test whether all ideals were primary and in general |
---|
84 | // position |
---|
85 | if(prim==0) |
---|
86 | { |
---|
87 | rest[size(rest)+1]=i+invphi(p); |
---|
88 | } |
---|
89 | else |
---|
90 | { |
---|
91 | result[size(result)+1]= |
---|
92 | list(std(i+invphi(p)),std(invphi(prim))); |
---|
93 | } |
---|
94 | } |
---|
95 | //-------treat the bad cases collected in the rest again |
---|
96 | for(k=1;k<=size(rest);k++) |
---|
97 | { |
---|
98 | result=result+zeroDecomp(rest[k]); |
---|
99 | } |
---|
100 | option(noredSB); |
---|
101 | setring BAS; |
---|
102 | list result=imap(R,result); |
---|
103 | kill R; |
---|
104 | return(result); |
---|
105 | } |
---|
106 | |
---|
107 | ring r = 32003,(x,y,z),dp; |
---|
108 | poly p = z2+1; |
---|
109 | poly q = z4+2; |
---|
110 | ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4; |
---|
111 | list pr= zeroDecomp(i); |
---|
112 | pr; |
---|
113 | kill r; |
---|
114 | |
---|
115 | proc prepareQuotientring(ideal i) |
---|
116 | "USAGE: prepareQuotientring(i); i standard basis |
---|
117 | RETURN: a list l of two strings: |
---|
118 | l[1] to define K[x\u,u ], u a maximal independent |
---|
119 | set for i |
---|
120 | l[2] to define K(u)[x\u ], u a maximal independent |
---|
121 | set for i |
---|
122 | both rings with lexicographical ordering |
---|
123 | " |
---|
124 | { |
---|
125 | string va,pa; |
---|
126 | //v describes the independent set u: var(j) is in |
---|
127 | //u iff v[j]!=0 |
---|
128 | intvec v=indepSet(i); |
---|
129 | int k; |
---|
130 | |
---|
131 | for(k=1;k<=size(v);k++) |
---|
132 | { |
---|
133 | if(v[k]!=0) |
---|
134 | { |
---|
135 | pa=pa+"var("+string(k)+"),"; |
---|
136 | } |
---|
137 | else |
---|
138 | { |
---|
139 | va=va+"var("+string(k)+"),"; |
---|
140 | } |
---|
141 | } |
---|
142 | |
---|
143 | pa=pa[1..size(pa)-1]; |
---|
144 | va=va[1..size(va)-1]; |
---|
145 | |
---|
146 | string newring |
---|
147 | ="ring nring=("+charstr(basering)+"),("+va+","+pa+"),lp;"; |
---|
148 | string quotring |
---|
149 | ="ring quring=("+charstr(basering)+","+pa+"),("+va+"),lp;"; |
---|
150 | return(newring,quotring); |
---|
151 | } |
---|
152 | |
---|
153 | ring s=(0,x),(a,b,c,d,e,f,g),dp; |
---|
154 | ideal i=x*b*c,d^2,f-g; |
---|
155 | i=std(i); |
---|
156 | def Q=basering; |
---|
157 | list l= prepareQuotientring(i); |
---|
158 | l; |
---|
159 | execute (l[1]); |
---|
160 | basering; |
---|
161 | execute (l[2]); |
---|
162 | basering; |
---|
163 | setring Q; |
---|
164 | |
---|
165 | proc prepareSat(ideal i) |
---|
166 | { |
---|
167 | int k; |
---|
168 | poly p=leadcoef(i[1]); |
---|
169 | for(k=2;k<=size(i);k++) |
---|
170 | { |
---|
171 | p=p*leadcoef(i[k]); |
---|
172 | } |
---|
173 | return(p); |
---|
174 | } |
---|
175 | |
---|
176 | //LIB"elim.lib"; |
---|
177 | proc decomp(ideal i) |
---|
178 | "USAGE: decomp(i); i ideal |
---|
179 | RETURN: list l of lists of two ideals such that the |
---|
180 | intersection(l[j][1], j=1..)=i, the l[i][1] are |
---|
181 | primary and the l[i][2] their radicals |
---|
182 | NOTE: algorithm of Gianni/Trager/Zacharias |
---|
183 | " |
---|
184 | { |
---|
185 | if(i==0) |
---|
186 | { |
---|
187 | return(list(i,i)); |
---|
188 | } |
---|
189 | def BAS = basering; |
---|
190 | ideal j; |
---|
191 | int n=nvars(BAS); |
---|
192 | int k; |
---|
193 | |
---|
194 | ideal SBi=std(i); |
---|
195 | int d=dim(SBi); |
---|
196 | //---the trivial case and the zero-dimensional case |
---|
197 | if ((d==0)||(d==-1)) |
---|
198 | { |
---|
199 | return(zeroDecomp(i)); |
---|
200 | } |
---|
201 | //---prepare the quotientring with respect to a maximal |
---|
202 | // independent set |
---|
203 | list quotring=prepareQuotientring(SBi); |
---|
204 | execute (quotring[1]); |
---|
205 | //---we use this ring to compute a standardbasis of i*quring |
---|
206 | // which is in i |
---|
207 | ideal i=std(imap(BAS,i)); |
---|
208 | //---pass to the quotientring with respect to a maximal |
---|
209 | // independent set |
---|
210 | execute(quotring[2]); |
---|
211 | ideal i=imap(nring,i); |
---|
212 | kill nring; |
---|
213 | //---computation of the zerodimensional decomposition |
---|
214 | list ra=zeroDecomp(i); |
---|
215 | //---preparation for saturation |
---|
216 | list p; |
---|
217 | for(k=1;k<=size(ra);k++) |
---|
218 | { |
---|
219 | p[k]=list(prepareSat(ra[k][1]),prepareSat(ra[k][2])); |
---|
220 | } |
---|
221 | poly q=prepareSat(i); |
---|
222 | //---back to the original ring |
---|
223 | setring BAS; |
---|
224 | list p=imap(quring,p); |
---|
225 | list ra=imap(quring,ra); |
---|
226 | |
---|
227 | poly q=imap(quring,q); |
---|
228 | |
---|
229 | kill quring; |
---|
230 | //---compute the intersection of ra with BAS |
---|
231 | for(k=1;k<=size(ra);k++) |
---|
232 | { |
---|
233 | ra[k]=list(sat(ra[k][1],p[k][1])[1], |
---|
234 | sat(ra[k][2],p[k][2])[1]); |
---|
235 | } |
---|
236 | q=q^sat(i,q)[2]; |
---|
237 | |
---|
238 | //---i=intersection((i:q),(i,q)) and ra is the primary |
---|
239 | // decomposition of i:q |
---|
240 | if(deg(q)>0) |
---|
241 | { |
---|
242 | ra=ra+decomp(i+q); |
---|
243 | } |
---|
244 | return(ra); |
---|
245 | } |
---|
246 | |
---|
247 | ring r = 0,(x,y,z),dp; |
---|
248 | ideal i = intersect(ideal(x,y,z)^3,ideal(x-y-z)^2, |
---|
249 | ideal(x-y,x-z)^2); |
---|
250 | list pr= decomp(i); |
---|
251 | pr; |
---|
252 | kill r; |
---|
253 | |
---|
254 | proc equidimensional(ideal i) |
---|
255 | "USAGE: equidimensional(i); i ideal |
---|
256 | RETURN: list l of two ideals such that intersetion(l[1], |
---|
257 | l[2])=i |
---|
258 | l[1] is equidimensional and dim(l[1])>dim(l[2]) |
---|
259 | " |
---|
260 | { |
---|
261 | def BAS = basering; |
---|
262 | |
---|
263 | ideal SBi=std(i); |
---|
264 | int d=dim(SBi); |
---|
265 | int n=nvars(BAS); |
---|
266 | int k; |
---|
267 | list result; |
---|
268 | |
---|
269 | //-----the trivial cases |
---|
270 | if ((d==-1)||(n==d)||(n==1)||(d==0)) |
---|
271 | { |
---|
272 | result=i,ideal(1); |
---|
273 | return(result); |
---|
274 | } |
---|
275 | //-----prepare the quotientring with respect to a maximal |
---|
276 | // independent set |
---|
277 | list quotring=prepareQuotientring(SBi); |
---|
278 | execute (quotring[1]); |
---|
279 | //-----we use this ring to compute a standardbasis of |
---|
280 | // i*quring which is in i |
---|
281 | ideal eq=std(imap(BAS,i)); |
---|
282 | //-----pass to the quotientring with respect to a maximal |
---|
283 | // independent set |
---|
284 | execute (quotring[2]); |
---|
285 | ideal eq=imap(nring,eq); |
---|
286 | kill nring; |
---|
287 | //-----preparation for saturation |
---|
288 | poly p=prepareSat(eq); |
---|
289 | //-----back to the original ring |
---|
290 | setring BAS; |
---|
291 | poly p=imap(quring,p); |
---|
292 | ideal eq=imap(quring,eq); |
---|
293 | kill quring; |
---|
294 | //-----compute the intersection of eq with BAS |
---|
295 | eq=sat(eq,p)[1]; |
---|
296 | |
---|
297 | SBi=std(quotient(i,eq)); |
---|
298 | |
---|
299 | if(d>dim(SBi)) |
---|
300 | { |
---|
301 | result=eq,SBi; |
---|
302 | return(result); |
---|
303 | } |
---|
304 | result=equidimensional(i); |
---|
305 | result=intersect(result[1],eq),result[2]; |
---|
306 | return(result); |
---|
307 | } |
---|
308 | |
---|
309 | ring r = 0,(x,y,z),dp; |
---|
310 | ideal i = intersect(ideal(x,y,z)^3,ideal(x-y-z)^2, |
---|
311 | ideal(x-y,x-z)^2); |
---|
312 | list pr= equidimensional(i); |
---|
313 | pr; |
---|
314 | dim(std(pr[1])); |
---|
315 | dim(std(pr[2])); |
---|
316 | option(redSB); |
---|
317 | std(i); |
---|
318 | std(intersect(pr[1],pr[2])); |
---|
319 | kill r; |
---|
320 | |
---|
321 | proc squarefree(poly f, int i) |
---|
322 | { |
---|
323 | poly h=gcd(f,diff(f,var(i))); |
---|
324 | poly g=lift(h,f)[1][1]; // f/h |
---|
325 | return(g); |
---|
326 | } |
---|
327 | |
---|
328 | proc radical(ideal i) |
---|
329 | "USAGE: radical(i); i ideal |
---|
330 | RETURN: ideal = the radical of i |
---|
331 | NOTE: algorithm of Krick/Logar |
---|
332 | " |
---|
333 | { |
---|
334 | def BAS = basering; |
---|
335 | ideal j; |
---|
336 | int n=nvars(BAS); |
---|
337 | int k; |
---|
338 | |
---|
339 | option(redSB); |
---|
340 | ideal SBi=std(i); |
---|
341 | option(noredSB); |
---|
342 | int d=dim(SBi); |
---|
343 | |
---|
344 | //-----the trivial cases |
---|
345 | if ((d==-1)||(n==d)||(n==1)) |
---|
346 | { |
---|
347 | return(ideal(squarefree(SBi[1],1))); |
---|
348 | } |
---|
349 | //-----the zero-dimensional case |
---|
350 | if (d==0) |
---|
351 | { |
---|
352 | j=finduni(SBi); |
---|
353 | for(k=1;k<=size(j);k++) |
---|
354 | { |
---|
355 | i=i,squarefree(cleardenom(j[k]),k); |
---|
356 | } |
---|
357 | return(std(i)); |
---|
358 | } |
---|
359 | //-----prepare the quotientring with respect to a maximal |
---|
360 | // independent set |
---|
361 | list quotring=prepareQuotientring(SBi); |
---|
362 | execute (quotring[1]); |
---|
363 | //-----we use this ring to compute a standardbasis of |
---|
364 | // i*quring which is in i |
---|
365 | ideal i=std(imap(BAS,i)); |
---|
366 | //-----pass to the quotientring with respect to a maximal |
---|
367 | // independent set |
---|
368 | execute( quotring[2]); |
---|
369 | ideal i=imap(nring,i); |
---|
370 | kill nring; |
---|
371 | //-----computation of the zerodimensional radical |
---|
372 | ideal ra=radical(i); |
---|
373 | //-----preparation for saturation |
---|
374 | poly p=prepareSat(ra); |
---|
375 | poly q=prepareSat(i); |
---|
376 | //-----back to the original ring |
---|
377 | setring BAS; |
---|
378 | poly p=imap(quring,p); |
---|
379 | poly q=imap(quring,q); |
---|
380 | ideal ra=imap(quring,ra); |
---|
381 | kill quring; |
---|
382 | //-----compute the intersection of ra with BAS |
---|
383 | ra=sat(ra,p)[1]; |
---|
384 | //----now we have radical(i)=intersection(ra,radical((i,q))) |
---|
385 | return(intersect(ra,radical(i+q))); |
---|
386 | } |
---|
387 | |
---|
388 | ring r = 0,(x,y,z),dp; |
---|
389 | ideal i = |
---|
390 | intersect(ideal(x,y,z)^3,ideal(x-y-z)^2,ideal(x-y,x-z)^2); |
---|
391 | ideal j = radical(i); |
---|
392 | j; |
---|
393 | |
---|
394 | tst_status(1);$ |
---|