Context Navigation

primdec.lib @ 0ae4ce

spielwiese

Last change on this file since 0ae4ce was 0ae4ce, checked in by Anne Frühbis-Krüger <anne@…>, 23 years ago
*anne: added category to libraries for "Commutative Algebra" git-svn-id: file:///usr/local/Singular/svn/trunk@4941 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 120.7 KB

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

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

Download in other formats:

Original Format