Context Navigation

primdec.lib @ 8942a5

spielwiese

Last change on this file since 8942a5 was 8942a5, checked in by Gert-Martin Greuel <greuel@…>, 23 years ago
* GMG: Kosmetik git-svn-id: file:///usr/local/Singular/svn/trunk@4982 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 120.2 KB

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

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

Download in other formats:

Original Format