Context Navigation

primdec.lib @ 8ef575

spielwiese

Last change on this file since 8ef575 was 8ef575, checked in by Hans Schönemann <hannes@…>, 25 years ago
*hannes/pfister: fixes for small characteristics git-svn-id: file:///usr/local/Singular/svn/trunk@3115 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 102.1 KB

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

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

Download in other formats:

Original Format