Context Navigation

primdec.lib @ 333b889

spielwiese

Last change on this file since 333b889 was 333b889, checked in by Hans Schönemann <hannes@…>, 24 years ago
*hannes: fixed primdecSY(ideal(1)); git-svn-id: file:///usr/local/Singular/svn/trunk@4232 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 103.7 KB

Line
1	// $Id: primdec.lib,v 1.54 2000-03-29 15:26:52 Singular 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.54 2000-03-29 15:26:52 Singular 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
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	///////////////////////////////////////////////////////////////////////////////
1764	proc decomp(ideal i,list #)
1765	"USAGE: decomp(i); i ideal (for primary decomposition) (resp.
1766	decomp(i,1); (for the minimal associated primes) )
1767	RETURN: list = list of primary ideals and their associated primes
1768	(at even positions in the list)
1769	(resp. a list of the minimal associated primes)
1770	NOTE: Algorithm of Gianni, Traeger, Zacharias
1771	EXAMPLE: example decomp; shows an example
1772	"
1773	{
1774	def @P = basering;
1775	list primary,indep,ltras;
1776	intvec @vh,isat;
1777	int @wr,@k,@n,@m,@n1,@n2,@n3,homo,seri,keepdi;
1778	ideal peek=i;
1779	ideal ser,tras;
1780
1781	if(size(#)>0)
1782	{
1783	if((#[1]==1)\|\|(#[1]==2))
1784	{
1785	@wr=#[1];
1786	if(size(#)>1)
1787	{
1788	seri=1;
1789	peek=#[2];
1790	ser=#[3];
1791	}
1792	}
1793	else
1794	{
1795	seri=1;
1796	peek=#[1];
1797	ser=#[2];
1798	}
1799	}
1800
1801	homo=homog(i);
1802	if((find(ordstr(basering),"w")!=0)\|\|(find(ordstr(basering),"W")!=0))
1803	{
1804	homo=0;
1805	}
1806
1807	if(homo==1)
1808	{
1809	if(attrib(i,"isSB")!=1)
1810	{
1811	ltras=mstd(i);
1812	attrib(ltras[1],"isSB",1);
1813	}
1814	else
1815	{
1816	ltras=i,i;
1817	}
1818	tras=ltras[1];
1819	if(dim(tras)==0)
1820	{
1821	primary[1]=ltras[2];
1822	primary[2]=maxideal(1);
1823	if(@wr>0)
1824	{
1825	list l;
1826	l[1]=maxideal(1);
1827	l[2]=maxideal(1);
1828	return(l);
1829	}
1830	return(primary);
1831	}
1832	intvec @hilb=hilb(tras,1);
1833	intvec keephilb=@hilb;
1834	}
1835
1836	//----------------------------------------------------------------
1837	//i is the zero-ideal
1838	//----------------------------------------------------------------
1839
1840	if(size(i)==0)
1841	{
1842	primary=i,i;
1843	return(primary);
1844	}
1845
1846	//----------------------------------------------------------------
1847	//pass to the lexicographical ordering and compute a standardbasis
1848	//----------------------------------------------------------------
1849
1850	execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),(C,lp);";
1851	option(redSB);
1852
1853	ideal ser=fetch(@P,ser);
1854
1855	if(homo==1)
1856	{
1857	if((ordstr(@P)[1]!="(C,lp)")&&(ordstr(@P)[3]!="(C,lp)"))
1858	{
1859	ideal @j=std(fetch(@P,i),@hilb);
1860	}
1861	else
1862	{
1863	ideal @j=fetch(@P,tras);
1864	attrib(@j,"isSB",1);
1865	}
1866	}
1867	else
1868	{
1869	ideal @j=groebner(fetch(@P,i));
1870	}
1871	option(noredSB);
1872	if(seri==1)
1873	{
1874	ideal peek=fetch(@P,peek);
1875	attrib(peek,"isSB",1);
1876	}
1877	else
1878	{
1879	ideal peek=@j;
1880	}
1881	if(size(ser)==0)
1882	{
1883	ideal fried;
1884	@n=size(@j);
1885	for(@k=1;@k<=@n;@k++)
1886	{
1887	if(deg(lead(@j[@k]))==1)
1888	{
1889	fried[size(fried)+1]=@j[@k];
1890	@j[@k]=0;
1891	}
1892	}
1893	if(size(fried)>0)
1894	{
1895	@j=simplify(@j,2);
1896	attrib(@j,"isSB",1);
1897	list pr=decomp(@j);
1898	for(@k=1;@k<=size(pr);@k++)
1899	{
1900	@j=pr[@k]+fried;
1901	pr[@k]=@j;
1902	}
1903	setring @P;
1904	return(fetch(gnir,pr));
1905	}
1906	}
1907
1908	//----------------------------------------------------------------
1909	//j is the ring
1910	//----------------------------------------------------------------
1911
1912	if (dim(@j)==-1)
1913	{
1914	setring @P;
1915	return(ideal(0));
1916	}
1917
1918	//----------------------------------------------------------------
1919	// the case of one variable
1920	//----------------------------------------------------------------
1921
1922	if(nvars(basering)==1)
1923	{
1924
1925	list fac=factor(@j[1]);
1926	list gprimary;
1927	for(@k=1;@k<=size(fac[1]);@k++)
1928	{
1929	if(@wr==0)
1930	{
1931	gprimary[2*@k-1]=ideal(fac[1][@k]^fac[2][@k]);
1932	gprimary[2*@k]=ideal(fac[1][@k]);
1933	}
1934	else
1935	{
1936	gprimary[2*@k-1]=ideal(fac[1][@k]);
1937	gprimary[2*@k]=ideal(fac[1][@k]);
1938	}
1939	}
1940	setring @P;
1941	primary=fetch(gnir,gprimary);
1942
1943	return(primary);
1944	}
1945
1946	//------------------------------------------------------------------
1947	//the zero-dimensional case
1948	//------------------------------------------------------------------
1949	if (dim(@j)==0)
1950	{
1951	option(redSB);
1952	list gprimary= zero_decomp(@j,ser,@wr);
1953	option(noredSB);
1954	setring @P;
1955	primary=fetch(gnir,gprimary);
1956	if(size(ser)>0)
1957	{
1958	primary=cleanPrimary(primary);
1959	}
1960	return(primary);
1961	}
1962
1963	poly @gs,@gh,@p;
1964	string @va,quotring;
1965	list quprimary,htprimary,collectprimary,lsau,lnew,allindep,restindep;
1966	ideal @h;
1967	int jdim=dim(@j);
1968	list fett;
1969	int lauf,di,newtest;
1970	//------------------------------------------------------------------
1971	//search for a maximal independent set indep,i.e.
1972	//look for subring such that the intersection with the ideal is zero
1973	//j intersected with K[var(indep[3]+1),...,var(nvar] is zero,
1974	//indep[1] is the new varstring and indep[2] the string for block-ordering
1975	//------------------------------------------------------------------
1976
1977	if(@wr!=1)
1978	{
1979	allindep=independSet(@j);
1980	for(@m=1;@m<=size(allindep);@m++)
1981	{
1982	if(allindep[@m][3]==jdim)
1983	{
1984	di++;
1985	indep[di]=allindep[@m];
1986	}
1987	else
1988	{
1989	lauf++;
1990	restindep[lauf]=allindep[@m];
1991	}
1992	}
1993	}
1994	else
1995	{
1996	indep=maxIndependSet(@j);
1997	}
1998
1999	ideal jkeep=@j;
2000	if(ordstr(@P)[1]=="w")
2001	{
2002	execute "ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),("+ordstr(@P)+");";
2003	}
2004	else
2005	{
2006	execute "ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),(C,dp);";
2007	}
2008
2009	if(homo==1)
2010	{
2011	if((ordstr(@P)[3]=="d")\|\|(ordstr(@P)[1]=="d")\|\|(ordstr(@P)[1]=="w")
2012	\|\|(ordstr(@P)[3]=="w"))
2013	{
2014	ideal jwork=imap(@P,tras);
2015	attrib(jwork,"isSB",1);
2016	}
2017	else
2018	{
2019	ideal jwork=std(imap(gnir,@j),@hilb);
2020	}
2021
2022	}
2023	else
2024	{
2025	ideal jwork=groebner(imap(gnir,@j));
2026	}
2027	list hquprimary;
2028	poly @p,@q;
2029	ideal @h,fac,ser;
2030	di=dim(jwork);
2031	keepdi=di;
2032
2033	setring gnir;
2034	for(@m=1;@m<=size(indep);@m++)
2035	{
2036	isat=0;
2037	@n2=0;
2038	option(redSB);
2039	if((indep[@m][1]==varstr(basering))&&(@m==1))
2040	//this is the good case, nothing to do, just to have the same notations
2041	//change the ring
2042	{
2043	execute "ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
2044	+ordstr(basering)+");";
2045	ideal @j=fetch(gnir,@j);
2046	attrib(@j,"isSB",1);
2047	ideal ser=fetch(gnir,ser);
2048
2049	}
2050	else
2051	{
2052	@va=string(maxideal(1));
2053	execute "ring gnir1 = ("+charstr(basering)+"),("+indep[@m][1]+"),("
2054	+indep[@m][2]+");";
2055	execute "map phi=gnir,"+@va+";";
2056	if(homo==1)
2057	{
2058	ideal @j=std(phi(@j),@hilb);
2059	}
2060	else
2061	{
2062	ideal @j=groebner(phi(@j));
2063	}
2064	ideal ser=phi(ser);
2065
2066	}
2067	option(noredSB);
2068	if((deg(@j[1])==0)\|\|(dim(@j)<jdim))
2069	{
2070	setring gnir;
2071	break;
2072	}
2073	for (lauf=1;lauf<=size(@j);lauf++)
2074	{
2075	fett[lauf]=size(@j[lauf]);
2076	}
2077	//------------------------------------------------------------------------
2078	//we have now the following situation:
2079	//j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
2080	//to this quotientring, j is their still a standardbasis, the
2081	//leading coefficients of the polynomials there (polynomials in
2082	//K[var(nnp+1),..,var(nva)]) are collected in the list h,
2083	//we need their ggt, gh, because of the following: let
2084	//(j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..the rest..]
2085	//intersected with K[var(1),...,var(nva)] is (j:gh^n)
2086	//on the other hand j=(j,gh^n) intersected with (j:gh^n)
2087
2088	//------------------------------------------------------------------------
2089
2090	//arrangement for quotientring K(var(nnp+1),..,var(nva))[..the rest..] and
2091	//map phi:K[var(1),...,var(nva)] --->K(var(nnpr+1),..,var(nva))[..rest..]
2092	//------------------------------------------------------------------------
2093
2094	quotring=prepareQuotientring(nvars(basering)-indep[@m][3]);
2095
2096	//---------------------------------------------------------------------
2097	//we pass to the quotientring K(var(nnp+1),..,var(nva))[..the rest..]
2098	//---------------------------------------------------------------------
2099
2100	execute quotring;
2101
2102	// @j considered in the quotientring
2103	ideal @j=imap(gnir1,@j);
2104	ideal ser=imap(gnir1,ser);
2105
2106	kill gnir1;
2107
2108	//j is a standardbasis in the quotientring but usually not minimal
2109	//here it becomes minimal
2110
2111	@j=clearSB(@j,fett);
2112	attrib(@j,"isSB",1);
2113
2114	//we need later ggt(h[1],...)=gh for saturation
2115	ideal @h;
2116	if(deg(@j[1])>0)
2117	{
2118	for(@n=1;@n<=size(@j);@n++)
2119	{
2120	@h[@n]=leadcoef(@j[@n]);
2121	}
2122	//the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..]
2123	option(redSB);
2124	list uprimary= zero_decomp(@j,ser,@wr);
2125	option(noredSB);
2126	}
2127	else
2128	{
2129	list uprimary;
2130	uprimary[1]=ideal(1);
2131	uprimary[2]=ideal(1);
2132	}
2133	//we need the intersection of the ideals in the list quprimary with the
2134	//polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
2135	//but fi polynomials, then the intersection of q with the polynomialring
2136	//is the saturation of the ideal generated by f1,...,fr with respect to
2137	//h which is the lcm of the leading coefficients of the fi considered in
2138	//in the quotientring: this is coded in saturn
2139
2140	list saturn;
2141	ideal hpl;
2142
2143	for(@n=1;@n<=size(uprimary);@n++)
2144	{
2145	uprimary[@n]=interred(uprimary[@n]); // temporary fix
2146	hpl=0;
2147	for(@n1=1;@n1<=size(uprimary[@n]);@n1++)
2148	{
2149	hpl=hpl,leadcoef(uprimary[@n][@n1]);
2150	}
2151	saturn[@n]=hpl;
2152	}
2153
2154	//--------------------------------------------------------------------
2155	//we leave the quotientring K(var(nnp+1),..,var(nva))[..the rest..]
2156	//back to the polynomialring
2157	//---------------------------------------------------------------------
2158	setring gnir;
2159
2160	collectprimary=imap(quring,uprimary);
2161	lsau=imap(quring,saturn);
2162	@h=imap(quring,@h);
2163
2164	kill quring;
2165
2166
2167	@n2=size(quprimary);
2168	@n3=@n2;
2169
2170	for(@n1=1;@n1<=size(collectprimary)/2;@n1++)
2171	{
2172	if(deg(collectprimary[2*@n1][1])>0)
2173	{
2174	@n2++;
2175	quprimary[@n2]=collectprimary[2*@n1-1];
2176	lnew[@n2]=lsau[2*@n1-1];
2177	@n2++;
2178	lnew[@n2]=lsau[2*@n1];
2179	quprimary[@n2]=collectprimary[2*@n1];
2180	}
2181	}
2182
2183	//here the intersection with the polynomialring
2184	//mentioned above is really computed
2185	for(@n=@n3/2+1;@n<=@n2/2;@n++)
2186	{
2187	if(specialIdealsEqual(quprimary[2@n-1],quprimary[2@n]))
2188	{
2189	quprimary[2@n-1]=sat2(quprimary[2@n-1],lnew[2*@n-1])[1];
2190	quprimary[2@n]=quprimary[2@n-1];
2191	}
2192	else
2193	{
2194	if(@wr==0)
2195	{
2196	quprimary[2@n-1]=sat2(quprimary[2@n-1],lnew[2*@n-1])[1];
2197	}
2198	quprimary[2@n]=sat2(quprimary[2@n],lnew[2*@n])[1];
2199	}
2200	}
2201
2202	if(size(@h)>0)
2203	{
2204	//---------------------------------------------------------------
2205	//we change to @Phelp to have the ordering dp for saturation
2206	//---------------------------------------------------------------
2207	setring @Phelp;
2208	@h=imap(gnir,@h);
2209	if(@wr!=1)
2210	{
2211	@q=minSat(jwork,@h)[2];
2212	}
2213	else
2214	{
2215	fac=ideal(0);
2216	for(lauf=1;lauf<=ncols(@h);lauf++)
2217	{
2218	if(deg(@h[lauf])>0)
2219	{
2220	fac=fac+factorize(@h[lauf],1);
2221	}
2222	}
2223	fac=simplify(fac,4);
2224	@q=1;
2225	for(lauf=1;lauf<=size(fac);lauf++)
2226	{
2227	@q=@q*fac[lauf];
2228	}
2229	}
2230	jwork=std(jwork,@q);
2231	keepdi=dim(jwork);
2232	if(keepdi<di)
2233	{
2234	setring gnir;
2235	@j=imap(@Phelp,jwork);
2236	break;
2237	}
2238	if(homo==1)
2239	{
2240	@hilb=hilb(jwork,1);
2241	}
2242
2243	setring gnir;
2244	@j=imap(@Phelp,jwork);
2245	}
2246	}
2247	if((size(quprimary)==0)&&(@wr>0))
2248	{
2249	@j=ideal(1);
2250	quprimary[1]=ideal(1);
2251	quprimary[2]=ideal(1);
2252	}
2253	if((size(quprimary)==0))
2254	{
2255	keepdi=di-1;
2256	}
2257	//---------------------------------------------------------------
2258	//notice that j=sat(j,gh) intersected with (j,gh^n)
2259	//we finished with sat(j,gh) and have to start with (j,gh^n)
2260	//---------------------------------------------------------------
2261	if((deg(@j[1])!=0)&&(@wr!=1))
2262	{
2263	if(size(quprimary)>0)
2264	{
2265	setring @Phelp;
2266	ser=imap(gnir,ser);
2267	hquprimary=imap(gnir,quprimary);
2268	if(@wr==0)
2269	{
2270	ideal htest=hquprimary[1];
2271	for (@n1=2;@n1<=size(hquprimary)/2;@n1++)
2272	{
2273	htest=intersect(htest,hquprimary[2*@n1-1]);
2274	}
2275	}
2276	else
2277	{
2278	ideal htest=hquprimary[2];
2279
2280	for (@n1=2;@n1<=size(hquprimary)/2;@n1++)
2281	{
2282	htest=intersect(htest,hquprimary[2*@n1]);
2283	}
2284	}
2285
2286	if(size(ser)>0)
2287	{
2288	ser=intersect(htest,ser);
2289	}
2290	else
2291	{
2292	ser=htest;
2293	}
2294	setring gnir;
2295	ser=imap(@Phelp,ser);
2296	}
2297	if(size(reduce(ser,peek,1))!=0)
2298	{
2299	for(@m=1;@m<=size(restindep);@m++)
2300	{
2301	// if(restindep[@m][3]>=keepdi)
2302	// {
2303	isat=0;
2304	@n2=0;
2305	option(redSB);
2306
2307	if(restindep[@m][1]==varstr(basering))
2308	//the good case, nothing to do, just to have the same notations
2309	//change the ring
2310	{
2311	execute "ring gnir1 = ("+charstr(basering)+"),("+
2312	varstr(basering)+"),("+ordstr(basering)+");";
2313	ideal @j=fetch(gnir,jkeep);
2314	attrib(@j,"isSB",1);
2315	}
2316	else
2317	{
2318	@va=string(maxideal(1));
2319	execute "ring gnir1 = ("+charstr(basering)+"),("+
2320	restindep[@m][1]+"),(" +restindep[@m][2]+");";
2321	execute "map phi=gnir,"+@va+";";
2322	if(homo==1)
2323	{
2324	ideal @j=std(phi(jkeep),keephilb);
2325	}
2326	else
2327	{
2328	ideal @j=groebner(phi(jkeep));
2329	}
2330	ideal ser=phi(ser);
2331	}
2332	option(noredSB);
2333
2334	for (lauf=1;lauf<=size(@j);lauf++)
2335	{
2336	fett[lauf]=size(@j[lauf]);
2337	}
2338	//------------------------------------------------------------------
2339	//we have now the following situation:
2340	//j intersected with K[var(nnp+1),..,var(nva)] is zero so we may
2341	//pass to this quotientring, j is their still a standardbasis, the
2342	//leading coefficients of the polynomials there (polynomials in
2343	//K[var(nnp+1),..,var(nva)]) are collected in the list h,
2344	//we need their ggt, gh, because of the following:
2345	//let (j:gh^n)=(j:gh^infinity) then
2346	//j*K(var(nnp+1),..,var(nva))[..the rest..]
2347	//intersected with K[var(1),...,var(nva)] is (j:gh^n)
2348	//on the other hand j=(j,gh^n) intersected with (j:gh^n)
2349
2350	//------------------------------------------------------------------
2351
2352	//the arrangement for the quotientring
2353	// K(var(nnp+1),..,var(nva))[..the rest..]
2354	//and the map phi:K[var(1),...,var(nva)] ---->
2355	//--->K(var(nnpr+1),..,var(nva))[..the rest..]
2356	//------------------------------------------------------------------
2357
2358	quotring=prepareQuotientring(nvars(basering)-restindep[@m][3]);
2359
2360	//------------------------------------------------------------------
2361	//we pass to the quotientring K(var(nnp+1),..,var(nva))[..rest..]
2362	//------------------------------------------------------------------
2363
2364	execute quotring;
2365
2366	// @j considered in the quotientring
2367	ideal @j=imap(gnir1,@j);
2368	ideal ser=imap(gnir1,ser);
2369
2370	kill gnir1;
2371
2372	//j is a standardbasis in the quotientring but usually not minimal
2373	//here it becomes minimal
2374	@j=clearSB(@j,fett);
2375	attrib(@j,"isSB",1);
2376
2377	//we need later ggt(h[1],...)=gh for saturation
2378	ideal @h;
2379
2380	for(@n=1;@n<=size(@j);@n++)
2381	{
2382	@h[@n]=leadcoef(@j[@n]);
2383	}
2384	//the primary decomposition of j*K(var(nnp+1),..,var(nva))[..rest..]
2385
2386	option(redSB);
2387	list uprimary= zero_decomp(@j,ser,@wr);
2388	option(noredSB);
2389
2390
2391	//we need the intersection of the ideals in the list quprimary with
2392	//the polynomialring, i.e. let q=(f1,...,fr) in the quotientring
2393	//such an ideal but fi polynomials, then the intersection of q with
2394	//the polynomialring is the saturation of the ideal generated by
2395	//f1,...,fr with respect toh which is the lcm of the leading
2396	//coefficients of the fi considered in the quotientring:
2397	//this is coded in saturn
2398
2399	list saturn;
2400	ideal hpl;
2401
2402	for(@n=1;@n<=size(uprimary);@n++)
2403	{
2404	hpl=0;
2405	for(@n1=1;@n1<=size(uprimary[@n]);@n1++)
2406	{
2407	hpl=hpl,leadcoef(uprimary[@n][@n1]);
2408	}
2409	saturn[@n]=hpl;
2410	}
2411	//------------------------------------------------------------------
2412	//we leave the quotientring K(var(nnp+1),..,var(nva))[..rest..]
2413	//back to the polynomialring
2414	//------------------------------------------------------------------
2415	setring gnir;
2416
2417	collectprimary=imap(quring,uprimary);
2418	lsau=imap(quring,saturn);
2419	@h=imap(quring,@h);
2420
2421	kill quring;
2422
2423
2424	@n2=size(quprimary);
2425	@n3=@n2;
2426
2427	for(@n1=1;@n1<=size(collectprimary)/2;@n1++)
2428	{
2429	if(deg(collectprimary[2*@n1][1])>0)
2430	{
2431	@n2++;
2432	quprimary[@n2]=collectprimary[2*@n1-1];
2433	lnew[@n2]=lsau[2*@n1-1];
2434	@n2++;
2435	lnew[@n2]=lsau[2*@n1];
2436	quprimary[@n2]=collectprimary[2*@n1];
2437	}
2438	}
2439
2440
2441	//here the intersection with the polynomialring
2442	//mentioned above is really computed
2443
2444	for(@n=@n3/2+1;@n<=@n2/2;@n++)
2445	{
2446	if(specialIdealsEqual(quprimary[2@n-1],quprimary[2@n]))
2447	{
2448	quprimary[2@n-1]=sat2(quprimary[2@n-1],lnew[2*@n-1])[1];
2449	quprimary[2@n]=quprimary[2@n-1];
2450	}
2451	else
2452	{
2453	if(@wr==0)
2454	{
2455	quprimary[2@n-1]=sat2(quprimary[2@n-1],lnew[2*@n-1])[1];
2456	}
2457	quprimary[2@n]=sat2(quprimary[2@n],lnew[2*@n])[1];
2458	}
2459	}
2460	if(@n2>=@n3+2)
2461	{
2462	setring @Phelp;
2463	ser=imap(gnir,ser);
2464	hquprimary=imap(gnir,quprimary);
2465	for(@n=@n3/2+1;@n<=@n2/2;@n++)
2466	{
2467	if(@wr==0)
2468	{
2469	ser=intersect(ser,hquprimary[2*@n-1]);
2470	}
2471	else
2472	{
2473	ser=intersect(ser,hquprimary[2*@n]);
2474	}
2475	}
2476	setring gnir;
2477	ser=imap(@Phelp,ser);
2478	}
2479
2480	// }
2481	}
2482	if(size(reduce(ser,peek,1))!=0)
2483	{
2484	if(@wr>0)
2485	{
2486	htprimary=decomp(@j,@wr,peek,ser);
2487	}
2488	else
2489	{
2490	htprimary=decomp(@j,peek,ser);
2491	}
2492	// here we collect now both results primary(sat(j,gh))
2493	// and primary(j,gh^n)
2494	@n=size(quprimary);
2495	for (@k=1;@k<=size(htprimary);@k++)
2496	{
2497	quprimary[@n+@k]=htprimary[@k];
2498	}
2499	}
2500	}
2501
2502	}
2503	//---------------------------------------------------------------------------
2504	//back to the ring we started with
2505	//the final result: primary
2506	//---------------------------------------------------------------------------
2507
2508	setring @P;
2509	primary=imap(gnir,quprimary);
2510	return(primary);
2511	}
2512
2513
2514	example
2515	{ "EXAMPLE:"; echo = 2;
2516	ring r = 32003,(x,y,z),lp;
2517	poly p = z2+1;
2518	poly q = z4+2;
2519	ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
2520	list pr= decomp(i);
2521	pr;
2522	testPrimary( pr, i);
2523	}
2524
2525	///////////////////////////////////////////////////////////////////////////////
2526	proc radicalKL (list m,ideal ser,list #)
2527	{
2528	ideal i=m[2];
2529	//--------------------------------------------------------------------------
2530	//i is the zero-ideal
2531	//-------------------------------------------------------------------------
2532
2533	if(size(i)==0)
2534	{
2535	return(ideal(0));
2536	}
2537
2538	def @P = basering;
2539	list indep,allindep,restindep,fett,@mu;
2540	intvec @vh,isat;
2541	int @wr,@k,@n,@m,@n1,@n2,@n3,lauf,di;
2542	ideal @j,@j1,fac,@h,collectrad,htrad,lsau;
2543	ideal rad=ideal(1);
2544	ideal te=ser;
2545
2546	poly @p,@q;
2547	string @va,quotring;
2548	int homo=homog(i);
2549	if((find(ordstr(basering),"w")!=0)\|\|(find(ordstr(basering),"W")!=0)\|\|(find(ordstr(basering),"a")!=0))
2550	{
2551	homo=0;
2552	}
2553	if(size(#)>0)
2554	{
2555	@wr=#[1];
2556	}
2557	@j=m[1];
2558	@j1=m[2];
2559	int jdim=dim(@j);
2560	if(size(reduce(ser,@j,1))==0)
2561	{
2562	return(ser);
2563	}
2564	if(homo==1)
2565	{
2566	if(jdim==0)
2567	{
2568	option(noredSB);
2569	return(maxideal(1));
2570	}
2571	intvec @hilb=hilb(@j,1);
2572	}
2573
2574
2575	//---------------------------------------------------------------------------
2576	//j is the ring
2577	//---------------------------------------------------------------------------
2578
2579	if (jdim==-1)
2580	{
2581	option(noredSB);
2582	return(ideal(0));
2583	}
2584
2585	//---------------------------------------------------------------------------
2586	// the case of one variable
2587	//---------------------------------------------------------------------------
2588
2589	if(nvars(basering)==1)
2590	{
2591	fac=factorize(@j[1],1);
2592	@p=1;
2593	for(@k=1;@k<=size(fac);@k++)
2594	{
2595	@p=@p*fac[@k];
2596	}
2597	option(noredSB);
2598
2599	return(ideal(@p));
2600	}
2601	//---------------------------------------------------------------------------
2602	//the case of a complete intersection
2603	//---------------------------------------------------------------------------
2604	if(jdim+size(@j1)==nvars(basering))
2605	{
2606	// ideal jac=minor(jacob(@j1),size(@j1));
2607	// return(quotient(@j1,jac));
2608	}
2609
2610	//---------------------------------------------------------------------------
2611	//the zero-dimensional case
2612	//---------------------------------------------------------------------------
2613
2614	if (jdim==0)
2615	{
2616	@j1=finduni(@j);
2617	for(@k=1;@k<=size(@j1);@k++)
2618	{
2619	fac=factorize(cleardenom(@j1[@k]),1);
2620	@p=fac[1];
2621	for(@n=2;@n<=size(fac);@n++)
2622	{
2623	@p=@p*fac[@n];
2624	}
2625	@j=@j,@p;
2626	}
2627	@j=std(@j);
2628	option(noredSB);
2629	return(@j);
2630	}
2631
2632	//-------------------------------------------------------------------------
2633	//search for a maximal independent set indep,i.e.
2634	//look for subring such that the intersection with the ideal is zero
2635	//j intersected with K[var(indep[3]+1),...,var(nvar] is zero,
2636	//indep[1] is the new varstring, indep[2] the string for the block-ordering
2637	//-------------------------------------------------------------------------
2638
2639	indep=maxIndependSet(@j);
2640
2641	di=dim(@j);
2642
2643	for(@m=1;@m<=size(indep);@m++)
2644	{
2645	if((indep[@m][1]==varstr(basering))&&(@m==1))
2646	//this is the good case, nothing to do, just to have the same notations
2647	//change the ring
2648	{
2649	execute "ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
2650	+ordstr(basering)+");";
2651	ideal @j=fetch(@P,@j);
2652	attrib(@j,"isSB",1);
2653	}
2654	else
2655	{
2656	@va=string(maxideal(1));
2657	execute "ring gnir1 = ("+charstr(basering)+"),("+indep[@m][1]+"),("
2658	+indep[@m][2]+");";
2659	execute "map phi=@P,"+@va+";";
2660	if(homo==1)
2661	{
2662	ideal @j=std(phi(@j),@hilb);
2663	}
2664	else
2665	{
2666	ideal @j=groebner(phi(@j));
2667	}
2668	}
2669	if((deg(@j[1])==0)\|\|(dim(@j)<jdim))
2670	{
2671	setring @P;
2672	break;
2673	}
2674	for (lauf=1;lauf<=size(@j);lauf++)
2675	{
2676	fett[lauf]=size(@j[lauf]);
2677	}
2678	//------------------------------------------------------------------------
2679	//we have now the following situation:
2680	//j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
2681	//to this quotientring, j is their still a standardbasis, the
2682	//leading coefficients of the polynomials there (polynomials in
2683	//K[var(nnp+1),..,var(nva)]) are collected in the list h,
2684	//we need their ggt, gh, because of the following:
2685	//let (j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..rest..]
2686	//intersected with K[var(1),...,var(nva)] is (j:gh^n)
2687	//on the other hand j=(j,gh^n) intersected with (j:gh^n)
2688
2689	//------------------------------------------------------------------------
2690
2691	//the arrangement for the quotientring K(var(nnp+1),..,var(nva))[..rest..]
2692	//and the map phi:K[var(1),...,var(nva)] ----->
2693	//K(var(nnpr+1),..,var(nva))[..the rest..]
2694	//------------------------------------------------------------------------
2695	quotring=prepareQuotientring(nvars(basering)-indep[@m][3]);
2696
2697	//------------------------------------------------------------------------
2698	//we pass to the quotientring K(var(nnp+1),..,var(nva))[..the rest..]
2699	//------------------------------------------------------------------------
2700
2701	execute quotring;
2702
2703	// @j considered in the quotientring
2704	ideal @j=imap(gnir1,@j);
2705
2706	kill gnir1;
2707
2708	//j is a standardbasis in the quotientring but usually not minimal
2709	//here it becomes minimal
2710
2711	@j=clearSB(@j,fett);
2712	attrib(@j,"isSB",1);
2713
2714	//we need later ggt(h[1],...)=gh for saturation
2715	ideal @h;
2716	if(deg(@j[1])>0)
2717	{
2718	for(@n=1;@n<=size(@j);@n++)
2719	{
2720	@h[@n]=leadcoef(@j[@n]);
2721	}
2722	//the primary decomposition of j*K(var(nnp+1),..,var(nva))[..rest..]
2723	option(redSB);
2724
2725	@j=interred(@j);
2726
2727	attrib(@j,"isSB",1);
2728	list @mo=@j,@j;
2729	ideal zero_rad= radicalKL(@mo,ideal(1));
2730	}
2731	else
2732	{
2733	ideal zero_rad=ideal(1);
2734	}
2735
2736	//we need the intersection of the ideals in the list quprimary with the
2737	//polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
2738	//but fi polynomials, then the intersection of q with the polynomialring
2739	//is the saturation of the ideal generated by f1,...,fr with respect to
2740	//h which is the lcm of the leading coefficients of the fi considered in
2741	//the quotientring: this is coded in saturn
2742
2743	ideal hpl;
2744
2745	for(@n=1;@n<=size(zero_rad);@n++)
2746	{
2747	hpl=hpl,leadcoef(zero_rad[@n]);
2748	}
2749
2750	//------------------------------------------------------------------------
2751	//we leave the quotientring K(var(nnp+1),..,var(nva))[..the rest..]
2752	//back to the polynomialring
2753	//------------------------------------------------------------------------
2754	setring @P;
2755
2756	collectrad=imap(quring,zero_rad);
2757	lsau=simplify(imap(quring,hpl),2);
2758	@h=imap(quring,@h);
2759
2760	kill quring;
2761
2762
2763	//here the intersection with the polynomialring
2764	//mentioned above is really computed
2765
2766	collectrad=sat2(collectrad,lsau)[1];
2767
2768	if(deg(@h[1])>=0)
2769	{
2770	fac=ideal(0);
2771	for(lauf=1;lauf<=ncols(@h);lauf++)
2772	{
2773	if(deg(@h[lauf])>0)
2774	{
2775	fac=fac+factorize(@h[lauf],1);
2776	}
2777	}
2778	fac=simplify(fac,4);
2779	@q=1;
2780	for(lauf=1;lauf<=size(fac);lauf++)
2781	{
2782	@q=@q*fac[lauf];
2783	}
2784
2785
2786	@mu=mstd(quotient(@j+ideal(@q),rad));
2787	@j=@mu[1];
2788	attrib(@j,"isSB",1);
2789
2790	}
2791	if((deg(rad[1])>0)&&(deg(collectrad[1])>0))
2792	{
2793	rad=intersect(rad,collectrad);
2794	}
2795	else
2796	{
2797	if(deg(collectrad[1])>0)
2798	{
2799	rad=collectrad;
2800	}
2801	}
2802
2803	te=simplify(reduce(te*rad,@j),2);
2804
2805	if((dim(@j)<di)\|\|(size(te)==0))
2806	{
2807	break;
2808	}
2809	if(homo==1)
2810	{
2811	@hilb=hilb(@j,1);
2812	}
2813	}
2814
2815	if(((@wr==1)&&(dim(@j)<di))\|\|(deg(@j[1])==0)\|\|(size(te)==0))
2816	{
2817	return(rad);
2818	}
2819	// rad=intersect(rad,radicalKL(@mu,rad,@wr));
2820	rad=intersect(rad,radicalKL(@mu,ideal(1),@wr));
2821
2822
2823	option(noredSB);
2824	return(rad);
2825	}
2826
2827	///////////////////////////////////////////////////////////////////////////////
2828
2829	proc radicalEHV(ideal i,ideal re,list #)
2830	{
2831	ideal J,I,I0,radI0,L,radI1,I2,radI2;
2832	int l,il;
2833	if(size(#)>0)
2834	{
2835	il=#[1];
2836	}
2837
2838	option(redSB);
2839	list m=mstd(i);
2840	I=m[2];
2841	option(noredSB);
2842	if(size(reduce(re,m[1],1))==0)
2843	{
2844	return(re);
2845	}
2846	int cod=nvars(basering)-dim(m[1]);
2847	if((nvars(basering)<=5)&&(size(m[2])<=5))
2848	{
2849	if(cod==size(m[2]))
2850	{
2851	J=minor(jacob(I),cod);
2852	return(quotient(I,J));
2853	}
2854
2855	for(l=1;l<=cod;l++)
2856	{
2857	I0[l]=I[l];
2858	}
2859	if(dim(std(I0))+cod==nvars(basering))
2860	{
2861	J=minor(jacob(I0),cod);
2862	radI0=quotient(I0,J);
2863	L=quotient(radI0,I);
2864	radI1=quotient(radI0,L);
2865
2866	if(size(reduce(radI1,m[1],1))==0)
2867	{
2868	return(I);
2869	}
2870	if(il==1)
2871	{
2872
2873	return(radI1);
2874	}
2875
2876	I2=sat(I,radI1)[1];
2877
2878	if(deg(I2[1])<=0)
2879	{
2880	return(radI1);
2881	}
2882	return(intersect(radI1,radicalEHV(I2,re,il)));
2883	}
2884	}
2885	return(radicalKL(m,re,il));
2886	}
2887	///////////////////////////////////////////////////////////////////////////////
2888
2889	proc Ann(module M)
2890	{
2891	M=prune(M); //to obtain a small embedding
2892	ideal ann=quotient1(M,freemodule(nrows(M)));
2893	return(ann);
2894	}
2895	///////////////////////////////////////////////////////////////////////////////
2896
2897	//computes the equidimensional part of the ideal i of codimension e
2898	proc int_ass_primary_e(ideal i, int e)
2899	{
2900	if(homog(i)!=1)
2901	{
2902	i=std(i);
2903	}
2904	list re=sres(i,0); //the resolution
2905	re=minres(re); //minimized resolution
2906	ideal ann=AnnExt_R(e,re);
2907	if(nvars(basering)-dim(std(ann))!=e)
2908	{
2909	return(ideal(1));
2910	}
2911	return(ann);
2912	}
2913
2914	///////////////////////////////////////////////////////////////////////////////
2915
2916	//computes the annihilator of Ext^n(R/i,R) with given resolution re
2917	//n is not necessarily the number of variables
2918	proc AnnExt_R(int n,list re)
2919	{
2920	if(n<nvars(basering))
2921	{
2922	matrix f=transpose(re[n+1]); //Hom(_,R)
2923	module k=nres(f,2)[2]; //the kernel
2924	matrix g=transpose(re[n]); //the image of Hom(_,R)
2925
2926	ideal ann=quotient1(g,k); //the anihilator
2927	}
2928	else
2929	{
2930	ideal ann=Ann(transpose(re[n]));
2931	}
2932	return(ann);
2933	}
2934	///////////////////////////////////////////////////////////////////////////////
2935
2936	proc analyze(list pr)
2937	{
2938	int ii,jj;
2939	for(ii=1;ii<=size(pr)/2;ii++)
2940	{
2941	dim(std(pr[2*ii]));
2942	idealsEqual(pr[2ii-1],pr[2ii]);
2943	"===========================";
2944	}
2945
2946	for(ii=size(pr)/2;ii>1;ii--)
2947	{
2948	for(jj=1;jj<ii;jj++)
2949	{
2950	if(size(reduce(pr[2jj],std(pr[2ii],1)))==0)
2951	{
2952	"eingebette Komponente";
2953	jj;
2954	ii;
2955	}
2956	}
2957	}
2958	}
2959
2960	///////////////////////////////////////////////////////////////////////////////
2961	//
2962	// Shimoyama-Yokoyama
2963	//
2964	///////////////////////////////////////////////////////////////////////////////
2965
2966	proc simplifyIdeal(ideal i)
2967	{
2968	def r=basering;
2969
2970	int j,k;
2971	map phi;
2972	poly p;
2973
2974	ideal iwork=i;
2975	ideal imap1=maxideal(1);
2976	ideal imap2=maxideal(1);
2977
2978
2979	for(j=1;j<=nvars(basering);j++)
2980	{
2981	for(k=1;k<=size(i);k++)
2982	{
2983	if(deg(iwork[k]/var(j))==0)
2984	{
2985	p=-1/leadcoef(iwork[k]/var(j))*iwork[k];
2986	imap1[j]=p+2*var(j);
2987	phi=r,imap1;
2988	iwork=phi(iwork);
2989	iwork=subst(iwork,var(j),0);
2990	iwork[k]=var(j);
2991	imap1=maxideal(1);
2992	imap2[j]=-p;
2993	break;
2994	}
2995	}
2996	}
2997	return(iwork,imap2);
2998	}
2999
3000
3001	///////////////////////////////////////////////////////
3002	// ini_mod
3003	// input: a polynomial p
3004	// output: the initial term of p as needed
3005	// in the context of characteristic sets
3006	//////////////////////////////////////////////////////
3007
3008	proc ini_mod(poly p)
3009	{
3010	if (p==0)
3011	{
3012	return(0);
3013	}
3014	int n; matrix m;
3015	for( n=nvars(basering); n>0; n=n-1)
3016	{
3017	m=coef(p,var(n));
3018	if(m[1,1]!=1)
3019	{
3020	p=m[2,1];
3021	break;
3022	}
3023	}
3024	if(deg(p)==0)
3025	{
3026	p=0;
3027	}
3028	return(p);
3029	}
3030	///////////////////////////////////////////////////////
3031	// min_ass_prim_charsets
3032	// input: generators of an ideal PS and an integer cho
3033	// If cho=0, the given ordering of the variables is used.
3034	// Otherwise, the system tries to find an "optimal ordering",
3035	// which in some cases may considerably speed up the algorithm
3036	// output: the minimal associated primes of PS
3037	// algorithm: via characteriostic sets
3038	//////////////////////////////////////////////////////
3039
3040
3041	proc min_ass_prim_charsets (ideal PS, int cho)
3042	{
3043	if((cho<0) and (cho>1))
3044	{
3045	"ERROR: <int> must be 0 or 1"
3046	return();
3047	}
3048	if(system("version")>933)
3049	{
3050	option(notWarnSB);
3051	}
3052	if(cho==0)
3053	{
3054	return(min_ass_prim_charsets0(PS));
3055	}
3056	else
3057	{
3058	return(min_ass_prim_charsets1(PS));
3059	}
3060	}
3061	///////////////////////////////////////////////////////
3062	// min_ass_prim_charsets0
3063	// input: generators of an ideal PS
3064	// output: the minimal associated primes of PS
3065	// algorithm: via characteristic sets
3066	// the given ordering of the variables is used
3067	//////////////////////////////////////////////////////
3068
3069
3070	proc min_ass_prim_charsets0 (ideal PS)
3071	{
3072
3073	matrix m=char_series(PS); // We compute an irreducible
3074	// characteristic series
3075	int i,j,k;
3076	list PSI;
3077	list PHI; // the ideals given by the characteristic series
3078	for(i=nrows(m);i>=1; i--)
3079	{
3080	PHI[i]=ideal(m[i,1..ncols(m)]);
3081	}
3082	// We compute the radical of each ideal in PHI
3083	ideal I,JS,II;
3084	int sizeJS, sizeII;
3085	for(i=size(PHI);i>=1; i--)
3086	{
3087	I=0;
3088	for(j=size(PHI[i]);j>0;j=j-1)
3089	{
3090	I=I+ini_mod(PHI[i][j]);
3091	}
3092	JS=std(PHI[i]);
3093	sizeJS=size(JS);
3094	for(j=size(I);j>0;j=j-1)
3095	{
3096	II=0;
3097	sizeII=0;
3098	k=0;
3099	while(k<=sizeII) // successive saturation
3100	{
3101	option(returnSB);
3102	II=quotient(JS,I[j]);
3103	option(noreturnSB);
3104	//std
3105	// II=std(II);
3106	sizeII=size(II);
3107	if(sizeII==sizeJS)
3108	{
3109	for(k=1;k<=sizeII;k++)
3110	{
3111	if(leadexp(II[k])!=leadexp(JS[k])) break;
3112	}
3113	}
3114	JS=II;
3115	sizeJS=sizeII;
3116	}
3117	}
3118	PSI=insert(PSI,JS);
3119	}
3120	int sizePSI=size(PSI);
3121	// We eliminate redundant ideals
3122	for(i=1;i<sizePSI;i++)
3123	{
3124	for(j=i+1;j<=sizePSI;j++)
3125	{
3126	if(size(PSI[i])!=0)
3127	{
3128	if(size(PSI[j])!=0)
3129	{
3130	if(size(NF(PSI[i],PSI[j],1))==0)
3131	{
3132	PSI[j]=ideal(0);
3133	}
3134	else
3135	{
3136	if(size(NF(PSI[j],PSI[i],1))==0)
3137	{
3138	PSI[i]=ideal(0);
3139	}
3140	}
3141	}
3142	}
3143	}
3144	}
3145	for(i=sizePSI;i>=1;i--)
3146	{
3147	if(size(PSI[i])==0)
3148	{
3149	PSI=delete(PSI,i);
3150	}
3151	}
3152	return (PSI);
3153	}
3154
3155	///////////////////////////////////////////////////////
3156	// min_ass_prim_charsets1
3157	// input: generators of an ideal PS
3158	// output: the minimal associated primes of PS
3159	// algorithm: via characteristic sets
3160	// input: generators of an ideal PS and an integer i
3161	// The system tries to find an "optimal ordering" of
3162	// the variables
3163	//////////////////////////////////////////////////////
3164
3165
3166	proc min_ass_prim_charsets1 (ideal PS)
3167	{
3168	def oldring=basering;
3169	string n=system("neworder",PS);
3170	execute "ring r=("+charstr(oldring)+"),("+n+"),dp;";
3171	ideal PS=imap(oldring,PS);
3172	matrix m=char_series(PS); // We compute an irreducible
3173	// characteristic series
3174	int i,j,k;
3175	ideal I;
3176	list PSI;
3177	list PHI; // the ideals given by the characteristic series
3178	list ITPHI; // their initial terms
3179	for(i=nrows(m);i>=1; i--)
3180	{
3181	PHI[i]=ideal(m[i,1..ncols(m)]);
3182	I=0;
3183	for(j=size(PHI[i]);j>0;j=j-1)
3184	{
3185	I=I,ini_mod(PHI[i][j]);
3186	}
3187	I=I[2..ncols(I)];
3188	ITPHI[i]=I;
3189	}
3190	setring oldring;
3191	matrix m=imap(r,m);
3192	list PHI=imap(r,PHI);
3193	list ITPHI=imap(r,ITPHI);
3194	// We compute the radical of each ideal in PHI
3195	ideal I,JS,II;
3196	int sizeJS, sizeII;
3197	for(i=size(PHI);i>=1; i--)
3198	{
3199	I=0;
3200	for(j=size(PHI[i]);j>0;j=j-1)
3201	{
3202	I=I+ITPHI[i][j];
3203	}
3204	JS=std(PHI[i]);
3205	sizeJS=size(JS);
3206	for(j=size(I);j>0;j=j-1)
3207	{
3208	II=0;
3209	sizeII=0;
3210	k=0;
3211	while(k<=sizeII) // successive iteration
3212	{
3213	option(returnSB);
3214	II=quotient(JS,I[j]);
3215	option(noreturnSB);
3216	//std
3217	// II=std(II);
3218	sizeII=size(II);
3219	if(sizeII==sizeJS)
3220	{
3221	for(k=1;k<=sizeII;k++)
3222	{
3223	if(leadexp(II[k])!=leadexp(JS[k])) break;
3224	}
3225	}
3226	JS=II;
3227	sizeJS=sizeII;
3228	}
3229	}
3230	PSI=insert(PSI,JS);
3231	}
3232	int sizePSI=size(PSI);
3233	// We eliminate redundant ideals
3234	for(i=1;i<sizePSI;i++)
3235	{
3236	for(j=i+1;j<=sizePSI;j++)
3237	{
3238	if(size(PSI[i])!=0)
3239	{
3240	if(size(PSI[j])!=0)
3241	{
3242	if(size(NF(PSI[i],PSI[j],1))==0)
3243	{
3244	PSI[j]=ideal(0);
3245	}
3246	else
3247	{
3248	if(size(NF(PSI[j],PSI[i],1))==0)
3249	{
3250	PSI[i]=ideal(0);
3251	}
3252	}
3253	}
3254	}
3255	}
3256	}
3257	for(i=sizePSI;i>=1;i--)
3258	{
3259	if(size(PSI[i])==0)
3260	{
3261	PSI=delete(PSI,i);
3262	}
3263	}
3264	return (PSI);
3265	}
3266
3267
3268	/////////////////////////////////////////////////////
3269	// proc prim_dec
3270	// input: generators of an ideal I and an integer choose
3271	// If choose=0, min_ass_prim_charsets with the given
3272	// ordering of the variables is used.
3273	// If choose=1, min_ass_prim_charsets with the "optimized"
3274	// ordering of the variables is used.
3275	// If choose=2, minAssPrimes from primdec.lib is used
3276	// If choose=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
3277	// output: a primary decomposition of I, i.e., a list
3278	// of pairs consisting of a standard basis of a primary component
3279	// of I and a standard basis of the corresponding associated prime.
3280	// To compute the minimal associated primes of a given ideal
3281	// min_ass_prim_l is called, i.e., the minimal associated primes
3282	// are computed via characteristic sets.
3283	// In the homogeneous case, the performance of the procedure
3284	// will be improved if I is already given by a minimal set of
3285	// generators. Apply minbase if necessary.
3286	//////////////////////////////////////////////////////////
3287
3288
3289	proc prim_dec(ideal I, int choose)
3290	{
3291	if((choose<0) or (choose>3))
3292	{
3293	"ERROR: <int> must be 0 or 1 or 2 or 3";
3294	return();
3295	}
3296	if(system("version")>933)
3297	{
3298	option(notWarnSB);
3299	}
3300	ideal H=1; // The intersection of the primary components
3301	list U; // the leaves of the decomposition tree, i.e.,
3302	// pairs consisting of a primary component of I
3303	// and the corresponding associated prime
3304	list W; // the non-leaf vertices in the decomposition tree.
3305	// every entry has 6 components:
3306	// 1- the vertex itself , i.e., a standard bais of the
3307	// given ideal I (type 1), or a standard basis of a
3308	// pseudo-primary component arising from
3309	// pseudo-primary decomposition (type 2), or a
3310	// standard basis of a remaining component arising from
3311	// pseudo-primary decomposition or extraction (type 3)
3312	// 2- the type of the vertex as indicated above
3313	// 3- the weighted_tree_depth of the vertex
3314	// 4- the tester of the vertex
3315	// 5- a standard basis of the associated prime
3316	// of a vertex of type 2, or 0 otherwise
3317	// 6- a list of pairs consisting of a standard
3318	// basis of a minimal associated prime ideal
3319	// of the father of the vertex and the
3320	// irreducible factors of the "minimal
3321	// divisor" of the seperator or extractor
3322	// corresponding to the prime ideal
3323	// as computed by the procedure minsat,
3324	// if the vertex is of type 3, or
3325	// the empty list otherwise
3326	ideal SI=std(I);
3327	if(SI[1]==1) // primdecSY(ideal(1))
3328	{
3329	return(list());
3330	}
3331	int ncolsSI=ncols(SI);
3332	int ncolsH=1;
3333	W[1]=list(I,1,0,poly(1),ideal(0),list()); // The root of the tree
3334	int weighted_tree_depth;
3335	int i,j;
3336	int check;
3337	list V; // current vertex
3338	list VV; // new vertex
3339	list QQ;
3340	list WI;
3341	ideal Qi,SQ,SRest,fac;
3342	poly tester;
3343
3344	while(1)
3345	{
3346	i=1;
3347	while(1)
3348	{
3349	while(i<=size(W)) // find vertex V of smallest weighted tree-depth
3350	{
3351	if (W[i][3]<=weighted_tree_depth) break;
3352	i++;
3353	}
3354	if (i<=size(W)) break;
3355	i=1;
3356	weighted_tree_depth++;
3357	}
3358	V=W[i];
3359	W=delete(W,i); // delete V from W
3360
3361	// now proceed by type of vertex V
3362
3363	if (V[2]==2) // extraction needed
3364	{
3365	SQ,SRest,fac=extraction(V[1],V[5]);
3366	// standard basis of primary component,
3367	// standard basis of remaining component,
3368	// irreducible factors of
3369	// the "minimal divisor" of the extractor
3370	// as computed by the procedure minsat,
3371	check=0;
3372	for(j=1;j<=ncolsH;j++)
3373	{
3374	if (NF(H[j],SQ,1)!=0) // Q is not redundant
3375	{
3376	check=1;
3377	break;
3378	}
3379	}
3380	if(check==1) // Q is not redundant
3381	{
3382	QQ=list();
3383	QQ[1]=list(SQ,V[5]); // primary component, associated prime,
3384	// i.e., standard bases thereof
3385	U=U+QQ;
3386	H=intersect(H,SQ);
3387	H=std(H);
3388	ncolsH=ncols(H);
3389	check=0;
3390	if(ncolsH==ncolsSI)
3391	{
3392	for(j=1;j<=ncolsSI;j++)
3393	{
3394	if(leadexp(H[j])!=leadexp(SI[j]))
3395	{
3396	check=1;
3397	break;
3398	}
3399	}
3400	}
3401	else
3402	{
3403	check=1;
3404	}
3405	if(check==0) // H==I => U is a primary decomposition
3406	{
3407	return(U);
3408	}
3409	}
3410	if (SRest[1]!=1) // the remaining component is not
3411	// the whole ring
3412	{
3413	if (rad_con(V[4],SRest)==0) // the new vertex is not the
3414	// root of a redundant subtree
3415	{
3416	VV[1]=SRest; // remaining component
3417	VV[2]=3; // pseudoprimdec_special
3418	VV[3]=V[3]+1; // weighted depth
3419	VV[4]=V[4]; // the tester did not change
3420	VV[5]=ideal(0);
3421	VV[6]=list(list(V[5],fac));
3422	W=insert(W,VV,size(W));
3423	}
3424	}
3425	}
3426	else
3427	{
3428	if (V[2]==3) // pseudo_prim_dec_special is needed
3429	{
3430	QQ,SRest=pseudo_prim_dec_special_charsets(V[1],V[6],choose);
3431	// QQ = quadruples:
3432	// standard basis of pseudo-primary component,
3433	// standard basis of corresponding prime,
3434	// seperator, irreducible factors of
3435	// the "minimal divisor" of the seperator
3436	// as computed by the procedure minsat,
3437	// SRest=standard basis of remaining component
3438	}
3439	else // V is the root, pseudo_prim_dec is needed
3440	{
3441	QQ,SRest=pseudo_prim_dec_charsets(I,SI,choose);
3442	// QQ = quadruples:
3443	// standard basis of pseudo-primary component,
3444	// standard basis of corresponding prime,
3445	// seperator, irreducible factors of
3446	// the "minimal divisor" of the seperator
3447	// as computed by the procedure minsat,
3448	// SRest=standard basis of remaining component
3449
3450	}
3451	//check
3452	for(i=size(QQ);i>=1;i--)
3453	//for(i=1;i<=size(QQ);i++)
3454	{
3455	tester=QQ[i][3]*V[4];
3456	Qi=QQ[i][2];
3457	if(NF(tester,Qi,1)!=0) // the new vertex is not the
3458	// root of a redundant subtree
3459	{
3460	VV[1]=QQ[i][1];
3461	VV[2]=2;
3462	VV[3]=V[3]+1;
3463	VV[4]=tester; // the new tester as computed above
3464	VV[5]=Qi; // QQ[i][2];
3465	VV[6]=list();
3466	W=insert(W,VV,size(W));
3467	}
3468	}
3469	if (SRest[1]!=1) // the remaining component is not
3470	// the whole ring
3471	{
3472	if (rad_con(V[4],SRest)==0) // the vertex is not the root
3473	// of a redundant subtree
3474	{
3475	VV[1]=SRest;
3476	VV[2]=3;
3477	VV[3]=V[3]+2;
3478	VV[4]=V[4]; // the tester did not change
3479	VV[5]=ideal(0);
3480	WI=list();
3481	for(i=1;i<=size(QQ);i++)
3482	{
3483	WI=insert(WI,list(QQ[i][2],QQ[i][4]));
3484	}
3485	VV[6]=WI;
3486	W=insert(W,VV,size(W));
3487	}
3488	}
3489	}
3490	}
3491	}
3492
3493	//////////////////////////////////////////////////////////////////////////
3494	// proc pseudo_prim_dec_charsets
3495	// input: Generators of an arbitrary ideal I, a standard basis SI of I,
3496	// and an integer choo
3497	// If choo=0, min_ass_prim_charsets with the given
3498	// ordering of the variables is used.
3499	// If choo=1, min_ass_prim_charsets with the "optimized"
3500	// ordering of the variables is used.
3501	// If choo=2, minAssPrimes from primdec.lib is used
3502	// If choo=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
3503	// output: a pseudo primary decomposition of I, i.e., a list
3504	// of pseudo primary components together with a standard basis of the
3505	// remaining component. Each pseudo primary component is
3506	// represented by a quadrupel: A standard basis of the component,
3507	// a standard basis of the corresponding associated prime, the
3508	// seperator of the component, and the irreducible factors of the
3509	// "minimal divisor" of the seperator as computed by the procedure minsat,
3510	// calls proc pseudo_prim_dec_i
3511	//////////////////////////////////////////////////////////////////////////
3512
3513
3514	proc pseudo_prim_dec_charsets (ideal I, ideal SI, int choo)
3515	{
3516	list L; // The list of minimal associated primes,
3517	// each one given by a standard basis
3518	if((choo==0) or (choo==1))
3519	{
3520	L=min_ass_prim_charsets(I,choo);
3521	}
3522	else
3523	{
3524	if(choo==2)
3525	{
3526	L=minAssPrimes(I);
3527	}
3528	else
3529	{
3530	L=minAssPrimes(I,1);
3531	}
3532	for(int i=size(L);i>=1;i=i-1)
3533	{
3534	L[i]=std(L[i]);
3535	}
3536	}
3537	return (pseudo_prim_dec_i(SI,L));
3538	}
3539
3540	////////////////////////////////////////////////////////////////
3541	// proc pseudo_prim_dec_special_charsets
3542	// input: a standard basis of an ideal I whose radical is the
3543	// intersection of the radicals of ideals generated by one prime ideal
3544	// P_i together with one polynomial f_i, the list V6 must be the list of
3545	// pairs (standard basis of P_i, irreducible factors of f_i),
3546	// and an integer choo
3547	// If choo=0, min_ass_prim_charsets with the given
3548	// ordering of the variables is used.
3549	// If choo=1, min_ass_prim_charsets with the "optimized"
3550	// ordering of the variables is used.
3551	// If choo=2, minAssPrimes from primdec.lib is used
3552	// If choo=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
3553	// output: a pseudo primary decomposition of I, i.e., a list
3554	// of pseudo primary components together with a standard basis of the
3555	// remaining component. Each pseudo primary component is
3556	// represented by a quadrupel: A standard basis of the component,
3557	// a standard basis of the corresponding associated prime, the
3558	// seperator of the component, and the irreducible factors of the
3559	// "minimal divisor" of the seperator as computed by the procedure minsat,
3560	// calls proc pseudo_prim_dec_i
3561	////////////////////////////////////////////////////////////////
3562
3563
3564	proc pseudo_prim_dec_special_charsets (ideal SI,list V6, int choo)
3565	{
3566	int i,j,l;
3567	list m;
3568	list L;
3569	int sizeL;
3570	ideal P,SP; ideal fac;
3571	int dimSP;
3572	for(l=size(V6);l>=1;l--) // creates a list of associated primes
3573	// of I, possibly redundant
3574	{
3575	P=V6[l][1];
3576	fac=V6[l][2];
3577	for(i=ncols(fac);i>=1;i--)
3578	{
3579	SP=P+fac[i];
3580	SP=std(SP);
3581	if(SP[1]!=1)
3582	{
3583	if((choo==0) or (choo==1))
3584	{
3585	m=min_ass_prim_charsets(SP,choo); // a list of SB
3586	}
3587	else
3588	{
3589	if(choo==2)
3590	{
3591	m=minAssPrimes(SP);
3592	}
3593	else
3594	{
3595	m=minAssPrimes(SP,1);
3596	}
3597	for(j=size(m);j>=1;j=j-1)
3598	{
3599	m[j]=std(m[j]);
3600	}
3601	}
3602	dimSP=dim(SP);
3603	for(j=size(m);j>=1; j--)
3604	{
3605	if(dim(m[j])==dimSP)
3606	{
3607	L=insert(L,m[j],size(L));
3608	}
3609	}
3610	}
3611	}
3612	}
3613	sizeL=size(L);
3614	for(i=1;i<sizeL;i++) // get rid of redundant primes
3615	{
3616	for(j=i+1;j<=sizeL;j++)
3617	{
3618	if(size(L[i])!=0)
3619	{
3620	if(size(L[j])!=0)
3621	{
3622	if(size(NF(L[i],L[j],1))==0)
3623	{
3624	L[j]=ideal(0);
3625	}
3626	else
3627	{
3628	if(size(NF(L[j],L[i],1))==0)
3629	{
3630	L[i]=ideal(0);
3631	}
3632	}
3633	}
3634	}
3635	}
3636	}
3637	for(i=sizeL;i>=1;i--)
3638	{
3639	if(size(L[i])==0)
3640	{
3641	L=delete(L,i);
3642	}
3643	}
3644	return (pseudo_prim_dec_i(SI,L));
3645	}
3646
3647
3648	////////////////////////////////////////////////////////////////
3649	// proc pseudo_prim_dec_i
3650	// input: A standard basis of an arbitrary ideal I, and standard bases
3651	// of the minimal associated primes of I
3652	// output: a pseudo primary decomposition of I, i.e., a list
3653	// of pseudo primary components together with a standard basis of the
3654	// remaining component. Each pseudo primary component is
3655	// represented by a quadrupel: A standard basis of the component Q_i,
3656	// a standard basis of the corresponding associated prime P_i, the
3657	// seperator of the component, and the irreducible factors of the
3658	// "minimal divisor" of the seperator as computed by the procedure minsat,
3659	////////////////////////////////////////////////////////////////
3660
3661
3662	proc pseudo_prim_dec_i (ideal SI, list L)
3663	{
3664	list Q;
3665	if (size(L)==1) // one minimal associated prime only
3666	// the ideal is already pseudo primary
3667	{
3668	Q=SI,L[1],1;
3669	list QQ;
3670	QQ[1]=Q;
3671	return (QQ,ideal(1));
3672	}
3673
3674	poly f0,f,g;
3675	ideal fac;
3676	int i,j,k,l;
3677	ideal SQi;
3678	ideal I'=SI;
3679	list QP;
3680	int sizeL=size(L);
3681	for(i=1;i<=sizeL;i++)
3682	{
3683	fac=0;
3684	for(j=1;j<=sizeL;j++) // compute the seperator sep_i
3685	// of the i-th component
3686	{
3687	if (i!=j) // search g not in L[i], but L[j]
3688	{
3689	for(k=1;k<=ncols(L[j]);k++)
3690	{
3691	if(NF(L[j][k],L[i],1)!=0)
3692	{
3693	break;
3694	}
3695	}
3696	fac=fac+L[j][k];
3697	}
3698	}
3699	// delete superfluous polynomials
3700	fac=simplify(fac,8);
3701	// saturation
3702	SQi,f0,f,fac=minsat_ppd(SI,fac);
3703	I'=I',f;
3704	QP=SQi,L[i],f0,fac;
3705	// the quadrupel:
3706	// a standard basis of Q_i,
3707	// a standard basis of P_i,
3708	// sep_i,
3709	// irreducible factors of
3710	// the "minimal divisor" of the seperator
3711	// as computed by the procedure minsat,
3712	Q[i]=QP;
3713	}
3714	I'=std(I');
3715	return (Q, I');
3716	// I' = remaining component
3717	}
3718
3719
3720	////////////////////////////////////////////////////////////////
3721	// proc extraction
3722	// input: A standard basis of a pseudo primary ideal I, and a standard
3723	// basis of the unique minimal associated prime P of I
3724	// output: an extraction of I, i.e., a standard basis of the primary
3725	// component Q of I with associated prime P, a standard basis of the
3726	// remaining component, and the irreducible factors of the
3727	// "minimal divisor" of the extractor as computed by the procedure minsat
3728	////////////////////////////////////////////////////////////////
3729
3730
3731	proc extraction (ideal SI, ideal SP)
3732	{
3733	list indsets=indepSet(SP,0);
3734	poly f;
3735	if(size(indsets)!=0) //check, whether dim P != 0
3736	{
3737	intvec v; // a maximal independent set of variables
3738	// modulo P
3739	string U; // the independent variables
3740	string A; // the dependent variables
3741	int j,k;
3742	int a; // the size of A
3743	int degf;
3744	ideal g;
3745	list polys;
3746	int sizepolys;
3747	list newpoly;
3748	def R=basering;
3749	//intvec hv=hilb(SI,1);
3750	for (k=1;k<=size(indsets);k++)
3751	{
3752	v=indsets[k];
3753	for (j=1;j<=nvars(R);j++)
3754	{
3755	if (v[j]==1)
3756	{
3757	U=U+varstr(j)+",";
3758	}
3759	else
3760	{
3761	A=A+varstr(j)+",";
3762	a++;
3763	}
3764	}
3765
3766	U[size(U)]=")"; // we compute the extractor of I (w.r.t. U)
3767	execute "ring RAU="+charstr(basering)+",("+A+U+",(dp("+string(a)+"),dp);";
3768	ideal I=imap(R,SI);
3769	//I=std(I,hv); // the standard basis in (R[U])[A]
3770	I=std(I); // the standard basis in (R[U])[A]
3771	A[size(A)]=")";
3772	execute "ring Rloc=("+charstr(basering)+","+U+",("+A+",dp;";
3773	ideal I=imap(RAU,I);
3774	//"std in lokalisierung:"+newline,I;
3775	ideal h;
3776	for(j=ncols(I);j>=1;j--)
3777	{
3778	h[j]=leadcoef(I[j]); // consider I in (R(U))[A]
3779	}
3780	setring R;
3781	g=imap(Rloc,h);
3782	kill RAU,Rloc;
3783	U="";
3784	A="";
3785	a=0;
3786	f=lcm(g);
3787	newpoly[1]=f;
3788	polys=polys+newpoly;
3789	newpoly=list();
3790	}
3791	f=polys[1];
3792	degf=deg(f);
3793	sizepolys=size(polys);
3794	for (k=2;k<=sizepolys;k++)
3795	{
3796	if (deg(polys[k])<degf)
3797	{
3798	f=polys[k];
3799	degf=deg(f);
3800	}
3801	}
3802	}
3803	else
3804	{
3805	f=1;
3806	}
3807	poly f0,h0; ideal SQ; ideal fac;
3808	if(f!=1)
3809	{
3810	SQ,f0,h0,fac=minsat(SI,f);
3811	return(SQ,std(SI+h0),fac);
3812	// the tripel
3813	// a standard basis of Q,
3814	// a standard basis of remaining component,
3815	// irreducible factors of
3816	// the "minimal divisor" of the extractor
3817	// as computed by the procedure minsat
3818	}
3819	else
3820	{
3821	return(SI,ideal(1),ideal(1));
3822	}
3823	}
3824
3825	/////////////////////////////////////////////////////
3826	// proc minsat
3827	// input: a standard basis of an ideal I and a polynomial p
3828	// output: a standard basis IS of the saturation of I w.r. to p,
3829	// the maximal squarefree factor f0 of p,
3830	// the "minimal divisor" f of f0 such that the saturation of
3831	// I w.r. to f equals the saturation of I w.r. to f0 (which is IS),
3832	// the irreducible factors of f
3833	//////////////////////////////////////////////////////////
3834
3835
3836	proc minsat(ideal SI, poly p)
3837	{
3838	ideal fac=factorize(p,1); //the irreducible factors of p
3839	fac=sort(fac)[1];
3840	int i,k;
3841	poly f0=1;
3842	for(i=ncols(fac);i>=1;i--)
3843	{
3844	f0=f0*fac[i];
3845	}
3846	poly f=1;
3847	ideal iold;
3848	list quotM;
3849	quotM[1]=SI;
3850	quotM[2]=fac;
3851	quotM[3]=f0;
3852	// we deal seperately with the first quotient;
3853	// factors, which do not contribute to this one,
3854	// are omitted
3855	iold=quotM[1];
3856	quotM=minquot(quotM);
3857	fac=quotM[2];
3858	if(quotM[3]==1)
3859	{
3860	return(quotM[1],f0,f,fac);
3861	}
3862	while(special_ideals_equal(iold,quotM[1])==0)
3863	{
3864	f=f*quotM[3];
3865	iold=quotM[1];
3866	quotM=minquot(quotM);
3867	}
3868	return(quotM[1],f0,f,fac); // the quadrupel ((I:p),f0,f, irr. factors of f)
3869	}
3870
3871	/////////////////////////////////////////////////////
3872	// proc minsat_ppd
3873	// input: a standard basis of an ideal I and a polynomial p
3874	// output: a standard basis IS of the saturation of I w.r. to p,
3875	// the maximal squarefree factor f0 of p,
3876	// the "minimal divisor" f of f0 such that the saturation of
3877	// I w.r. to f equals the saturation of I w.r. to f0 (which is IS),
3878	// the irreducible factors of f
3879	//////////////////////////////////////////////////////////
3880
3881
3882	proc minsat_ppd(ideal SI, ideal fac)
3883	{
3884	fac=sort(fac)[1];
3885	int i,k;
3886	poly f0=1;
3887	for(i=ncols(fac);i>=1;i--)
3888	{
3889	f0=f0*fac[i];
3890	}
3891	poly f=1;
3892	ideal iold;
3893	list quotM;
3894	quotM[1]=SI;
3895	quotM[2]=fac;
3896	quotM[3]=f0;
3897	// we deal seperately with the first quotient;
3898	// factors, which do not contribute to this one,
3899	// are omitted
3900	iold=quotM[1];
3901	quotM=minquot(quotM);
3902	fac=quotM[2];
3903	if(quotM[3]==1)
3904	{
3905	return(quotM[1],f0,f,fac);
3906	}
3907	while(special_ideals_equal(iold,quotM[1])==0)
3908	{
3909	f=f*quotM[3];
3910	iold=quotM[1];
3911	quotM=minquot(quotM);
3912	k++;
3913	}
3914	return(quotM[1],f0,f,fac); // the quadrupel ((I:p),f0,f, irr. factors of f)
3915	}
3916	/////////////////////////////////////////////////////////////////
3917	// proc minquot
3918	// input: a list with 3 components: a standard basis
3919	// of an ideal I, a set of irreducible polynomials, and
3920	// there product f0
3921	// output: a standard basis of the ideal (I:f0), the irreducible
3922	// factors of the "minimal divisor" f of f0 with (I:f0) = (I:f),
3923	// the "minimal divisor" f
3924	/////////////////////////////////////////////////////////////////
3925
3926	proc minquot(list tsil)
3927	{
3928	int i,j,k,action;
3929	ideal verg;
3930	list l;
3931	poly g;
3932	ideal laedi=tsil[1];
3933	ideal fac=tsil[2];
3934	poly f=tsil[3];
3935
3936	//std
3937	// ideal star=quotient(laedi,f);
3938	// star=std(star);
3939	option(returnSB);
3940	ideal star=quotient(laedi,f);
3941	option(noreturnSB);
3942	if(special_ideals_equal(laedi,star)==1)
3943	{
3944	return(laedi,ideal(1),1);
3945	}
3946	action=1;
3947	while(action==1)
3948	{
3949	if(size(fac)==1)
3950	{
3951	action=0;
3952	break;
3953	}
3954	for(i=1;i<=size(fac);i++)
3955	{
3956	g=1;
3957	for(j=1;j<=size(fac);j++)
3958	{
3959	if(i!=j)
3960	{
3961	g=g*fac[j];
3962	}
3963	}
3964	//std
3965	// verg=quotient(laedi,g);
3966	// verg=std(verg);
3967	option(returnSB);
3968	verg=quotient(laedi,g);
3969	option(noreturnSB);
3970	if(special_ideals_equal(verg,star)==1)
3971	{
3972	f=g;
3973	fac[i]=0;
3974	fac=simplify(fac,2);
3975	break;
3976	}
3977	if(i==size(fac))
3978	{
3979	action=0;
3980	}
3981	}
3982	}
3983	l=star,fac,f;
3984	return(l);
3985	}
3986	/////////////////////////////////////////////////
3987	// proc special_ideals_equal
3988	// input: standard bases of ideal k1 and k2 such that
3989	// k1 is contained in k2, or k2 is contained ink1
3990	// output: 1, if k1 equals k2, 0 otherwise
3991	//////////////////////////////////////////////////
3992
3993	proc special_ideals_equal( ideal k1, ideal k2)
3994	{
3995	int j;
3996	if(size(k1)==size(k2))
3997	{
3998	for(j=1;j<=size(k1);j++)
3999	{
4000	if(leadexp(k1[j])!=leadexp(k2[j]))
4001	{
4002	return(0);
4003	}
4004	}
4005	return(1);
4006	}
4007	return(0);
4008	}
4009
4010
4011	///////////////////////////////////////////////////////////////////////////////
4012
4013	proc convList(list l)
4014	{
4015	int i;
4016	list re,he;
4017	for(i=1;i<=size(l)/2;i++)
4018	{
4019	he=l[2i-1],l[2i];
4020	re[i]=he;
4021	}
4022	return(re);
4023	}
4024	///////////////////////////////////////////////////////////////////////////////
4025
4026	proc reconvList(list l)
4027	{
4028	int i;
4029	list re;
4030	for(i=1;i<=size(l);i++)
4031	{
4032	re[2*i-1]=l[i][1];
4033	re[2*i]=l[i][2];
4034	}
4035	return(re);
4036	}
4037
4038	///////////////////////////////////////////////////////////////////////////////
4039	//
4040	// The main procedures
4041	//
4042	///////////////////////////////////////////////////////////////////////////////
4043
4044	proc primdecGTZ(ideal i)
4045	"USAGE: primdecGTZ(i); i ideal
4046	RETURN: a list, say pr, of primary ideals and their associated primes
4047	pr[i][1], resp. pr[i][2] is the i-th primary resp. prime component
4048	NOTE: Algorithm of Gianni, Traeger, Zacharias
4049	designed for characteristic 0, works also in char k > 0, if it
4050	terminates (may result in an infinite loop in small characteristic!)
4051	EXAMPLE: example primdecGTZ; shows an example
4052	"
4053	{
4054	return(convList(decomp(i)));
4055	}
4056	example
4057	{ "EXAMPLE:"; echo = 2;
4058	ring r = 32003,(x,y,z),lp;
4059	poly p = z2+1;
4060	poly q = z4+2;
4061	ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
4062	list pr = primdecGTZ(i);
4063	pr;
4064	}
4065	///////////////////////////////////////////////////////////////////////////////
4066
4067	proc primdecSY(ideal i, list #))
4068	"USAGE: primdecSY(i); i ideal, c int
4069	if c=0, the given ordering of the variables is used.
4070	if c=1, minAssChar tries to use an optimal ordering,
4071	if c=2, minAssGTZ is used
4072	if c=3, minAssGTZ and facstd is used
4073	RETURN: a list, say pr, of primary ideals and their associated primes
4074	pr[i][1], resp. pr[i][2] is the i-th primary resp. prime component
4075	NOTE: Algorithm of Shimoyama-Yokoyama
4076	implemented for characteristic 0, works also in char k > 0,
4077	the result may be not completely decomposed in small characteristic
4078	EXAMPLE: example primdecSY; shows an example
4079	"
4080	{
4081	if (size(#)==1)
4082	{ return(prim_dec(i,#[1])); }
4083	else
4084	{ return(prim_dec(i,1)); }
4085	}
4086	example
4087	{ "EXAMPLE:"; echo = 2;
4088	ring r = 32003,(x,y,z),lp;
4089	poly p = z2+1;
4090	poly q = z4+2;
4091	ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
4092	list pr = primdecSY(i);
4093	pr;
4094	}
4095	///////////////////////////////////////////////////////////////////////////////
4096	proc minAssGTZ(ideal i)
4097	"USAGE: minAssGTZ(i); i ideal
4098	RETURN: list = the minimal associated prime ideals of i
4099	NOTE: designed for characteristic 0, works also in char k > 0 if it
4100	terminates, may result in an infinite loop in small characteristic
4101	EXAMPLE: example minAssGTZ; shows an example
4102	"
4103	{
4104	return(minAssPrimes(i,1));
4105	}
4106	example
4107	{ "EXAMPLE:"; echo = 2;
4108	ring r = 32003,(x,y,z),dp;
4109	poly p = z2+1;
4110	poly q = z4+2;
4111	ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
4112	list pr= minAssGTZ(i);
4113	pr;
4114	}
4115
4116	///////////////////////////////////////////////////////////////////////////////
4117	proc minAssChar(ideal i, list #)
4118	"USAGE: minAssChar(i[,c]); i ideal,
4119	if c=0, the given ordering of the variables is used.
4120	Otherwise, the system tries to find an optimal ordering,
4121	which in some cases may considerably speed up the algorithm
4122	RETURN: list = the minimal associated prime ideals of i
4123	NOTE: implemented for characteristic 0, works also in char k >> 0,
4124	the result may be not compltely decomposed in small characteristic
4125	EXAMPLE: example minAssChar; shows an example
4126	"
4127	{
4128	if (size(#)==1)
4129	{ return(min_ass_prim_charsets(i,#[1])); }
4130	else
4131	{ return(min_ass_prim_charsets(i,1)); }
4132	}
4133	example
4134	{ "EXAMPLE:"; echo = 2;
4135	ring r = 32003,(x,y,z),dp;
4136	poly p = z2+1;
4137	poly q = z4+2;
4138	ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
4139	list pr= minAssChar(i);
4140	pr;
4141	}
4142	///////////////////////////////////////////////////////////////////////////////
4143	proc equiRadical(ideal i)
4144	"USAGE: equiRadical(i); i ideal
4145	RETURN: ideal, intersection of associated primes of i of maximal dimension
4146	NOTE: designed for characteristic 0, works also in char k > 0 if it
4147	terminates, may result in an infinite loop in small characteristic
4148	EXAMPLE: example equiRadical; shows an example
4149	"
4150	{
4151	return(radical(i,1));
4152	}
4153	example
4154	{ "EXAMPLE:"; echo = 2;
4155	ring r = 32003,(x,y,z),dp;
4156	poly p = z2+1;
4157	poly q = z4+2;
4158	ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
4159	ideal pr= equiRadical(i);
4160	pr;
4161	}
4162	///////////////////////////////////////////////////////////////////////////////
4163	proc radical(ideal i,list #)
4164	"USAGE: radical(i); i ideal
4165	RETURN: ideal = the radical of i
4166	NOTE: a combination of the algorithms of Krick/Logar and
4167	Eisenbud/Huneke/Vasconcelos
4168	designed for characteristic 0, works also in char k > 0 if it
4169	terminates, may result in an infinite loop in small characteristic
4170	EXAMPLE: example radical; shows an example
4171	"
4172	{
4173	def @P=basering;
4174	int j,il;
4175	if(size(#)>0)
4176	{
4177	il=#[1];
4178	}
4179	ideal re=1;
4180	option(redSB);
4181	list qr=simplifyIdeal(i);
4182
4183	map phi=@P,qr[2];
4184	i=qr[1];
4185
4186	list pr=facstd(i);
4187	if(size(pr)==1)
4188	{
4189	attrib(pr[1],"isSB",1);
4190	if((dim(pr[1])==0)&&(homog(pr[1])==1))
4191	{
4192	ideal @res=maxideal(1);
4193	return(phi(@res));
4194	}
4195	if(dim(pr[1])>1)
4196	{
4197	execute "ring gnir = ("+charstr(basering)+"),
4198	("+varstr(basering)+"),(C,lp);";
4199	ideal i=fetch(@P,i);
4200	list @pr=facstd(i);
4201	setring @P;
4202	pr=fetch(gnir,@pr);
4203	}
4204	}
4205	option(noredSB);
4206	int s=size(pr);
4207	if(s==1)
4208	{
4209	i=radicalEHV(i,ideal(1),il);
4210	return(phi(i));
4211	}
4212	intvec pos;
4213	pos[s]=0;
4214	if(il==1)
4215	{
4216	int ndim,k;
4217	attrib(pr[1],"isSB",1);
4218	int odim=dim(pr[1]);
4219	int count=1;
4220
4221	for(j=2;j<=s;j++)
4222	{
4223	attrib(pr[j],"isSB",1);
4224	ndim=dim(pr[j]);
4225	if(ndim>odim)
4226	{
4227	for(k=count;k<=j-1;k++)
4228	{
4229	pos[k]=1;
4230	}
4231	count=j;
4232	odim=ndim;
4233	}
4234	if(ndim<odim)
4235	{
4236	pos[j]=1;
4237	}
4238	}
4239	}
4240	for(j=1;j<=s;j++)
4241	{
4242	if(pos[s+1-j]==0)
4243	{
4244	re=intersect(re,radicalEHV(pr[s+1-j],re,il));
4245	}
4246	}
4247	re=interred(re);
4248	return(phi(re));
4249	}
4250	example
4251	{ "EXAMPLE:"; echo = 2;
4252	ring r = 32003,(x,y,z),dp;
4253	poly p = z2+1;
4254	poly q = z4+2;
4255	ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
4256	ideal pr= radical(i);
4257	pr;
4258	}
4259	///////////////////////////////////////////////////////////////////////////////
4260	proc prepareAss(ideal i)
4261	"USAGE: prepareAss(i); i ideal
4262	RETURN: list = the radicals of the maximal dimensional components of i
4263	NOTE: uses algorithm of Eisenbud, Huneke and Vasconcelos
4264	EXAMPLE: example prepareAss; shows an example
4265	"
4266	{
4267	ideal j=std(i);
4268	int cod=nvars(basering)-dim(j);
4269	int e;
4270	list er;
4271	ideal ann;
4272	if(homog(i)==1)
4273	{
4274	list re=sres(i,0); //the resolution
4275	re=minres(re); //minimized resolution
4276	}
4277	else
4278	{
4279	list re=mres(i,0);
4280	}
4281	for(e=cod;e<=nvars(basering);e++)
4282	{
4283	ann=AnnExt_R(e,re);
4284
4285	if(nvars(basering)-dim(std(ann))==e)
4286	{
4287	er[size(er)+1]=equiRadical(ann);
4288	}
4289	}
4290	return(er);
4291	}
4292	example
4293	{ "EXAMPLE:"; echo = 2;
4294	ring r = 32003,(x,y,z),dp;
4295	poly p = z2+1;
4296	poly q = z4+2;
4297	ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
4298	list pr= prepareAss(i);
4299	pr;
4300	}
4301
4302	proc testPrimary(list pr, ideal k)
4303	"USAGE: testPrimary(pr,k); pr a list, result of primdecGTZ(k) or primdecSY(k)
4304	RETURN: int, 1 if intersection of the primary ideals in pr is k, 0 if not
4305	EXAMPLE: example testPrimary; shows an example
4306	"
4307	{
4308	int i;
4309	pr=reconvList(pr);
4310	ideal j=pr[1];
4311	for (i=2;i<=size(pr)/2;i++)
4312	{
4313	j=intersect(j,pr[2*i-1]);
4314	}
4315	return(idealsEqual(j,k));
4316	}
4317	example
4318	{ "EXAMPLE:"; echo = 2;
4319	ring r = 32003,(x,y,z),dp;
4320	poly p = z2+1;
4321	poly q = z4+2;
4322	ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
4323	list pr = primdecGTZ(i);
4324	testPrimary(pr,i);
4325	}
4326

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

Download in other formats:

Original Format