Context Navigation

primdec.lib @ f0daaa2

spielwiese

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

Line
1	///////////////////////////////////////////////////////////////////////////////
2	version="$Id: primdec.lib,v 1.115 2006-05-15 10:59:17 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	if (deg(quprimary[1][1])<=0) { @n=0; }
3511	for (@k=1;@k<=size(htprimary);@k++)
3512	{
3513	quprimary[@n+@k]=htprimary[@k];
3514	}
3515	}
3516	}
3517
3518	}
3519	else
3520	{
3521	if(abspri)
3522	{
3523	list resu,tempo;
3524	for(ab=1;ab<=size(quprimary)/2;ab++)
3525	{
3526	tempo=quprimary[2ab-1],quprimary[2ab],
3527	absprimary[ab],abskeep[ab];
3528	resu[ab]=tempo;
3529	}
3530	quprimary=resu;
3531	}
3532	}
3533	//---------------------------------------------------------------------------
3534	//back to the ring we started with
3535	//the final result: primary
3536	//---------------------------------------------------------------------------
3537
3538	setring @P;
3539	primary=imap(gnir,quprimary);
3540	return(primary);
3541	}
3542
3543
3544	example
3545	{ "EXAMPLE:"; echo = 2;
3546	ring r = 32003,(x,y,z),lp;
3547	poly p = z2+1;
3548	poly q = z4+2;
3549	ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
3550	list pr= decomp(i);
3551	pr;
3552	testPrimary( pr, i);
3553	}
3554
3555	///////////////////////////////////////////////////////////////////////////////
3556	static proc powerCoeffs(poly f,int e)
3557	//computes a polynomial with the same monomials as f but coefficients
3558	//the p^e th power of the coefficients of f
3559	{
3560	int i;
3561	poly g;
3562	int ex=char(basering)^e;
3563	for(i=1;i<=size(f);i++)
3564	{
3565	g=g+leadcoef(f[i])^ex*leadmonom(f[i]);
3566	}
3567	return(g);
3568	}
3569	///////////////////////////////////////////////////////////////////////////////
3570
3571	proc sep(poly f,int i, list #)
3572	"USAGE: input: a polynomial f depending on the i-th variable and optional
3573	an integer k considering the polynomial f defined over Fp(t1,...,tm)
3574	as polynomial over Fp(t(1)^(p^-k),...,t(m)^(p^-k))
3575	RETURN: the separabel part of f as polynomial in Fp(t1,...,tm)
3576	and an integer k to indicate that f should be considerd
3577	as polynomial over Fp(t(1)^(p^-k),...,t(m)^(p^-k))
3578	EXAMPLE: example sep; shows an example
3579	{
3580	def R=basering;
3581	int k;
3582	if(size(#)>0){k=#[1];}
3583
3584
3585	poly h=gcd(f,diff(f,var(i)));
3586	if((reduce(f,std(h))!=0)\|\|(reduce(diff(f,var(i)),std(h))!=0))
3587	{
3588	ERROR("FEHLER IN GCD");
3589	}
3590	poly g1=lift(h,f)[1][1]; // f/h
3591	poly h1;
3592
3593	while(h!=h1)
3594	{
3595	h1=h;
3596	h=gcd(h,diff(h,var(i)));
3597	}
3598
3599	if(deg(h1)==0){return(list(g1,k));} //in characteristic 0 we return here
3600
3601	k++;
3602
3603	ideal ma=maxideal(1);
3604	ma[i]=var(i)^char(R);
3605	map phi=R,ma;
3606	ideal hh=h; //this is technical because preimage works only for ideals
3607
3608	poly u=preimage(R,phi,hh)[1]; //h=u(x(i)^p)
3609
3610	list g2=sep(u,i,k); //we consider u(t(1)^(p^-1),...,t(m)^(p^-1))
3611	g1=powerCoeffs(g1,g2[2]-k+1); //to have g1 over the same field as g2[1]
3612
3613	list g3=sep(g1*g2[1],i,g2[2]);
3614	return(g3);
3615	}
3616	example
3617	{ "EXAMPLE:"; echo = 2;
3618	ring R=(5,t,s),(x,y,z),dp;
3619	poly f=(x^25-tx^5+t)(x^3+s);
3620	sep(f,1);
3621	}
3622
3623	///////////////////////////////////////////////////////////////////////////////
3624	proc zeroRad(ideal I,list #)
3625	"USAGE: zeroRad(I) , I a zero-dimensional ideal
3626	RETURN: the radical of I
3627	NOTE: Algorithm of Kemper
3628	EXAMPLE: example zeroRad; shows an example
3629	{
3630	if(homog(I)==1){return(maxideal(1));}
3631	//I needs to be a reduced standard basis
3632	def R=basering;
3633	int m=npars(R);
3634	int n=nvars(R);
3635	int p=char(R);
3636	int d=vdim(I);
3637	int i,k;
3638	list l;
3639	if(((p==0)\|\|(p>d))&&(d==deg(I[1])))
3640	{
3641	intvec e=leadexp(I[1]);
3642	for(i=1;i<=nvars(basering);i++)
3643	{
3644	if(e[i]!=0) break;
3645	}
3646	I[1]=sep(I[1],i)[1];
3647	return(interred(I));
3648	}
3649	intvec op=option(get);
3650
3651	option(redSB);
3652	ideal F=finduni(I);//F[i] generates I intersected with K[var(i)]
3653
3654	option(set,op);
3655	if(size(#)>0){I=#[1];}
3656
3657	for(i=1;i<=n;i++)
3658	{
3659	l[i]=sep(F[i],i);
3660	F[i]=l[i][1];
3661	if(l[i][2]>k){k=l[i][2];} //computation of the maximal k
3662	}
3663
3664	if((k==0)\|\|(m==0)){return(interred(I+F));} //the separable case
3665
3666	for(i=1;i<=n;i++) //consider all polynomials over
3667	{ //Fp(t(1)^(p^-k),...,t(m)^(p^-k))
3668	F[i]=powerCoeffs(F[i],k-l[i][2]);
3669	}
3670
3671	string cR="ring @R="+string(p)+",("+parstr(R)+","+varstr(R)+"),dp;";
3672	execute(cR);
3673	ideal F=imap(R,F);
3674
3675	string nR="ring @S="+string(p)+",(y(1..m),"+varstr(R)+","+parstr(R)+"),dp;";
3676	execute(nR);
3677
3678	ideal G=fetch(@R,F); //G[i](t(1)^(p^-k),...,t(m)^(p^-k),x(i))=sep(F[i])
3679
3680	ideal I=imap(R,I);
3681	ideal J=I+G;
3682	poly el=1;
3683	k=p^k;
3684	for(i=1;i<=m;i++)
3685	{
3686	J=J,var(i)^k-var(m+n+i);
3687	el=el*y(i);
3688	}
3689
3690	J=eliminate(J,el);
3691	setring R;
3692	ideal J=imap(@S,J);
3693	return(J);
3694	}
3695	example
3696	{ "EXAMPLE:"; echo = 2;
3697	ring R=(5,t),(x,y),dp;
3698	ideal I=x^5-t,y^5-t;
3699	zeroRad(I);
3700	}
3701
3702	///////////////////////////////////////////////////////////////////////////////
3703	static proc radicalKL (ideal i,ideal ser,list #)
3704	{
3705	attrib(i,"isSB",1); // i needs to be a reduced standard basis
3706	list indep,fett;
3707	intvec @w,@hilb,op;
3708	int @wr,@n,@m,lauf,di;
3709	ideal fac,@h,collectrad,lsau;
3710	poly @q;
3711	string @va,quotring;
3712
3713	def @P = basering;
3714	int jdim=dim(i);
3715	int homo=homog(i);
3716	ideal rad=ideal(1);
3717	ideal te=ser;
3718	if(size(#)>0)
3719	{
3720	@wr=#[1];
3721	}
3722	if(homo==1)
3723	{
3724	for(@n=1;@n<=nvars(basering);@n++)
3725	{
3726	@w[@n]=ord(var(@n));
3727	}
3728	@hilb=hilb(i,1,@w);
3729	}
3730
3731
3732	//---------------------------------------------------------------------------
3733	//j is the ring
3734	//---------------------------------------------------------------------------
3735
3736	if (jdim==-1)
3737	{
3738
3739	return(ideal(1));
3740	}
3741
3742	//---------------------------------------------------------------------------
3743	//the zero-dimensional case
3744	//---------------------------------------------------------------------------
3745
3746	if (jdim==0)
3747	{
3748	return(zeroRad(i));
3749	}
3750	//-------------------------------------------------------------------------
3751	//search for a maximal independent set indep,i.e.
3752	//look for subring such that the intersection with the ideal is zero
3753	//j intersected with K[var(indep[3]+1),...,var(nvar] is zero,
3754	//indep[1] is the new varstring, indep[2] the string for the block-ordering
3755	//-------------------------------------------------------------------------
3756
3757	indep=maxIndependSet(i);
3758
3759	for(@m=1;@m<=size(indep);@m++)
3760	{
3761	if((indep[@m][1]==varstr(basering))&&(@m==1))
3762	//this is the good case, nothing to do, just to have the same notations
3763	//change the ring
3764	{
3765	execute("ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
3766	+ordstr(basering)+");");
3767	ideal @j=fetch(@P,i);
3768	attrib(@j,"isSB",1);
3769	}
3770	else
3771	{
3772	@va=string(maxideal(1));
3773	execute("ring gnir1 = ("+charstr(basering)+"),("+indep[@m][1]+"),("
3774	+indep[@m][2]+");");
3775	execute("map phi=@P,"+@va+";");
3776	if(homo==1)
3777	{
3778	ideal @j=std(phi(i),@hilb,@w);
3779	}
3780	else
3781	{
3782	ideal @j=groebner(phi(i));
3783	}
3784	}
3785	if((deg(@j[1])==0)\|\|(dim(@j)<jdim))
3786	{
3787	setring @P;
3788	break;
3789	}
3790	for (lauf=1;lauf<=size(@j);lauf++)
3791	{
3792	fett[lauf]=size(@j[lauf]);
3793	}
3794	//------------------------------------------------------------------------
3795	//we have now the following situation:
3796	//j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
3797	//to this quotientring, j is their still a standardbasis, the
3798	//leading coefficients of the polynomials there (polynomials in
3799	//K[var(nnp+1),..,var(nva)]) are collected in the list h,
3800	//we need their ggt, gh, because of the following:
3801	//let (j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..rest..]
3802	//intersected with K[var(1),...,var(nva)] is (j:gh^n)
3803	//on the other hand j=(j,gh^n) intersected with (j:gh^n)
3804
3805	//------------------------------------------------------------------------
3806
3807	//the arrangement for the quotientring K(var(nnp+1),..,var(nva))[..rest..]
3808	//and the map phi:K[var(1),...,var(nva)] ----->
3809	//K(var(nnpr+1),..,var(nva))[..the rest..]
3810	//------------------------------------------------------------------------
3811	quotring=prepareQuotientring(nvars(basering)-indep[@m][3]);
3812
3813	//------------------------------------------------------------------------
3814	//we pass to the quotientring K(var(nnp+1),..,var(nva))[..the rest..]
3815	//------------------------------------------------------------------------
3816
3817	execute(quotring);
3818
3819	// @j considered in the quotientring
3820	ideal @j=imap(gnir1,@j);
3821
3822	kill gnir1;
3823
3824	//j is a standardbasis in the quotientring but usually not minimal
3825	//here it becomes minimal
3826
3827	@j=clearSB(@j,fett);
3828
3829	//we need later ggt(h[1],...)=gh for saturation
3830	ideal @h;
3831	if(deg(@j[1])>0)
3832	{
3833	for(@n=1;@n<=size(@j);@n++)
3834	{
3835	@h[@n]=leadcoef(@j[@n]);
3836	}
3837	op=option(get);
3838	option(redSB);
3839	@j=interred(@j); //to obtain a reduced standardbasis
3840	attrib(@j,"isSB",1);
3841	option(set,op);
3842
3843	ideal zero_rad= zeroRad(@j);
3844	}
3845	else
3846	{
3847	ideal zero_rad=ideal(1);
3848	}
3849
3850	//we need the intersection of the ideals in the list quprimary with the
3851	//polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
3852	//but fi polynomials, then the intersection of q with the polynomialring
3853	//is the saturation of the ideal generated by f1,...,fr with respect to
3854	//h which is the lcm of the leading coefficients of the fi considered in
3855	//the quotientring: this is coded in saturn
3856
3857	zero_rad=std(zero_rad);
3858
3859	ideal hpl;
3860
3861	for(@n=1;@n<=size(zero_rad);@n++)
3862	{
3863	hpl=hpl,leadcoef(zero_rad[@n]);
3864	}
3865
3866	//------------------------------------------------------------------------
3867	//we leave the quotientring K(var(nnp+1),..,var(nva))[..the rest..]
3868	//back to the polynomialring
3869	//------------------------------------------------------------------------
3870	setring @P;
3871
3872	collectrad=imap(quring,zero_rad);
3873	lsau=simplify(imap(quring,hpl),2);
3874	@h=imap(quring,@h);
3875
3876	kill quring;
3877
3878
3879	//here the intersection with the polynomialring
3880	//mentioned above is really computed
3881
3882	collectrad=sat2(collectrad,lsau)[1];
3883	if(deg(@h[1])>=0)
3884	{
3885	fac=ideal(0);
3886	for(lauf=1;lauf<=ncols(@h);lauf++)
3887	{
3888	if(deg(@h[lauf])>0)
3889	{
3890	fac=fac+factorize(@h[lauf],1);
3891	}
3892	}
3893	fac=simplify(fac,4);
3894	@q=1;
3895	for(lauf=1;lauf<=size(fac);lauf++)
3896	{
3897	@q=@q*fac[lauf];
3898	}
3899	op=option(get);
3900	option(returnSB);
3901	option(redSB);
3902	i=quotient(i+ideal(@q),rad);
3903	attrib(i,"isSB",1);
3904	option(set,op);
3905
3906	}
3907	if((deg(rad[1])>0)&&(deg(collectrad[1])>0))
3908	{
3909	rad=intersect(rad,collectrad);
3910	te=intersect(te,collectrad);
3911	te=simplify(reduce(te,i,1),2);
3912	}
3913	else
3914	{
3915	if(deg(collectrad[1])>0)
3916	{
3917	rad=collectrad;
3918	te=intersect(te,collectrad);
3919	te=simplify(reduce(te,i,1),2);
3920	}
3921	}
3922
3923	if((dim(i)<jdim)\|\|(size(te)==0))
3924	{
3925	break;
3926	}
3927	if(homo==1)
3928	{
3929	@hilb=hilb(i,1,@w);
3930	}
3931	}
3932	if(((@wr==1)&&(dim(i)<jdim))\|\|(deg(i[1])==0)\|\|(size(te)==0))
3933	{
3934	return(rad);
3935	}
3936	rad=intersect(rad,radicalKL(i,ideal(1),@wr));
3937	return(rad);
3938	}
3939
3940	///////////////////////////////////////////////////////////////////////////////
3941
3942	proc radicalEHV(ideal i)
3943	"USAGE: radicalEHV(i); i ideal.
3944	RETURN: ideal, the radical of i.
3945	NOTE: Uses the algorithm of Eisenbud/Huneke/Vasconcelos, which
3946	reduces the computation to the complete intersection case,
3947	by taking, in the general case, a generic linear combination
3948	of the input.
3949	Works only in characteristic 0 or p large.
3950	EXAMPLE: example radicalEHV; shows an example
3951	"
3952	{
3953	if(ord_test(basering)!=1)
3954	{
3955	ERROR(
3956	"// Not implemented for this ordering, please change to global ordering."
3957	);
3958	}
3959	if((char(basering)<100)&&(char(basering)!=0))
3960	{
3961	"WARNING: The characteristic is too small, the result may be wrong";
3962	}
3963	ideal J,I,I0,radI0,L,radI1,I2,radI2;
3964	int l,n;
3965	intvec op=option(get);
3966	matrix M;
3967
3968	option(redSB);
3969	list m=mstd(i);
3970	I=m[2];
3971	option(set,op);
3972
3973	int cod=nvars(basering)-dim(m[1]);
3974	//-------------------complete intersection case:----------------------
3975	if(cod==size(m[2]))
3976	{
3977	J=minor(jacob(I),cod);
3978	return(quotient(I,J));
3979	}
3980	//-----first codim elements of I are a complete intersection:---------
3981	for(l=1;l<=cod;l++)
3982	{
3983	I0[l]=I[l];
3984	}
3985	n=dim(std(I0))+cod-nvars(basering);
3986	//-----last codim elements of I are a complete intersection:----------
3987	if(n!=0)
3988	{
3989	for(l=1;l<=cod;l++)
3990	{
3991	I0[l]=I[size(I)-l+1];
3992	}
3993	n=dim(std(I0))+cod-nvars(basering);
3994	}
3995	//-----taking a generic linear combination of the input:--------------
3996	if(n!=0)
3997	{
3998	M=transpose(sparsetriag(size(m[2]),cod,95,1));
3999	I0=ideal(M*transpose(I));
4000	n=dim(std(I0))+cod-nvars(basering);
4001	}
4002	//-----taking a more generic linear combination of the input:---------
4003	if(n!=0)
4004	{
4005	M=transpose(sparsetriag(size(m[2]),cod,0,100));
4006	I0=ideal(M*transpose(I));
4007	n=dim(std(I0))+cod-nvars(basering);
4008	}
4009	if(n==0)
4010	{
4011	J=minor(jacob(I0),cod);
4012	radI0=quotient(I0,J);
4013	L=quotient(radI0,I);
4014	radI1=quotient(radI0,L);
4015
4016	if(size(reduce(radI1,m[1],1))==0)
4017	{
4018	return(I);
4019	}
4020
4021	I2=sat(I,radI1)[1];
4022
4023	if(deg(I2[1])<=0)
4024	{
4025	return(radI1);
4026	}
4027	return(intersect(radI1,radicalEHV(I2)));
4028	}
4029	//---------------------general case-------------------------------------
4030	return(radical(I));
4031	}
4032	example
4033	{ "EXAMPLE:"; echo = 2;
4034	ring r = 0,(x,y,z),dp;
4035	poly p = z2+1;
4036	poly q = z3+2;
4037	ideal i = p*q^2,y-z2;
4038	ideal pr= radicalEHV(i);
4039	pr;
4040	}
4041
4042	///////////////////////////////////////////////////////////////////////////////
4043
4044	proc Ann(module M)
4045	"USAGE: Ann(M); M module
4046	RETURN: ideal, the annihilator of coker(M)
4047	NOTE: The output is the ideal of all elements a of the basering R such that
4048	a * R^m is contained in M (m=number of rows of M).
4049	EXAMPLE: example Ann; shows an example
4050	"
4051	{
4052	M=prune(M); //to obtain a small embedding
4053	ideal ann=quotient1(M,freemodule(nrows(M)));
4054	return(ann);
4055	}
4056	example
4057	{ "EXAMPLE:"; echo = 2;
4058	ring r = 0,(x,y,z),lp;
4059	module M = x2-y2,z3;
4060	Ann(M);
4061	M = [1,x2],[y,x];
4062	Ann(M);
4063	qring Q=std(xy-1);
4064	module M=imap(r,M);
4065	Ann(M);
4066	}
4067
4068	///////////////////////////////////////////////////////////////////////////////
4069
4070	//computes the equidimensional part of the ideal i of codimension e
4071	static proc int_ass_primary_e(ideal i, int e)
4072	{
4073	if(homog(i)!=1)
4074	{
4075	i=std(i);
4076	}
4077	list re=sres(i,0); //the resolution
4078	re=minres(re); //minimized resolution
4079	ideal ann=AnnExt_R(e,re);
4080	if(nvars(basering)-dim(std(ann))!=e)
4081	{
4082	return(ideal(1));
4083	}
4084	return(ann);
4085	}
4086
4087	///////////////////////////////////////////////////////////////////////////////
4088
4089	//computes the annihilator of Ext^n(R/i,R) with given resolution re
4090	//n is not necessarily the number of variables
4091	static proc AnnExt_R(int n,list re)
4092	{
4093	if(n<nvars(basering))
4094	{
4095	matrix f=transpose(re[n+1]); //Hom(_,R)
4096	module k=nres(f,2)[2]; //the kernel
4097	matrix g=transpose(re[n]); //the image of Hom(_,R)
4098
4099	ideal ann=quotient1(g,k); //the anihilator
4100	}
4101	else
4102	{
4103	ideal ann=Ann(transpose(re[n]));
4104	}
4105	return(ann);
4106	}
4107	///////////////////////////////////////////////////////////////////////////////
4108
4109	static proc analyze(list pr)
4110	{
4111	int ii,jj;
4112	for(ii=1;ii<=size(pr)/2;ii++)
4113	{
4114	dim(std(pr[2*ii]));
4115	idealsEqual(pr[2ii-1],pr[2ii]);
4116	"===========================";
4117	}
4118
4119	for(ii=size(pr)/2;ii>1;ii--)
4120	{
4121	for(jj=1;jj<ii;jj++)
4122	{
4123	if(size(reduce(pr[2jj],std(pr[2ii],1)))==0)
4124	{
4125	"eingebette Komponente";
4126	jj;
4127	ii;
4128	}
4129	}
4130	}
4131	}
4132
4133	///////////////////////////////////////////////////////////////////////////////
4134	//
4135	// Shimoyama-Yokoyama
4136	//
4137	///////////////////////////////////////////////////////////////////////////////
4138
4139	static proc simplifyIdeal(ideal i)
4140	{
4141	def r=basering;
4142
4143	int j,k;
4144	map phi;
4145	poly p;
4146
4147	ideal iwork=i;
4148	ideal imap1=maxideal(1);
4149	ideal imap2=maxideal(1);
4150
4151
4152	for(j=1;j<=nvars(basering);j++)
4153	{
4154	for(k=1;k<=size(i);k++)
4155	{
4156	if(deg(iwork[k]/var(j))==0)
4157	{
4158	p=-1/leadcoef(iwork[k]/var(j))*iwork[k];
4159	imap1[j]=p+2*var(j);
4160	phi=r,imap1;
4161	iwork=phi(iwork);
4162	iwork=subst(iwork,var(j),0);
4163	iwork[k]=var(j);
4164	imap1=maxideal(1);
4165	imap2[j]=-p;
4166	break;
4167	}
4168	}
4169	}
4170	return(iwork,imap2);
4171	}
4172
4173
4174	///////////////////////////////////////////////////////
4175	// ini_mod
4176	// input: a polynomial p
4177	// output: the initial term of p as needed
4178	// in the context of characteristic sets
4179	//////////////////////////////////////////////////////
4180
4181	static proc ini_mod(poly p)
4182	{
4183	if (p==0)
4184	{
4185	return(0);
4186	}
4187	int n; matrix m;
4188	for( n=nvars(basering); n>0; n=n-1)
4189	{
4190	m=coef(p,var(n));
4191	if(m[1,1]!=1)
4192	{
4193	p=m[2,1];
4194	break;
4195	}
4196	}
4197	if(deg(p)==0)
4198	{
4199	p=0;
4200	}
4201	return(p);
4202	}
4203	///////////////////////////////////////////////////////
4204	// min_ass_prim_charsets
4205	// input: generators of an ideal PS and an integer cho
4206	// If cho=0, the given ordering of the variables is used.
4207	// Otherwise, the system tries to find an "optimal ordering",
4208	// which in some cases may considerably speed up the algorithm
4209	// output: the minimal associated primes of PS
4210	// algorithm: via characteriostic sets
4211	//////////////////////////////////////////////////////
4212
4213
4214	static proc min_ass_prim_charsets (ideal PS, int cho)
4215	{
4216	if((cho<0) and (cho>1))
4217	{
4218	"ERROR: <int> must be 0 or 1"
4219	return();
4220	}
4221	if(system("version")>933)
4222	{
4223	option(notWarnSB);
4224	}
4225	if(cho==0)
4226	{
4227	return(min_ass_prim_charsets0(PS));
4228	}
4229	else
4230	{
4231	return(min_ass_prim_charsets1(PS));
4232	}
4233	}
4234	///////////////////////////////////////////////////////
4235	// min_ass_prim_charsets0
4236	// input: generators of an ideal PS
4237	// output: the minimal associated primes of PS
4238	// algorithm: via characteristic sets
4239	// the given ordering of the variables is used
4240	//////////////////////////////////////////////////////
4241
4242
4243	static proc min_ass_prim_charsets0 (ideal PS)
4244	{
4245	intvec op;
4246	matrix m=char_series(PS); // We compute an irreducible
4247	// characteristic series
4248	int i,j,k;
4249	list PSI;
4250	list PHI; // the ideals given by the characteristic series
4251	for(i=nrows(m);i>=1; i--)
4252	{
4253	PHI[i]=ideal(m[i,1..ncols(m)]);
4254	}
4255	// We compute the radical of each ideal in PHI
4256	ideal I,JS,II;
4257	int sizeJS, sizeII;
4258	for(i=size(PHI);i>=1; i--)
4259	{
4260	I=0;
4261	for(j=size(PHI[i]);j>0;j=j-1)
4262	{
4263	I=I+ini_mod(PHI[i][j]);
4264	}
4265	JS=std(PHI[i]);
4266	sizeJS=size(JS);
4267	for(j=size(I);j>0;j=j-1)
4268	{
4269	II=0;
4270	sizeII=0;
4271	k=0;
4272	while(k<=sizeII) // successive saturation
4273	{
4274	op=option(get);
4275	option(returnSB);
4276	II=quotient(JS,I[j]);
4277	option(set,op);
4278	sizeII=size(II);
4279	if(sizeII==sizeJS)
4280	{
4281	for(k=1;k<=sizeII;k++)
4282	{
4283	if(leadexp(II[k])!=leadexp(JS[k])) break;
4284	}
4285	}
4286	JS=II;
4287	sizeJS=sizeII;
4288	}
4289	}
4290	PSI=insert(PSI,JS);
4291	}
4292	int sizePSI=size(PSI);
4293	// We eliminate redundant ideals
4294	for(i=1;i<sizePSI;i++)
4295	{
4296	for(j=i+1;j<=sizePSI;j++)
4297	{
4298	if(size(PSI[i])!=0)
4299	{
4300	if(size(PSI[j])!=0)
4301	{
4302	if(size(NF(PSI[i],PSI[j],1))==0)
4303	{
4304	PSI[j]=ideal(0);
4305	}
4306	else
4307	{
4308	if(size(NF(PSI[j],PSI[i],1))==0)
4309	{
4310	PSI[i]=ideal(0);
4311	}
4312	}
4313	}
4314	}
4315	}
4316	}
4317	for(i=sizePSI;i>=1;i--)
4318	{
4319	if(size(PSI[i])==0)
4320	{
4321	PSI=delete(PSI,i);
4322	}
4323	}
4324	return (PSI);
4325	}
4326
4327	///////////////////////////////////////////////////////
4328	// min_ass_prim_charsets1
4329	// input: generators of an ideal PS
4330	// output: the minimal associated primes of PS
4331	// algorithm: via characteristic sets
4332	// input: generators of an ideal PS and an integer i
4333	// The system tries to find an "optimal ordering" of
4334	// the variables
4335	//////////////////////////////////////////////////////
4336
4337
4338	static proc min_ass_prim_charsets1 (ideal PS)
4339	{
4340	intvec op;
4341	def oldring=basering;
4342	string n=system("neworder",PS);
4343	execute("ring r=("+charstr(oldring)+"),("+n+"),dp;");
4344	ideal PS=imap(oldring,PS);
4345	matrix m=char_series(PS); // We compute an irreducible
4346	// characteristic series
4347	int i,j,k;
4348	ideal I;
4349	list PSI;
4350	list PHI; // the ideals given by the characteristic series
4351	list ITPHI; // their initial terms
4352	for(i=nrows(m);i>=1; i--)
4353	{
4354	PHI[i]=ideal(m[i,1..ncols(m)]);
4355	I=0;
4356	for(j=size(PHI[i]);j>0;j=j-1)
4357	{
4358	I=I,ini_mod(PHI[i][j]);
4359	}
4360	I=I[2..ncols(I)];
4361	ITPHI[i]=I;
4362	}
4363	setring oldring;
4364	matrix m=imap(r,m);
4365	list PHI=imap(r,PHI);
4366	list ITPHI=imap(r,ITPHI);
4367	// We compute the radical of each ideal in PHI
4368	ideal I,JS,II;
4369	int sizeJS, sizeII;
4370	for(i=size(PHI);i>=1; i--)
4371	{
4372	I=0;
4373	for(j=size(PHI[i]);j>0;j=j-1)
4374	{
4375	I=I+ITPHI[i][j];
4376	}
4377	JS=std(PHI[i]);
4378	sizeJS=size(JS);
4379	for(j=size(I);j>0;j=j-1)
4380	{
4381	II=0;
4382	sizeII=0;
4383	k=0;
4384	while(k<=sizeII) // successive iteration
4385	{
4386	op=option(get);
4387	option(returnSB);
4388	II=quotient(JS,I[j]);
4389	option(set,op);
4390	//std
4391	// II=std(II);
4392	sizeII=size(II);
4393	if(sizeII==sizeJS)
4394	{
4395	for(k=1;k<=sizeII;k++)
4396	{
4397	if(leadexp(II[k])!=leadexp(JS[k])) break;
4398	}
4399	}
4400	JS=II;
4401	sizeJS=sizeII;
4402	}
4403	}
4404	PSI=insert(PSI,JS);
4405	}
4406	int sizePSI=size(PSI);
4407	// We eliminate redundant ideals
4408	for(i=1;i<sizePSI;i++)
4409	{
4410	for(j=i+1;j<=sizePSI;j++)
4411	{
4412	if(size(PSI[i])!=0)
4413	{
4414	if(size(PSI[j])!=0)
4415	{
4416	if(size(NF(PSI[i],PSI[j],1))==0)
4417	{
4418	PSI[j]=ideal(0);
4419	}
4420	else
4421	{
4422	if(size(NF(PSI[j],PSI[i],1))==0)
4423	{
4424	PSI[i]=ideal(0);
4425	}
4426	}
4427	}
4428	}
4429	}
4430	}
4431	for(i=sizePSI;i>=1;i--)
4432	{
4433	if(size(PSI[i])==0)
4434	{
4435	PSI=delete(PSI,i);
4436	}
4437	}
4438	return (PSI);
4439	}
4440
4441
4442	/////////////////////////////////////////////////////
4443	// proc prim_dec
4444	// input: generators of an ideal I and an integer choose
4445	// If choose=0, min_ass_prim_charsets with the given
4446	// ordering of the variables is used.
4447	// If choose=1, min_ass_prim_charsets with the "optimized"
4448	// ordering of the variables is used.
4449	// If choose=2, minAssPrimes from primdec.lib is used
4450	// If choose=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
4451	// output: a primary decomposition of I, i.e., a list
4452	// of pairs consisting of a standard basis of a primary component
4453	// of I and a standard basis of the corresponding associated prime.
4454	// To compute the minimal associated primes of a given ideal
4455	// min_ass_prim_l is called, i.e., the minimal associated primes
4456	// are computed via characteristic sets.
4457	// In the homogeneous case, the performance of the procedure
4458	// will be improved if I is already given by a minimal set of
4459	// generators. Apply minbase if necessary.
4460	//////////////////////////////////////////////////////////
4461
4462
4463	static proc prim_dec(ideal I, int choose)
4464	{
4465	if((choose<0) or (choose>3))
4466	{
4467	"ERROR: <int> must be 0 or 1 or 2 or 3";
4468	return();
4469	}
4470	if(system("version")>933)
4471	{
4472	option(notWarnSB);
4473	}
4474	ideal H=1; // The intersection of the primary components
4475	list U; // the leaves of the decomposition tree, i.e.,
4476	// pairs consisting of a primary component of I
4477	// and the corresponding associated prime
4478	list W; // the non-leaf vertices in the decomposition tree.
4479	// every entry has 6 components:
4480	// 1- the vertex itself , i.e., a standard bais of the
4481	// given ideal I (type 1), or a standard basis of a
4482	// pseudo-primary component arising from
4483	// pseudo-primary decomposition (type 2), or a
4484	// standard basis of a remaining component arising from
4485	// pseudo-primary decomposition or extraction (type 3)
4486	// 2- the type of the vertex as indicated above
4487	// 3- the weighted_tree_depth of the vertex
4488	// 4- the tester of the vertex
4489	// 5- a standard basis of the associated prime
4490	// of a vertex of type 2, or 0 otherwise
4491	// 6- a list of pairs consisting of a standard
4492	// basis of a minimal associated prime ideal
4493	// of the father of the vertex and the
4494	// irreducible factors of the "minimal
4495	// divisor" of the seperator or extractor
4496	// corresponding to the prime ideal
4497	// as computed by the procedure minsat,
4498	// if the vertex is of type 3, or
4499	// the empty list otherwise
4500	ideal SI=std(I);
4501	if(SI[1]==1) // primdecSY(ideal(1))
4502	{
4503	return(list());
4504	}
4505	int ncolsSI=ncols(SI);
4506	int ncolsH=1;
4507	W[1]=list(I,1,0,poly(1),ideal(0),list()); // The root of the tree
4508	int weighted_tree_depth;
4509	int i,j;
4510	int check;
4511	list V; // current vertex
4512	list VV; // new vertex
4513	list QQ;
4514	list WI;
4515	ideal Qi,SQ,SRest,fac;
4516	poly tester;
4517
4518	while(1)
4519	{
4520	i=1;
4521	while(1)
4522	{
4523	while(i<=size(W)) // find vertex V of smallest weighted tree-depth
4524	{
4525	if (W[i][3]<=weighted_tree_depth) break;
4526	i++;
4527	}
4528	if (i<=size(W)) break;
4529	i=1;
4530	weighted_tree_depth++;
4531	}
4532	V=W[i];
4533	W=delete(W,i); // delete V from W
4534
4535	// now proceed by type of vertex V
4536
4537	if (V[2]==2) // extraction needed
4538	{
4539	SQ,SRest,fac=extraction(V[1],V[5]);
4540	// standard basis of primary component,
4541	// standard basis of remaining component,
4542	// irreducible factors of
4543	// the "minimal divisor" of the extractor
4544	// as computed by the procedure minsat,
4545	check=0;
4546	for(j=1;j<=ncolsH;j++)
4547	{
4548	if (NF(H[j],SQ,1)!=0) // Q is not redundant
4549	{
4550	check=1;
4551	break;
4552	}
4553	}
4554	if(check==1) // Q is not redundant
4555	{
4556	QQ=list();
4557	QQ[1]=list(SQ,V[5]); // primary component, associated prime,
4558	// i.e., standard bases thereof
4559	U=U+QQ;
4560	H=intersect(H,SQ);
4561	H=std(H);
4562	ncolsH=ncols(H);
4563	check=0;
4564	if(ncolsH==ncolsSI)
4565	{
4566	for(j=1;j<=ncolsSI;j++)
4567	{
4568	if(leadexp(H[j])!=leadexp(SI[j]))
4569	{
4570	check=1;
4571	break;
4572	}
4573	}
4574	}
4575	else
4576	{
4577	check=1;
4578	}
4579	if(check==0) // H==I => U is a primary decomposition
4580	{
4581	return(U);
4582	}
4583	}
4584	if (SRest[1]!=1) // the remaining component is not
4585	// the whole ring
4586	{
4587	if (rad_con(V[4],SRest)==0) // the new vertex is not the
4588	// root of a redundant subtree
4589	{
4590	VV[1]=SRest; // remaining component
4591	VV[2]=3; // pseudoprimdec_special
4592	VV[3]=V[3]+1; // weighted depth
4593	VV[4]=V[4]; // the tester did not change
4594	VV[5]=ideal(0);
4595	VV[6]=list(list(V[5],fac));
4596	W=insert(W,VV,size(W));
4597	}
4598	}
4599	}
4600	else
4601	{
4602	if (V[2]==3) // pseudo_prim_dec_special is needed
4603	{
4604	QQ,SRest=pseudo_prim_dec_special_charsets(V[1],V[6],choose);
4605	// QQ = quadruples:
4606	// standard basis of pseudo-primary component,
4607	// standard basis of corresponding prime,
4608	// seperator, irreducible factors of
4609	// the "minimal divisor" of the seperator
4610	// as computed by the procedure minsat,
4611	// SRest=standard basis of remaining component
4612	}
4613	else // V is the root, pseudo_prim_dec is needed
4614	{
4615	QQ,SRest=pseudo_prim_dec_charsets(I,SI,choose);
4616	// QQ = quadruples:
4617	// standard basis of pseudo-primary component,
4618	// standard basis of corresponding prime,
4619	// seperator, irreducible factors of
4620	// the "minimal divisor" of the seperator
4621	// as computed by the procedure minsat,
4622	// SRest=standard basis of remaining component
4623
4624	}
4625	//check
4626	for(i=size(QQ);i>=1;i--)
4627	//for(i=1;i<=size(QQ);i++)
4628	{
4629	tester=QQ[i][3]*V[4];
4630	Qi=QQ[i][2];
4631	if(NF(tester,Qi,1)!=0) // the new vertex is not the
4632	// root of a redundant subtree
4633	{
4634	VV[1]=QQ[i][1];
4635	VV[2]=2;
4636	VV[3]=V[3]+1;
4637	VV[4]=tester; // the new tester as computed above
4638	VV[5]=Qi; // QQ[i][2];
4639	VV[6]=list();
4640	W=insert(W,VV,size(W));
4641	}
4642	}
4643	if (SRest[1]!=1) // the remaining component is not
4644	// the whole ring
4645	{
4646	if (rad_con(V[4],SRest)==0) // the vertex is not the root
4647	// of a redundant subtree
4648	{
4649	VV[1]=SRest;
4650	VV[2]=3;
4651	VV[3]=V[3]+2;
4652	VV[4]=V[4]; // the tester did not change
4653	VV[5]=ideal(0);
4654	WI=list();
4655	for(i=1;i<=size(QQ);i++)
4656	{
4657	WI=insert(WI,list(QQ[i][2],QQ[i][4]));
4658	}
4659	VV[6]=WI;
4660	W=insert(W,VV,size(W));
4661	}
4662	}
4663	}
4664	}
4665	}
4666
4667	//////////////////////////////////////////////////////////////////////////
4668	// proc pseudo_prim_dec_charsets
4669	// input: Generators of an arbitrary ideal I, a standard basis SI of I,
4670	// and an integer choo
4671	// If choo=0, min_ass_prim_charsets with the given
4672	// ordering of the variables is used.
4673	// If choo=1, min_ass_prim_charsets with the "optimized"
4674	// ordering of the variables is used.
4675	// If choo=2, minAssPrimes from primdec.lib is used
4676	// If choo=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
4677	// output: a pseudo primary decomposition of I, i.e., a list
4678	// of pseudo primary components together with a standard basis of the
4679	// remaining component. Each pseudo primary component is
4680	// represented by a quadrupel: A standard basis of the component,
4681	// a standard basis of the corresponding associated prime, the
4682	// seperator of the component, and the irreducible factors of the
4683	// "minimal divisor" of the seperator as computed by the procedure minsat,
4684	// calls proc pseudo_prim_dec_i
4685	//////////////////////////////////////////////////////////////////////////
4686
4687
4688	static proc pseudo_prim_dec_charsets (ideal I, ideal SI, int choo)
4689	{
4690	list L; // The list of minimal associated primes,
4691	// each one given by a standard basis
4692	if((choo==0) or (choo==1))
4693	{
4694	L=min_ass_prim_charsets(I,choo);
4695	}
4696	else
4697	{
4698	if(choo==2)
4699	{
4700	L=minAssPrimes(I);
4701	}
4702	else
4703	{
4704	L=minAssPrimes(I,1);
4705	}
4706	for(int i=size(L);i>=1;i=i-1)
4707	{
4708	L[i]=std(L[i]);
4709	}
4710	}
4711	return (pseudo_prim_dec_i(SI,L));
4712	}
4713
4714	////////////////////////////////////////////////////////////////
4715	// proc pseudo_prim_dec_special_charsets
4716	// input: a standard basis of an ideal I whose radical is the
4717	// intersection of the radicals of ideals generated by one prime ideal
4718	// P_i together with one polynomial f_i, the list V6 must be the list of
4719	// pairs (standard basis of P_i, irreducible factors of f_i),
4720	// and an integer choo
4721	// If choo=0, min_ass_prim_charsets with the given
4722	// ordering of the variables is used.
4723	// If choo=1, min_ass_prim_charsets with the "optimized"
4724	// ordering of the variables is used.
4725	// If choo=2, minAssPrimes from primdec.lib is used
4726	// If choo=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
4727	// output: a pseudo primary decomposition of I, i.e., a list
4728	// of pseudo primary components together with a standard basis of the
4729	// remaining component. Each pseudo primary component is
4730	// represented by a quadrupel: A standard basis of the component,
4731	// a standard basis of the corresponding associated prime, the
4732	// seperator of the component, and the irreducible factors of the
4733	// "minimal divisor" of the seperator as computed by the procedure minsat,
4734	// calls proc pseudo_prim_dec_i
4735	////////////////////////////////////////////////////////////////
4736
4737
4738	static proc pseudo_prim_dec_special_charsets (ideal SI,list V6, int choo)
4739	{
4740	int i,j,l;
4741	list m;
4742	list L;
4743	int sizeL;
4744	ideal P,SP; ideal fac;
4745	int dimSP;
4746	for(l=size(V6);l>=1;l--) // creates a list of associated primes
4747	// of I, possibly redundant
4748	{
4749	P=V6[l][1];
4750	fac=V6[l][2];
4751	for(i=ncols(fac);i>=1;i--)
4752	{
4753	SP=P+fac[i];
4754	SP=std(SP);
4755	if(SP[1]!=1)
4756	{
4757	if((choo==0) or (choo==1))
4758	{
4759	m=min_ass_prim_charsets(SP,choo); // a list of SB
4760	}
4761	else
4762	{
4763	if(choo==2)
4764	{
4765	m=minAssPrimes(SP);
4766	}
4767	else
4768	{
4769	m=minAssPrimes(SP,1);
4770	}
4771	for(j=size(m);j>=1;j=j-1)
4772	{
4773	m[j]=std(m[j]);
4774	}
4775	}
4776	dimSP=dim(SP);
4777	for(j=size(m);j>=1; j--)
4778	{
4779	if(dim(m[j])==dimSP)
4780	{
4781	L=insert(L,m[j],size(L));
4782	}
4783	}
4784	}
4785	}
4786	}
4787	sizeL=size(L);
4788	for(i=1;i<sizeL;i++) // get rid of redundant primes
4789	{
4790	for(j=i+1;j<=sizeL;j++)
4791	{
4792	if(size(L[i])!=0)
4793	{
4794	if(size(L[j])!=0)
4795	{
4796	if(size(NF(L[i],L[j],1))==0)
4797	{
4798	L[j]=ideal(0);
4799	}
4800	else
4801	{
4802	if(size(NF(L[j],L[i],1))==0)
4803	{
4804	L[i]=ideal(0);
4805	}
4806	}
4807	}
4808	}
4809	}
4810	}
4811	for(i=sizeL;i>=1;i--)
4812	{
4813	if(size(L[i])==0)
4814	{
4815	L=delete(L,i);
4816	}
4817	}
4818	return (pseudo_prim_dec_i(SI,L));
4819	}
4820
4821
4822	////////////////////////////////////////////////////////////////
4823	// proc pseudo_prim_dec_i
4824	// input: A standard basis of an arbitrary ideal I, and standard bases
4825	// of the minimal associated primes of I
4826	// output: a pseudo primary decomposition of I, i.e., a list
4827	// of pseudo primary components together with a standard basis of the
4828	// remaining component. Each pseudo primary component is
4829	// represented by a quadrupel: A standard basis of the component Q_i,
4830	// a standard basis of the corresponding associated prime P_i, the
4831	// seperator of the component, and the irreducible factors of the
4832	// "minimal divisor" of the seperator as computed by the procedure minsat,
4833	////////////////////////////////////////////////////////////////
4834
4835
4836	static proc pseudo_prim_dec_i (ideal SI, list L)
4837	{
4838	list Q;
4839	if (size(L)==1) // one minimal associated prime only
4840	// the ideal is already pseudo primary
4841	{
4842	Q=SI,L[1],1;
4843	list QQ;
4844	QQ[1]=Q;
4845	return (QQ,ideal(1));
4846	}
4847
4848	poly f0,f,g;
4849	ideal fac;
4850	int i,j,k,l;
4851	ideal SQi;
4852	ideal I'=SI;
4853	list QP;
4854	int sizeL=size(L);
4855	for(i=1;i<=sizeL;i++)
4856	{
4857	fac=0;
4858	for(j=1;j<=sizeL;j++) // compute the seperator sep_i
4859	// of the i-th component
4860	{
4861	if (i!=j) // search g not in L[i], but L[j]
4862	{
4863	for(k=1;k<=ncols(L[j]);k++)
4864	{
4865	if(NF(L[j][k],L[i],1)!=0)
4866	{
4867	break;
4868	}
4869	}
4870	fac=fac+L[j][k];
4871	}
4872	}
4873	// delete superfluous polynomials
4874	fac=simplify(fac,8);
4875	// saturation
4876	SQi,f0,f,fac=minsat_ppd(SI,fac);
4877	I'=I',f;
4878	QP=SQi,L[i],f0,fac;
4879	// the quadrupel:
4880	// a standard basis of Q_i,
4881	// a standard basis of P_i,
4882	// sep_i,
4883	// irreducible factors of
4884	// the "minimal divisor" of the seperator
4885	// as computed by the procedure minsat,
4886	Q[i]=QP;
4887	}
4888	I'=std(I');
4889	return (Q, I');
4890	// I' = remaining component
4891	}
4892
4893
4894	////////////////////////////////////////////////////////////////
4895	// proc extraction
4896	// input: A standard basis of a pseudo primary ideal I, and a standard
4897	// basis of the unique minimal associated prime P of I
4898	// output: an extraction of I, i.e., a standard basis of the primary
4899	// component Q of I with associated prime P, a standard basis of the
4900	// remaining component, and the irreducible factors of the
4901	// "minimal divisor" of the extractor as computed by the procedure minsat
4902	////////////////////////////////////////////////////////////////
4903
4904
4905	static proc extraction (ideal SI, ideal SP)
4906	{
4907	list indsets=indepSet(SP,0);
4908	poly f;
4909	if(size(indsets)!=0) //check, whether dim P != 0
4910	{
4911	intvec v; // a maximal independent set of variables
4912	// modulo P
4913	string U; // the independent variables
4914	string A; // the dependent variables
4915	int j,k;
4916	int a; // the size of A
4917	int degf;
4918	ideal g;
4919	list polys;
4920	int sizepolys;
4921	list newpoly;
4922	def R=basering;
4923	//intvec hv=hilb(SI,1);
4924	for (k=1;k<=size(indsets);k++)
4925	{
4926	v=indsets[k];
4927	for (j=1;j<=nvars(R);j++)
4928	{
4929	if (v[j]==1)
4930	{
4931	U=U+varstr(j)+",";
4932	}
4933	else
4934	{
4935	A=A+varstr(j)+",";
4936	a++;
4937	}
4938	}
4939
4940	U[size(U)]=")"; // we compute the extractor of I (w.r.t. U)
4941	execute("ring RAU=("+charstr(basering)+"),("+A+U+",(dp("+string(a)+"),dp);");
4942	ideal I=imap(R,SI);
4943	//I=std(I,hv); // the standard basis in (R[U])[A]
4944	I=std(I); // the standard basis in (R[U])[A]
4945	A[size(A)]=")";
4946	execute("ring Rloc=("+charstr(basering)+","+U+",("+A+",dp;");
4947	ideal I=imap(RAU,I);
4948	//"std in lokalisierung:"+newline,I;
4949	ideal h;
4950	for(j=ncols(I);j>=1;j--)
4951	{
4952	h[j]=leadcoef(I[j]); // consider I in (R(U))[A]
4953	}
4954	setring R;
4955	g=imap(Rloc,h);
4956	kill RAU,Rloc;
4957	U="";
4958	A="";
4959	a=0;
4960	f=lcm(g);
4961	newpoly[1]=f;
4962	polys=polys+newpoly;
4963	newpoly=list();
4964	}
4965	f=polys[1];
4966	degf=deg(f);
4967	sizepolys=size(polys);
4968	for (k=2;k<=sizepolys;k++)
4969	{
4970	if (deg(polys[k])<degf)
4971	{
4972	f=polys[k];
4973	degf=deg(f);
4974	}
4975	}
4976	}
4977	else
4978	{
4979	f=1;
4980	}
4981	poly f0,h0; ideal SQ; ideal fac;
4982	if(f!=1)
4983	{
4984	SQ,f0,h0,fac=minsat(SI,f);
4985	return(SQ,std(SI+h0),fac);
4986	// the tripel
4987	// a standard basis of Q,
4988	// a standard basis of remaining component,
4989	// irreducible factors of
4990	// the "minimal divisor" of the extractor
4991	// as computed by the procedure minsat
4992	}
4993	else
4994	{
4995	return(SI,ideal(1),ideal(1));
4996	}
4997	}
4998
4999	/////////////////////////////////////////////////////
5000	// proc minsat
5001	// input: a standard basis of an ideal I and a polynomial p
5002	// output: a standard basis IS of the saturation of I w.r. to p,
5003	// the maximal squarefree factor f0 of p,
5004	// the "minimal divisor" f of f0 such that the saturation of
5005	// I w.r. to f equals the saturation of I w.r. to f0 (which is IS),
5006	// the irreducible factors of f
5007	//////////////////////////////////////////////////////////
5008
5009
5010	static proc minsat(ideal SI, poly p)
5011	{
5012	ideal fac=factorize(p,1); //the irreducible factors of p
5013	fac=sort(fac)[1];
5014	int i,k;
5015	poly f0=1;
5016	for(i=ncols(fac);i>=1;i--)
5017	{
5018	f0=f0*fac[i];
5019	}
5020	poly f=1;
5021	ideal iold;
5022	list quotM;
5023	quotM[1]=SI;
5024	quotM[2]=fac;
5025	quotM[3]=f0;
5026	// we deal seperately with the first quotient;
5027	// factors, which do not contribute to this one,
5028	// are omitted
5029	iold=quotM[1];
5030	quotM=minquot(quotM);
5031	fac=quotM[2];
5032	if(quotM[3]==1)
5033	{
5034	return(quotM[1],f0,f,fac);
5035	}
5036	while(special_ideals_equal(iold,quotM[1])==0)
5037	{
5038	f=f*quotM[3];
5039	iold=quotM[1];
5040	quotM=minquot(quotM);
5041	}
5042	return(quotM[1],f0,f,fac); // the quadrupel ((I:p),f0,f, irr. factors of f)
5043	}
5044
5045	/////////////////////////////////////////////////////
5046	// proc minsat_ppd
5047	// input: a standard basis of an ideal I and a polynomial p
5048	// output: a standard basis IS of the saturation of I w.r. to p,
5049	// the maximal squarefree factor f0 of p,
5050	// the "minimal divisor" f of f0 such that the saturation of
5051	// I w.r. to f equals the saturation of I w.r. to f0 (which is IS),
5052	// the irreducible factors of f
5053	//////////////////////////////////////////////////////////
5054
5055
5056	static proc minsat_ppd(ideal SI, ideal fac)
5057	{
5058	fac=sort(fac)[1];
5059	int i,k;
5060	poly f0=1;
5061	for(i=ncols(fac);i>=1;i--)
5062	{
5063	f0=f0*fac[i];
5064	}
5065	poly f=1;
5066	ideal iold;
5067	list quotM;
5068	quotM[1]=SI;
5069	quotM[2]=fac;
5070	quotM[3]=f0;
5071	// we deal seperately with the first quotient;
5072	// factors, which do not contribute to this one,
5073	// are omitted
5074	iold=quotM[1];
5075	quotM=minquot(quotM);
5076	fac=quotM[2];
5077	if(quotM[3]==1)
5078	{
5079	return(quotM[1],f0,f,fac);
5080	}
5081	while(special_ideals_equal(iold,quotM[1])==0)
5082	{
5083	f=f*quotM[3];
5084	iold=quotM[1];
5085	quotM=minquot(quotM);
5086	k++;
5087	}
5088	return(quotM[1],f0,f,fac); // the quadrupel ((I:p),f0,f, irr. factors of f)
5089	}
5090	/////////////////////////////////////////////////////////////////
5091	// proc minquot
5092	// input: a list with 3 components: a standard basis
5093	// of an ideal I, a set of irreducible polynomials, and
5094	// there product f0
5095	// output: a standard basis of the ideal (I:f0), the irreducible
5096	// factors of the "minimal divisor" f of f0 with (I:f0) = (I:f),
5097	// the "minimal divisor" f
5098	/////////////////////////////////////////////////////////////////
5099
5100	static proc minquot(list tsil)
5101	{
5102	intvec op;
5103	int i,j,k,action;
5104	ideal verg;
5105	list l;
5106	poly g;
5107	ideal laedi=tsil[1];
5108	ideal fac=tsil[2];
5109	poly f=tsil[3];
5110
5111	//std
5112	// ideal star=quotient(laedi,f);
5113	// star=std(star);
5114	op=option(get);
5115	option(returnSB);
5116	ideal star=quotient(laedi,f);
5117	option(set,op);
5118	if(special_ideals_equal(laedi,star)==1)
5119	{
5120	return(laedi,ideal(1),1);
5121	}
5122	action=1;
5123	while(action==1)
5124	{
5125	if(size(fac)==1)
5126	{
5127	action=0;
5128	break;
5129	}
5130	for(i=1;i<=size(fac);i++)
5131	{
5132	g=1;
5133	for(j=1;j<=size(fac);j++)
5134	{
5135	if(i!=j)
5136	{
5137	g=g*fac[j];
5138	}
5139	}
5140	//std
5141	// verg=quotient(laedi,g);
5142	// verg=std(verg);
5143	op=option(get);
5144	option(returnSB);
5145	verg=quotient(laedi,g);
5146	option(set,op);
5147	if(special_ideals_equal(verg,star)==1)
5148	{
5149	f=g;
5150	fac[i]=0;
5151	fac=simplify(fac,2);
5152	break;
5153	}
5154	if(i==size(fac))
5155	{
5156	action=0;
5157	}
5158	}
5159	}
5160	l=star,fac,f;
5161	return(l);
5162	}
5163	/////////////////////////////////////////////////
5164	// proc special_ideals_equal
5165	// input: standard bases of ideal k1 and k2 such that
5166	// k1 is contained in k2, or k2 is contained ink1
5167	// output: 1, if k1 equals k2, 0 otherwise
5168	//////////////////////////////////////////////////
5169
5170	static proc special_ideals_equal( ideal k1, ideal k2)
5171	{
5172	int j;
5173	if(size(k1)==size(k2))
5174	{
5175	for(j=1;j<=size(k1);j++)
5176	{
5177	if(leadexp(k1[j])!=leadexp(k2[j]))
5178	{
5179	return(0);
5180	}
5181	}
5182	return(1);
5183	}
5184	return(0);
5185	}
5186
5187
5188	///////////////////////////////////////////////////////////////////////////////
5189
5190	static proc convList(list l)
5191	{
5192	int i;
5193	list re,he;
5194	for(i=1;i<=size(l)/2;i++)
5195	{
5196	he=l[2i-1],l[2i];
5197	re[i]=he;
5198	}
5199	return(re);
5200	}
5201	///////////////////////////////////////////////////////////////////////////////
5202
5203	static proc reconvList(list l)
5204	{
5205	int i;
5206	list re;
5207	for(i=1;i<=size(l);i++)
5208	{
5209	re[2*i-1]=l[i][1];
5210	re[2*i]=l[i][2];
5211	}
5212	return(re);
5213	}
5214
5215	///////////////////////////////////////////////////////////////////////////////
5216	//
5217	// The main procedures
5218	//
5219	///////////////////////////////////////////////////////////////////////////////
5220
5221	proc primdecGTZ(ideal i)
5222	"USAGE: primdecGTZ(i); i ideal
5223	RETURN: a list pr of primary ideals and their associated primes:
5224	@format
5225	pr[i][1] the i-th primary component,
5226	pr[i][2] the i-th prime component.
5227	@end format
5228	NOTE: Algorithm of Gianni/Trager/Zacharias.
5229	Designed for characteristic 0, works also in char k > 0, if it
5230	terminates (may result in an infinite loop in small characteristic!)
5231	EXAMPLE: example primdecGTZ; shows an example
5232	"
5233	{
5234	if(ord_test(basering)!=1)
5235	{
5236	ERROR(
5237	"// Not implemented for this ordering, please change to global ordering."
5238	);
5239	}
5240	if(minpoly!=0)
5241	{
5242	return(algeDeco(i,0));
5243	ERROR(
5244	"// Not implemented yet for algebraic extensions.Simulate the ring extension by adding the minpoly to the ideal"
5245	);
5246	}
5247	return(convList(decomp(i)));
5248	}
5249	example
5250	{ "EXAMPLE:"; echo = 2;
5251	ring r = 0,(x,y,z),lp;
5252	poly p = z2+1;
5253	poly q = z3+2;
5254	ideal i = p*q^2,y-z2;
5255	list pr = primdecGTZ(i);
5256	pr;
5257	}
5258	///////////////////////////////////////////////////////////////////////////////
5259
5260	proc absPrimdecGTZ(ideal I)
5261	"USAGE: absPrimdecGTZ(I); I ideal
5262	ASSUME: Ground field has characteristic 0.
5263	RETURN: a ring containing two lists: @code{absolute_primes} (the absolute
5264	prime components of I) and @code{primary_decomp} (the output of
5265	@code{primdecGTZ(I)}).
5266	The list absolute_primes has to be interpreted as follows:
5267	each entry describes a class of conjugated absolute primes,
5268	@format
5269	absolute_primes[i][1] the absolute prime component,
5270	absolute_primes[i][2] the number of conjugates.
5271	@end format
5272	The first entry of @code{absolute_primes[i][1]} is the minimal
5273	polynomial of a minimal finite field extension over which the
5274	absolute prime component is defined.
5275	NOTE: Algorithm of Gianni/Trager/Zacharias combined with the
5276	@code{absFactorize} command.
5277	SEE ALSO: primdecGTZ; absFactorize
5278	EXAMPLE: example absPrimdecGTZ; shows an example
5279	"
5280	{
5281	if (char(basering) != 0) {
5282	ERROR("primdec.lib::absPrimdecGTZ is only implemented for "+
5283	+"characteristic 0");
5284	}
5285
5286	if(ord_test(basering)!=1)
5287	{
5288	ERROR(
5289	"// Not implemented for this ordering, please change to global ordering."
5290	);
5291	}
5292	if(minpoly!=0)
5293	{
5294	//return(algeDeco(i,0));
5295	ERROR(
5296	"// Not implemented yet for algebraic extensions.Simulate the ring extension by adding the minpoly to the ideal"
5297	);
5298	}
5299	def R=basering;
5300	int n=nvars(R);
5301	list L=decomp(I,3);
5302	string newvar=L[1][3];
5303	int k=find(newvar,",",find(newvar,",")+1);
5304	newvar=newvar[k+1..size(newvar)];
5305	list lR=ringlist(R);
5306	int i,d;
5307	intvec vv;
5308	for(i=1;i<=n;i++){vv[i]=1;}
5309
5310	list orst;
5311	orst[1]=list("dp",vv);
5312	orst[2]=list("dp",intvec(1));
5313	orst[3]=list("C",0);
5314	lR[3]=orst;
5315	lR[2][n+1] = newvar;
5316	def Rz = ring(lR);
5317	setring Rz;
5318	list L=imap(R,L);
5319	list absolute_primes,primary_decomp;
5320	ideal I,M,N,K;
5321	M=maxideal(1);
5322	N=maxideal(1);
5323	poly p,q,f,g;
5324	map phi,psi;
5325	for(i=1;i<=size(L);i++)
5326	{
5327	I=L[i][2];
5328	execute("K="+L[i][3]+";");
5329	p=K[1];
5330	q=K[2];
5331	execute("f="+L[i][4]+";");
5332	g=2*var(n)-f;
5333	M[n]=f;
5334	N[n]=g;
5335	d=deg(p);
5336	phi=Rz,M;
5337	psi=Rz,N;
5338	I=phi(I),p,q;
5339	I=std(I);
5340	absolute_primes[i]=list(psi(I),d);
5341	primary_decomp[i]=list(L[i][1],L[i][2]);
5342	}
5343	export(primary_decomp);
5344	export(absolute_primes);
5345	setring R;
5346	dbprint( printlevel-voice+3,"
5347	// 'absPrimdecGTZ' created a ring, in which two lists absolute_primes (the
5348	// absolute prime components) and primary_decomp (the primary and prime
5349	// components over the current basering) are stored.
5350	// To access the list of absolute prime components, type (if the name S was
5351	// assigned to the return value):
5352	setring S; absolute_primes; ");
5353
5354	return(Rz);
5355	}
5356	example
5357	{ "EXAMPLE:"; echo = 2;
5358	ring r = 0,(x,y,z),lp;
5359	poly p = z2+1;
5360	poly q = z3+2;
5361	ideal i = p*q^2,y-z2;
5362	def S = absPrimdecGTZ(i);
5363	setring S;
5364	absolute_primes;
5365	}
5366	///////////////////////////////////////////////////////////////////////////////
5367
5368	proc primdecSY(ideal i, list #)
5369	"USAGE: primdecSY(i); i ideal, c int
5370	RETURN: a list pr of primary ideals and their associated primes:
5371	@format
5372	pr[i][1] the i-th primary component,
5373	pr[i][2] the i-th prime component.
5374	@end format
5375	NOTE: Algorithm of Shimoyama/Yokoyama.
5376	@format
5377	if c=0, the given ordering of the variables is used,
5378	if c=1, minAssChar tries to use an optimal ordering,
5379	if c=2, minAssGTZ is used,
5380	if c=3, minAssGTZ and facstd are used.
5381	@end format
5382	EXAMPLE: example primdecSY; shows an example
5383	"
5384	{
5385	if(ord_test(basering)!=1)
5386	{
5387	ERROR(
5388	"// Not implemented for this ordering, please change to global ordering."
5389	);
5390	}
5391	i=simplify(i,2);
5392	if ((i[1]==0)\|\|(i[1]==1))
5393	{
5394	list L=list(ideal(i[1]),ideal(i[1]));
5395	return(list(L));
5396	}
5397	if(minpoly!=0)
5398	{
5399	return(algeDeco(i,1));
5400	}
5401	if (size(#)==1)
5402	{ return(prim_dec(i,#[1])); }
5403	else
5404	{ return(prim_dec(i,1)); }
5405	}
5406	example
5407	{ "EXAMPLE:"; echo = 2;
5408	ring r = 0,(x,y,z),lp;
5409	poly p = z2+1;
5410	poly q = z3+2;
5411	ideal i = p*q^2,y-z2;
5412	list pr = primdecSY(i);
5413	pr;
5414	}
5415	///////////////////////////////////////////////////////////////////////////////
5416	proc minAssGTZ(ideal i,list #)
5417	"USAGE: minAssGTZ(i); i ideal
5418	minAssGTZ(i,1); i ideal does not use the factorizing Groebner
5419	RETURN: a list, the minimal associated prime ideals of i.
5420	NOTE: Designed for characteristic 0, works also in char k > 0 based
5421	on an algorithm of Yokoyama
5422	EXAMPLE: example minAssGTZ; shows an example
5423	"
5424	{
5425	if(ord_test(basering)!=1)
5426	{
5427	ERROR(
5428	"// Not implemented for this ordering, please change to global ordering."
5429	);
5430	}
5431	if(minpoly!=0)
5432	{
5433	return(algeDeco(i,2));
5434	}
5435	if(size(#)==0)
5436	{
5437	return(minAssPrimes(i,1));
5438	}
5439	return(minAssPrimes(i));
5440	}
5441	example
5442	{ "EXAMPLE:"; echo = 2;
5443	ring r = 0,(x,y,z),dp;
5444	poly p = z2+1;
5445	poly q = z3+2;
5446	ideal i = p*q^2,y-z2;
5447	list pr = minAssGTZ(i);
5448	pr;
5449	}
5450
5451	///////////////////////////////////////////////////////////////////////////////
5452	proc minAssChar(ideal i, list #)
5453	"USAGE: minAssChar(i[,c]); i ideal, c int.
5454	RETURN: list, the minimal associated prime ideals of i.
5455	NOTE: If c=0, the given ordering of the variables is used. @*
5456	Otherwise, the system tries to find an optimal ordering,
5457	which in some cases may considerably speed up the algorithm. @*
5458	Due to a bug in the factorization, the result may be not completely
5459	decomposed in small characteristic.
5460	EXAMPLE: example minAssChar; shows an example
5461	"
5462	{
5463	if(ord_test(basering)!=1)
5464	{
5465	ERROR(
5466	"// Not implemented for this ordering, please change to global ordering."
5467	);
5468	}
5469	if (size(#)==1)
5470	{ return(min_ass_prim_charsets(i,#[1])); }
5471	else
5472	{ return(min_ass_prim_charsets(i,1)); }
5473	}
5474	example
5475	{ "EXAMPLE:"; echo = 2;
5476	ring r = 0,(x,y,z),dp;
5477	poly p = z2+1;
5478	poly q = z3+2;
5479	ideal i = p*q^2,y-z2;
5480	list pr = minAssChar(i);
5481	pr;
5482	}
5483	///////////////////////////////////////////////////////////////////////////////
5484	proc equiRadical(ideal i)
5485	"USAGE: equiRadical(i); i ideal
5486	RETURN: ideal, intersection of associated primes of i of maximal dimension.
5487	NOTE: A combination of the algorithms of Krick/Logar and Kemper is used.
5488	Works also in positive characteristic (Kempers algorithm).
5489	EXAMPLE: example equiRadical; shows an example
5490	"
5491	{
5492	if(ord_test(basering)!=1)
5493	{
5494	ERROR(
5495	"// Not implemented for this ordering, please change to global ordering."
5496	);
5497	}
5498	return(radical(i,1));
5499	}
5500	example
5501	{ "EXAMPLE:"; echo = 2;
5502	ring r = 0,(x,y,z),dp;
5503	poly p = z2+1;
5504	poly q = z3+2;
5505	ideal i = p*q^2,y-z2;
5506	ideal pr= equiRadical(i);
5507	pr;
5508	}
5509
5510	///////////////////////////////////////////////////////////////////////////////
5511	proc radical(ideal i, list #)
5512	"USAGE: radical(i); i ideal.
5513	Optional parameters in list #: (can be entered in any order)
5514	1, "equiRad" -> equiRadical is computed
5515	0, "fullRad" -> full radical is computed (default)
5516	"SL" -> the new algorithm is used (default)
5517	"KL" -> the old algorithm is used
5518	"KLdp" -> the old algorithm, with modifications
5519	"facstd" -> uses facstd to first decompose the ideal (default for non homogeneous ideals)
5520	"noFacstd" -> does not use facstd (default for homogeneous ideals)
5521	RETURN: ideal, the radical of i.
5522	NOTE: A combination of the algorithms of Krick/Logar and Kemper is used.
5523	Works also in positive characteristic (Kempers algorithm).
5524	EXAMPLE: example radical; shows an example
5525	"
5526	{
5527	dbprint(printlevel - voice, "Radical, version 2006.05.08");
5528	if(ord_test(basering)!=1)
5529	{
5530	ERROR(
5531	"// Not implemented for this ordering, please change to global ordering."
5532	);
5533	}
5534	def @P=basering;
5535	int j;
5536	int il;
5537	string algorithm;
5538	int useFac;
5539
5540	// Set input parameters
5541	algorithm = "SL"; // Default: SL algorithm
5542	il = 0; // Default: Full radical (not only equidim part)
5543	if (homog(i) == 1) { // Default: facStd is used, except if the ideal is homogeneous.
5544	useFac = 0;
5545	} else {
5546	useFac = 1;
5547	};
5548	if(size(#) > 0){
5549	int valid;
5550	for(j = 1; j <= size(#); j++){
5551	valid = 0;
5552	if((typeof(#[j]) == "int") or (typeof(#[j]) == "number")) {
5553	il = #[j]; // If il == 1, equiRadical is computed
5554	valid = 1;
5555	};
5556	if(typeof(#[j]) == "string"){
5557	if(#[j] == "KL") {
5558	algorithm = "KL";
5559	valid = 1};
5560	if(#[j] == "KLdp") {
5561	algorithm = "KLdp";
5562	valid = 1};
5563	if(#[j] == "SL") {
5564	algorithm = "SL";
5565	valid = 1};
5566	if(#[j] == "noFacstd") {
5567	useFac = 0;
5568	valid = 1};
5569	if(#[j] == "facstd") {
5570	useFac = 1;
5571	valid = 1};
5572	if(#[j] == "equiRad") {
5573	il = 1;
5574	valid = 1};
5575	if(#[j] == "fullRad") {
5576	il = 0;
5577	valid = 1};
5578	};
5579	if(valid == 0) {
5580	dbprint(1, "Warning! The following input parameter was not recognized:", #[j]);
5581	};
5582	};
5583	};
5584
5585	if(size(i) == 0){return(ideal(0));}
5586	ideal rad = 1;
5587	intvec op = option(get);
5588	list qr = simplifyIdeal(i);
5589	// SL - Line removed. The ideal isave was never used.
5590	//ideal isave=i;
5591	map phi = @P, qr[2];
5592
5593	option(redSB);
5594	i = groebner(qr[1]);
5595	option(set, op);
5596	int di = dim(i);
5597
5598	if(di == 0)
5599	{
5600	i = zeroRad(i, qr[1]);
5601	return(interred(phi(i)));
5602	}
5603
5604	option(redSB);
5605	list pr;
5606	if(useFac == 1)
5607	{
5608	pr = facstd(i);
5609	} else {
5610	pr = i
5611	};
5612	option(set, op);
5613	int s = size(pr);
5614	if(useFac == 1) {
5615	dbprint(printlevel - voice, "Number of components returned by facstd: ", s);
5616	};
5617	for(j = 1; j <= s; j++)
5618	{
5619	attrib(pr[s + 1 - j], "isSB", 1);
5620	if((size(reduce(rad, pr[s + 1 - j], 1)) != 0) && ((dim(pr[s + 1 - j]) == di) \|\| !il))
5621	{
5622	// SL Debug messages
5623	dbprint(printlevel-voice, "We shall compute the radical of ", pr[s + 1 - j]);
5624	dbprint(printlevel-voice, "The dimension is: ", dim(pr[s+1-j]));
5625
5626	if(algorithm == "KL") {
5627	rad = intersect(rad, radicalKL(pr[s + 1 - j], rad, il));
5628	}
5629	if(algorithm == "KLdp") {
5630	rad = intersect(rad, newRadicalKL(pr[s + 1 - j], rad, il));
5631	}
5632	if(algorithm == "SL") {
5633	rad = intersect(rad, newRadicalSL(pr[s + 1 - j], il));
5634	};
5635	} else
5636	{
5637	// SL Debug
5638	dbprint(printlevel-voice, "The radical of this component is not needed.");
5639	dbprint(printlevel-voice, "size(reduce(rad, pr[s + 1 - j], 1))", size(reduce(rad, pr[s + 1 - j], 1)));
5640	dbprint(printlevel-voice, "dim(pr[s + 1 - j])", dim(pr[s + 1 - j]));
5641	dbprint(printlevel-voice, "il", il);
5642	};
5643	}
5644	return(interred(phi(rad)));
5645	}
5646	example
5647	{ "EXAMPLE:"; echo = 2;
5648	ring r = 0,(x,y,z),dp;
5649	poly p = z2+1;
5650	poly q = z3+2;
5651	ideal i = p*q^2,y-z2;
5652	ideal pr = radical(i);
5653	pr;
5654	}
5655
5656	///////////////////////////////////////////////////////////////////////////////
5657	//
5658	// Computes the radical of I using KL algorithm.
5659	// The only difference with the the previous implementation of KL algorithm is
5660	// that now it uses block dp instead of lp ordering for the reduction to the
5661	// zerodimensional case.
5662	// The reduction step has been moved to the new routine newRadicalReduction, so that it can be
5663	// used also by the new algorithm.
5664	// This algorithm should be left in the final version of the library,
5665	// and radicalKL removed.
5666	//
5667	static proc newRadicalKL(ideal I, ideal ser, list #) {
5668	// I The ideal for which the radical is computed
5669	// ideal ser Same as in radicalKL (used to reduce components already obtained)
5670	// list # If #[1] = 1, equiradical is computed.
5671
5672	// I needs to be a Groebner basis.
5673	if (attrib(I, "isSB") != 1) {
5674	I = groebner(I);
5675	};
5676
5677	ideal rad; // The radical
5678	int allIndep = 1; // All max independent sets are used
5679
5680	list result = newRadicalReduction(I, ser, allIndep, #);
5681	int done = result[3];
5682	rad = result[1];
5683	if (done == 0) {
5684	rad = intersect(rad, newRadicalKL(result[2], ideal(1), #));
5685	};
5686	return(rad);
5687	};
5688
5689
5690	///////////////////////////////////////////////////////////////////////////////
5691	//
5692	// Computes the radical of I via the new algorithm, using zerodimensional radical in
5693	// the zero dimensional case.
5694	// For the reduction to the zerodimensional case, it uses the procedure
5695	// radical of primdec.lib, with some modifications to avoid the recursion.
5696	//
5697	static proc newRadicalSL(ideal I, list #)
5698	// Input = I, ideal
5699	// #, list. If #[1] = 1, then computes only the equiradical.
5700	// Output = (P, primaryDec) where P = rad(I) and primaryDec is the list of the radicals
5701	// obtained in intermediate steps.
5702	{
5703	ideal rad = 1;
5704	ideal equiRad = 1;
5705	list primes;
5706	int k; // Counter
5707	int il; // If il = 1, only the equiradical is required.
5708	int iDim; // The dimension of I
5709	int stop = 0; // Checks if the radical has been obtained
5710
5711	if (attrib(I, "isSB") != 1) {
5712	I = groebner(I);
5713	};
5714	iDim = dim(I);
5715
5716	// Checks if only equiradical is required
5717	if (size(#) > 0) {
5718	il = #[1];
5719	};
5720
5721	while(stop == 0) {
5722	dbprint (printlevel-voice, "// We call radLoopR to find new prime ideals.");
5723	primes = newRadicalSLIteration(I, rad); // A list of primes or intersections of primes, not included in P
5724	dbprint (printlevel - voice, "// Output of Iteration Step:");
5725	dbprint (printlevel - voice, primes);
5726	if (size(primes) > 0) {
5727	dbprint (printlevel - voice, "// We intersect P with the ideal just obtained.");
5728	for(k = 1; k <= size(primes); k++) {
5729	rad = intersect(rad, primes[k]);
5730	if (il == 1){
5731	if (attrib(primes[k], "isSB") != 1) {
5732	primes[k] = groebner(primes[k]);
5733	};
5734	if (iDim == dim(primes[k])) {
5735	equiRad = intersect(equiRad, primes[k]);
5736	};
5737	};
5738	};
5739	} else {
5740	stop = 1
5741	};
5742	};
5743	if (il == 0) {
5744	return(rad);
5745	} else {
5746	return(equiRad);
5747	};
5748	};
5749
5750	//////////////////////////////////////////////////////////////////////////
5751	// Based on radicalKL.
5752	// It contains all of proc radicalKL except the recursion call.
5753
5754	// Output:
5755	// #1 -> output ideal, the part of the radical that has been computed
5756	// #2 -> complementary ideal, the part of the ideal I whose radical remains to be computed
5757	// = (I, h) in KL algorithm
5758	// This is not used in the new algorithm. It is part of KL algorithm
5759	// #3 -> done, 1: output = radical, there is no need to continue
5760	// 0: radical = output \cap \sqrt{complementary ideal}
5761	// This is not used in the new algorithm. It is part of KL algorithm
5762
5763	static proc newRadicalReduction(ideal I, ideal ser, int allIndep, list #) {
5764	// allMaximal 1 -> Indicates that the reduction to the zerodim case
5765	// must be done for all indep set of the leading terms ideal
5766	// 0 -> Otherwise
5767	// ideal ser Only for radicalKL. (Same as in radicalKL)
5768	// list # Only for radicalKL (If #[1] = 1, only equiradical is required.
5769	// It is used to set the value of done.)
5770
5771	attrib(I, "isSB", 1); // I needs to be a reduced standard basis
5772	list indep, fett;
5773	intvec @w, @hilb, op;
5774	int @wr, @n, @m, lauf, di;
5775	ideal fac, @h, collectrad, lsau;
5776	poly @q;
5777	string @va, quotring;
5778
5779	def @P = basering;
5780	int jdim = dim(I); // Computes the dimension of I
5781	int homo = homog(I); // Finds out if I is homogeneous
5782	ideal rad = ideal(1); // The unit ideal
5783	ideal te = ser;
5784	if(size(#) > 0)
5785	{
5786	@wr = #[1];
5787	}
5788	if(homo == 1)
5789	{
5790	for(@n = 1; @n <= nvars(basering); @n++)
5791	{
5792	@w[@n] = ord(var(@n));
5793	}
5794	@hilb = hilb(I, 1, @w);
5795	}
5796
5797	// SL 2006.04.11 1 Debug messages
5798	dbprint(printlevel-voice, "//Computes the radical of the ideal:", I);
5799	// SL 2006.04.11 2 Debug messages
5800
5801	//---------------------------------------------------------------------------
5802	//j is the ring
5803	//---------------------------------------------------------------------------
5804
5805	if (jdim==-1)
5806	{
5807	return(ideal(1), ideal(1), 1);
5808	}
5809
5810	//---------------------------------------------------------------------------
5811	//the zero-dimensional case
5812	//---------------------------------------------------------------------------
5813
5814	if (jdim==0)
5815	{
5816	return(zeroRad(I), ideal(1), 1);
5817	}
5818
5819	//-------------------------------------------------------------------------
5820	//search for a maximal independent set indep,i.e.
5821	//look for subring such that the intersection with the ideal is zero
5822	//j intersected with K[var(indep[3]+1),...,var(nvar)] is zero,
5823	//indep[1] is the new varstring, indep[2] the string for the block-ordering
5824	//-------------------------------------------------------------------------
5825
5826	// SL 2006-04-24 1 If allIndep = 0, then it only computes one maximal
5827	// independent set.
5828	// This looks better for the new algorithm but not for KL
5829	// algorithm
5830	list parameters = allIndep;
5831	indep = newMaxIndependSet(I, parameters);
5832	// SL 2006-04-24 2
5833
5834	for(@m = 1; @m <= size(indep); @m++)
5835	{
5836	if((indep[@m][1] == varstr(basering)) && (@m == 1))
5837	//this is the good case, nothing to do, just to have the same notations
5838	//change the ring
5839	{
5840	execute("ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
5841	+ordstr(basering)+");");
5842	ideal @j = fetch(@P, I);
5843	attrib(@j, "isSB", 1);
5844	}
5845	else
5846	{
5847	@va = string(maxideal(1));
5848
5849	execute("ring gnir1 = (" + charstr(basering) + "), (" + indep[@m][1] + "),("
5850	+ indep[@m][2] + ");");
5851	execute("map phi = @P," + @va + ";");
5852	if(homo == 1)
5853	{
5854	ideal @j = std(phi(I), @hilb, @w);
5855	}
5856	else
5857	{
5858	ideal @j = groebner(phi(I));
5859	}
5860	}
5861	if((deg(@j[1]) == 0) \|\| (dim(@j) < jdim))
5862	{
5863	setring @P;
5864	break;
5865	}
5866	for (lauf = 1; lauf <= size(@j); lauf++)
5867	{
5868	fett[lauf] = size(@j[lauf]);
5869	}
5870	//------------------------------------------------------------------------
5871	// We have now the following situation:
5872	// j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
5873	// to this quotientring, j is there still a standardbasis, the
5874	// leading coefficients of the polynomials there (polynomials in
5875	// K[var(nnp+1),..,var(nva)]) are collected in the list h,
5876	// we need their LCM, gh, because of the following:
5877	// let (j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..rest..]
5878	// intersected with K[var(1),...,var(nva)] is (j:gh^n)
5879	// on the other hand j = ((j, gh^n) intersected with (j : gh^n))
5880
5881	//------------------------------------------------------------------------
5882	// The arrangement for the quotientring K(var(nnp+1),..,var(nva))[..rest..]
5883	// and the map phi:K[var(1),...,var(nva)] ----->
5884	// K(var(nnpr+1),..,var(nva))[..the rest..]
5885	//------------------------------------------------------------------------
5886	quotring = newPrepareQuotientRing(nvars(basering) - indep[@m][3]);
5887	//------------------------------------------------------------------------
5888	// We pass to the quotientring K(var(nnp+1),..,var(nva))[..the rest..]
5889	//------------------------------------------------------------------------
5890
5891	execute(quotring);
5892
5893	// @j considered in the quotientring
5894	ideal @j = imap(gnir1, @j);
5895
5896	kill gnir1;
5897
5898	// j is a standardbasis in the quotientring but usually not minimal
5899	// here it becomes minimal
5900
5901	@j = clearSB(@j, fett);
5902
5903	// We need later LCM(h[1],...) = gh for saturation
5904	ideal @h;
5905	if(deg(@j[1]) > 0)
5906	{
5907	for(@n = 1; @n <= size(@j); @n++)
5908	{
5909	@h[@n] = leadcoef(@j[@n]);
5910	}
5911	op = option(get);
5912	option(redSB);
5913	@j = interred(@j); //to obtain a reduced standardbasis
5914	attrib(@j, "isSB", 1);
5915	option(set, op);
5916
5917	// SL 1 Debug messages
5918	dbprint(printlevel - voice, "zero_rad", basering, @j, dim(groebner(@j)));
5919	ideal zero_rad = zeroRad(@j);
5920	dbprint(printlevel - voice, "zero_rad passed");
5921	// SL 2
5922	}
5923	else
5924	{
5925	ideal zero_rad = ideal(1);
5926	}
5927
5928	// We need the intersection of the ideals in the list quprimary with the
5929	// polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
5930	// but fi polynomials, then the intersection of q with the polynomialring
5931	// is the saturation of the ideal generated by f1,...,fr with respect to
5932	// h which is the lcm of the leading coefficients of the fi considered in
5933	// the quotientring: this is coded in saturn
5934
5935	zero_rad = std(zero_rad);
5936
5937	ideal hpl;
5938
5939	for(@n = 1; @n <= size(zero_rad); @n++)
5940	{
5941	hpl = hpl, leadcoef(zero_rad[@n]);
5942	}
5943
5944	//------------------------------------------------------------------------
5945	// We leave the quotientring K(var(nnp+1),..,var(nva))[..the rest..]
5946	// back to the polynomialring
5947	//------------------------------------------------------------------------
5948	setring @P;
5949
5950	collectrad = imap(quring, zero_rad);
5951	lsau = simplify(imap(quring, hpl), 2);
5952	@h = imap(quring, @h);
5953
5954	kill quring;
5955
5956	// Here the intersection with the polynomialring
5957	// mentioned above is really computed
5958
5959	collectrad = sat2(collectrad, lsau)[1];
5960	if(deg(@h[1])>=0)
5961	{
5962	fac = ideal(0);
5963	for(lauf = 1; lauf <= ncols(@h); lauf++)
5964	{
5965	if(deg(@h[lauf]) > 0)
5966	{
5967	fac = fac + factorize(@h[lauf], 1);
5968	}
5969	}
5970	fac = simplify(fac, 4);
5971	@q = 1;
5972	for(lauf = 1; lauf <= size(fac); lauf++)
5973	{
5974	@q = @q * fac[lauf];
5975	}
5976	op = option(get);
5977	option(returnSB);
5978	option(redSB);
5979	I = quotient(I + ideal(@q), rad);
5980	attrib(I, "isSB", 1);
5981	option(set, op);
5982	}
5983	if((deg(rad[1]) > 0) && (deg(collectrad[1]) > 0))
5984	{
5985	rad = intersect(rad, collectrad);
5986	te = intersect(te, collectrad);
5987	te = simplify(reduce(te, I, 1), 2);
5988	}
5989	else
5990	{
5991	if(deg(collectrad[1]) > 0)
5992	{
5993	rad = collectrad;
5994	te = intersect(te, collectrad);
5995	te = simplify(reduce(te, I, 1), 2);
5996	}
5997	}
5998
5999	if((dim(I) < jdim)\|\|(size(te) == 0))
6000	{
6001	break;
6002	}
6003	if(homo==1)
6004	{
6005	@hilb = hilb(I, 1, @w);
6006	}
6007	}
6008
6009	// SL 2006.04.11 1 Debug messages
6010	dbprint (printlevel-voice, "// Part of the Radical already computed:", rad);
6011	dbprint (printlevel-voice, "// Dimension:", dim(groebner(rad)));
6012	// SL 2006.04.11 2 Debug messages
6013
6014	// SL 2006.04.21 1 New variable "done".
6015	// It tells if the radical is already computed or
6016	// if it still has to be computed the radical of the new ideal I
6017	int done;
6018	if(((@wr == 1) && (dim(I)<jdim)) \|\| (deg(I[1])==0) \|\| (size(te) == 0))
6019	{
6020	done = 1;
6021	} else {
6022	done = 0;
6023	};
6024	// SL 2006.04.21 2
6025
6026	// SL 2006.04.21 1 See details of the output at the beggining of this proc.
6027	list result = rad, I, done;
6028	return(result);
6029	// SL 2006.04.21 2
6030	};
6031
6032
6033
6034	///////////////////////////////////////////////////////////////////////////////
6035	// Given an ideal I and an ideal P (intersection of some minimal prime ideals
6036	// associated to I), it calculates the intersection of new minimal prime ideals
6037	// associated to I which where not used to calculate P.
6038	// This version uses ZD Radical in the zerodimensional case.
6039	static proc newRadicalSLIteration (ideal I, ideal P);
6040	// Input: I, ideal. The ideal from which new prime components will be obtained.
6041	// P, ideal. Intersection of some prime ideals of I.
6042	// Output: ideal. Intersection of some primes of I different from the ones in P.
6043	{
6044	int k = 1; // Counter
6045	int good = 0; // Checks if an element of P is in rad(I)
6046
6047	dbprint (printlevel-voice, "// We search for an element in P - sqrt(I).");
6048	while ((k <= size(P)) and (good == 0)) {
6049	dbprint (printlevel-voice, "// We try with:", P[k]);
6050	good = 1 - rad_con(P[k], I);
6051	k++;
6052	};
6053	k--;
6054	if (good == 0) {
6055	dbprint (printlevel-voice, "// No element was found, P = sqrt(I).");
6056	list emptyList = list();
6057	return (emptyList);
6058	};
6059	dbprint(printlevel - voice, "// That one was good!");
6060	dbprint(printlevel - voice, "// We saturate I with respect to this element.");
6061	if (P[k] != 1) {
6062	ideal J = sat(I, P[k])[1];
6063	} else {
6064	dbprint(printlevel - voice, "// The polynomial is 1, the saturation in not actually computed.");
6065	ideal J = I;
6066	};
6067
6068	// We now call proc radicalNew;
6069	dbprint(printlevel - voice, "// We do the reduction to the zerodimensional case, via radical.");
6070	dbprint(printlevel - voice, "// The ideal is ", J);
6071	dbprint(printlevel - voice, "// The dimension is ", dim(groebner(J)));
6072
6073	int allMaximal = 0; // Compute the zerodim reduction for only one indep set.
6074	ideal re = 1; // No reduction is need, there are not redundant components.
6075	list emptyList = list(); // Look for primes of any dimension, not only of max dimension.
6076	list result = newRadicalReduction(J, re, allMaximal, emptyList);
6077
6078	return(result[1]);
6079	};
6080
6081
6082
6083	///////////////////////////////////////////////////////////////////////////////////
6084	// Based on maxIndependSet
6085	// Added list # as parameter
6086	// If the first element of # is 0, the output is only 1 max indep set.
6087	// If no list is specified or #[1] = 1, the output is all the max indep set of the
6088	// leading terms ideal. This is the original output of maxIndependSet
6089
6090	// The ordering given in the output has been changed to block dp instead of lp.
6091
6092	static proc newMaxIndependSet(ideal j, list #)
6093	// #[1]=0 --> computes only one max indep set
6094	"USAGE: maxIndependentSet(i); i ideal
6095	RETURN: list = #1. new varstring with the maximal independent set at the end,
6096	#2. ordstring with the corresponding block ordering,
6097	#3. the number of independent variables
6098	// SL #3 explanation modified,
6099	// it was: integer where the independent set starts in the varstring
6100	NOTE:
6101	EXAMPLE: example maxIndependentSet; shows an example
6102	"
6103	{
6104	int n, k, di;
6105	list resu, hilf;
6106	string var1, var2;
6107	list v = indepSet(j, 0);
6108
6109	// SL 2006.04.21 1 Lines modified to use only one independent Set
6110	int allMaximal;
6111	if (size(#) > 0) {
6112	allMaximal = #[1];
6113	} else {
6114	allMaximal = 1;
6115	};
6116
6117	int nMax;
6118	if (allMaximal == 1) {
6119	nMax = size(v);
6120	} else {
6121	nMax = 1;
6122	};
6123
6124	for(n = 1; n <= nMax; n++)
6125	// SL 2006.04.21 2
6126	{
6127	di = 0;
6128	var1 = "";
6129	var2 = "";
6130	for(k = 1; k <= size(v[n]); k++)
6131	{
6132	if(v[n][k] != 0)
6133	{
6134	di++;
6135	var2 = var2 + "var(" + string(k) + "), ";
6136	}
6137	else
6138	{
6139	var1 = var1 + "var(" + string(k) + "), ";
6140	}
6141	}
6142	if(di > 0)
6143	{
6144	var1 = var1 + var2;
6145	var1 = var1[1..size(var1) - 2]; // The "- 2" removes the trailer comma
6146	hilf[1] = var1;
6147	// SL 2006.21.04 1 The order is now block dp instead of lp
6148	hilf[2] = "dp(" + string(nvars(basering) - di) + "), dp(" + string(di) + ")";
6149	// SL 2006.21.04 2
6150	hilf[3] = di;
6151	resu[n] = hilf;
6152	}
6153	else
6154	{
6155	resu[n] = varstr(basering), ordstr(basering), 0;
6156	}
6157	}
6158	return(resu);
6159	}
6160	example
6161	{ "EXAMPLE:"; echo = 2;
6162	ring s1 = (0, x, y), (a, b, c, d, e, f, g), lp;
6163	ideal i = ea - fbg, fa + be, ec - fdg, fc + de;
6164	i = std(i);
6165	list l = newMaxIndependSet(i);
6166	l;
6167	i = i, g;
6168	l = newMaxIndependSet(i);
6169	l;
6170
6171	ring s = 0, (x, y, z), lp;
6172	ideal i = z, yx;
6173	list l = newMaxIndependSet(i);
6174	l;
6175	}
6176
6177
6178	///////////////////////////////////////////////////////////////////////////////
6179	// based on prepareQuotientring
6180	// The order returned is now (C, dp) instead of (C, lp)
6181
6182	static proc newPrepareQuotientRing (int nnp)
6183	"USAGE: newPrepareQuotientRing(nnp); nnp int
6184	RETURN: string = to define Kvar(nnp+1),...,var(nvars)[..rest ]
6185	NOTE:
6186	EXAMPLE: example newPrepareQuotientRing; shows an example
6187	"
6188	{
6189	ideal @ih,@jh;
6190	int npar=npars(basering);
6191	int @n;
6192
6193	string quotring= "ring quring = ("+charstr(basering);
6194	for(@n = nnp + 1; @n <= nvars(basering); @n++)
6195	{
6196	quotring = quotring + ", var(" + string(@n) + ")";
6197	@ih = @ih + var(@n);
6198	}
6199
6200	quotring = quotring+"),(var(1)";
6201	@jh = @jh + var(1);
6202	for(@n = 2; @n <= nnp; @n++)
6203	{
6204	quotring = quotring + ", var(" + string(@n) + ")";
6205	@jh = @jh + var(@n);
6206	}
6207	// SL 2006-04-21 1 The order returned is now (C, dp) instead of (C, lp)
6208	quotring = quotring + "), (C, dp);";
6209	// SL 2006-04-21 2
6210
6211	return(quotring);
6212	}
6213	example
6214	{ "EXAMPLE:"; echo = 2;
6215	ring s1=(0,x),(a,b,c,d,e,f,g),lp;
6216	def @Q=basering;
6217	list l= newPrepareQuotientRing(3);
6218	l;
6219	execute(l[1]);
6220	execute(l[2]);
6221	basering;
6222	phi;
6223	setring @Q;
6224
6225	}
6226
6227	///////////////////////////////////////////////////////////////////////////////
6228	proc prepareAss(ideal i)
6229	"USAGE: prepareAss(i); i ideal
6230	RETURN: list, the radicals of the maximal dimensional components of i.
6231	NOTE: Uses algorithm of Eisenbud/Huneke/Vasconcelos.
6232	EXAMPLE: example prepareAss; shows an example
6233	"
6234	{
6235	if(ord_test(basering)!=1)
6236	{
6237	ERROR(
6238	"// Not implemented for this ordering, please change to global ordering."
6239	);
6240	}
6241	ideal j=std(i);
6242	int cod=nvars(basering)-dim(j);
6243	int e;
6244	list er;
6245	ideal ann;
6246	if(homog(i)==1)
6247	{
6248	list re=sres(j,0); //the resolution
6249	re=minres(re); //minimized resolution
6250	}
6251	else
6252	{
6253	list re=mres(i,0);
6254	}
6255	for(e=cod;e<=nvars(basering);e++)
6256	{
6257	ann=AnnExt_R(e,re);
6258
6259	if(nvars(basering)-dim(std(ann))==e)
6260	{
6261	er[size(er)+1]=equiRadical(ann);
6262	}
6263	}
6264	return(er);
6265	}
6266	example
6267	{ "EXAMPLE:"; echo = 2;
6268	ring r = 0,(x,y,z),dp;
6269	poly p = z2+1;
6270	poly q = z3+2;
6271	ideal i = p*q^2,y-z2;
6272	list pr = prepareAss(i);
6273	pr;
6274	}
6275	///////////////////////////////////////////////////////////////////////////////
6276	proc equidimMaxEHV(ideal i)
6277	"USAGE: equidimMaxEHV(i); i ideal
6278	RETURN: ideal, the equidimensional component (of maximal dimension) of i.
6279	NOTE: Uses algorithm of Eisenbud, Huneke and Vasconcelos.
6280	EXAMPLE: example equidimMaxEHV; shows an example
6281	"
6282	{
6283	if(ord_test(basering)!=1)
6284	{
6285	ERROR(
6286	"// Not implemented for this ordering, please change to global ordering."
6287	);
6288	}
6289	ideal j=groebner(i);
6290	int cod=nvars(basering)-dim(j);
6291	int e;
6292	ideal ann;
6293	if(homog(i)==1)
6294	{
6295	list re=sres(j,0); //the resolution
6296	re=minres(re); //minimized resolution
6297	}
6298	else
6299	{
6300	list re=mres(i,0);
6301	}
6302	ann=AnnExt_R(cod,re);
6303	return(ann);
6304	}
6305	example
6306	{ "EXAMPLE:"; echo = 2;
6307	ring r = 0,(x,y,z),dp;
6308	ideal i=intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5));
6309	equidimMaxEHV(i);
6310	}
6311
6312	proc testPrimary(list pr, ideal k)
6313	"USAGE: testPrimary(pr,k); pr a list, k an ideal.
6314	ASSUME: pr is the result of primdecGTZ(k) or primdecSY(k).
6315	RETURN: int, 1 if the intersection of the ideals in pr is k, 0 if not
6316	EXAMPLE: example testPrimary; shows an example
6317	"
6318	{
6319	int i;
6320	pr=reconvList(pr);
6321	ideal j=pr[1];
6322	for (i=2;i<=size(pr)/2;i++)
6323	{
6324	j=intersect(j,pr[2*i-1]);
6325	}
6326	return(idealsEqual(j,k));
6327	}
6328	example
6329	{ "EXAMPLE:"; echo = 2;
6330	ring r = 32003,(x,y,z),dp;
6331	poly p = z2+1;
6332	poly q = z4+2;
6333	ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
6334	list pr = primdecGTZ(i);
6335	testPrimary(pr,i);
6336	}
6337
6338	///////////////////////////////////////////////////////////////////////////////
6339	proc zerodec(ideal I)
6340	"USAGE: zerodec(I); I ideal
6341	ASSUME: I is zero-dimensional, the characteristic of the ground field is 0
6342	RETURN: list of primary ideals, the zero-dimensional decomposition of I
6343	NOTE: The algorithm (of Monico), works well only for a small total number
6344	of solutions (@code{vdim(std(I))} should be < 100) and without
6345	parameters. In practice, it works also in large characteristic p>0
6346	but may fail for small p.
6347	@* If printlevel > 0 (default = 0) additional information is displayed.
6348	EXAMPLE: example zerodec; shows an example
6349	"
6350	{
6351	if(ord_test(basering)!=1)
6352	{
6353	ERROR(
6354	"// Not implemented for this ordering, please change to global ordering."
6355	);
6356	}
6357	def R=basering;
6358	poly q;
6359	int j,time;
6360	matrix m;
6361	list re;
6362	poly va=var(1);
6363	ideal J=groebner(I);
6364	ideal ba=kbase(J);
6365	int d=vdim(J);
6366	dbprint(printlevel-voice+2,"// multiplicity of ideal : "+ string(d));
6367	//------ compute matrix of multiplication on R/I with generic element p -----
6368	int e=nvars(basering);
6369	poly p=randomLast(100)[e]+random(-50,50); //the generic element
6370	matrix n[d][d];
6371	time = timer;
6372	for(j=2;j<=e;j++)
6373	{
6374	va=va*var(j);
6375	}
6376	for(j=1;j<=d;j++)
6377	{
6378	q=reduce(p*ba[j],J);
6379	m=coeffs(q,ba,va);
6380	n[j,1..d]=m[1..d,1];
6381	}
6382	dbprint(printlevel-voice+2,
6383	"// time for computing multiplication matrix (with generic element) : "+
6384	string(timer-time));
6385	//---------------- compute characteristic polynomial of matrix --------------
6386	execute("ring P1=("+charstr(R)+"),T,dp;");
6387	matrix n=imap(R,n);
6388	time = timer;
6389	poly charpol=det(n-T*freemodule(d));
6390	dbprint(printlevel-voice+2,"// time for computing char poly: "+
6391	string(timer-time));
6392	//------------------- factorize characteristic polynomial -------------------
6393	//check first if constant term of charpoly is != 0 (which is true for
6394	//sufficiently generic element)
6395	if(charpol[size(charpol)]!=0)
6396	{
6397	time = timer;
6398	list fac=factor(charpol);
6399	testFactor(fac,charpol);
6400	dbprint(printlevel-voice+2,"// time for factorizing char poly: "+
6401	string(timer-time));
6402	int f=size(fac[1]);
6403	//--------------------------- the irreducible case --------------------------
6404	if(f==1)
6405	{
6406	setring R;
6407	re=I;
6408	return(re);
6409	}
6410	//---------------------------- the reducible case ---------------------------
6411	//if f_i are the irreducible factors of charpoly, mult=ri, then <I,g_i^ri>
6412	//are the primary components where g_i = f_i(p). However, substituting p in
6413	//f_i may result in a huge object although the final result may be small.
6414	//Hence it is better to simultaneously reduce with I. For this we need a new
6415	//ring.
6416	execute("ring P=("+charstr(R)+"),(T,"+varstr(R)+"),(dp(1),dp);");
6417	list rfac=imap(P1,fac);
6418	intvec ov=option(get);;
6419	option(redSB);
6420	list re1;
6421	ideal new = T-imap(R,p),imap(R,J);
6422	attrib(new, "isSB",1); //we know that new is a standard basis
6423	for(j=1;j<=f;j++)
6424	{
6425	re1[j]=reduce(rfac[1][j]^rfac[2][j],new);
6426	}
6427	setring R;
6428	re = imap(P,re1);
6429	for(j=1;j<=f;j++)
6430	{
6431	J=I,re[j];
6432	re[j]=interred(J);
6433	}
6434	option(set,ov);
6435	return(re);
6436	}
6437	else
6438	//------------------- choice of generic element failed -------------------
6439	{
6440	dbprint(printlevel-voice+2,"// try new generic element!");
6441	setring R;
6442	return(zerodec(I));
6443	}
6444	}
6445	example
6446	{ "EXAMPLE:"; echo = 2;
6447	ring r = 0,(x,y),dp;
6448	ideal i = x2-2,y2-2;
6449	list pr = zerodec(i);
6450	pr;
6451	}
6452	////////////////////////////////////////////////////////////////////////////
6453	/*
6454	//Beispiele Wenk-Dipl (in ~/Texfiles/Diplom/Wenk/Examples/)
6455	//Zeiten: Singular/Singular/Singular -r123456789 -v :wilde13 (PentiumPro200)
6456	//Singular for HPUX-9 version 1-3-8 (2000060214) Jun 2 2000 15:31:26
6457	//(wilde13)
6458
6459	//1. vdim=20, 3 Komponenten
6460	//zerodec-time:2(1) (matrix:1 charpoly:0 factor:1)
6461	//primdecGTZ-time: 1(0)
6462	ring rs= 0,(a,b,c),dp;
6463	poly f1= a^2bc + ab^2c + abc^2 + abc + ab + ac + b*c;
6464	poly f2= a^2b^2c + ab^2c^2 + a^2bc + abc + b*c + a + c;
6465	poly f3= a^2b^2c^2 + a^2b^2c + ab^2c + abc + a*c + c + 1;
6466	ideal gls=f1,f2,f3;
6467	int time=timer;
6468	printlevel =1;
6469	time=timer; list pr1=zerodec(gls); timer-time;size(pr1);
6470	time=timer; list pr =primdecGTZ(gls); timer-time;size(pr);
6471	time=timer; ideal ra =radical(gls); timer-time;size(pr);
6472
6473	//2.cyclic5 vdim=70, 20 Komponenten
6474	//zerodec-time:36(28) (matrix:1(0) charpoly:18(19) factor:17(9)
6475	//primdecGTZ-time: 28(5)
6476	//radical : 0
6477	ring rs= 0,(a,b,c,d,e),dp;
6478	poly f0= a + b + c + d + e + 1;
6479	poly f1= a + b + c + d + e;
6480	poly f2= ab + bc + cd + ae + d*e;
6481	poly f3= abc + bcd + abe + ade + cde;
6482	poly f4= abcd + abce + abde + acde + bcd*e;
6483	poly f5= abcde - 1;
6484	ideal gls= f1,f2,f3,f4,f5;
6485
6486	//3. random vdim=40, 1 Komponente
6487	//zerodec-time:126(304) (matrix:1 charpoly:115(298) factor:10(5))
6488	//primdecGTZ-time:17 (11)
6489	ring rs=0,(x,y,z),dp;
6490	poly f1=2x^2 + 4x + 3y^2 + 7xz + 9yz + 5z^2;
6491	poly f2=7x^3 + 8xy + 12y^2 + 18xz + 3y^4z + 10*z^3 + 12;
6492	poly f3=3x^4 + 1xyz + 6y^3 + 3xz^2 + 2yz^2 + 4z^2 + 5;
6493	ideal gls=f1,f2,f3;
6494
6495	//4. introduction into resultants, sturmfels, vdim=28, 1 Komponente
6496	//zerodec-time:4 (matrix:0 charpoly:0 factor:4)
6497	//primdecGTZ-time:1
6498	ring rs=0,(x,y),dp;
6499	poly f1= x4+y4-1;
6500	poly f2= x5y2-4x3y3+x2y5-1;
6501	ideal gls=f1,f2;
6502
6503	//5. 3 quadratic equations with random coeffs, vdim=8, 1 Komponente
6504	//zerodec-time:0(0) (matrix:0 charpoly:0 factor:0)
6505	//primdecGTZ-time:1(0)
6506	ring rs=0,(x,y,z),dp;
6507	poly f1=2x^2 + 4xy + 3y^2 + 7xz + 9yz + 5*z^2 + 2;
6508	poly f2=7x^2 + 8xy + 12y^2 + 18xz + 3yz + 10*z^2 + 12;
6509	poly f3=3x^2 + 1xy + 6y^2 + 3xz + 2yz + 4*z^2 + 5;
6510	ideal gls=f1,f2,f3;
6511
6512	//6. 3 polys vdim=24, 1 Komponente
6513	// run("ex14",2);
6514	//zerodec-time:5(4) (matrix:0 charpoly:3(3) factor:2(1))
6515	//primdecGTZ-time:4 (2)
6516	ring rs=0,(x1,x2,x3,x4),dp;
6517	poly f1=16x1^2 + 3x2^2 + 5x3^4 - 1 - 4x4 + x4^3;
6518	poly f2=5x1^3 + 3x2^2 + 4x3^2x4 + 2x1x4 - 1 + x4 + 4*x1 + x2 + x3 + x4;
6519	poly f3=-4x1^2 + x2^2 + x3^2 - 3 + x4^2 + 4x1 + x2 + x3 + x4;
6520	poly f4=-4*x1 + x2 + x3 + x4;
6521	ideal gls=f1,f2,f3,f4;
6522
6523	//7. ex43, PoSSo, caprasse vdim=56, 16 Komponenten
6524	//zerodec-time:23(15) (matrix:0 charpoly:16(13) factor:3(2))
6525	//primdecGTZ-time:3 (2)
6526	ring rs= 0,(y,z,x,t),dp;
6527	ideal gls=y^2z+2yxt-z-2*x,
6528	4y^2zx-zx^3+2y^3t+4yx^2t-10y^2+4zx+4x^2-10y*t+2,
6529	2yzt+xt^2-2*z-x,
6530	-z^3x+4yz^2t+4zxt^2+2yt^3+4z^2+4zx-10yt-10*t^2+2;
6531
6532	//8. Arnborg-System, n=6 (II), vdim=156, 90 Komponenten
6533	//zerodec-time (char32003):127(45)(matrix:2(0) charpoly:106(37) factor:16(7))
6534	//primdecGTZ-time(char32003) :81 (18)
6535	//ring rs= 0,(a,b,c,d,x,f),dp;
6536	ring rs= 32003,(a,b,c,d,x,f),dp;
6537	ideal gls=a+b+c+d+x+f, ab+bc+cd+dx+xf+af, abc+bcd+cdx+d*xf+axf+abf,
6538	abcd+bcdx+cdxf+adxf+abxf+abcf, abcdx+bcdxf+acdxf+abd*xf+abcxf+abcdf,
6539	abcd*xf-1;
6540
6541	//9. ex42, PoSSo, Methan6_1, vdim=27, 2 Komponenten
6542	//zerodec-time:610 (matrix:10 charpoly:557 factor:26)
6543	//primdecGTZ-time: 118
6544	//zerodec-time(char32003):2
6545	//primdecGTZ-time(char32003):4
6546	//ring rs= 0,(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10),dp;
6547	ring rs= 32003,(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10),dp;
6548	ideal gls=64x2x7-10x1x8+10x7x9+11x7x10-320000*x1,
6549	-32x2x7-5x2x8-5x2x10+160000x1-5000x2,
6550	-x3x8+x6x8+x9x10+210x6+1300000,
6551	-x4*x8+700000,
6552	x10^2-2*x5,
6553	-x6x8+x7x9-210*x6,
6554	-64x2x7-10x7x9-11x7x10+320000x1-16x7+7000000,
6555	-10x1x8-10x2x8-10x3x8-10x4x8-10x6x8+10x2x10+11x7x10
6556	+20000x2+14x5,
6557	x4x8-x7x9-x9x10-410x9,
6558	10x2x8+10x3x8+10x6x8+10x7x9-10x2x10-11x7x10-10x9x10
6559	-10x10^2+1400x6-4200*x10;
6560
6561	//10. ex33, walk-s7, Diplomarbeit von Tim, vdim=114
6562	//zerfaellt in unterschiedlich viele Komponenten in versch. Charkteristiken:
6563	//char32003:30, char0:3(2xdeg1,1xdeg112!), char181:4(2xdeg1,1xdeg28,1xdeg84)
6564	//char 0: zerodec-time:10075 (ca 3h) (matrix:3 charpoly:9367, factor:680
6565	// + 24 sec fuer Normalform (anstatt einsetzen), total [29623k])
6566	// primdecGTZ-time: 214
6567	//char 32003:zerodec-time:197(68) (matrix:2(1) charpoly:173(60) factor:15(6))
6568	// primdecGTZ-time:14 (5)
6569	//char 181:zerodec-time:(87) (matrix:(1) charpoly:(58) factor:(25))
6570	// primdecGTZ-time:(2)
6571	//in char181 stimmen Ergebnisse von zerodec und primdecGTZ ueberein (gecheckt)
6572
6573	//ring rs= 0,(a,b,c,d,e,f,g),dp;
6574	ring rs= 32003,(a,b,c,d,e,f,g),dp;
6575	poly f1= 2gb + 2fc + 2ed + a2 + a;
6576	poly f2= 2gc + 2fd + e2 + 2ba + b;
6577	poly f3= 2gd + 2fe + 2ca + c + b2;
6578	poly f4= 2ge + f2 + 2da + d + 2cb;
6579	poly f5= 2fg + 2ea + e + 2db + c2;
6580	poly f6= g2 + 2fa + f + 2eb + 2dc;
6581	poly f7= 2ga + g + 2fb + 2ec + d2;
6582	ideal gls= f1,f2,f3,f4,f5,f6,f7;
6583
6584	~/Singular/Singular/Singular -r123456789 -v
6585	LIB"./primdec.lib";
6586	timer=1;
6587	int time=timer;
6588	printlevel =1;
6589	option(prot,mem);
6590	time=timer; list pr1=zerodec(gls); timer-time;
6591
6592	time=timer; list pr =primdecGTZ(gls); timer-time;
6593	time=timer; list pr =primdecSY(gls); timer-time;
6594	time=timer; ideal ra =radical(gls); timer-time;size(pr);
6595	LIB"all.lib";
6596
6597	ring R=0,(a,b,c,d,e,f),dp;
6598	ideal I=cyclic(6);
6599	minAssGTZ(I);
6600
6601
6602	ring S=(2,a,b),(x,y),lp;
6603	ideal I=x8-b,y4+a;
6604	minAssGTZ(I);
6605
6606	ring S1=2,(x,y,a,b),lp;
6607	ideal I=x8-b,y4+a;
6608	minAssGTZ(I);
6609
6610
6611	ring S2=(2,z),(x,y),dp;
6612	minpoly=z2+z+1;
6613	ideal I=y3+y+1,x4+x+1;
6614	primdecGTZ(I);
6615	minAssGTZ(I);
6616
6617	ring S3=2,(x,y,z),dp;
6618	ideal I=y3+y+1,x4+x+1,z2+z+1;
6619	primdecGTZ(I);
6620	minAssGTZ(I);
6621
6622
6623	ring R1=2,(x,y,z),lp;
6624	ideal I=y6+y5+y3+y2+1,x4+x+1,z2+z+1;
6625	primdecGTZ(I);
6626	minAssGTZ(I);
6627
6628
6629	ring R2=(2,z),(x,y),lp;
6630	minpoly=z3+z+1;
6631	ideal I=y2+y+(z2+z+1),x4+x+1;
6632	primdecGTZ(I);
6633	minAssGTZ(I);
6634
6635	*/

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

Download in other formats:

Original Format