Context Navigation

primdec.lib @ bf241b

spielwiese

Last change on this file since bf241b was bf241b, checked in by Hans Schönemann <hannes@…>, 22 years ago
*hannes: handle primdecSY(0) git-svn-id: file:///usr/local/Singular/svn/trunk@5776 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 123.6 KB

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

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

Download in other formats:

Original Format