Context Navigation

primdec.lib @ c2868b

spielwiese

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

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

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

Download in other formats:

Original Format