Context Navigation

primdec.lib @ 68e324

spielwiese

Last change on this file since 68e324 was 82716e, checked in by Hans Schönemann <hannes@…>, 26 years ago
*hannes: typos in the info-help-string git-svn-id: file:///usr/local/Singular/svn/trunk@1773 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 100.6 KB

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

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

Download in other formats:

Original Format