Context Navigation

primdec.lib @ 5480da

spielwiese

Last change on this file since 5480da was 5480da, checked in by Kai Krüger <krueger@…>, 26 years ago
Added new help for libraries git-svn-id: file:///usr/local/Singular/svn/trunk@1318 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 100.6 KB

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

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

Download in other formats:

Original Format