Context Navigation

primdec.lib @ 63be42

spielwiese

Last change on this file since 63be42 was a3432c, checked in by Hans Schönemann <hannes@…>, 26 years ago
minor optimizations git-svn-id: file:///usr/local/Singular/svn/trunk@2519 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 100.9 KB

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

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

Download in other formats:

Original Format