Context Navigation

primdec.lib @ 9491d6

spielwiese

Last change on this file since 9491d6 was 9491d6, checked in by Hans Schönemann <hannes@…>, 18 years ago
*pfister/hannes: bug_42 fixed git-svn-id: file:///usr/local/Singular/svn/trunk@9125 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 144.8 KB

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

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

Download in other formats:

Original Format