Context Navigation

primdec.lib @ 9050ca

spielwiese

Last change on this file since 9050ca was 9050ca, checked in by Hans Schönemann <hannes@…>, 25 years ago
*hannes: GMG's fixes git-svn-id: file:///usr/local/Singular/svn/trunk@3360 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 102.0 KB

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

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

Download in other formats:

Original Format