Context Navigation

primdec.lib @ 7e1114d

fieker-DuValspielwiese

Last change on this file since 7e1114d was 3e23be, checked in by Hans Schönemann <hannes@…>, 25 years ago
*hannes: added () for parametric rings git-svn-id: file:///usr/local/Singular/svn/trunk@2940 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 101.1 KB

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

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

Download in other formats:

Original Format