Context Navigation

primdec.lib @ 55fcff

spielwiese

Last change on this file since 55fcff was 55fcff, checked in by Gert-Martin Greuel <greuel@…>, 23 years ago
* Greuel: modification of decomp git-svn-id: file:///usr/local/Singular/svn/trunk@4441 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 120.7 KB

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

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

Download in other formats:

Original Format