Context Navigation

primdec.lib @ ba94539

spielwiese

Last change on this file since ba94539 was ba94539, checked in by Gerhard Pfister <pfister@…>, 24 years ago
neue Prozedur equidim eingebaut. git-svn-id: file:///usr/local/Singular/svn/trunk@4291 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 106.2 KB

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

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

Download in other formats:

Original Format