Context Navigation

primdec.lib @ d950c5

spielwiese

Last change on this file since d950c5 was d950c5, checked in by Hans Schönemann <hannes@…>, 26 years ago
* hannes: version 1.24 computes >1 h, bachk to 1.22 (primdec.lib) git-svn-id: file:///usr/local/Singular/svn/trunk@2428 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 100.9 KB

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

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

Download in other formats:

Original Format