Context Navigation

primdec.lib @ 80b3cd

spielwiese

Last change on this file since 80b3cd was 80b3cd, checked in by Olaf Bachmann <obachman@…>, 25 years ago
* merged changes of gmg manually git-svn-id: file:///usr/local/Singular/svn/trunk@3390 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 102.6 KB

Line
1	// $Id: primdec.lib,v 1.44 1999-08-03 13:01:02 obachman Exp $
2	///////////////////////////////////////////////////////////////////////////////
3	// primdec.lib //
4	// algorithms for primary decomposition based on //
5	// the ideas of Gianni,Trager,Zacharias //
6	// written by Gerhard Pfister //
7	// //
8	// algorithms for primary decomposition based on //
9	// the ideas of Shimoyama/Yokoyama //
10	// written by Wolfram Decker and Hans Schoenemann //
11	///////////////////////////////////////////////////////////////////////////////
12
13	version="$Id: primdec.lib,v 1.44 1999-08-03 13:01:02 obachman Exp $";
14	info="
15	LIBRARY: primdec.lib: PROCEDURES FOR PRIMARY DECOMPOSITION
16	AUTHORS: Gerhard Pfister, email: pfister@mathematik.uni-kl.de (GTZ)
17	Wolfram Decker, email: decker@math.uni-sb.de (SY)
18	Hans Schoenemann, email: hannes@mathematik.uni-kl.de (SY)
19
20	PROCEDURES:
21	primdecGTZ(I); complete primary decomposition via Gianni,Trager,Zacharias
22	primdecSY(I...); complete primary decomposition via Shimoyama-Yokoyama
23	minAssGTZ(I); the minimal associated primes via Gianni,Trager,Zacharias
24	minAssChar(I...); the minimal associated primes using characteristic sets
25	testPrimary(L,k); tests the result of the primary decomposition
26	radical(I); computes the radical of the ideal I
27	equiRadical(I); the radical of the equidimensional part of the ideal I
28	prepareAss(I); list of radicals of the equidimensional components of I
29
30	REMARK:
31	These procedures are implemented to be used in characteristic 0.
32	They work also in positive characteristic >> 0.
33	In small characteristic primdecGTZ, minAssGTZ, radical and equiRadical may not
34	terminate and primdecSY and minAssChar may not give a complete decomposition.
35	";
36
37	LIB "general.lib";
38	LIB "elim.lib";
39	LIB "poly.lib";
40	LIB "random.lib";
41	LIB "inout.lib";
42	///////////////////////////////////////////////////////////////////////////////
43	//
44	// Gianni/Trager/Zacharias
45	//
46	///////////////////////////////////////////////////////////////////////////////
47
48	proc sat1 (ideal id, poly p)
49	"USAGE: sat1(id,j); id ideal, j polynomial
50	RETURN: saturation of id with respect to j (= union_(k=1...) of id:j^k)
51	NOTE: result is a std basis in the basering
52	EXAMPLE: example sat; shows an example
53	"
54	{
55	int @k;
56	ideal inew=std(id);
57	ideal iold;
58	option(returnSB);
59	while(specialIdealsEqual(iold,inew)==0 )
60	{
61	iold=inew;
62	inew=quotient(iold,p);
63	@k++;
64	}
65	@k--;
66	option(noreturnSB);
67	list L =inew,p^@k;
68	return (L);
69	}
70
71	///////////////////////////////////////////////////////////////////////////////
72
73	proc sat2 (ideal id, ideal h)
74	"USAGE: sat2(id,j); id ideal, j polynomial
75	RETURN: saturation of id with respect to j (= union_(k=1...) of id:j^k)
76	NOTE: result is a std basis in the basering
77	EXAMPLE: example sat2; shows an example
78	"
79	{
80	int @k,@i;
81	def @P= basering;
82	if(ordstr(basering)[1,2]!="dp")
83	{
84	execute "ring @Phelp=("+charstr(@P)+"),("+varstr(@P)+"),(C,dp);";
85	ideal inew=std(imap(@P,id));
86	ideal @h=imap(@P,h);
87	}
88	else
89	{
90	ideal @h=h;
91	ideal inew=std(id);
92	}
93	ideal fac;
94
95	for(@i=1;@i<=ncols(@h);@i++)
96	{
97	if(deg(@h[@i])>0)
98	{
99	fac=fac+factorize(@h[@i],1);
100	}
101	}
102	fac=simplify(fac,4);
103	poly @f=1;
104	if(deg(fac[1])>0)
105	{
106	ideal iold;
107
108	for(@i=1;@i<=size(fac);@i++)
109	{
110	@f=@f*fac[@i];
111	}
112	option(returnSB);
113	while(specialIdealsEqual(iold,inew)==0 )
114	{
115	iold=inew;
116	if(deg(iold[size(iold)])!=1)
117	{
118	inew=quotient(iold,@f);
119	}
120	else
121	{
122	inew=iold;
123	}
124	@k++;
125	}
126	option(noreturnSB);
127	@k--;
128	}
129
130	if(ordstr(@P)[1,2]!="dp")
131	{
132	setring @P;
133	ideal inew=std(imap(@Phelp,inew));
134	poly @f=imap(@Phelp,@f);
135	}
136	list L =inew,@f^@k;
137	return (L);
138	}
139
140	///////////////////////////////////////////////////////////////////////////////
141
142	proc minSat(ideal inew, ideal h)
143	{
144	int i,k;
145	poly f=1;
146	ideal iold,fac;
147	list quotM,l;
148
149	for(i=1;i<=ncols(h);i++)
150	{
151	if(deg(h[i])>0)
152	{
153	fac=fac+factorize(h[i],1);
154	}
155	}
156	fac=simplify(fac,4);
157	if(size(fac)==0)
158	{
159	l=inew,1;
160	return(l);
161	}
162	fac=sort(fac)[1];
163	for(i=1;i<=size(fac);i++)
164	{
165	f=f*fac[i];
166	}
167	quotM[1]=inew;
168	quotM[2]=fac;
169	quotM[3]=f;
170	f=1;
171	option(returnSB);
172	while(specialIdealsEqual(iold,quotM[1])==0)
173	{
174	if(k>0)
175	{
176	f=f*quotM[3];
177	}
178	iold=quotM[1];
179	quotM=quotMin(quotM);
180	k++;
181	}
182	option(noreturnSB);
183	l=quotM[1],f;
184	return(l);
185	}
186
187	proc quotMin(list tsil)
188	{
189	int i,j,k,action;
190	ideal verg;
191	list l;
192	poly g;
193
194	ideal laedi=tsil[1];
195	ideal fac=tsil[2];
196	poly f=tsil[3];
197
198	ideal star=quotient(laedi,f);
199	action=1;
200
201	while(action==1)
202	{
203	if(size(fac)==1)
204	{
205	action=0;
206	break;
207	}
208	for(i=1;i<=size(fac);i++)
209	{
210	g=1;
211	verg=laedi;
212
213	for(j=1;j<=size(fac);j++)
214	{
215	if(i!=j)
216	{
217	g=g*fac[j];
218	}
219	}
220	verg=quotient(laedi,g);
221
222	if(specialIdealsEqual(verg,star)==1)
223	{
224	f=g;
225	fac[i]=0;
226	fac=simplify(fac,2);
227	break;
228	}
229	if(i==size(fac))
230	{
231	action=0;
232	}
233	}
234	}
235	l=star,fac,f;
236	return(l);
237	}
238
239	///////////////////////////////////////////////////////////////////////////////
240
241	proc testFactor(list act,poly p)
242	{
243	poly keep=p;
244
245	int i;
246	poly q=act[1][1]^act[2][1];
247	for(i=2;i<=size(act[1]);i++)
248	{
249	q=q*act[1][i]^act[2][i];
250	}
251	q=1/leadcoef(q)*q;
252	p=1/leadcoef(p)*p;
253	if(p-q!=0)
254	{
255	"ERROR IN FACTOR";
256	basering;
257
258	act;
259	keep;
260	pause();
261
262	p;
263	q;
264	pause();
265	}
266	}
267	///////////////////////////////////////////////////////////////////////////////
268
269	proc factor(poly p)
270	"USAGE: factor(p) p poly
271	RETURN: list=;
272	NOTE:
273	EXAMPLE: example factor; shows an example
274	"
275	{
276
277	ideal @i;
278	list @l;
279	intvec @v,@w;
280	int @j,@k,@n;
281
282	if(deg(p)<=1)
283	{
284	@i=ideal(p);
285	@v=1;
286	@l[1]=@i;
287	@l[2]=@v;
288	return(@l);
289	}
290	if (size(p)==1)
291	{
292	@w=leadexp(p);
293	for(@j=1;@j<=nvars(basering);@j++)
294	{
295	if(@w[@j]!=0)
296	{
297	@k++;
298	@v[@k]=@w[@j];
299	@i=@i+ideal(var(@j));
300	}
301	}
302	@l[1]=@i;
303	@l[2]=@v;
304	return(@l);
305	}
306	// @l=factorize(p,2);
307	@l=factorize(p);
308	// if(npars(basering)>0)
309	// {
310	for(@j=1;@j<=size(@l[1]);@j++)
311	{
312	if(deg(@l[1][@j])==0)
313	{
314	@n++;
315	}
316	}
317	if(@n>0)
318	{
319	if(@n==size(@l[1]))
320	{
321	@l[1]=ideal(1);
322	@v=1;
323	@l[2]=@v;
324	}
325	else
326	{
327	@k=0;
328	int pleh;
329	for(@j=1;@j<=size(@l[1]);@j++)
330	{
331	if(deg(@l[1][@j])!=0)
332	{
333	@k++;
334	@i=@i+ideal(@l[1][@j]);
335	if(size(@i)==pleh)
336	{
337	"//factorization error";
338	@l;
339	@k--;
340	@v[@k]=@v[@k]+@l[2][@j];
341	}
342	else
343	{
344	pleh++;
345	@v[@k]=@l[2][@j];
346	}
347	}
348	}
349	@l[1]=@i;
350	@l[2]=@v;
351	}
352	}
353	// }
354	return(@l);
355	}
356	example
357	{ "EXAMPLE:"; echo = 2;
358	ring r = 0,(x,y,z),lp;
359	poly p = (x+y)^2*(y-z)^3;
360	list l = factor(p);
361	l;
362	ring r1 =(0,b,d,f,g),(a,c,e),lp;
363	poly p =(1d)e^2+(1df^2*g);
364	list l = factor(p);
365	l;
366	ring r2 =(0,b,f,g),(a,c,e,d),lp;
367	poly p =(1d)e^2+(1df^2*g);
368	list l = factor(p);
369	l;
370
371	}
372
373
374
375	///////////////////////////////////////////////////////////////////////////////
376
377	proc idealsEqual( ideal k, ideal j)
378	{
379	return(stdIdealsEqual(std(k),std(j)));
380	}
381
382	proc specialIdealsEqual( ideal k1, ideal k2)
383	{
384	int j;
385
386	if(size(k1)==size(k2))
387	{
388	for(j=1;j<=size(k1);j++)
389	{
390	if(leadexp(k1[j])!=leadexp(k2[j]))
391	{
392	return(0);
393	}
394	}
395	return(1);
396	}
397	return(0);
398	}
399
400	proc stdIdealsEqual( ideal k1, ideal k2)
401	{
402	int j;
403
404	if(size(k1)==size(k2))
405	{
406	for(j=1;j<=size(k1);j++)
407	{
408	if(leadexp(k1[j])!=leadexp(k2[j]))
409	{
410	return(0);
411	}
412	}
413	attrib(k2,"isSB",1);
414	if(size(reduce(k1,k2,1))==0)
415	{
416	return(1);
417	}
418	}
419	return(0);
420	}
421	///////////////////////////////////////////////////////////////////////////////
422
423
424	proc primaryTest (ideal i, poly p)
425	{
426	int m=1;
427	int n=nvars(basering);
428	int e;
429	poly t;
430	ideal h;
431
432	ideal prm=p;
433	attrib(prm,"isSB",1);
434
435	while (n>1)
436	{
437	n=n-1;
438	m=m+1;
439
440	//search for i[m] which has a power of var(n) as leading term
441	if (n==1)
442	{
443	m=size(i);
444	}
445	else
446	{
447	while (lead(i[m])/var(n-1)==0)
448	{
449	m=m+1;
450	}
451	m=m-1;
452	}
453	//check whether i[m] =(c*var(n)+h)^e modulo prm for some
454	//h in K[var(n+1),...,var(nvars(basering))], c in K
455	//if not (0) is returned, else var(n)+h is added to prm
456
457	e=deg(lead(i[m]));
458	t=leadcoef(i[m])evar(n)+(i[m]-lead(i[m]))/var(n)^(e-1);
459
460	i[m]=poly(e)^eleadcoef(i[m])^(e-1)i[m];
461
462	if(i[m]==0)
463	{
464	if ( voice >=15 )
465	{
466	"// WARNING: The characteristic is perhaps too small to use";
467	"// the algorithm of Gianni/Trager/Zacharias.";
468	"// This may result in an infinte loop";
469	"// current nesting level in primaryTest:",voice;"";
470	}
471
472	}
473	if (reduce(i[m]-t^e,prm,1) !=0)
474	{
475	return(ideal(0));
476	}
477	h=interred(t);
478	t=h[1];
479
480	prm = prm,t;
481	attrib(prm,"isSB",1);
482	}
483	return(prm);
484	}
485
486
487	///////////////////////////////////////////////////////////////////////////////
488	proc splitPrimary(list l,ideal ser,int @wr,list sact)
489	{
490	int i,j,k,s,r,w;
491	list keepresult,act,keepprime;
492	poly @f;
493	int sl=size(l);
494
495	for(i=1;i<=sl/2;i++)
496	{
497	if(sact[2][i]>1)
498	{
499	keepprime[i]=l[2*i-1]+ideal(sact[1][i]);
500	}
501	else
502	{
503	keepprime[i]=l[2*i-1];
504	}
505	}
506	i=0;
507	while(i<size(l)/2)
508	{
509	i++;
510	if((size(ser)>0)&&(size(reduce(ser,l[2*i-1],1))==0))
511	{
512	l[2*i-1]=ideal(1);
513	l[2*i]=ideal(1);
514	continue;
515	}
516
517	if(size(l[2*i])==0)
518	{
519	if(homog(l[2*i-1])==1)
520	{
521	l[2*i]=maxideal(1);
522	continue;
523	}
524	j=0;
525	if(i<=sl/2)
526	{
527	j=1;
528	}
529	while(j<size(l[2*i-1]))
530	{
531	j++;
532	act=factor(l[2*i-1][j]);
533	r=size(act[1]);
534	attrib(l[2*i-1],"isSB",1);
535	if((r==1)&&(vdim(l[2i-1])==deg(l[2i-1][j])))
536	{
537	l[2i]=std(l[2i-1],act[1][1]);
538	break;
539	}
540	if((r==1)&&(act[2][1]>1))
541	{
542	keepprime[i]=interred(keepprime[i]+ideal(act[1][1]));
543	if(homog(keepprime[i])==1)
544	{
545	l[2*i]=maxideal(1);
546	break;
547	}
548	}
549	if(gcdTest(act[1])==1)
550	{
551	for(k=2;k<=r;k++)
552	{
553	keepprime[size(l)/2+k-1]=interred(keepprime[i]+ideal(act[1][k]));
554	}
555	keepprime[i]=interred(keepprime[i]+ideal(act[1][1]));
556	for(k=1;k<=r;k++)
557	{
558	if(@wr==0)
559	{
560	keepresult[k]=std(l[2*i-1],act[1][k]^act[2][k]);
561	}
562	else
563	{
564	keepresult[k]=std(l[2*i-1],act[1][k]);
565	}
566	}
567	l[2*i-1]=keepresult[1];
568	if(vdim(keepresult[1])==deg(act[1][1]))
569	{
570	l[2*i]=keepresult[1];
571	}
572	if((homog(keepresult[1])==1)\|\|(homog(keepprime[i])==1))
573	{
574	l[2*i]=maxideal(1);
575	}
576	s=size(l)-2;
577	for(k=2;k<=r;k++)
578	{
579	l[s+2*k-1]=keepresult[k];
580	keepprime[s/2+k]=interred(keepresult[k]+ideal(act[1][k]));
581	if(vdim(keepresult[k])==deg(act[1][k]))
582	{
583	l[s+2*k]=keepresult[k];
584	}
585	else
586	{
587	l[s+2*k]=ideal(0);
588	}
589	if((homog(keepresult[k])==1)\|\|(homog(keepprime[s/2+k])==1))
590	{
591	l[s+2*k]=maxideal(1);
592	}
593	}
594	i--;
595	break;
596	}
597	if(r>=2)
598	{
599	s=size(l);
600	@f=act[1][1];
601	act=sat1(l[2*i-1],act[1][1]);
602	if(deg(act[1][1])>0)
603	{
604	l[s+1]=std(l[2*i-1],act[2]);
605	if(homog(l[s+1])==1)
606	{
607	l[s+2]=maxideal(1);
608	}
609	else
610	{
611	l[s+2]=ideal(0);
612	}
613	keepprime[s/2+1]=interred(keepprime[i]+ideal(@f));
614	if(homog(keepprime[s/2+1])==1)
615	{
616	l[s+2]=maxideal(1);
617	}
618	keepprime[i]=act[1];
619	l[2*i-1]=act[1];
620	attrib(l[2*i-1],"isSB",1);
621	if(homog(l[2*i-1])==1)
622	{
623	l[2*i]=maxideal(1);
624	}
625
626	i--;
627	break;
628	}
629	}
630	}
631	}
632	}
633	if(sl==size(l))
634	{
635	return(l);
636	}
637	for(i=1;i<=size(l)/2;i++)
638	{
639	attrib(l[2*i-1],"isSB",1);
640
641	if((size(ser)>0)&&(size(reduce(ser,l[2i-1],1))==0)&&(deg(l[2i-1][1])>0))
642	{
643	"Achtung in split";
644
645	l[2*i-1]=ideal(1);
646	l[2*i]=ideal(1);
647	}
648	if((size(l[2i])==0)&&(specialIdealsEqual(keepprime[i],l[2i-1])!=1))
649	{
650	keepprime[i]=std(keepprime[i]);
651	if(homog(keepprime[i])==1)
652	{
653	l[2*i]=maxideal(1);
654	}
655	else
656	{
657	act=zero_decomp(keepprime[i],ideal(0),@wr,1);
658	if(size(act)==2)
659	{
660	l[2*i]=act[2];
661	}
662	}
663	}
664	}
665	return(l);
666	}
667	example
668	{ "EXAMPLE:"; echo=2;
669	ring r = 32003,(x,y,z),lp;
670	ideal i1=x(x+1),yz,(z+1)(z-1);
671	ideal i2=xy,yz,(x-2)*(x+3);
672	list l=i1,ideal(0),i2,ideal(0),i2,ideal(1);
673	list l1=splitPrimary(l,ideal(0),0);
674	l1;
675	}
676
677	///////////////////////////////////////////////////////////////////////////////
678
679	proc zero_decomp (ideal j,ideal ser,int @wr,list #)
680	"USAGE: zero_decomp(j,ser,@wr); j,ser ideals, @wr=0 or 1
681	(@wr=0 for primary decomposition, @wr=1 for computaion of associated
682	primes)
683	RETURN: list = list of primary ideals and their radicals (at even positions
684	in the list) if the input is zero-dimensional and a standardbases
685	with respect to lex-ordering
686	If ser!=(0) and ser is contained in j or if j is not zero-dimen-
687	sional then ideal(1),ideal(1) is returned
688	NOTE: Algorithm of Gianni, Traeger, Zacharias
689	EXAMPLE: example zero_decomp; shows an example
690	"
691	{
692	def @P = basering;
693	int uytrewq;
694	int nva = nvars(basering);
695	int @k,@s,@n,@k1,zz;
696	list primary,lres0,lres1,act,@lh,@wh;
697	map phi,psi,phi1,psi1;
698	ideal jmap,jmap1,jmap2,helpprim,@qh,@qht,ser1;
699	intvec @vh,@hilb;
700	string @ri;
701	poly @f;
702
703	if (dim(j)>0)
704	{
705	primary[1]=ideal(1);
706	primary[2]=ideal(1);
707	return(primary);
708	}
709
710	j=interred(j);
711
712	attrib(j,"isSB",1);
713	if(vdim(j)==deg(j[1]))
714	{
715	act=factor(j[1]);
716	for(@k=1;@k<=size(act[1]);@k++)
717	{
718	@qh=j;
719	if(@wr==0)
720	{
721	@qh[1]=act[1][@k]^act[2][@k];
722	}
723	else
724	{
725	@qh[1]=act[1][@k];
726	}
727	primary[2*@k-1]=interred(@qh);
728	@qh=j;
729	@qh[1]=act[1][@k];
730	primary[2*@k]=interred(@qh);
731	attrib( primary[2*@k-1],"isSB",1);
732
733	if((size(ser)>0)&&(size(reduce(ser,primary[2*@k-1],1))==0))
734	{
735	primary[2*@k-1]=ideal(1);
736	primary[2*@k]=ideal(1);
737	}
738	}
739	return(primary);
740	}
741
742	if(homog(j)==1)
743	{
744	primary[1]=j;
745	if((size(ser)>0)&&(size(reduce(ser,j,1))==0))
746	{
747	primary[1]=ideal(1);
748	primary[2]=ideal(1);
749	return(primary);
750	}
751	if(dim(j)==-1)
752	{
753	primary[1]=ideal(1);
754	primary[2]=ideal(1);
755	}
756	else
757	{
758	primary[2]=maxideal(1);
759	}
760	return(primary);
761	}
762
763	//the first element in the standardbase is factorized
764	if(deg(j[1])>0)
765	{
766	act=factor(j[1]);
767	testFactor(act,j[1]);
768	}
769	else
770	{
771	primary[1]=ideal(1);
772	primary[2]=ideal(1);
773	return(primary);
774	}
775
776	//with the factors new ideals (hopefully the primary decomposition)
777	//are created
778
779	if(size(act[1])>1)
780	{
781	if(size(#)>1)
782	{
783	primary[1]=ideal(1);
784	primary[2]=ideal(1);
785	primary[3]=ideal(1);
786	primary[4]=ideal(1);
787	return(primary);
788	}
789	for(@k=1;@k<=size(act[1]);@k++)
790	{
791	if(@wr==0)
792	{
793	primary[2*@k-1]=std(j,act[1][@k]^act[2][@k]);
794	}
795	else
796	{
797	primary[2*@k-1]=std(j,act[1][@k]);
798	}
799	if((act[2][@k]==1)&&(vdim(primary[2*@k-1])==deg(act[1][@k])))
800	{
801	primary[2@k] = primary[2@k-1];
802	}
803	else
804	{
805	primary[2@k] = primaryTest(primary[2@k-1],act[1][@k]);
806	}
807	}
808	}
809	else
810	{
811	primary[1]=j;
812	if((size(#)>0)&&(act[2][1]>1))
813	{
814	act[2]=1;
815	primary[1]=std(primary[1],act[1][1]);
816	}
817
818	if((act[2][1]==1)&&(vdim(primary[1])==deg(act[1][1])))
819	{
820	primary[2]=primary[1];
821	}
822	else
823	{
824	primary[2]=primaryTest(primary[1],act[1][1]);
825	}
826	}
827	if(size(#)==0)
828	{
829
830	primary=splitPrimary(primary,ser,@wr,act);
831
832	}
833
834	//test whether all ideals in the decomposition are primary and
835	//in general position
836	//if not after a random coordinate transformation of the last
837	//variable the corresponding ideal is decomposed again.
838
839	@k=0;
840	while(@k<(size(primary)/2))
841	{
842	@k++;
843	if (size(primary[2*@k])==0)
844	{
845	for(zz=1;zz<size(primary[2*@k-1])-1;zz++)
846	{
847	if(vdim(primary[2@k-1])==deg(primary[2@k-1][zz]))
848	{
849	primary[2@k]=primary[2@k-1];
850	}
851	}
852	}
853	}
854
855	@k=0;
856	ideal keep;
857	while(@k<(size(primary)/2))
858	{
859	@k++;
860	if (size(primary[2*@k])==0)
861	{
862
863	jmap=randomLast(100);
864	jmap1=maxideal(1);
865	jmap2=maxideal(1);
866	@qht=primary[2*@k-1];
867
868	for(@n=2;@n<=size(primary[2*@k-1]);@n++)
869	{
870	if(deg(lead(primary[2*@k-1][@n]))==1)
871	{
872	for(zz=1;zz<=nva;zz++)
873	{
874	if(lead(primary[2*@k-1][@n])/var(zz)!=0)
875	{
876	jmap1[zz]=-1/leadcoef(primary[2@k-1][@n])primary[2*@k-1][@n]
877	+2/leadcoef(primary[2@k-1][@n])lead(primary[2*@k-1][@n]);
878	jmap2[zz]=primary[2*@k-1][@n];
879	@qht[@n]=var(zz);
880
881	}
882	}
883	jmap[nva]=subst(jmap[nva],lead(primary[2*@k-1][@n]),0);
884	}
885	}
886	if(size(subst(jmap[nva],var(1),0)-var(nva))!=0)
887	{
888	jmap[nva]=subst(jmap[nva],var(1),0);
889	}
890	phi1=@P,jmap1;
891	phi=@P,jmap;
892
893	for(@n=1;@n<=nva;@n++)
894	{
895	jmap[@n]=-(jmap[@n]-2*var(@n));
896	}
897	psi=@P,jmap;
898	psi1=@P,jmap2;
899
900	@qh=phi(@qht);
901
902	if(npars(@P)>0)
903	{
904	@ri= "ring @Phelp ="
905	+string(char(@P))+",
906	("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);";
907	}
908	else
909	{
910	@ri= "ring @Phelp ="
911	+string(char(@P))+",("+varstr(@P)+",@t),(C,dp);";
912	}
913	execute(@ri);
914	ideal @qh=homog(imap(@P,@qht),@t);
915
916	ideal @qh1=std(@qh);
917	@hilb=hilb(@qh1,1);
918	@ri= "ring @Phelp1 ="
919	+string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);";
920	execute(@ri);
921	ideal @qh=homog(imap(@P,@qh),@t);
922	kill @Phelp;
923	@qh=std(@qh,@hilb);
924	@qh=subst(@qh,@t,1);
925	setring @P;
926	@qh=imap(@Phelp1,@qh);
927	kill @Phelp1;
928	@qh=clearSB(@qh);
929	attrib(@qh,"isSB",1);
930	ser1=phi1(ser);
931
932
933	@lh=zero_decomp (@qh,phi(ser1),@wr);
934	// @lh=zero_decomp (@qh,psi(ser),@wr);
935
936
937	kill lres0;
938	list lres0;
939	if(size(@lh)==2)
940	{
941	helpprim=@lh[2];
942	lres0[1]=primary[2*@k-1];
943	ser1=psi(helpprim);
944	lres0[2]=psi1(ser1);
945	if(size(reduce(lres0[2],lres0[1],1))==0)
946	{
947	primary[2@k]=primary[2@k-1];
948	continue;
949	}
950	}
951	else
952	{
953	//act=factor(@qh[1]);
954	//if(2*size(act[1])==size(@lh))
955	//{
956
957	// for(@n=1;@n<=size(act[1]);@n++)
958	// {
959	// @f=act[1][@n]^act[2][@n];
960	// ser1=psi(@f);
961	// lres0[2@n-1]=interred(primary[2@k-1]+psi1(ser1));
962	// helpprim=@lh[2*@n];
963	// ser1=psi(helpprim);
964	// lres0[2*@n]=psi1(ser1);
965	// }
966	// }
967	// else
968	// {
969	lres1=psi(@lh);
970	lres0=psi1(lres1);
971	//}
972	}
973	if(npars(@P)>0)
974	{
975	@ri= "ring @Phelp ="
976	+string(char(@P))+",
977	("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);";
978	}
979	else
980	{
981	@ri= "ring @Phelp ="
982	+string(char(@P))+",("+varstr(@P)+",@t),(C,dp);";
983	}
984	execute(@ri);
985	list @lvec;
986	list @lr=imap(@P,lres0);
987	ideal @lr1;
988
989	if(size(@lr)==2)
990	{
991	@lr[2]=homog(@lr[2],@t);
992	@lr1=std(@lr[2]);
993	@lvec[2]=hilb(@lr1,1);
994	}
995	else
996	{
997	for(@n=1;@n<=size(@lr)/2;@n++)
998	{
999	if(specialIdealsEqual(@lr[2@n-1],@lr[2@n])==1)
1000	{
1001	@lr[2@n-1]=homog(@lr[2@n-1],@t);
1002	@lr1=std(@lr[2*@n-1]);
1003	@lvec[2*@n-1]=hilb(@lr1,1);
1004	@lvec[2@n]=@lvec[2@n-1];
1005	}
1006	else
1007	{
1008	@lr[2@n-1]=homog(@lr[2@n-1],@t);
1009	@lr1=std(@lr[2*@n-1]);
1010	@lvec[2*@n-1]=hilb(@lr1,1);
1011	@lr[2@n]=homog(@lr[2@n],@t);
1012	@lr1=std(@lr[2*@n]);
1013	@lvec[2*@n]=hilb(@lr1,1);
1014
1015	}
1016	}
1017	}
1018	@ri= "ring @Phelp1 ="
1019	+string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);";
1020	execute(@ri);
1021	list @lr=imap(@Phelp,@lr);
1022
1023	kill @Phelp;
1024	if(size(@lr)==2)
1025	{
1026	@lr[2]=std(@lr[2],@lvec[2]);
1027	@lr[2]=subst(@lr[2],@t,1);
1028
1029	}
1030	else
1031	{
1032	for(@n=1;@n<=size(@lr)/2;@n++)
1033	{
1034	if(specialIdealsEqual(@lr[2@n-1],@lr[2@n])==1)
1035	{
1036	@lr[2@n-1]=std(@lr[2@n-1],@lvec[2*@n-1]);
1037	@lr[2@n-1]=subst(@lr[2@n-1],@t,1);
1038	@lr[2@n]=@lr[2@n-1];
1039	attrib(@lr[2*@n],"isSB",1);
1040	}
1041	else
1042	{
1043	@lr[2@n-1]=std(@lr[2@n-1],@lvec[2*@n-1]);
1044	@lr[2@n-1]=subst(@lr[2@n-1],@t,1);
1045	@lr[2@n]=std(@lr[2@n],@lvec[2*@n]);
1046	@lr[2@n]=subst(@lr[2@n],@t,1);
1047	}
1048	}
1049	}
1050	kill @lvec;
1051	setring @P;
1052	lres0=imap(@Phelp1,@lr);
1053	kill @Phelp1;
1054	for(@n=1;@n<=size(lres0);@n++)
1055	{
1056	lres0[@n]=clearSB(lres0[@n]);
1057	attrib(lres0[@n],"isSB",1);
1058	}
1059
1060	primary[2*@k-1]=lres0[1];
1061	primary[2*@k]=lres0[2];
1062	@s=size(primary)/2;
1063	for(@n=1;@n<=size(lres0)/2-1;@n++)
1064	{
1065	primary[2@s+2@n-1]=lres0[2*@n+1];
1066	primary[2@s+2@n]=lres0[2*@n+2];
1067	}
1068	@k--;
1069	}
1070	}
1071	return(primary);
1072	}
1073	example
1074	{ "EXAMPLE:"; echo = 2;
1075	ring r = 0,(x,y,z),lp;
1076	poly p = z2+1;
1077	poly q = z4+2;
1078	ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
1079	i=std(i);
1080	list pr= zero_decomp(i,ideal(0),0);
1081	pr;
1082	}
1083
1084	///////////////////////////////////////////////////////////////////////////////
1085
1086	proc ggt (ideal i)
1087	"USAGE: ggt(i); i list of polynomials
1088	RETURN: poly = ggt(i[1],...,i[size(i)])
1089	NOTE:
1090	EXAMPLE: example ggt; shows an example
1091	"
1092	{
1093	int k;
1094	poly p=i[1];
1095	if(deg(p)==0)
1096	{
1097	return(1);
1098	}
1099
1100
1101	for (k=2;k<=size(i);k++)
1102	{
1103	if(deg(i[k])==0)
1104	{
1105	return(1)
1106	}
1107	p=GCD(p,i[k]);
1108	if(deg(p)==0)
1109	{
1110	return(1);
1111	}
1112	}
1113	return(p);
1114	}
1115	example
1116	{ "EXAMPLE:"; echo = 2;
1117	ring r = 0,(x,y,z),lp;
1118	poly p = (x+y)*(y+z);
1119	poly q = (z4+2)*(y+z);
1120	ideal l=p,q;
1121	poly pr= ggt(l);
1122	pr;
1123	}
1124	///////////////////////////////////////////////////////////////////////////////
1125	proc gcdTest(ideal act)
1126	{
1127	int i,j;
1128	if(size(act)<=1)
1129	{
1130	return(0);
1131	}
1132	for (i=1;i<=size(act)-1;i++)
1133	{
1134	for(j=i+1;j<=size(act);j++)
1135	{
1136	if(deg(std(ideal(act[i],act[j]))[1])>0)
1137	{
1138	return(0);
1139	}
1140	}
1141	}
1142	return(1);
1143	}
1144
1145	///////////////////////////////////////////////////////////////////////////////
1146	proc coeffLcm(ideal h)
1147	{
1148	string @pa=parstr(basering);
1149	if(size(@pa)==0)
1150	{
1151	return(lcmP(h));
1152	}
1153	def bsr= basering;
1154	string @id=string(h);
1155	execute "ring @r=0,("+@pa+","+varstr(bsr)+"),(C,dp);";
1156	execute "ideal @i="+@id+";";
1157	poly @p=lcmP(@i);
1158	string @ps=string(@p);
1159	setring bsr;
1160	execute "poly @p="+@ps+";";
1161	return(@p);
1162	}
1163	example
1164	{
1165	"EXAMPLE:"; echo = 2;
1166	ring r =( 0,a,b),(x,y,z),lp;
1167	poly p = (a+b)*(y-z);
1168	poly q = (a+b)*(y+z);
1169	ideal l=p,q;
1170	poly pr= coeffLcm(l);
1171	pr;
1172	}
1173
1174	///////////////////////////////////////////////////////////////////////////////
1175
1176	proc lcmP(ideal i)
1177	"USAGE: lcm(i); i list of polynomials
1178	RETURN: poly = lcm(i[1],...,i[size(i)])
1179	NOTE:
1180	EXAMPLE: example lcm; shows an example
1181	"
1182	{
1183	int k,j;
1184	poly p,q;
1185	i=simplify(i,10);
1186	for(j=1;j<=size(i);j++)
1187	{
1188	if(deg(i[j])>0)
1189	{
1190	p=i[j];
1191	break;
1192	}
1193	}
1194	if(deg(p)==-1)
1195	{
1196	return(1);
1197	}
1198	for (k=j+1;k<=size(i);k++)
1199	{
1200	if(deg(i[k])!=0)
1201	{
1202	q=GCD(p,i[k]);
1203	if(deg(q)==0)
1204	{
1205	p=p*i[k];
1206	}
1207	else
1208	{
1209	p=p/q;
1210	p=p*i[k];
1211	}
1212	}
1213	}
1214	return(p);
1215	}
1216	example
1217	{ "EXAMPLE:"; echo = 2;
1218	ring r = 0,(x,y,z),lp;
1219	poly p = (x+y)*(y+z);
1220	poly q = (z4+2)*(y+z);
1221	ideal l=p,q;
1222	poly pr= lcmP(l);
1223	pr;
1224	l=1,-1,p,1,-1,q,1;
1225	pr=lcmP(l);
1226	pr;
1227	}
1228
1229	///////////////////////////////////////////////////////////////////////////////
1230	proc clearSB (ideal i,list #)
1231	"USAGE: clearSB(i); i ideal which is SB ordered by monomial ordering
1232	RETURN: ideal = minimal SB
1233	NOTE:
1234	EXAMPLE: example clearSB; shows an example
1235	"
1236	{
1237	int k,j;
1238	poly m;
1239	int c=size(i);
1240
1241	if(size(#)==0)
1242	{
1243	for(j=1;j<c;j++)
1244	{
1245	if(deg(i[j])==0)
1246	{
1247	i=ideal(1);
1248	return(i);
1249	}
1250	if(deg(i[j])>0)
1251	{
1252	m=lead(i[j]);
1253	for(k=j+1;k<=c;k++)
1254	{
1255	if(size(lead(i[k])/m)>0)
1256	{
1257	i[k]=0;
1258	}
1259	}
1260	}
1261	}
1262	}
1263	else
1264	{
1265	j=0;
1266	while(j<c-1)
1267	{
1268	j++;
1269	if(deg(i[j])==0)
1270	{
1271	i=ideal(1);
1272	return(i);
1273	}
1274	if(deg(i[j])>0)
1275	{
1276	m=lead(i[j]);
1277	for(k=j+1;k<=c;k++)
1278	{
1279	if(size(lead(i[k])/m)>0)
1280	{
1281	if((leadexp(m)!=leadexp(i[k]))\|\|(#[j]<=#[k]))
1282	{
1283	i[k]=0;
1284	}
1285	else
1286	{
1287	i[j]=0;
1288	break;
1289	}
1290	}
1291	}
1292	}
1293	}
1294	}
1295	return(simplify(i,2));
1296	}
1297	example
1298	{ "EXAMPLE:"; echo = 2;
1299	ring r = (0,a,b),(x,y,z),dp;
1300	ideal i=ax2+y,a2x+y,bx;
1301	list l=1,2,1;
1302	ideal j=clearSB(i,l);
1303	j;
1304	}
1305
1306	///////////////////////////////////////////////////////////////////////////////
1307
1308	proc independSet (ideal j)
1309	"USAGE: independentSet(i); i ideal
1310	RETURN: list = new varstring with the independent set at the end,
1311	ordstring with the corresponding block ordering,
1312	the integer where the independent set starts in the varstring
1313	NOTE:
1314	EXAMPLE: example independentSet; shows an example
1315	"
1316	{
1317	int n,k,di;
1318	list resu,hilf;
1319	string var1,var2;
1320	list v=indepSet(j,1);
1321
1322	for(n=1;n<=size(v);n++)
1323	{
1324	di=0;
1325	var1="";
1326	var2="";
1327	for(k=1;k<=size(v[n]);k++)
1328	{
1329	if(v[n][k]!=0)
1330	{
1331	di++;
1332	var2=var2+"var("+string(k)+"),";
1333	}
1334	else
1335	{
1336	var1=var1+"var("+string(k)+"),";
1337	}
1338	}
1339	if(di>0)
1340	{
1341	var1=var1+var2;
1342	var1=var1[1..size(var1)-1];
1343	hilf[1]=var1;
1344	hilf[2]="lp";
1345	//"lp("+string(nvars(basering)-di)+"),dp("+string(di)+")";
1346	hilf[3]=di;
1347	resu[n]=hilf;
1348	}
1349	else
1350	{
1351	resu[n]=varstr(basering),ordstr(basering),0;
1352	}
1353	}
1354	return(resu);
1355	}
1356	example
1357	{ "EXAMPLE:"; echo = 2;
1358	ring s1=(0,x,y),(a,b,c,d,e,f,g),lp;
1359	ideal i=ea-fbg,fa+be,ec-fdg,fc+de;
1360	i=std(i);
1361	list l=independSet(i);
1362	l;
1363	i=i,g;
1364	l=independSet(i);
1365	l;
1366
1367	ring s=0,(x,y,z),lp;
1368	ideal i=z,yx;
1369	list l=independSet(i);
1370	l;
1371
1372
1373	}
1374	///////////////////////////////////////////////////////////////////////////////
1375
1376	proc maxIndependSet (ideal j)
1377	"USAGE: maxIndependentSet(i); i ideal
1378	RETURN: list = new varstring with the maximal independent set at the end,
1379	ordstring with the corresponding block ordering,
1380	the integer where the independent set starts in the varstring
1381	NOTE:
1382	EXAMPLE: example maxIndependentSet; shows an example
1383	"
1384	{
1385	int n,k,di;
1386	list resu,hilf;
1387	string var1,var2;
1388	list v=indepSet(j,0);
1389
1390	for(n=1;n<=size(v);n++)
1391	{
1392	di=0;
1393	var1="";
1394	var2="";
1395	for(k=1;k<=size(v[n]);k++)
1396	{
1397	if(v[n][k]!=0)
1398	{
1399	di++;
1400	var2=var2+"var("+string(k)+"),";
1401	}
1402	else
1403	{
1404	var1=var1+"var("+string(k)+"),";
1405	}
1406	}
1407	if(di>0)
1408	{
1409	var1=var1+var2;
1410	var1=var1[1..size(var1)-1];
1411	hilf[1]=var1;
1412	hilf[2]="lp";
1413	hilf[3]=di;
1414	resu[n]=hilf;
1415	}
1416	else
1417	{
1418	resu[n]=varstr(basering),ordstr(basering),0;
1419	}
1420	}
1421	return(resu);
1422	}
1423	example
1424	{ "EXAMPLE:"; echo = 2;
1425	ring s1=(0,x,y),(a,b,c,d,e,f,g),lp;
1426	ideal i=ea-fbg,fa+be,ec-fdg,fc+de;
1427	i=std(i);
1428	list l=maxIndependSet(i);
1429	l;
1430	i=i,g;
1431	l=maxIndependSet(i);
1432	l;
1433
1434	ring s=0,(x,y,z),lp;
1435	ideal i=z,yx;
1436	list l=maxIndependSet(i);
1437	l;
1438
1439
1440	}
1441
1442	///////////////////////////////////////////////////////////////////////////////
1443
1444	proc prepareQuotientring (int nnp)
1445	"USAGE: prepareQuotientring(nnp); nnp int
1446	RETURN: string = to define Kvar(nnp+1),...,var(nvars)[..rest ]
1447	NOTE:
1448	EXAMPLE: example independentSet; shows an example
1449	"
1450	{
1451	ideal @ih,@jh;
1452	int npar=npars(basering);
1453	int @n;
1454
1455	string quotring= "ring quring = ("+charstr(basering);
1456	for(@n=nnp+1;@n<=nvars(basering);@n++)
1457	{
1458	quotring=quotring+",var("+string(@n)+")";
1459	@ih=@ih+var(@n);
1460	}
1461
1462	quotring=quotring+"),(var(1)";
1463	@jh=@jh+var(1);
1464	for(@n=2;@n<=nnp;@n++)
1465	{
1466	quotring=quotring+",var("+string(@n)+")";
1467	@jh=@jh+var(@n);
1468	}
1469	quotring=quotring+"),(C,lp);";
1470
1471	return(quotring);
1472
1473	}
1474	example
1475	{ "EXAMPLE:"; echo = 2;
1476	ring s1=(0,x),(a,b,c,d,e,f,g),lp;
1477	def @Q=basering;
1478	list l= prepareQuotientring(3);
1479	l;
1480	execute l[1];
1481	execute l[2];
1482	basering;
1483	phi;
1484	setring @Q;
1485
1486	}
1487
1488	///////////////////////////////////////////////////////////////////////////////
1489
1490	proc projdim(list l)
1491	{
1492	int i=size(l)+1;
1493
1494	while(i>2)
1495	{
1496	i--;
1497	if((size(l[i])>0)&&(deg(l[i][1])>0))
1498	{
1499	return(i);
1500	}
1501	}
1502	}
1503
1504	///////////////////////////////////////////////////////////////////////////////
1505	proc cleanPrimary(list l)
1506	{
1507	int i,j;
1508	list lh;
1509	for(i=1;i<=size(l)/2;i++)
1510	{
1511	if(deg(l[2*i-1][1])>0)
1512	{
1513	j++;
1514	lh[j]=l[2*i-1];
1515	j++;
1516	lh[j]=l[2*i];
1517	}
1518	}
1519	return(lh);
1520	}
1521	///////////////////////////////////////////////////////////////////////////////
1522
1523	proc minAssPrimes(ideal i, list #)
1524	"USAGE: minAssPrimes(i); i ideal
1525	minAssPrimes(i,1); i ideal (to use also the factorizing Groebner)
1526	RETURN: list = the minimal associated prime ideals of i
1527	EXAMPLE: example minAssPrimes; shows an example
1528	"
1529	{
1530	def @P=basering;
1531	list qr=simplifyIdeal(i);
1532	map phi=@P,qr[2];
1533	i=qr[1];
1534
1535	execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),("
1536	+ordstr(basering)+");";
1537
1538
1539	ideal i=fetch(@P,i);
1540	if(size(#)==0)
1541	{
1542	int @wr;
1543	list tluser,@res;
1544	list primary=decomp(i,2);
1545
1546	@res[1]=primary;
1547
1548	tluser=union(@res);
1549	setring @P;
1550	list @res=imap(gnir,tluser);
1551	return(phi(@res));
1552	}
1553	list @res,empty;
1554	ideal ser;
1555	option(redSB);
1556	list @pr=facstd(i);
1557	if(size(@pr)==1)
1558	{
1559	attrib(@pr[1],"isSB",1);
1560	if((dim(@pr[1])==0)&&(homog(@pr[1])==1))
1561	{
1562	setring @P;
1563	list @res=maxideal(1);
1564	return(phi(@res));
1565	}
1566	if(dim(@pr[1])>1)
1567	{
1568	setring @P;
1569	// kill gnir;
1570	execute "ring gnir1 = ("+charstr(basering)+"),
1571	("+varstr(basering)+"),(C,lp);";
1572	ideal i=fetch(@P,i);
1573	list @pr=facstd(i);
1574	// ideal ser;
1575	setring gnir;
1576	@pr=fetch(gnir1,@pr);
1577	kill gnir1;
1578	}
1579	}
1580	option(noredSB);
1581	int j,k,odim,ndim,count;
1582	attrib(@pr[1],"isSB",1);
1583	if(#[1]==77)
1584	{
1585	odim=dim(@pr[1]);
1586	count=1;
1587	intvec pos;
1588	pos[size(@pr)]=0;
1589	for(j=2;j<=size(@pr);j++)
1590	{
1591	attrib(@pr[j],"isSB",1);
1592	ndim=dim(@pr[j]);
1593	if(ndim>odim)
1594	{
1595	for(k=count;k<=j-1;k++)
1596	{
1597	pos[k]=1;
1598	}
1599	count=j;
1600	odim=ndim;
1601	}
1602	if(ndim<odim)
1603	{
1604	pos[j]=1;
1605	}
1606	}
1607	for(j=1;j<=size(@pr);j++)
1608	{
1609	if(pos[j]!=1)
1610	{
1611	@res[j]=decomp(@pr[j],2);
1612	}
1613	else
1614	{
1615	@res[j]=empty;
1616	}
1617	}
1618	}
1619	else
1620	{
1621	ser=ideal(1);
1622	for(j=1;j<=size(@pr);j++)
1623	{
1624	//@pr[j];
1625	//pause();
1626	@res[j]=decomp(@pr[j],2);
1627	// @res[j]=decomp(@pr[j],2,@pr[j],ser);
1628	// for(k=1;k<=size(@res[j]);k++)
1629	// {
1630	// ser=intersect(ser,@res[j][k]);
1631	// }
1632	}
1633	}
1634
1635	@res=union(@res);
1636	setring @P;
1637	list @res=imap(gnir,@res);
1638	return(phi(@res));
1639	}
1640	example
1641	{ "EXAMPLE:"; echo = 2;
1642	ring r = 32003,(x,y,z),lp;
1643	poly p = z2+1;
1644	poly q = z4+2;
1645	ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
1646	list pr= minAssPrimes(i); pr;
1647
1648	minAssPrimes(i,1);
1649	}
1650
1651	proc union(list li)
1652	{
1653	int i,j,k;
1654
1655	def P=basering;
1656
1657	execute "ring ir = ("+charstr(basering)+"),("+varstr(basering)+"),(C,lp);";
1658	list l=fetch(P,li);
1659	list @erg;
1660
1661	for(k=1;k<=size(l);k++)
1662	{
1663	for(j=1;j<=size(l[k])/2;j++)
1664	{
1665	if(deg(l[k][2*j][1])!=0)
1666	{
1667	i++;
1668	@erg[i]=l[k][2*j];
1669	}
1670	}
1671	}
1672
1673	list @wos;
1674	i=0;
1675	ideal i1,i2;
1676	while(i<size(@erg)-1)
1677	{
1678	i++;
1679	k=i+1;
1680	i1=lead(@erg[i]);
1681	attrib(i1,"isSB",1);
1682	attrib(@erg[i],"isSB",1);
1683
1684	while(k<=size(@erg))
1685	{
1686	if(deg(@erg[i][1])==0)
1687	{
1688	break;
1689	}
1690	i2=lead(@erg[k]);
1691	attrib(@erg[k],"isSB",1);
1692	attrib(i2,"isSB",1);
1693
1694	if(size(reduce(i1,i2,1))==0)
1695	{
1696	if(size(reduce(@erg[i],@erg[k],1))==0)
1697	{
1698	@erg[k]=ideal(1);
1699	i2=ideal(1);
1700	}
1701	}
1702	if(size(reduce(i2,i1,1))==0)
1703	{
1704	if(size(reduce(@erg[k],@erg[i],1))==0)
1705	{
1706	break;
1707	}
1708	}
1709	k++;
1710	if(k>size(@erg))
1711	{
1712	@wos[size(@wos)+1]=@erg[i];
1713	}
1714	}
1715	}
1716	if(deg(@erg[size(@erg)][1])!=0)
1717	{
1718	@wos[size(@wos)+1]=@erg[size(@erg)];
1719	}
1720	setring P;
1721	list @ser=fetch(ir,@wos);
1722	return(@ser);
1723	}
1724	///////////////////////////////////////////////////////////////////////////////
1725	proc radicalOld(ideal i)
1726	{
1727	list pr=minAssPrimes(i,1);
1728	int j;
1729	ideal k=pr[1];
1730	for(j=2;j<=size(pr);j++)
1731	{
1732	k=intersect(k,pr[j]);
1733	}
1734	return(k);
1735	}
1736
1737	///////////////////////////////////////////////////////////////////////////////
1738	proc decomp(ideal i,list #)
1739	"USAGE: decomp(i); i ideal (for primary decomposition) (resp.
1740	decomp(i,1); (for the minimal associated primes) )
1741	RETURN: list = list of primary ideals and their associated primes
1742	(at even positions in the list)
1743	(resp. a list of the minimal associated primes)
1744	NOTE: Algorithm of Gianni, Traeger, Zacharias
1745	EXAMPLE: example decomp; shows an example
1746	"
1747	{
1748	def @P = basering;
1749	list primary,indep,ltras;
1750	intvec @vh,isat;
1751	int @wr,@k,@n,@m,@n1,@n2,@n3,homo,seri,keepdi;
1752	ideal peek=i;
1753	ideal ser,tras;
1754
1755	if(size(#)>0)
1756	{
1757	if((#[1]==1)\|\|(#[1]==2))
1758	{
1759	@wr=#[1];
1760	if(size(#)>1)
1761	{
1762	seri=1;
1763	peek=#[2];
1764	ser=#[3];
1765	}
1766	}
1767	else
1768	{
1769	seri=1;
1770	peek=#[1];
1771	ser=#[2];
1772	}
1773	}
1774
1775	homo=homog(i);
1776
1777	if(homo==1)
1778	{
1779	if(attrib(i,"isSB")!=1)
1780	{
1781	ltras=mstd(i);
1782	attrib(ltras[1],"isSB",1);
1783	}
1784	else
1785	{
1786	ltras=i,i;
1787	}
1788	tras=ltras[1];
1789	if(dim(tras)==0)
1790	{
1791	primary[1]=ltras[2];
1792	primary[2]=maxideal(1);
1793	if(@wr>0)
1794	{
1795	list l;
1796	l[1]=maxideal(1);
1797	l[2]=maxideal(1);
1798	return(l);
1799	}
1800	return(primary);
1801	}
1802	intvec @hilb=hilb(tras,1);
1803	intvec keephilb=@hilb;
1804	}
1805
1806	//----------------------------------------------------------------
1807	//i is the zero-ideal
1808	//----------------------------------------------------------------
1809
1810	if(size(i)==0)
1811	{
1812	primary=i,i;
1813	return(primary);
1814	}
1815
1816	//----------------------------------------------------------------
1817	//pass to the lexicographical ordering and compute a standardbasis
1818	//----------------------------------------------------------------
1819
1820	execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),(C,lp);";
1821	option(redSB);
1822
1823	ideal ser=fetch(@P,ser);
1824
1825	if(homo==1)
1826	{
1827	if((ordstr(@P)[1]!="(C,lp)")&&(ordstr(@P)[3]!="(C,lp)"))
1828	{
1829	ideal @j=std(fetch(@P,i),@hilb);
1830	}
1831	else
1832	{
1833	ideal @j=fetch(@P,tras);
1834	attrib(@j,"isSB",1);
1835	}
1836	}
1837	else
1838	{
1839	ideal @j=std(fetch(@P,i));
1840	}
1841	option(noredSB);
1842	if(seri==1)
1843	{
1844	ideal peek=fetch(@P,peek);
1845	attrib(peek,"isSB",1);
1846	}
1847	else
1848	{
1849	ideal peek=@j;
1850	}
1851	if(size(ser)==0)
1852	{
1853	ideal fried;
1854	@n=size(@j);
1855	for(@k=1;@k<=@n;@k++)
1856	{
1857	if(deg(lead(@j[@k]))==1)
1858	{
1859	fried[size(fried)+1]=@j[@k];
1860	@j[@k]=0;
1861	}
1862	}
1863	if(size(fried)>0)
1864	{
1865	@j=simplify(@j,2);
1866	attrib(@j,"isSB",1);
1867	list pr=decomp(@j);
1868	for(@k=1;@k<=size(pr);@k++)
1869	{
1870	@j=pr[@k]+fried;
1871	pr[@k]=@j;
1872	}
1873	setring @P;
1874	return(fetch(gnir,pr));
1875	}
1876	}
1877
1878	//----------------------------------------------------------------
1879	//j is the ring
1880	//----------------------------------------------------------------
1881
1882	if (dim(@j)==-1)
1883	{
1884	setring @P;
1885	return(ideal(0));
1886	}
1887
1888	//----------------------------------------------------------------
1889	// the case of one variable
1890	//----------------------------------------------------------------
1891
1892	if(nvars(basering)==1)
1893	{
1894
1895	list fac=factor(@j[1]);
1896	list gprimary;
1897	for(@k=1;@k<=size(fac[1]);@k++)
1898	{
1899	if(@wr==0)
1900	{
1901	gprimary[2*@k-1]=ideal(fac[1][@k]^fac[2][@k]);
1902	gprimary[2*@k]=ideal(fac[1][@k]);
1903	}
1904	else
1905	{
1906	gprimary[2*@k-1]=ideal(fac[1][@k]);
1907	gprimary[2*@k]=ideal(fac[1][@k]);
1908	}
1909	}
1910	setring @P;
1911	primary=fetch(gnir,gprimary);
1912
1913	return(primary);
1914	}
1915
1916	//------------------------------------------------------------------
1917	//the zero-dimensional case
1918	//------------------------------------------------------------------
1919
1920	if (dim(@j)==0)
1921	{
1922	option(redSB);
1923	list gprimary= zero_decomp(@j,ser,@wr);
1924	option(noredSB);
1925	setring @P;
1926	primary=fetch(gnir,gprimary);
1927	if(size(ser)>0)
1928	{
1929	primary=cleanPrimary(primary);
1930	}
1931	return(primary);
1932	}
1933
1934	poly @gs,@gh,@p;
1935	string @va,quotring;
1936	list quprimary,htprimary,collectprimary,lsau,lnew,allindep,restindep;
1937	ideal @h;
1938	int jdim=dim(@j);
1939	list fett;
1940	int lauf,di,newtest;
1941	//------------------------------------------------------------------
1942	//search for a maximal independent set indep,i.e.
1943	//look for subring such that the intersection with the ideal is zero
1944	//j intersected with K[var(indep[3]+1),...,var(nvar] is zero,
1945	//indep[1] is the new varstring and indep[2] the string for block-ordering
1946	//------------------------------------------------------------------
1947
1948	if(@wr!=1)
1949	{
1950	allindep=independSet(@j);
1951	for(@m=1;@m<=size(allindep);@m++)
1952	{
1953	if(allindep[@m][3]==jdim)
1954	{
1955	di++;
1956	indep[di]=allindep[@m];
1957	}
1958	else
1959	{
1960	lauf++;
1961	restindep[lauf]=allindep[@m];
1962	}
1963	}
1964	}
1965	else
1966	{
1967	indep=maxIndependSet(@j);
1968	}
1969
1970	ideal jkeep=@j;
1971
1972	if(ordstr(@P)[1]=="w")
1973	{
1974	execute "ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),("+ordstr(@P)+");";
1975	}
1976	else
1977	{
1978	execute "ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),(C,dp);";
1979	}
1980
1981	if(homo==1)
1982	{
1983	if((ordstr(@P)[3]=="d")\|\|(ordstr(@P)[1]=="d")\|\|(ordstr(@P)[1]=="w")
1984	\|\|(ordstr(@P)[3]=="w"))
1985	{
1986	ideal jwork=imap(@P,tras);
1987	attrib(jwork,"isSB",1);
1988	}
1989	else
1990	{
1991	ideal jwork=std(imap(gnir,@j),@hilb);
1992	}
1993
1994	}
1995	else
1996	{
1997	ideal jwork=std(imap(gnir,@j));
1998	}
1999	list hquprimary;
2000	poly @p,@q;
2001	ideal @h,fac,ser;
2002	di=dim(jwork);
2003	keepdi=di;
2004
2005	setring gnir;
2006	for(@m=1;@m<=size(indep);@m++)
2007	{
2008	isat=0;
2009	@n2=0;
2010	option(redSB);
2011	if((indep[@m][1]==varstr(basering))&&(@m==1))
2012	//this is the good case, nothing to do, just to have the same notations
2013	//change the ring
2014	{
2015	execute "ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
2016	+ordstr(basering)+");";
2017	ideal @j=fetch(gnir,@j);
2018	attrib(@j,"isSB",1);
2019	ideal ser=fetch(gnir,ser);
2020
2021	}
2022	else
2023	{
2024	@va=string(maxideal(1));
2025	execute "ring gnir1 = ("+charstr(basering)+"),("+indep[@m][1]+"),("
2026	+indep[@m][2]+");";
2027	execute "map phi=gnir,"+@va+";";
2028	if(homo==1)
2029	{
2030	ideal @j=std(phi(@j),@hilb);
2031	}
2032	else
2033	{
2034	ideal @j=std(phi(@j));
2035	}
2036	ideal ser=phi(ser);
2037
2038	}
2039	option(noredSB);
2040	if((deg(@j[1])==0)\|\|(dim(@j)<jdim))
2041	{
2042	setring gnir;
2043	break;
2044	}
2045	for (lauf=1;lauf<=size(@j);lauf++)
2046	{
2047	fett[lauf]=size(@j[lauf]);
2048	}
2049	//------------------------------------------------------------------------
2050	//we have now the following situation:
2051	//j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
2052	//to this quotientring, j is their still a standardbasis, the
2053	//leading coefficients of the polynomials there (polynomials in
2054	//K[var(nnp+1),..,var(nva)]) are collected in the list h,
2055	//we need their ggt, gh, because of the following: let
2056	//(j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..the rest..]
2057	//intersected with K[var(1),...,var(nva)] is (j:gh^n)
2058	//on the other hand j=(j,gh^n) intersected with (j:gh^n)
2059
2060	//------------------------------------------------------------------------
2061
2062	//arrangement for quotientring K(var(nnp+1),..,var(nva))[..the rest..] and
2063	//map phi:K[var(1),...,var(nva)] --->K(var(nnpr+1),..,var(nva))[..rest..]
2064	//------------------------------------------------------------------------
2065
2066	quotring=prepareQuotientring(nvars(basering)-indep[@m][3]);
2067
2068	//---------------------------------------------------------------------
2069	//we pass to the quotientring K(var(nnp+1),..,var(nva))[..the rest..]
2070	//---------------------------------------------------------------------
2071
2072	execute quotring;
2073
2074	// @j considered in the quotientring
2075	ideal @j=imap(gnir1,@j);
2076	ideal ser=imap(gnir1,ser);
2077
2078	kill gnir1;
2079
2080	//j is a standardbasis in the quotientring but usually not minimal
2081	//here it becomes minimal
2082
2083	@j=clearSB(@j,fett);
2084	attrib(@j,"isSB",1);
2085
2086	//we need later ggt(h[1],...)=gh for saturation
2087	ideal @h;
2088	if(deg(@j[1])>0)
2089	{
2090	for(@n=1;@n<=size(@j);@n++)
2091	{
2092	@h[@n]=leadcoef(@j[@n]);
2093	}
2094	//the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..]
2095	option(redSB);
2096	list uprimary= zero_decomp(@j,ser,@wr);
2097	option(noredSB);
2098	}
2099	else
2100	{
2101	list uprimary;
2102	uprimary[1]=ideal(1);
2103	uprimary[2]=ideal(1);
2104	}
2105
2106	//we need the intersection of the ideals in the list quprimary with the
2107	//polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
2108	//but fi polynomials, then the intersection of q with the polynomialring
2109	//is the saturation of the ideal generated by f1,...,fr with respect to
2110	//h which is the lcm of the leading coefficients of the fi considered in
2111	//in the quotientring: this is coded in saturn
2112
2113	list saturn;
2114	ideal hpl;
2115
2116	for(@n=1;@n<=size(uprimary);@n++)
2117	{
2118	uprimary[@n]=interred(uprimary[@n]); // temporary fix
2119	hpl=0;
2120	for(@n1=1;@n1<=size(uprimary[@n]);@n1++)
2121	{
2122	hpl=hpl,leadcoef(uprimary[@n][@n1]);
2123	}
2124	saturn[@n]=hpl;
2125	}
2126
2127	//--------------------------------------------------------------------
2128	//we leave the quotientring K(var(nnp+1),..,var(nva))[..the rest..]
2129	//back to the polynomialring
2130	//---------------------------------------------------------------------
2131	setring gnir;
2132
2133	collectprimary=imap(quring,uprimary);
2134	lsau=imap(quring,saturn);
2135	@h=imap(quring,@h);
2136
2137	kill quring;
2138
2139
2140	@n2=size(quprimary);
2141	@n3=@n2;
2142
2143	for(@n1=1;@n1<=size(collectprimary)/2;@n1++)
2144	{
2145	if(deg(collectprimary[2*@n1][1])>0)
2146	{
2147	@n2++;
2148	quprimary[@n2]=collectprimary[2*@n1-1];
2149	lnew[@n2]=lsau[2*@n1-1];
2150	@n2++;
2151	lnew[@n2]=lsau[2*@n1];
2152	quprimary[@n2]=collectprimary[2*@n1];
2153	}
2154	}
2155
2156	//here the intersection with the polynomialring
2157	//mentioned above is really computed
2158	for(@n=@n3/2+1;@n<=@n2/2;@n++)
2159	{
2160	if(specialIdealsEqual(quprimary[2@n-1],quprimary[2@n]))
2161	{
2162	quprimary[2@n-1]=sat2(quprimary[2@n-1],lnew[2*@n-1])[1];
2163	quprimary[2@n]=quprimary[2@n-1];
2164	}
2165	else
2166	{
2167	if(@wr==0)
2168	{
2169	quprimary[2@n-1]=sat2(quprimary[2@n-1],lnew[2*@n-1])[1];
2170	}
2171	quprimary[2@n]=sat2(quprimary[2@n],lnew[2*@n])[1];
2172	}
2173	}
2174
2175	if(size(@h)>0)
2176	{
2177	//---------------------------------------------------------------
2178	//we change to @Phelp to have the ordering dp for saturation
2179	//---------------------------------------------------------------
2180	setring @Phelp;
2181	@h=imap(gnir,@h);
2182	if(@wr!=1)
2183	{
2184	@q=minSat(jwork,@h)[2];
2185	}
2186	else
2187	{
2188	fac=ideal(0);
2189	for(lauf=1;lauf<=ncols(@h);lauf++)
2190	{
2191	if(deg(@h[lauf])>0)
2192	{
2193	fac=fac+factorize(@h[lauf],1);
2194	}
2195	}
2196	fac=simplify(fac,4);
2197	@q=1;
2198	for(lauf=1;lauf<=size(fac);lauf++)
2199	{
2200	@q=@q*fac[lauf];
2201	}
2202	}
2203	jwork=std(jwork,@q);
2204	keepdi=dim(jwork);
2205	if(keepdi<di)
2206	{
2207	setring gnir;
2208	@j=imap(@Phelp,jwork);
2209	break;
2210	}
2211	if(homo==1)
2212	{
2213	@hilb=hilb(jwork,1);
2214	}
2215
2216	setring gnir;
2217	@j=imap(@Phelp,jwork);
2218	}
2219	}
2220	if((size(quprimary)==0)&&(@wr>0))
2221	{
2222	@j=ideal(1);
2223	quprimary[1]=ideal(1);
2224	quprimary[2]=ideal(1);
2225	}
2226	if((size(quprimary)==0))
2227	{
2228	keepdi=di-1;
2229	}
2230	//---------------------------------------------------------------
2231	//notice that j=sat(j,gh) intersected with (j,gh^n)
2232	//we finished with sat(j,gh) and have to start with (j,gh^n)
2233	//---------------------------------------------------------------
2234	if((deg(@j[1])!=0)&&(@wr!=1))
2235	{
2236	if(size(quprimary)>0)
2237	{
2238	setring @Phelp;
2239	ser=imap(gnir,ser);
2240	hquprimary=imap(gnir,quprimary);
2241	if(@wr==0)
2242	{
2243	ideal htest=hquprimary[1];
2244	for (@n1=2;@n1<=size(hquprimary)/2;@n1++)
2245	{
2246	htest=intersect(htest,hquprimary[2*@n1-1]);
2247	}
2248	}
2249	else
2250	{
2251	ideal htest=hquprimary[2];
2252
2253	for (@n1=2;@n1<=size(hquprimary)/2;@n1++)
2254	{
2255	htest=intersect(htest,hquprimary[2*@n1]);
2256	}
2257	}
2258
2259	if(size(ser)>0)
2260	{
2261	ser=intersect(htest,ser);
2262	}
2263	else
2264	{
2265	ser=htest;
2266	}
2267	setring gnir;
2268	ser=imap(@Phelp,ser);
2269	}
2270	if(size(reduce(ser,peek,1))!=0)
2271	{
2272	for(@m=1;@m<=size(restindep);@m++)
2273	{
2274	// if(restindep[@m][3]>=keepdi)
2275	// {
2276	isat=0;
2277	@n2=0;
2278	option(redSB);
2279
2280	if(restindep[@m][1]==varstr(basering))
2281	//the good case, nothing to do, just to have the same notations
2282	//change the ring
2283	{
2284	execute "ring gnir1 = ("+charstr(basering)+"),("+
2285	varstr(basering)+"),("+ordstr(basering)+");";
2286	ideal @j=fetch(gnir,jkeep);
2287	attrib(@j,"isSB",1);
2288	}
2289	else
2290	{
2291	@va=string(maxideal(1));
2292	execute "ring gnir1 = ("+charstr(basering)+"),("+
2293	restindep[@m][1]+"),(" +restindep[@m][2]+");";
2294	execute "map phi=gnir,"+@va+";";
2295	if(homo==1)
2296	{
2297	ideal @j=std(phi(jkeep),keephilb);
2298	}
2299	else
2300	{
2301	ideal @j=std(phi(jkeep));
2302	}
2303	ideal ser=phi(ser);
2304	}
2305	option(noredSB);
2306
2307	for (lauf=1;lauf<=size(@j);lauf++)
2308	{
2309	fett[lauf]=size(@j[lauf]);
2310	}
2311	//------------------------------------------------------------------
2312	//we have now the following situation:
2313	//j intersected with K[var(nnp+1),..,var(nva)] is zero so we may
2314	//pass to this quotientring, j is their still a standardbasis, the
2315	//leading coefficients of the polynomials there (polynomials in
2316	//K[var(nnp+1),..,var(nva)]) are collected in the list h,
2317	//we need their ggt, gh, because of the following:
2318	//let (j:gh^n)=(j:gh^infinity) then
2319	//j*K(var(nnp+1),..,var(nva))[..the rest..]
2320	//intersected with K[var(1),...,var(nva)] is (j:gh^n)
2321	//on the other hand j=(j,gh^n) intersected with (j:gh^n)
2322
2323	//------------------------------------------------------------------
2324
2325	//the arrangement for the quotientring
2326	// K(var(nnp+1),..,var(nva))[..the rest..]
2327	//and the map phi:K[var(1),...,var(nva)] ---->
2328	//--->K(var(nnpr+1),..,var(nva))[..the rest..]
2329	//------------------------------------------------------------------
2330
2331	quotring=prepareQuotientring(nvars(basering)-restindep[@m][3]);
2332
2333	//------------------------------------------------------------------
2334	//we pass to the quotientring K(var(nnp+1),..,var(nva))[..rest..]
2335	//------------------------------------------------------------------
2336
2337	execute quotring;
2338
2339	// @j considered in the quotientring
2340	ideal @j=imap(gnir1,@j);
2341	ideal ser=imap(gnir1,ser);
2342
2343	kill gnir1;
2344
2345	//j is a standardbasis in the quotientring but usually not minimal
2346	//here it becomes minimal
2347	@j=clearSB(@j,fett);
2348	attrib(@j,"isSB",1);
2349
2350	//we need later ggt(h[1],...)=gh for saturation
2351	ideal @h;
2352
2353	for(@n=1;@n<=size(@j);@n++)
2354	{
2355	@h[@n]=leadcoef(@j[@n]);
2356	}
2357	//the primary decomposition of j*K(var(nnp+1),..,var(nva))[..rest..]
2358
2359	option(redSB);
2360	list uprimary= zero_decomp(@j,ser,@wr);
2361	option(noredSB);
2362
2363
2364	//we need the intersection of the ideals in the list quprimary with
2365	//the polynomialring, i.e. let q=(f1,...,fr) in the quotientring
2366	//such an ideal but fi polynomials, then the intersection of q with
2367	//the polynomialring is the saturation of the ideal generated by
2368	//f1,...,fr with respect toh which is the lcm of the leading
2369	//coefficients of the fi considered in the quotientring:
2370	//this is coded in saturn
2371
2372	list saturn;
2373	ideal hpl;
2374
2375	for(@n=1;@n<=size(uprimary);@n++)
2376	{
2377	hpl=0;
2378	for(@n1=1;@n1<=size(uprimary[@n]);@n1++)
2379	{
2380	hpl=hpl,leadcoef(uprimary[@n][@n1]);
2381	}
2382	saturn[@n]=hpl;
2383	}
2384	//------------------------------------------------------------------
2385	//we leave the quotientring K(var(nnp+1),..,var(nva))[..rest..]
2386	//back to the polynomialring
2387	//------------------------------------------------------------------
2388	setring gnir;
2389
2390	collectprimary=imap(quring,uprimary);
2391	lsau=imap(quring,saturn);
2392	@h=imap(quring,@h);
2393
2394	kill quring;
2395
2396
2397	@n2=size(quprimary);
2398	@n3=@n2;
2399
2400	for(@n1=1;@n1<=size(collectprimary)/2;@n1++)
2401	{
2402	if(deg(collectprimary[2*@n1][1])>0)
2403	{
2404	@n2++;
2405	quprimary[@n2]=collectprimary[2*@n1-1];
2406	lnew[@n2]=lsau[2*@n1-1];
2407	@n2++;
2408	lnew[@n2]=lsau[2*@n1];
2409	quprimary[@n2]=collectprimary[2*@n1];
2410	}
2411	}
2412
2413
2414	//here the intersection with the polynomialring
2415	//mentioned above is really computed
2416
2417	for(@n=@n3/2+1;@n<=@n2/2;@n++)
2418	{
2419	if(specialIdealsEqual(quprimary[2@n-1],quprimary[2@n]))
2420	{
2421	quprimary[2@n-1]=sat2(quprimary[2@n-1],lnew[2*@n-1])[1];
2422	quprimary[2@n]=quprimary[2@n-1];
2423	}
2424	else
2425	{
2426	if(@wr==0)
2427	{
2428	quprimary[2@n-1]=sat2(quprimary[2@n-1],lnew[2*@n-1])[1];
2429	}
2430	quprimary[2@n]=sat2(quprimary[2@n],lnew[2*@n])[1];
2431	}
2432	}
2433	if(@n2>=@n3+2)
2434	{
2435	setring @Phelp;
2436	ser=imap(gnir,ser);
2437	hquprimary=imap(gnir,quprimary);
2438	for(@n=@n3/2+1;@n<=@n2/2;@n++)
2439	{
2440	if(@wr==0)
2441	{
2442	ser=intersect(ser,hquprimary[2*@n-1]);
2443	}
2444	else
2445	{
2446	ser=intersect(ser,hquprimary[2*@n]);
2447	}
2448	}
2449	setring gnir;
2450	ser=imap(@Phelp,ser);
2451	}
2452
2453	// }
2454	}
2455	if(size(reduce(ser,peek,1))!=0)
2456	{
2457	if(@wr>0)
2458	{
2459	htprimary=decomp(@j,@wr,peek,ser);
2460	}
2461	else
2462	{
2463	htprimary=decomp(@j,peek,ser);
2464	}
2465	// here we collect now both results primary(sat(j,gh))
2466	// and primary(j,gh^n)
2467	@n=size(quprimary);
2468	for (@k=1;@k<=size(htprimary);@k++)
2469	{
2470	quprimary[@n+@k]=htprimary[@k];
2471	}
2472	}
2473	}
2474
2475	}
2476	//---------------------------------------------------------------------------
2477	//back to the ring we started with
2478	//the final result: primary
2479	//---------------------------------------------------------------------------
2480
2481	setring @P;
2482	primary=imap(gnir,quprimary);
2483	return(primary);
2484	}
2485
2486
2487	example
2488	{ "EXAMPLE:"; echo = 2;
2489	ring r = 32003,(x,y,z),lp;
2490	poly p = z2+1;
2491	poly q = z4+2;
2492	ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
2493	list pr= decomp(i);
2494	pr;
2495	testPrimary( pr, i);
2496	}
2497
2498	///////////////////////////////////////////////////////////////////////////////
2499	proc radicalKL (list m,ideal ser,list #)
2500	{
2501	ideal i=m[2];
2502	//--------------------------------------------------------------------------
2503	//i is the zero-ideal
2504	//-------------------------------------------------------------------------
2505
2506	if(size(i)==0)
2507	{
2508	return(ideal(0));
2509	}
2510
2511	def @P = basering;
2512	list indep,allindep,restindep,fett,@mu;
2513	intvec @vh,isat;
2514	int @wr,@k,@n,@m,@n1,@n2,@n3,lauf,di;
2515	ideal @j,@j1,fac,@h,collectrad,htrad,lsau;
2516	ideal rad=ideal(1);
2517	ideal te=ser;
2518
2519	poly @p,@q;
2520	string @va,quotring;
2521	int homo=homog(i);
2522
2523	if(size(#)>0)
2524	{
2525	@wr=#[1];
2526	}
2527	@j=m[1];
2528	@j1=m[2];
2529	int jdim=dim(@j);
2530	if(size(reduce(ser,@j,1))==0)
2531	{
2532	return(ser);
2533	}
2534	if(homo==1)
2535	{
2536	if(jdim==0)
2537	{
2538	option(noredSB);
2539	return(maxideal(1));
2540	}
2541	intvec @hilb=hilb(@j,1);
2542	}
2543
2544
2545	//---------------------------------------------------------------------------
2546	//j is the ring
2547	//---------------------------------------------------------------------------
2548
2549	if (jdim==-1)
2550	{
2551	option(noredSB);
2552	return(ideal(0));
2553	}
2554
2555	//---------------------------------------------------------------------------
2556	// the case of one variable
2557	//---------------------------------------------------------------------------
2558
2559	if(nvars(basering)==1)
2560	{
2561	fac=factorize(@j[1],1);
2562	@p=1;
2563	for(@k=1;@k<=size(fac);@k++)
2564	{
2565	@p=@p*fac[@k];
2566	}
2567	option(noredSB);
2568
2569	return(ideal(@p));
2570	}
2571	//---------------------------------------------------------------------------
2572	//the case of a complete intersection
2573	//---------------------------------------------------------------------------
2574	if(jdim+size(@j1)==nvars(basering))
2575	{
2576	// ideal jac=minor(jacob(@j1),size(@j1));
2577	// return(quotient(@j1,jac));
2578	}
2579
2580	//---------------------------------------------------------------------------
2581	//the zero-dimensional case
2582	//---------------------------------------------------------------------------
2583
2584	if (jdim==0)
2585	{
2586	@j1=finduni(@j);
2587	for(@k=1;@k<=size(@j1);@k++)
2588	{
2589	fac=factorize(cleardenom(@j1[@k]),1);
2590	@p=fac[1];
2591	for(@n=2;@n<=size(fac);@n++)
2592	{
2593	@p=@p*fac[@n];
2594	}
2595	@j=@j,@p;
2596	}
2597	@j=std(@j);
2598	option(noredSB);
2599	return(@j);
2600	}
2601
2602	//-------------------------------------------------------------------------
2603	//search for a maximal independent set indep,i.e.
2604	//look for subring such that the intersection with the ideal is zero
2605	//j intersected with K[var(indep[3]+1),...,var(nvar] is zero,
2606	//indep[1] is the new varstring, indep[2] the string for the block-ordering
2607	//-------------------------------------------------------------------------
2608
2609	indep=maxIndependSet(@j);
2610
2611	di=dim(@j);
2612
2613	for(@m=1;@m<=size(indep);@m++)
2614	{
2615	if((indep[@m][1]==varstr(basering))&&(@m==1))
2616	//this is the good case, nothing to do, just to have the same notations
2617	//change the ring
2618	{
2619	execute "ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
2620	+ordstr(basering)+");";
2621	ideal @j=fetch(@P,@j);
2622	attrib(@j,"isSB",1);
2623	}
2624	else
2625	{
2626	@va=string(maxideal(1));
2627	execute "ring gnir1 = ("+charstr(basering)+"),("+indep[@m][1]+"),("
2628	+indep[@m][2]+");";
2629	execute "map phi=@P,"+@va+";";
2630	if(homo==1)
2631	{
2632	ideal @j=std(phi(@j),@hilb);
2633	}
2634	else
2635	{
2636	ideal @j=std(phi(@j));
2637	}
2638	}
2639	if((deg(@j[1])==0)\|\|(dim(@j)<jdim))
2640	{
2641	setring @P;
2642	break;
2643	}
2644	for (lauf=1;lauf<=size(@j);lauf++)
2645	{
2646	fett[lauf]=size(@j[lauf]);
2647	}
2648	//------------------------------------------------------------------------
2649	//we have now the following situation:
2650	//j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
2651	//to this quotientring, j is their still a standardbasis, the
2652	//leading coefficients of the polynomials there (polynomials in
2653	//K[var(nnp+1),..,var(nva)]) are collected in the list h,
2654	//we need their ggt, gh, because of the following:
2655	//let (j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..rest..]
2656	//intersected with K[var(1),...,var(nva)] is (j:gh^n)
2657	//on the other hand j=(j,gh^n) intersected with (j:gh^n)
2658
2659	//------------------------------------------------------------------------
2660
2661	//the arrangement for the quotientring K(var(nnp+1),..,var(nva))[..rest..]
2662	//and the map phi:K[var(1),...,var(nva)] ----->
2663	//K(var(nnpr+1),..,var(nva))[..the rest..]
2664	//------------------------------------------------------------------------
2665
2666	quotring=prepareQuotientring(nvars(basering)-indep[@m][3]);
2667
2668	//------------------------------------------------------------------------
2669	//we pass to the quotientring K(var(nnp+1),..,var(nva))[..the rest..]
2670	//------------------------------------------------------------------------
2671
2672	execute quotring;
2673
2674	// @j considered in the quotientring
2675	ideal @j=imap(gnir1,@j);
2676
2677	kill gnir1;
2678
2679	//j is a standardbasis in the quotientring but usually not minimal
2680	//here it becomes minimal
2681
2682	@j=clearSB(@j,fett);
2683	attrib(@j,"isSB",1);
2684
2685	//we need later ggt(h[1],...)=gh for saturation
2686	ideal @h;
2687	if(deg(@j[1])>0)
2688	{
2689	for(@n=1;@n<=size(@j);@n++)
2690	{
2691	@h[@n]=leadcoef(@j[@n]);
2692	}
2693	//the primary decomposition of j*K(var(nnp+1),..,var(nva))[..rest..]
2694	option(redSB);
2695	@j=interred(@j);
2696	attrib(@j,"isSB",1);
2697	list @mo=@j,@j;
2698	ideal zero_rad= radicalKL(@mo,ideal(1));
2699	}
2700	else
2701	{
2702	ideal zero_rad=ideal(1);
2703	}
2704
2705	//we need the intersection of the ideals in the list quprimary with the
2706	//polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
2707	//but fi polynomials, then the intersection of q with the polynomialring
2708	//is the saturation of the ideal generated by f1,...,fr with respect to
2709	//h which is the lcm of the leading coefficients of the fi considered in
2710	//the quotientring: this is coded in saturn
2711
2712	ideal hpl;
2713
2714	for(@n=1;@n<=size(zero_rad);@n++)
2715	{
2716	hpl=hpl,leadcoef(zero_rad[@n]);
2717	}
2718
2719	//------------------------------------------------------------------------
2720	//we leave the quotientring K(var(nnp+1),..,var(nva))[..the rest..]
2721	//back to the polynomialring
2722	//------------------------------------------------------------------------
2723	setring @P;
2724
2725	collectrad=imap(quring,zero_rad);
2726	lsau=simplify(imap(quring,hpl),2);
2727	@h=imap(quring,@h);
2728
2729	kill quring;
2730
2731
2732	//here the intersection with the polynomialring
2733	//mentioned above is really computed
2734
2735	collectrad=sat2(collectrad,lsau)[1];
2736
2737	if(deg(@h[1])>=0)
2738	{
2739	fac=ideal(0);
2740	for(lauf=1;lauf<=ncols(@h);lauf++)
2741	{
2742	if(deg(@h[lauf])>0)
2743	{
2744	fac=fac+factorize(@h[lauf],1);
2745	}
2746	}
2747	fac=simplify(fac,4);
2748	@q=1;
2749	for(lauf=1;lauf<=size(fac);lauf++)
2750	{
2751	@q=@q*fac[lauf];
2752	}
2753
2754
2755	@mu=mstd(quotient(@j+ideal(@q),rad));
2756	@j=@mu[1];
2757	attrib(@j,"isSB",1);
2758
2759	}
2760	if((deg(rad[1])>0)&&(deg(collectrad[1])>0))
2761	{
2762	rad=intersect(rad,collectrad);
2763	}
2764	else
2765	{
2766	if(deg(collectrad[1])>0)
2767	{
2768	rad=collectrad;
2769	}
2770	}
2771
2772	te=simplify(reduce(te*rad,@j),2);
2773
2774	if((dim(@j)<di)\|\|(size(te)==0))
2775	{
2776	break;
2777	}
2778	if(homo==1)
2779	{
2780	@hilb=hilb(@j,1);
2781	}
2782	}
2783
2784	if(((@wr==1)&&(dim(@j)<di))\|\|(deg(@j[1])==0)\|\|(size(te)==0))
2785	{
2786	return(rad);
2787	}
2788	// rad=intersect(rad,radicalKL(@mu,rad,@wr));
2789	rad=intersect(rad,radicalKL(@mu,ideal(1),@wr));
2790
2791
2792	option(noredSB);
2793	return(rad);
2794	}
2795
2796	///////////////////////////////////////////////////////////////////////////////
2797
2798	proc radicalEHV(ideal i,ideal re,list #)
2799	{
2800	ideal J,I,I0,radI0,L,radI1,I2,radI2;
2801	int l,il;
2802	if(size(#)>0)
2803	{
2804	il=#[1];
2805	}
2806
2807	option(redSB);
2808	list m=mstd(i);
2809	I=m[2];
2810	option(noredSB);
2811	if(size(reduce(re,m[1],1))==0)
2812	{
2813	return(re);
2814	}
2815	int cod=nvars(basering)-dim(m[1]);
2816	if((nvars(basering)<=5)&&(size(m[2])<=5))
2817	{
2818	if(cod==size(m[2]))
2819	{
2820	J=minor(jacob(I),cod);
2821	return(quotient(I,J));
2822	}
2823
2824	for(l=1;l<=cod;l++)
2825	{
2826	I0[l]=I[l];
2827	}
2828	if(dim(std(I0))+cod==nvars(basering))
2829	{
2830	J=minor(jacob(I0),cod);
2831	radI0=quotient(I0,J);
2832	L=quotient(radI0,I);
2833	radI1=quotient(radI0,L);
2834
2835	if(size(reduce(radI1,m[1],1))==0)
2836	{
2837	return(I);
2838	}
2839	if(il==1)
2840	{
2841
2842	return(radI1);
2843	}
2844
2845	I2=sat(I,radI1)[1];
2846
2847	if(deg(I2[1])<=0)
2848	{
2849	return(radI1);
2850	}
2851	return(intersect(radI1,radicalEHV(I2,re,il)));
2852	}
2853	}
2854	return(radicalKL(m,re,il));
2855	}
2856	///////////////////////////////////////////////////////////////////////////////
2857
2858	proc Ann(module M)
2859	{
2860	M=prune(M); //to obtain a small embedding
2861	ideal ann=quotient1(M,freemodule(nrows(M)));
2862	return(ann);
2863	}
2864	///////////////////////////////////////////////////////////////////////////////
2865
2866	//computes the equidimensional part of the ideal i of codimension e
2867	proc int_ass_primary_e(ideal i, int e)
2868	{
2869	if(homog(i)!=1)
2870	{
2871	i=std(i);
2872	}
2873	list re=sres(i,0); //the resolution
2874	re=minres(re); //minimized resolution
2875	ideal ann=AnnExt_R(e,re);
2876	if(nvars(basering)-dim(std(ann))!=e)
2877	{
2878	return(ideal(1));
2879	}
2880	return(ann);
2881	}
2882
2883	///////////////////////////////////////////////////////////////////////////////
2884
2885	//computes the annihilator of Ext^n(R/i,R) with given resolution re
2886	//n is not necessarily the number of variables
2887	proc AnnExt_R(int n,list re)
2888	{
2889	if(n<nvars(basering))
2890	{
2891	matrix f=transpose(re[n+1]); //Hom(_,R)
2892	module k=nres(f,2)[2]; //the kernel
2893	matrix g=transpose(re[n]); //the image of Hom(_,R)
2894
2895	ideal ann=quotient1(g,k); //the anihilator
2896	}
2897	else
2898	{
2899	ideal ann=Ann(transpose(re[n]));
2900	}
2901	return(ann);
2902	}
2903	///////////////////////////////////////////////////////////////////////////////
2904
2905	proc analyze(list pr)
2906	{
2907	int ii,jj;
2908	for(ii=1;ii<=size(pr)/2;ii++)
2909	{
2910	dim(std(pr[2*ii]));
2911	idealsEqual(pr[2ii-1],pr[2ii]);
2912	"===========================";
2913	}
2914
2915	for(ii=size(pr)/2;ii>1;ii--)
2916	{
2917	for(jj=1;jj<ii;jj++)
2918	{
2919	if(size(reduce(pr[2jj],std(pr[2ii],1)))==0)
2920	{
2921	"eingebette Komponente";
2922	jj;
2923	ii;
2924	}
2925	}
2926	}
2927	}
2928
2929	///////////////////////////////////////////////////////////////////////////////
2930	//
2931	// Shimoyama-Yokoyama
2932	//
2933	///////////////////////////////////////////////////////////////////////////////
2934
2935	proc simplifyIdeal(ideal i)
2936	{
2937	def r=basering;
2938
2939	int j,k;
2940	map phi;
2941	poly p;
2942
2943	ideal iwork=i;
2944	ideal imap1=maxideal(1);
2945	ideal imap2=maxideal(1);
2946
2947
2948	for(j=1;j<=nvars(basering);j++)
2949	{
2950	for(k=1;k<=size(i);k++)
2951	{
2952	if(deg(iwork[k]/var(j))==0)
2953	{
2954	p=-1/leadcoef(iwork[k]/var(j))*iwork[k];
2955	imap1[j]=p+2*var(j);
2956	phi=r,imap1;
2957	iwork=phi(iwork);
2958	iwork=subst(iwork,var(j),0);
2959	iwork[k]=var(j);
2960	imap1=maxideal(1);
2961	imap2[j]=-p;
2962	break;
2963	}
2964	}
2965	}
2966	return(iwork,imap2);
2967	}
2968
2969
2970	///////////////////////////////////////////////////////
2971	// ini_mod
2972	// input: a polynomial p
2973	// output: the initial term of p as needed
2974	// in the context of characteristic sets
2975	//////////////////////////////////////////////////////
2976
2977	proc ini_mod(poly p)
2978	{
2979	if (p==0)
2980	{
2981	return(0);
2982	}
2983	int n; matrix m;
2984	for( n=nvars(basering); n>0; n=n-1)
2985	{
2986	m=coef(p,var(n));
2987	if(m[1,1]!=1)
2988	{
2989	p=m[2,1];
2990	break;
2991	}
2992	}
2993	if(deg(p)==0)
2994	{
2995	p=0;
2996	}
2997	return(p);
2998	}
2999	///////////////////////////////////////////////////////
3000	// min_ass_prim_charsets
3001	// input: generators of an ideal PS and an integer cho
3002	// If cho=0, the given ordering of the variables is used.
3003	// Otherwise, the system tries to find an "optimal ordering",
3004	// which in some cases may considerably speed up the algorithm
3005	// output: the minimal associated primes of PS
3006	// algorithm: via characteriostic sets
3007	//////////////////////////////////////////////////////
3008
3009
3010	proc min_ass_prim_charsets (ideal PS, int cho)
3011	{
3012	if((cho<0) and (cho>1))
3013	{
3014	"ERROR: <int> must be 0 or 1"
3015	return();
3016	}
3017	if(system("version")>933)
3018	{
3019	option(notWarnSB);
3020	}
3021	if(cho==0)
3022	{
3023	return(min_ass_prim_charsets0(PS));
3024	}
3025	else
3026	{
3027	return(min_ass_prim_charsets1(PS));
3028	}
3029	}
3030	///////////////////////////////////////////////////////
3031	// min_ass_prim_charsets0
3032	// input: generators of an ideal PS
3033	// output: the minimal associated primes of PS
3034	// algorithm: via characteristic sets
3035	// the given ordering of the variables is used
3036	//////////////////////////////////////////////////////
3037
3038
3039	proc min_ass_prim_charsets0 (ideal PS)
3040	{
3041
3042	matrix m=char_series(PS); // We compute an irreducible
3043	// characteristic series
3044	int i,j,k;
3045	list PSI;
3046	list PHI; // the ideals given by the characteristic series
3047	for(i=nrows(m);i>=1; i--)
3048	{
3049	PHI[i]=ideal(m[i,1..ncols(m)]);
3050	}
3051	// We compute the radical of each ideal in PHI
3052	ideal I,JS,II;
3053	int sizeJS, sizeII;
3054	for(i=size(PHI);i>=1; i--)
3055	{
3056	I=0;
3057	for(j=size(PHI[i]);j>0;j=j-1)
3058	{
3059	I=I+ini_mod(PHI[i][j]);
3060	}
3061	JS=std(PHI[i]);
3062	sizeJS=size(JS);
3063	for(j=size(I);j>0;j=j-1)
3064	{
3065	II=0;
3066	sizeII=0;
3067	k=0;
3068	while(k<=sizeII) // successive saturation
3069	{
3070	option(returnSB);
3071	II=quotient(JS,I[j]);
3072	option(noreturnSB);
3073	//std
3074	// II=std(II);
3075	sizeII=size(II);
3076	if(sizeII==sizeJS)
3077	{
3078	for(k=1;k<=sizeII;k++)
3079	{
3080	if(leadexp(II[k])!=leadexp(JS[k])) break;
3081	}
3082	}
3083	JS=II;
3084	sizeJS=sizeII;
3085	}
3086	}
3087	PSI=insert(PSI,JS);
3088	}
3089	int sizePSI=size(PSI);
3090	// We eliminate redundant ideals
3091	for(i=1;i<sizePSI;i++)
3092	{
3093	for(j=i+1;j<=sizePSI;j++)
3094	{
3095	if(size(PSI[i])!=0)
3096	{
3097	if(size(PSI[j])!=0)
3098	{
3099	if(size(NF(PSI[i],PSI[j],1))==0)
3100	{
3101	PSI[j]=ideal(0);
3102	}
3103	else
3104	{
3105	if(size(NF(PSI[j],PSI[i],1))==0)
3106	{
3107	PSI[i]=ideal(0);
3108	}
3109	}
3110	}
3111	}
3112	}
3113	}
3114	for(i=sizePSI;i>=1;i--)
3115	{
3116	if(size(PSI[i])==0)
3117	{
3118	PSI=delete(PSI,i);
3119	}
3120	}
3121	return (PSI);
3122	}
3123
3124	///////////////////////////////////////////////////////
3125	// min_ass_prim_charsets1
3126	// input: generators of an ideal PS
3127	// output: the minimal associated primes of PS
3128	// algorithm: via characteristic sets
3129	// input: generators of an ideal PS and an integer i
3130	// The system tries to find an "optimal ordering" of
3131	// the variables
3132	//////////////////////////////////////////////////////
3133
3134
3135	proc min_ass_prim_charsets1 (ideal PS)
3136	{
3137	def oldring=basering;
3138	string n=system("neworder",PS);
3139	execute "ring r=("+charstr(oldring)+"),("+n+"),dp;";
3140	ideal PS=imap(oldring,PS);
3141	matrix m=char_series(PS); // We compute an irreducible
3142	// characteristic series
3143	int i,j,k;
3144	ideal I;
3145	list PSI;
3146	list PHI; // the ideals given by the characteristic series
3147	list ITPHI; // their initial terms
3148	for(i=nrows(m);i>=1; i--)
3149	{
3150	PHI[i]=ideal(m[i,1..ncols(m)]);
3151	I=0;
3152	for(j=size(PHI[i]);j>0;j=j-1)
3153	{
3154	I=I,ini_mod(PHI[i][j]);
3155	}
3156	I=I[2..ncols(I)];
3157	ITPHI[i]=I;
3158	}
3159	setring oldring;
3160	matrix m=imap(r,m);
3161	list PHI=imap(r,PHI);
3162	list ITPHI=imap(r,ITPHI);
3163	// We compute the radical of each ideal in PHI
3164	ideal I,JS,II;
3165	int sizeJS, sizeII;
3166	for(i=size(PHI);i>=1; i--)
3167	{
3168	I=0;
3169	for(j=size(PHI[i]);j>0;j=j-1)
3170	{
3171	I=I+ITPHI[i][j];
3172	}
3173	JS=std(PHI[i]);
3174	sizeJS=size(JS);
3175	for(j=size(I);j>0;j=j-1)
3176	{
3177	II=0;
3178	sizeII=0;
3179	k=0;
3180	while(k<=sizeII) // successive iteration
3181	{
3182	option(returnSB);
3183	II=quotient(JS,I[j]);
3184	option(noreturnSB);
3185	//std
3186	// II=std(II);
3187	sizeII=size(II);
3188	if(sizeII==sizeJS)
3189	{
3190	for(k=1;k<=sizeII;k++)
3191	{
3192	if(leadexp(II[k])!=leadexp(JS[k])) break;
3193	}
3194	}
3195	JS=II;
3196	sizeJS=sizeII;
3197	}
3198	}
3199	PSI=insert(PSI,JS);
3200	}
3201	int sizePSI=size(PSI);
3202	// We eliminate redundant ideals
3203	for(i=1;i<sizePSI;i++)
3204	{
3205	for(j=i+1;j<=sizePSI;j++)
3206	{
3207	if(size(PSI[i])!=0)
3208	{
3209	if(size(PSI[j])!=0)
3210	{
3211	if(size(NF(PSI[i],PSI[j],1))==0)
3212	{
3213	PSI[j]=ideal(0);
3214	}
3215	else
3216	{
3217	if(size(NF(PSI[j],PSI[i],1))==0)
3218	{
3219	PSI[i]=ideal(0);
3220	}
3221	}
3222	}
3223	}
3224	}
3225	}
3226	for(i=sizePSI;i>=1;i--)
3227	{
3228	if(size(PSI[i])==0)
3229	{
3230	PSI=delete(PSI,i);
3231	}
3232	}
3233	return (PSI);
3234	}
3235
3236
3237	/////////////////////////////////////////////////////
3238	// proc prim_dec
3239	// input: generators of an ideal I and an integer choose
3240	// If choose=0, min_ass_prim_charsets with the given
3241	// ordering of the variables is used.
3242	// If choose=1, min_ass_prim_charsets with the "optimized"
3243	// ordering of the variables is used.
3244	// If choose=2, minAssPrimes from primdec.lib is used
3245	// If choose=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
3246	// output: a primary decomposition of I, i.e., a list
3247	// of pairs consisting of a standard basis of a primary component
3248	// of I and a standard basis of the corresponding associated prime.
3249	// To compute the minimal associated primes of a given ideal
3250	// min_ass_prim_l is called, i.e., the minimal associated primes
3251	// are computed via characteristic sets.
3252	// In the homogeneous case, the performance of the procedure
3253	// will be improved if I is already given by a minimal set of
3254	// generators. Apply minbase if necessary.
3255	//////////////////////////////////////////////////////////
3256
3257
3258	proc prim_dec(ideal I, int choose)
3259	{
3260	if((choose<0) or (choose>3))
3261	{
3262	"ERROR: <int> must be 0 or 1 or 2 or 3";
3263	return();
3264	}
3265	if(system("version")>933)
3266	{
3267	option(notWarnSB);
3268	}
3269	ideal H=1; // The intersection of the primary components
3270	list U; // the leaves of the decomposition tree, i.e.,
3271	// pairs consisting of a primary component of I
3272	// and the corresponding associated prime
3273	list W; // the non-leaf vertices in the decomposition tree.
3274	// every entry has 6 components:
3275	// 1- the vertex itself , i.e., a standard bais of the
3276	// given ideal I (type 1), or a standard basis of a
3277	// pseudo-primary component arising from
3278	// pseudo-primary decomposition (type 2), or a
3279	// standard basis of a remaining component arising from
3280	// pseudo-primary decomposition or extraction (type 3)
3281	// 2- the type of the vertex as indicated above
3282	// 3- the weighted_tree_depth of the vertex
3283	// 4- the tester of the vertex
3284	// 5- a standard basis of the associated prime
3285	// of a vertex of type 2, or 0 otherwise
3286	// 6- a list of pairs consisting of a standard
3287	// basis of a minimal associated prime ideal
3288	// of the father of the vertex and the
3289	// irreducible factors of the "minimal
3290	// divisor" of the seperator or extractor
3291	// corresponding to the prime ideal
3292	// as computed by the procedure minsat,
3293	// if the vertex is of type 3, or
3294	// the empty list otherwise
3295	ideal SI=std(I);
3296	int ncolsSI=ncols(SI);
3297	int ncolsH=1;
3298	W[1]=list(I,1,0,poly(1),ideal(0),list()); // The root of the tree
3299	int weighted_tree_depth;
3300	int i,j;
3301	int check;
3302	list V; // current vertex
3303	list VV; // new vertex
3304	list QQ;
3305	list WI;
3306	ideal Qi,SQ,SRest,fac;
3307	poly tester;
3308
3309	while(1)
3310	{
3311	i=1;
3312	while(1)
3313	{
3314	while(i<=size(W)) // find vertex V of smallest weighted tree-depth
3315	{
3316	if (W[i][3]<=weighted_tree_depth) break;
3317	i++;
3318	}
3319	if (i<=size(W)) break;
3320	i=1;
3321	weighted_tree_depth++;
3322	}
3323	V=W[i];
3324	W=delete(W,i); // delete V from W
3325
3326	// now proceed by type of vertex V
3327
3328	if (V[2]==2) // extraction needed
3329	{
3330	SQ,SRest,fac=extraction(V[1],V[5]);
3331	// standard basis of primary component,
3332	// standard basis of remaining component,
3333	// irreducible factors of
3334	// the "minimal divisor" of the extractor
3335	// as computed by the procedure minsat,
3336	check=0;
3337	for(j=1;j<=ncolsH;j++)
3338	{
3339	if (NF(H[j],SQ,1)!=0) // Q is not redundant
3340	{
3341	check=1;
3342	break;
3343	}
3344	}
3345	if(check==1) // Q is not redundant
3346	{
3347	QQ=list();
3348	QQ[1]=list(SQ,V[5]); // primary component, associated prime,
3349	// i.e., standard bases thereof
3350	U=U+QQ;
3351	H=intersect(H,SQ);
3352	H=std(H);
3353	ncolsH=ncols(H);
3354	check=0;
3355	if(ncolsH==ncolsSI)
3356	{
3357	for(j=1;j<=ncolsSI;j++)
3358	{
3359	if(leadexp(H[j])!=leadexp(SI[j]))
3360	{
3361	check=1;
3362	break;
3363	}
3364	}
3365	}
3366	else
3367	{
3368	check=1;
3369	}
3370	if(check==0) // H==I => U is a primary decomposition
3371	{
3372	return(U);
3373	}
3374	}
3375	if (SRest[1]!=1) // the remaining component is not
3376	// the whole ring
3377	{
3378	if (rad_con(V[4],SRest)==0) // the new vertex is not the
3379	// root of a redundant subtree
3380	{
3381	VV[1]=SRest; // remaining component
3382	VV[2]=3; // pseudoprimdec_special
3383	VV[3]=V[3]+1; // weighted depth
3384	VV[4]=V[4]; // the tester did not change
3385	VV[5]=ideal(0);
3386	VV[6]=list(list(V[5],fac));
3387	W=insert(W,VV,size(W));
3388	}
3389	}
3390	}
3391	else
3392	{
3393	if (V[2]==3) // pseudo_prim_dec_special is needed
3394	{
3395	QQ,SRest=pseudo_prim_dec_special_charsets(V[1],V[6],choose);
3396	// QQ = quadruples:
3397	// standard basis of pseudo-primary component,
3398	// standard basis of corresponding prime,
3399	// seperator, irreducible factors of
3400	// the "minimal divisor" of the seperator
3401	// as computed by the procedure minsat,
3402	// SRest=standard basis of remaining component
3403	}
3404	else // V is the root, pseudo_prim_dec is needed
3405	{
3406	QQ,SRest=pseudo_prim_dec_charsets(I,SI,choose);
3407	// QQ = quadruples:
3408	// standard basis of pseudo-primary component,
3409	// standard basis of corresponding prime,
3410	// seperator, irreducible factors of
3411	// the "minimal divisor" of the seperator
3412	// as computed by the procedure minsat,
3413	// SRest=standard basis of remaining component
3414
3415	}
3416	//check
3417	for(i=size(QQ);i>=1;i--)
3418	//for(i=1;i<=size(QQ);i++)
3419	{
3420	tester=QQ[i][3]*V[4];
3421	Qi=QQ[i][2];
3422	if(NF(tester,Qi,1)!=0) // the new vertex is not the
3423	// root of a redundant subtree
3424	{
3425	VV[1]=QQ[i][1];
3426	VV[2]=2;
3427	VV[3]=V[3]+1;
3428	VV[4]=tester; // the new tester as computed above
3429	VV[5]=Qi; // QQ[i][2];
3430	VV[6]=list();
3431	W=insert(W,VV,size(W));
3432	}
3433	}
3434	if (SRest[1]!=1) // the remaining component is not
3435	// the whole ring
3436	{
3437	if (rad_con(V[4],SRest)==0) // the vertex is not the root
3438	// of a redundant subtree
3439	{
3440	VV[1]=SRest;
3441	VV[2]=3;
3442	VV[3]=V[3]+2;
3443	VV[4]=V[4]; // the tester did not change
3444	VV[5]=ideal(0);
3445	WI=list();
3446	for(i=1;i<=size(QQ);i++)
3447	{
3448	WI=insert(WI,list(QQ[i][2],QQ[i][4]));
3449	}
3450	VV[6]=WI;
3451	W=insert(W,VV,size(W));
3452	}
3453	}
3454	}
3455	}
3456	}
3457
3458	//////////////////////////////////////////////////////////////////////////
3459	// proc pseudo_prim_dec_charsets
3460	// input: Generators of an arbitrary ideal I, a standard basis SI of I,
3461	// and an integer choo
3462	// If choo=0, min_ass_prim_charsets with the given
3463	// ordering of the variables is used.
3464	// If choo=1, min_ass_prim_charsets with the "optimized"
3465	// ordering of the variables is used.
3466	// If choo=2, minAssPrimes from primdec.lib is used
3467	// If choo=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
3468	// output: a pseudo primary decomposition of I, i.e., a list
3469	// of pseudo primary components together with a standard basis of the
3470	// remaining component. Each pseudo primary component is
3471	// represented by a quadrupel: A standard basis of the component,
3472	// a standard basis of the corresponding associated prime, the
3473	// seperator of the component, and the irreducible factors of the
3474	// "minimal divisor" of the seperator as computed by the procedure minsat,
3475	// calls proc pseudo_prim_dec_i
3476	//////////////////////////////////////////////////////////////////////////
3477
3478
3479	proc pseudo_prim_dec_charsets (ideal I, ideal SI, int choo)
3480	{
3481	list L; // The list of minimal associated primes,
3482	// each one given by a standard basis
3483	if((choo==0) or (choo==1))
3484	{
3485	L=min_ass_prim_charsets(I,choo);
3486	}
3487	else
3488	{
3489	if(choo==2)
3490	{
3491	L=minAssPrimes(I);
3492	}
3493	else
3494	{
3495	L=minAssPrimes(I,1);
3496	}
3497	for(int i=size(L);i>=1;i=i-1)
3498	{
3499	L[i]=std(L[i]);
3500	}
3501	}
3502	return (pseudo_prim_dec_i(SI,L));
3503	}
3504
3505	////////////////////////////////////////////////////////////////
3506	// proc pseudo_prim_dec_special_charsets
3507	// input: a standard basis of an ideal I whose radical is the
3508	// intersection of the radicals of ideals generated by one prime ideal
3509	// P_i together with one polynomial f_i, the list V6 must be the list of
3510	// pairs (standard basis of P_i, irreducible factors of f_i),
3511	// and an integer choo
3512	// If choo=0, min_ass_prim_charsets with the given
3513	// ordering of the variables is used.
3514	// If choo=1, min_ass_prim_charsets with the "optimized"
3515	// ordering of the variables is used.
3516	// If choo=2, minAssPrimes from primdec.lib is used
3517	// If choo=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
3518	// output: a pseudo primary decomposition of I, i.e., a list
3519	// of pseudo primary components together with a standard basis of the
3520	// remaining component. Each pseudo primary component is
3521	// represented by a quadrupel: A standard basis of the component,
3522	// a standard basis of the corresponding associated prime, the
3523	// seperator of the component, and the irreducible factors of the
3524	// "minimal divisor" of the seperator as computed by the procedure minsat,
3525	// calls proc pseudo_prim_dec_i
3526	////////////////////////////////////////////////////////////////
3527
3528
3529	proc pseudo_prim_dec_special_charsets (ideal SI,list V6, int choo)
3530	{
3531	int i,j,l;
3532	list m;
3533	list L;
3534	int sizeL;
3535	ideal P,SP; ideal fac;
3536	int dimSP;
3537	for(l=size(V6);l>=1;l--) // creates a list of associated primes
3538	// of I, possibly redundant
3539	{
3540	P=V6[l][1];
3541	fac=V6[l][2];
3542	for(i=ncols(fac);i>=1;i--)
3543	{
3544	SP=P+fac[i];
3545	SP=std(SP);
3546	if(SP[1]!=1)
3547	{
3548	if((choo==0) or (choo==1))
3549	{
3550	m=min_ass_prim_charsets(SP,choo); // a list of SB
3551	}
3552	else
3553	{
3554	if(choo==2)
3555	{
3556	m=minAssPrimes(SP);
3557	}
3558	else
3559	{
3560	m=minAssPrimes(SP,1);
3561	}
3562	for(j=size(m);j>=1;j=j-1)
3563	{
3564	m[j]=std(m[j]);
3565	}
3566	}
3567	dimSP=dim(SP);
3568	for(j=size(m);j>=1; j--)
3569	{
3570	if(dim(m[j])==dimSP)
3571	{
3572	L=insert(L,m[j],size(L));
3573	}
3574	}
3575	}
3576	}
3577	}
3578	sizeL=size(L);
3579	for(i=1;i<sizeL;i++) // get rid of redundant primes
3580	{
3581	for(j=i+1;j<=sizeL;j++)
3582	{
3583	if(size(L[i])!=0)
3584	{
3585	if(size(L[j])!=0)
3586	{
3587	if(size(NF(L[i],L[j],1))==0)
3588	{
3589	L[j]=ideal(0);
3590	}
3591	else
3592	{
3593	if(size(NF(L[j],L[i],1))==0)
3594	{
3595	L[i]=ideal(0);
3596	}
3597	}
3598	}
3599	}
3600	}
3601	}
3602	for(i=sizeL;i>=1;i--)
3603	{
3604	if(size(L[i])==0)
3605	{
3606	L=delete(L,i);
3607	}
3608	}
3609	return (pseudo_prim_dec_i(SI,L));
3610	}
3611
3612
3613	////////////////////////////////////////////////////////////////
3614	// proc pseudo_prim_dec_i
3615	// input: A standard basis of an arbitrary ideal I, and standard bases
3616	// of the minimal associated primes of I
3617	// output: a pseudo primary decomposition of I, i.e., a list
3618	// of pseudo primary components together with a standard basis of the
3619	// remaining component. Each pseudo primary component is
3620	// represented by a quadrupel: A standard basis of the component Q_i,
3621	// a standard basis of the corresponding associated prime P_i, the
3622	// seperator of the component, and the irreducible factors of the
3623	// "minimal divisor" of the seperator as computed by the procedure minsat,
3624	////////////////////////////////////////////////////////////////
3625
3626
3627	proc pseudo_prim_dec_i (ideal SI, list L)
3628	{
3629	list Q;
3630	if (size(L)==1) // one minimal associated prime only
3631	// the ideal is already pseudo primary
3632	{
3633	Q=SI,L[1],1;
3634	list QQ;
3635	QQ[1]=Q;
3636	return (QQ,ideal(1));
3637	}
3638
3639	poly f0,f,g;
3640	ideal fac;
3641	int i,j,k,l;
3642	ideal SQi;
3643	ideal I'=SI;
3644	list QP;
3645	int sizeL=size(L);
3646	for(i=1;i<=sizeL;i++)
3647	{
3648	fac=0;
3649	for(j=1;j<=sizeL;j++) // compute the seperator sep_i
3650	// of the i-th component
3651	{
3652	if (i!=j) // search g not in L[i], but L[j]
3653	{
3654	for(k=1;k<=ncols(L[j]);k++)
3655	{
3656	if(NF(L[j][k],L[i],1)!=0)
3657	{
3658	break;
3659	}
3660	}
3661	fac=fac+L[j][k];
3662	}
3663	}
3664	// delete superfluous polynomials
3665	fac=simplify(fac,8);
3666	// saturation
3667	SQi,f0,f,fac=minsat_ppd(SI,fac);
3668	I'=I',f;
3669	QP=SQi,L[i],f0,fac;
3670	// the quadrupel:
3671	// a standard basis of Q_i,
3672	// a standard basis of P_i,
3673	// sep_i,
3674	// irreducible factors of
3675	// the "minimal divisor" of the seperator
3676	// as computed by the procedure minsat,
3677	Q[i]=QP;
3678	}
3679	I'=std(I');
3680	return (Q, I');
3681	// I' = remaining component
3682	}
3683
3684
3685	////////////////////////////////////////////////////////////////
3686	// proc extraction
3687	// input: A standard basis of a pseudo primary ideal I, and a standard
3688	// basis of the unique minimal associated prime P of I
3689	// output: an extraction of I, i.e., a standard basis of the primary
3690	// component Q of I with associated prime P, a standard basis of the
3691	// remaining component, and the irreducible factors of the
3692	// "minimal divisor" of the extractor as computed by the procedure minsat
3693	////////////////////////////////////////////////////////////////
3694
3695
3696	proc extraction (ideal SI, ideal SP)
3697	{
3698	list indsets=indepSet(SP,0);
3699	poly f;
3700	if(size(indsets)!=0) //check, whether dim P != 0
3701	{
3702	intvec v; // a maximal independent set of variables
3703	// modulo P
3704	string U; // the independent variables
3705	string A; // the dependent variables
3706	int j,k;
3707	int a; // the size of A
3708	int degf;
3709	ideal g;
3710	list polys;
3711	int sizepolys;
3712	list newpoly;
3713	def R=basering;
3714	//intvec hv=hilb(SI,1);
3715	for (k=1;k<=size(indsets);k++)
3716	{
3717	v=indsets[k];
3718	for (j=1;j<=nvars(R);j++)
3719	{
3720	if (v[j]==1)
3721	{
3722	U=U+varstr(j)+",";
3723	}
3724	else
3725	{
3726	A=A+varstr(j)+",";
3727	a++;
3728	}
3729	}
3730
3731	U[size(U)]=")"; // we compute the extractor of I (w.r.t. U)
3732	execute "ring RAU="+charstr(basering)+",("+A+U+",(dp("+string(a)+"),dp);";
3733	ideal I=imap(R,SI);
3734	//I=std(I,hv); // the standard basis in (R[U])[A]
3735	I=std(I); // the standard basis in (R[U])[A]
3736	A[size(A)]=")";
3737	execute "ring Rloc=("+charstr(basering)+","+U+",("+A+",dp;";
3738	ideal I=imap(RAU,I);
3739	//"std in lokalisierung:"+newline,I;
3740	ideal h;
3741	for(j=ncols(I);j>=1;j--)
3742	{
3743	h[j]=leadcoef(I[j]); // consider I in (R(U))[A]
3744	}
3745	setring R;
3746	g=imap(Rloc,h);
3747	kill RAU,Rloc;
3748	U="";
3749	A="";
3750	a=0;
3751	f=lcm(g);
3752	newpoly[1]=f;
3753	polys=polys+newpoly;
3754	newpoly=list();
3755	}
3756	f=polys[1];
3757	degf=deg(f);
3758	sizepolys=size(polys);
3759	for (k=2;k<=sizepolys;k++)
3760	{
3761	if (deg(polys[k])<degf)
3762	{
3763	f=polys[k];
3764	degf=deg(f);
3765	}
3766	}
3767	}
3768	else
3769	{
3770	f=1;
3771	}
3772	poly f0,h0; ideal SQ; ideal fac;
3773	if(f!=1)
3774	{
3775	SQ,f0,h0,fac=minsat(SI,f);
3776	return(SQ,std(SI+h0),fac);
3777	// the tripel
3778	// a standard basis of Q,
3779	// a standard basis of remaining component,
3780	// irreducible factors of
3781	// the "minimal divisor" of the extractor
3782	// as computed by the procedure minsat
3783	}
3784	else
3785	{
3786	return(SI,ideal(1),ideal(1));
3787	}
3788	}
3789
3790	/////////////////////////////////////////////////////
3791	// proc minsat
3792	// input: a standard basis of an ideal I and a polynomial p
3793	// output: a standard basis IS of the saturation of I w.r. to p,
3794	// the maximal squarefree factor f0 of p,
3795	// the "minimal divisor" f of f0 such that the saturation of
3796	// I w.r. to f equals the saturation of I w.r. to f0 (which is IS),
3797	// the irreducible factors of f
3798	//////////////////////////////////////////////////////////
3799
3800
3801	proc minsat(ideal SI, poly p)
3802	{
3803	ideal fac=factorize(p,1); //the irreducible factors of p
3804	fac=sort(fac)[1];
3805	int i,k;
3806	poly f0=1;
3807	for(i=ncols(fac);i>=1;i--)
3808	{
3809	f0=f0*fac[i];
3810	}
3811	poly f=1;
3812	ideal iold;
3813	list quotM;
3814	quotM[1]=SI;
3815	quotM[2]=fac;
3816	quotM[3]=f0;
3817	// we deal seperately with the first quotient;
3818	// factors, which do not contribute to this one,
3819	// are omitted
3820	iold=quotM[1];
3821	quotM=minquot(quotM);
3822	fac=quotM[2];
3823	if(quotM[3]==1)
3824	{
3825	return(quotM[1],f0,f,fac);
3826	}
3827	while(special_ideals_equal(iold,quotM[1])==0)
3828	{
3829	f=f*quotM[3];
3830	iold=quotM[1];
3831	quotM=minquot(quotM);
3832	}
3833	return(quotM[1],f0,f,fac); // the quadrupel ((I:p),f0,f, irr. factors of f)
3834	}
3835
3836	/////////////////////////////////////////////////////
3837	// proc minsat_ppd
3838	// input: a standard basis of an ideal I and a polynomial p
3839	// output: a standard basis IS of the saturation of I w.r. to p,
3840	// the maximal squarefree factor f0 of p,
3841	// the "minimal divisor" f of f0 such that the saturation of
3842	// I w.r. to f equals the saturation of I w.r. to f0 (which is IS),
3843	// the irreducible factors of f
3844	//////////////////////////////////////////////////////////
3845
3846
3847	proc minsat_ppd(ideal SI, ideal fac)
3848	{
3849	fac=sort(fac)[1];
3850	int i,k;
3851	poly f0=1;
3852	for(i=ncols(fac);i>=1;i--)
3853	{
3854	f0=f0*fac[i];
3855	}
3856	poly f=1;
3857	ideal iold;
3858	list quotM;
3859	quotM[1]=SI;
3860	quotM[2]=fac;
3861	quotM[3]=f0;
3862	// we deal seperately with the first quotient;
3863	// factors, which do not contribute to this one,
3864	// are omitted
3865	iold=quotM[1];
3866	quotM=minquot(quotM);
3867	fac=quotM[2];
3868	if(quotM[3]==1)
3869	{
3870	return(quotM[1],f0,f,fac);
3871	}
3872	while(special_ideals_equal(iold,quotM[1])==0)
3873	{
3874	f=f*quotM[3];
3875	iold=quotM[1];
3876	quotM=minquot(quotM);
3877	k++;
3878	}
3879	return(quotM[1],f0,f,fac); // the quadrupel ((I:p),f0,f, irr. factors of f)
3880	}
3881	/////////////////////////////////////////////////////////////////
3882	// proc minquot
3883	// input: a list with 3 components: a standard basis
3884	// of an ideal I, a set of irreducible polynomials, and
3885	// there product f0
3886	// output: a standard basis of the ideal (I:f0), the irreducible
3887	// factors of the "minimal divisor" f of f0 with (I:f0) = (I:f),
3888	// the "minimal divisor" f
3889	/////////////////////////////////////////////////////////////////
3890
3891	proc minquot(list tsil)
3892	{
3893	int i,j,k,action;
3894	ideal verg;
3895	list l;
3896	poly g;
3897	ideal laedi=tsil[1];
3898	ideal fac=tsil[2];
3899	poly f=tsil[3];
3900
3901	//std
3902	// ideal star=quotient(laedi,f);
3903	// star=std(star);
3904	option(returnSB);
3905	ideal star=quotient(laedi,f);
3906	option(noreturnSB);
3907	if(special_ideals_equal(laedi,star)==1)
3908	{
3909	return(laedi,ideal(1),1);
3910	}
3911	action=1;
3912	while(action==1)
3913	{
3914	if(size(fac)==1)
3915	{
3916	action=0;
3917	break;
3918	}
3919	for(i=1;i<=size(fac);i++)
3920	{
3921	g=1;
3922	for(j=1;j<=size(fac);j++)
3923	{
3924	if(i!=j)
3925	{
3926	g=g*fac[j];
3927	}
3928	}
3929	//std
3930	// verg=quotient(laedi,g);
3931	// verg=std(verg);
3932	option(returnSB);
3933	verg=quotient(laedi,g);
3934	option(noreturnSB);
3935	if(special_ideals_equal(verg,star)==1)
3936	{
3937	f=g;
3938	fac[i]=0;
3939	fac=simplify(fac,2);
3940	break;
3941	}
3942	if(i==size(fac))
3943	{
3944	action=0;
3945	}
3946	}
3947	}
3948	l=star,fac,f;
3949	return(l);
3950	}
3951	/////////////////////////////////////////////////
3952	// proc special_ideals_equal
3953	// input: standard bases of ideal k1 and k2 such that
3954	// k1 is contained in k2, or k2 is contained ink1
3955	// output: 1, if k1 equals k2, 0 otherwise
3956	//////////////////////////////////////////////////
3957
3958	proc special_ideals_equal( ideal k1, ideal k2)
3959	{
3960	int j;
3961	if(size(k1)==size(k2))
3962	{
3963	for(j=1;j<=size(k1);j++)
3964	{
3965	if(leadexp(k1[j])!=leadexp(k2[j]))
3966	{
3967	return(0);
3968	}
3969	}
3970	return(1);
3971	}
3972	return(0);
3973	}
3974
3975
3976	///////////////////////////////////////////////////////////////////////////////
3977
3978	proc convList(list l)
3979	{
3980	int i;
3981	list re,he;
3982	for(i=1;i<=size(l)/2;i++)
3983	{
3984	he=l[2i-1],l[2i];
3985	re[i]=he;
3986	}
3987	return(re);
3988	}
3989	///////////////////////////////////////////////////////////////////////////////
3990
3991	proc reconvList(list l)
3992	{
3993	int i;
3994	list re;
3995	for(i=1;i<=size(l);i++)
3996	{
3997	re[2*i-1]=l[i][1];
3998	re[2*i]=l[i][2];
3999	}
4000	return(re);
4001	}
4002
4003	///////////////////////////////////////////////////////////////////////////////
4004	//
4005	// The main procedures
4006	//
4007	///////////////////////////////////////////////////////////////////////////////
4008
4009	proc primdecGTZ(ideal i)
4010	"USAGE: primdecGTZ(i); i ideal
4011	RETURN: a list, say pr, of primary ideals and their associated primes
4012	pr[i][1], resp. pr[i][2] is the i-th primary resp. prime component
4013	NOTE: Algorithm of Gianni, Traeger, Zacharias
4014	designed for characteristic 0, works also in char k > 0, if it
4015	terminates (may result in an infinite loop in small characteristic!)
4016	EXAMPLE: example primdecGTZ; shows an example
4017	"
4018	{
4019	return(convList(decomp(i)));
4020	}
4021	example
4022	{ "EXAMPLE:"; echo = 2;
4023	ring r = 32003,(x,y,z),lp;
4024	poly p = z2+1;
4025	poly q = z4+2;
4026	ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
4027	list pr = primdecGTZ(i);
4028	pr;
4029	}
4030	///////////////////////////////////////////////////////////////////////////////
4031
4032	proc primdecSY(ideal i, list #))
4033	"USAGE: primdecSY(i); i ideal, c int
4034	if c=0, the given ordering of the variables is used.
4035	if c=1, minAssChar tries to use an optimal ordering,
4036	if c=2, minAssGTZ is used
4037	if c=3, minAssGTZ and facstd is used
4038	RETURN: a list, say pr, of primary ideals and their associated primes
4039	pr[i][1], resp. pr[i][2] is the i-th primary resp. prime component
4040	NOTE: Algorithm of Shimoyama-Yokoyama
4041	implemented for characteristic 0, works also in char k > 0,
4042	the result may be not completely decomposed in small characteristic
4043	EXAMPLE: example primdecSY; shows an example
4044	"
4045	{
4046	if (size(#)==1)
4047	{ return(prim_dec(i,#[1])); }
4048	else
4049	{ return(prim_dec(i,1)); }
4050	}
4051	example
4052	{ "EXAMPLE:"; echo = 2;
4053	ring r = 32003,(x,y,z),lp;
4054	poly p = z2+1;
4055	poly q = z4+2;
4056	ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
4057	list pr = primdecSY(i);
4058	pr;
4059	}
4060	///////////////////////////////////////////////////////////////////////////////
4061	proc minAssGTZ(ideal i)
4062	"USAGE: minAssGTZ(i); i ideal
4063	RETURN: list = the minimal associated prime ideals of i
4064	NOTE: designed for characteristic 0, works also in char k > 0 if it
4065	terminates, may result in an infinite loop in small characteristic
4066	EXAMPLE: example minAssGTZ; shows an example
4067	"
4068	{
4069	return(minAssPrimes(i,1));
4070	}
4071	example
4072	{ "EXAMPLE:"; echo = 2;
4073	ring r = 32003,(x,y,z),dp;
4074	poly p = z2+1;
4075	poly q = z4+2;
4076	ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
4077	list pr= minAssGTZ(i);
4078	pr;
4079	}
4080
4081	///////////////////////////////////////////////////////////////////////////////
4082	proc minAssChar(ideal i, list #)
4083	"USAGE: minAssChar(i[,c]); i ideal,
4084	if c=0, the given ordering of the variables is used.
4085	Otherwise, the system tries to find an optimal ordering,
4086	which in some cases may considerably speed up the algorithm
4087	RETURN: list = the minimal associated prime ideals of i
4088	NOTE: implemented for characteristic 0, works also in char k >> 0,
4089	the result may be not compltely decomposed in small characteristic
4090	EXAMPLE: example minAssChar; shows an example
4091	"
4092	{
4093	if (size(#)==1)
4094	{ return(min_ass_prim_charsets(i,#[1])); }
4095	else
4096	{ return(min_ass_prim_charsets(i,1)); }
4097	}
4098	example
4099	{ "EXAMPLE:"; echo = 2;
4100	ring r = 32003,(x,y,z),dp;
4101	poly p = z2+1;
4102	poly q = z4+2;
4103	ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
4104	list pr= minAssChar(i);
4105	pr;
4106	}
4107	///////////////////////////////////////////////////////////////////////////////
4108	proc equiRadical(ideal i)
4109	"USAGE: equiRadical(i); i ideal
4110	RETURN: ideal, intersection of associated primes of i of maximal dimension
4111	NOTE: designed for characteristic 0, works also in char k > 0 if it
4112	terminates, may result in an infinite loop in small characteristic
4113	EXAMPLE: example equiRadical; shows an example
4114	"
4115	{
4116	return(radical(i,1));
4117	}
4118	example
4119	{ "EXAMPLE:"; echo = 2;
4120	ring r = 32003,(x,y,z),dp;
4121	poly p = z2+1;
4122	poly q = z4+2;
4123	ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
4124	ideal pr= equiRadical(i);
4125	pr;
4126	}
4127	///////////////////////////////////////////////////////////////////////////////
4128	proc radical(ideal i,list #)
4129	"USAGE: radical(i); i ideal
4130	RETURN: ideal = the radical of i
4131	NOTE: a combination of the algorithms of Krick/Logar and
4132	Eisenbud/Huneke/Vasconcelos
4133	designed for characteristic 0, works also in char k > 0 if it
4134	terminates, may result in an infinite loop in small characteristic
4135	EXAMPLE: example radical; shows an example
4136	"
4137	{
4138	def @P=basering;
4139	int j,il;
4140	if(size(#)>0)
4141	{
4142	il=#[1];
4143	}
4144	ideal re=1;
4145	option(redSB);
4146	list qr=simplifyIdeal(i);
4147	map phi=@P,qr[2];
4148	i=qr[1];
4149
4150	list pr=facstd(i);
4151
4152	if(size(pr)==1)
4153	{
4154	attrib(pr[1],"isSB",1);
4155	if((dim(pr[1])==0)&&(homog(pr[1])==1))
4156	{
4157	ideal @res=maxideal(1);
4158	return(phi(@res));
4159	}
4160	if(dim(pr[1])>1)
4161	{
4162	execute "ring gnir = ("+charstr(basering)+"),
4163	("+varstr(basering)+"),(C,lp);";
4164	ideal i=fetch(@P,i);
4165	list @pr=facstd(i);
4166	setring @P;
4167	pr=fetch(gnir,@pr);
4168	}
4169	}
4170	option(noredSB);
4171	int s=size(pr);
4172
4173	if(s==1)
4174	{
4175	i=radicalEHV(i,ideal(1),il);
4176	return(phi(i));
4177	}
4178	intvec pos;
4179	pos[s]=0;
4180	if(il==1)
4181	{
4182	int ndim,k;
4183	attrib(pr[1],"isSB",1);
4184	int odim=dim(pr[1]);
4185	int count=1;
4186
4187	for(j=2;j<=s;j++)
4188	{
4189	attrib(pr[j],"isSB",1);
4190	ndim=dim(pr[j]);
4191	if(ndim>odim)
4192	{
4193	for(k=count;k<=j-1;k++)
4194	{
4195	pos[k]=1;
4196	}
4197	count=j;
4198	odim=ndim;
4199	}
4200	if(ndim<odim)
4201	{
4202	pos[j]=1;
4203	}
4204	}
4205	}
4206	for(j=1;j<=s;j++)
4207	{
4208	if(pos[s+1-j]==0)
4209	{
4210	re=intersect(re,radicalEHV(pr[s+1-j],re,il));
4211	}
4212	}
4213	re=interred(re);
4214	return(phi(re));
4215	}
4216	example
4217	{ "EXAMPLE:"; echo = 2;
4218	ring r = 32003,(x,y,z),dp;
4219	poly p = z2+1;
4220	poly q = z4+2;
4221	ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
4222	ideal pr= radical(i);
4223	pr;
4224	}
4225	///////////////////////////////////////////////////////////////////////////////
4226	proc prepareAss(ideal i)
4227	"USAGE: prepareAss(i); i ideal
4228	RETURN: list = the radicals of the maximal dimensional components of i
4229	NOTE: uses algorithm of Eisenbud, Huneke and Vasconcelos
4230	EXAMPLE: example prepareAss; shows an example
4231	"
4232	{
4233	ideal j=std(i);
4234	int cod=nvars(basering)-dim(j);
4235	int e;
4236	list er;
4237	ideal ann;
4238	if(homog(i)==1)
4239	{
4240	list re=sres(i,0); //the resolution
4241	re=minres(re); //minimized resolution
4242	}
4243	else
4244	{
4245	list re=mres(i,0);
4246	}
4247	for(e=cod;e<=nvars(basering);e++)
4248	{
4249	ann=AnnExt_R(e,re);
4250
4251	if(nvars(basering)-dim(std(ann))==e)
4252	{
4253	er[size(er)+1]=equiRadical(ann);
4254	}
4255	}
4256	return(er);
4257	}
4258	example
4259	{ "EXAMPLE:"; echo = 2;
4260	ring r = 32003,(x,y,z),dp;
4261	poly p = z2+1;
4262	poly q = z4+2;
4263	ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
4264	list pr= prepareAss(i);
4265	pr;
4266	}
4267
4268	proc testPrimary(list pr, ideal k)
4269	"USAGE: testPrimary(pr,k); pr a list, result of primdecGTZ(k) or primdecSY(k)
4270	RETURN: int, 1 if intersection of the primary ideals in pr is k, 0 if not
4271	EXAMPLE: example testPrimary; shows an example
4272	"
4273	{
4274	int i;
4275	pr=reconvList(pr);
4276	ideal j=pr[1];
4277	for (i=2;i<=size(pr)/2;i++)
4278	{
4279	j=intersect(j,pr[2*i-1]);
4280	}
4281	return(idealsEqual(j,k));
4282	}
4283	example
4284	{ "EXAMPLE:"; echo = 2;
4285	ring r = 32003,(x,y,z),dp;
4286	poly p = z2+1;
4287	poly q = z4+2;
4288	ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
4289	list pr = primdecGTZ(i);
4290	testPrimary(pr,i);
4291	}
4292

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format