Context Navigation

primdec.lib @ 5d21375

fieker-DuValspielwiese

Last change on this file since 5d21375 was b9b906, checked in by Hans Schönemann <hannes@…>, 23 years ago
*hannes: lib format revisited git-svn-id: file:///usr/local/Singular/svn/trunk@5078 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 122.2 KB

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

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

Download in other formats:

Original Format