Context Navigation

primdec.lib @ 75089b

spielwiese

Last change on this file since 75089b was 75089b, checked in by Kai Krüger <krueger@…>, 26 years ago
made fixes requested by libparse (s.a. brackets, ....) git-svn-id: file:///usr/local/Singular/svn/trunk@1198 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 100.5 KB

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

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

Download in other formats:

Original Format