Context Navigation

primdec.lib @ 03f29c

spielwiese

Last change on this file since 03f29c was 03f29c, checked in by Gerhard Pfister <pfister@…>, 23 years ago
Prozeduren equidimMax, equidimMaxEHV eingebaut git-svn-id: file:///usr/local/Singular/svn/trunk@4292 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 109.2 KB

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

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

Download in other formats:

Original Format